INTERMATHS, VOL. 4, NO. 1 (2023), 25–47

https://doi.org/10.22481/intermaths.v4i1.12683

Article

cb licença creative commons

Simulation of the onset turbulent ﬂow around a

Isothermal Complex Geometries: an analysis of thermoﬂuid

dynamic ﬂow

Rômulo D. C. Santos

a,∗

, Quétila G. Silva

, and Samia R. L. Tananta

Instituto Federal de Educação, Ciência e Tecnologia do Acre (IFAC) (IFAC), Rio Branco - AC, Brasil

* Correspondence: lrccarvalho@gmail.com

Abstract: In this work, in the area of Computational Fluid Dynamics (CFD), more speciﬁcally in

the area of thermoﬂuid dynamics for two-dimensional ﬂows (2D), and also considering, the ﬂuid-body

interaction, allied to the phenomena of heat-transfer by mixed convection and the beginning of processes

of the turbulent ﬂow phenomenon in the ﬂuid-body interaction, a study is proposed that demonstrates

the eﬃciency in the analysis and simulation of these complex phenomena. We adopt an Eulerian

approach for a ﬁxed mesh, which is intended to represent the thermoﬂuid dynamic movement, working

together with a Lagrangian mesh, the latter being intended to discretize the immersed body. The

strategy, in this work, allows approaching complex isothermal geometries, which present a certain

aerodynamic degree on their surface, being popularly known as blunt body, where this, in turn, is

immersed in an incompressible Newtonian ﬂuid. One of the contributions of this work is the introduction

of a simple but eﬃcient method to calculate the Nusselt number. In relation to validating and modeling

the physical phenomena of interest, speciﬁcally the eﬀectiveness of the Immersed Boundary Method, we

conducted an implementation with low computational cost. This implementation was used to transfer

mixed convection heat and model turbulence using the Spalart-Allmaras model within the framework

of the URANS (Unsteady Reynolds Average Navier-Stokes) methodology. Numerical results showed

good convergence with data available in the literature, which conﬁrms the numerical precision and

reliability of the adopted model.

Keywords: Immersed Boundary Method; Mixed Convection; Onset Turbulence; Isothermal Bluﬀ

Body.

Classiﬁcation MSC: 76D05; 80A19

1 Introduction

The simulation of ﬂows around complex geometries presents a challenging task

for computational ﬂuid dynamics. Traditional body-ﬁtted numerical methods, which

strongly couple the solution of governing equations and implementation of boundary

conditions, usually need mesh generation, in general, with high computational cost. Even

considering simple geometries, the generating a high quality mesh, is not a trivial job.

To overcome the disadvantage of a strong coupling between the solution of governing

equations and the implementation of the boundary condition in the immersed boundary,

Peskin [

] presented the idea of Immersed Boundary Method (IBM) in 1971 to simulate

blood ﬂows in mitral valves of the human heart. The IBM is a numerical method that

Submitted 18 May 2023; Accepted 25 June 2023; Available online: 30 June 2023.

uses Cartesian Eulerian grid points to discretize the solution of governing equations

and Lagrangian points to represent the boundary of immersed body. Peskin’s idea

was to model the contour of the heart valve as a deformable elastic ﬁber with high

stiﬀness. A small deformation or movement of the boundary would produce a force

that tends to restore the boundary back to its original shape or position. The restoring

force at the Lagrangean point is then distributed in the surrounding ﬂuid (Eulerian

grid point) as a body force, and the Navier-Stokes equations with body strength are

resolved throughout the domain. The IBM was reviewed in [

], presenting there an

extensive list of bibliographic references. The easy computational implementation of

IBM has attracted many reseachers. Various modiﬁcations and reﬁnements have been

made. Among them, Lai and Peskin [

] presented a second-order accurate IBM and

applied it to simulate the ﬂow past a circular cylinder. The interaction between the ﬂuid

and immersed boundary was modeled by a suitable representative of the delta-function.

The work of Goldstein et al. [

] presents a method that uses control of the velocity on

the solid boundary to enforce the boundary conditions, and demonstrates its suitability

for a ﬂow past a rigid cylinder. In the process of obtaining the forcing momentum

term iteratively, two negative semi-empirically constants are required in the numerical

scheme, creating severe stability problems in the computation. The forcing momentum

and the parameters are used into nearby Eulerian grids [

],[

]–[

]. Instead of adopting

user-speciﬁed parameters, the forcing term imposing a speciﬁc boundary condition on

the body was numerically extracted, so this scheme can be classiﬁed as a discrete forcing

approach.

The work of Yang [

] proposes the use of a methodology called Virtual Boundary

Method (VBM): a bilinear interpolation procedure is combined to realize the data

transmission between grid and boundary points. The selection of constant values in the

forcing term is discussed based on the error of the velocity at the boundary. Then, the

VBM method is used to simulate the ﬂow around bluﬀ bodies. The VBM feasibility is

veriﬁed by the good agreement between the simulation’s results and other reference’s

outcomes. Mohd-Yusof [

] introduced a momentum forcing term that does not aﬀect

the stability of the discrete-time equation. The need to use a small computational

time step, which is regarded as an important advantage of this method over other

previous methods, is therefore avoided. Ye et al. [

] proposed the so-called Cartesian

grid method in non-staggered grids to simulate the unsteady and incompressible Navier-

Stokes equations in the physical domain with complex immersed boundary. In the

whole computational domain, a ﬁnite-volume method of second-order accuracy is applied

together with the two-step fractional-step procedure. Near the immersed boundary, an

interpolation procedure of second-order spatial accuracy is used.

The objectives of this study were to propose an immersed boundary method that

accounts for both momentum and heat transfer during mixed convection, taking into

consideration the onset of turbulence between a heated immersed body (referred to as

an isothermal immersed body) and the surrounding ﬂuid. Furthermore, the study aimed

to analyze the eﬀects resulting from the combination of these physical phenomena.. As

far as we know, heat-transfer by mixed convection combined with the phenomena of

turbulence around complex geometries is treated in a segregated way in many works,

some of them considering only forced convection for low Reynolds numbers. Santos

26 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

et al. [

] present the IBM for low and high Reynolds numbers coupled with the

Virtual Physical Model (VPM) proposed by Silva et al. [

] to calculate the Lagrangian

force ﬁeld, based only on the momentum equation. To simulate incompressible two-

dimensional ﬂows around heated square cylinders, a constant temperature was assumed

on the surface of the cylinders. A good numerical agreement was achieved, with a

margin of error of less than 3% compared to these previous works. The time evolution

of the drag and lift coeﬃcients, as well as, the Nusselt number, were obtained with this

methodology, being the parameters obtained from the Eulerian ﬁelds. The geometry

used in this work has singularities, which were taken into account in the construction of

the algorithm/code. This fact holds signiﬁcant importance as it enables its applicability

to various geometries. Across all simulations, the results consistently indicate that the

inﬂuence of the immersed heated body’s surface on the ﬂow intensiﬁes with higher

Reynolds numbers. A turbulence model was also used for transfer process between the

largest and smallest turbulent scales.

For this work, we will consider Newtonian and incompressible ﬂuids with constant

properties. The buoyancy term based on the Boussinesq approximation is presented,

considering the case of mixed convection. The momentum and heat transformations

between the Lagrangian and Eulerian coordinates were realized by using a smooth Dirac

delta-function, so that the nearby grid points are locally subject to momentum and

heat transfer. The stability regimes were validated for two empirical constants. The

governing equations were independently solved in each grid system coupled. The energy

generation term is disregarded/ignored in this context, because neither the eﬀects of

internal heat (for example, absorption or emission radiation) nor the humidity (which

could be responsible for the latent heat exchange) are considered. Coriolis force or

rotation eﬀects are also absent in this paper. The proposed methods were validated

for natural and forced (mixed) convection, with the isothermal and constant heat ﬂux

boundary conditions. Finally, the heat-transfer and the onset of turbulence around

complex geometry with the surrounding ﬂuid was examined. The numerical results

obtained are both accurate and stable, demonstrating a strong agreement with the

existing ﬁndings in the literature.

2 Mathematical methodology

2.1 Mathematical formulation for the ﬂuid motion and temperature

Consider a viscous incompressible ﬂow in a two-dimensional domain Ω containing an

immersed boundary in the form of closed curve Γ, as shown in Fig. (1). Like in previous

works, the forcing term is also introduced in the momentum equations and the heat

source term is incorporated in the energy equation. As a result, the governing equations

describing the forced/natural convection can be written in the primitive variables form

∂u

∂t

+ (u • ∇) u

= −∇p + µ∇

u + f , (1)

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 27

∂u

∂t

+ (u • ∇) u

= −∇p + µ∇

u + ρ

g (1 − β (T − T

∞

)) j + f , (2)

∂T

∂t

+ (u • ∇) T

= k∇

T + q , (3)

∇ • u = 0 , (4)

Fig. 1. A two-dimensional domain Ω containing an ishothermal bluff body (isothermal square

cylinder) Γ with non-uniform grid.

where Eqs. (1), (3) and (4) are for forced convection, while Eqs. (2), (3) and (4) are for

natural convection (in Eq. (2), the Boussinesq approximation is used). The symbols (i)

and

∞

denote velocity vector, pressure, temperature and reference temperature,

respectively, (ii)

and

are ﬂuid density at temperature

∞

, viscosity,

thermal diﬀusivity, thermal expansion coeﬁcient and speciﬁc heat at constant pressure,

respectively, (iii) g is a downward gravitational acceleration, (iv) the term

β (T − T

∞

)

accounts for the eﬀects of the ﬂuid temperature on the ﬂuid ﬂow, and (v) j is the unit

vector in the positive y-axis direction, respectively.

The governing equations are non-dimensionalized by using characteristic scales such

as the ﬂuid density

for the characteristic density, the far-ﬁeld velocity

∞

for

the characteristic velocity, and the body diameter

for the characteristic length.

The following characteristic scales are introduced:

D/U

∞

for the time,

∞

for the

pressure, and

∞

for the Eulerian momentum force f or simply forcing term. The

dimenionless incompressible Navier-Stokes equations and energy equation are

28 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

∂u

∂t

+ (u · ∇) u = −∇p +

∇

u +

gT + f , (5)

∂T

∂t

+ u · ∇T =

Re · P r

∇

T + q , (6)

where

is the Reynolds number, deﬁned mathematically by

∞

D/µ

, and

is the Grashof number, deﬁned here by

= g

β (T − T

∞

) D

/ν

P r

denotes the

Prandtl number, deﬁned by

P r

ν/α

, where

and

are the thermal diﬀusivity and

the kinematic viscosity, respectively. The non-dimensional continuity equation takes the

same form as the dimensional continuity equation.

The forcing term f in the momentum Eqs. (1) and (2) is the force density at the ﬂuid

(Eulerian grid) point distributed from the force density at the boundary (Lagrangian)

point, which can be expressed as

f (x, t) =

F (X (s) , t) δ (x − X (s, t)) ds , (7)

where x is the Eulerian coordinate. X and F stand for Lagrangian coordinate and

the boundary force density respectively.

δ (x − X (s, t))

is the Dirac delta function

representing the interaction between ﬂuid and the immersed boundary. Similarly, the

heat source/sink term

in the energy, Eq. (3), is the heat density transferred to the

ﬂuid from the heat ﬂux at the immersed boundary, which can be written as

q (x, t) =

Q (X (s) , t) δ (x − X (s, t)) ds . (8)

Here, Q (X (s) , t) is the virtual boundary heat ﬂux.

The boundary condition-enforced immersed boundary method developed in this paper

involves both velocity correction and temperature correction to satisfy the boundary

conditions. In the following, the velocity correction procedure will be describe ﬁrst,

followed by the temperature correction procedure.

2.2 Velocity correction procedure

With the help of fractional step technique, the velocity ﬁeld can be obtained in the

following predictor-corrector steps:

(1) Predictor step:

We start with the usual Navier-Stokes equations without the forcing term,

∂u

∂t

+ (u • ∇) u

= −∇p + µ∇

u , (9)

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 29

∂u

∂t

+ (u • ∇) u

= −∇p + µ∇

u + ρ

g (1 − β (T − T

∞

)) j , (10)

By solving the above Eq. (9) or (10) together with Eq. (4), using the projection

method developed by Chorin [

], we obtain the predicted velocity u

∗

. The basic

solution procedure is described as follows.

To advance the solution from time level

+ 1, the Crank-Nilson scheme is applied

for temporal discretization, and the immediate velocity

u can be given by

˜u − u

δt

= −

[(

u · ∇)

u + (u

· ∇) u

] +



∇

u + ∇



, (11)

˜u − u

δt

= −

[(

u · ∇)

u + (u

· ∇) u

] +



∇

u + ∇







(1 − β (T

− T

∞

)) j . (12)

Then, the immediate velocity

is corrected to the predicted velocity u

∗

by the

following form

∗

−

δt

= −∇p

n+1

. (13)

The Richardson’s number (

Gr/Re

), present in Eq. (12), expresses the relationship

between the buoyancy of a ﬂuid in relation to its viscosity. In this work, the Richardson’s

number is based on the distance between two isothermal plates,

= g

β (T − T

∞

) D

/ν

which represent the walls of a computing domain, separated by a characteristic distance

. Before using Eq. (11), the pressure ﬁeld should be obtained ﬁrst by solving the

following Poisson equation

∇

n+1

= ρ

∇ ·

δt

. (14)

In the implementation process, a staggered grid is used for the ﬁnite volume discretiza-

tion of Eqs. (9) and (10), i.e. a scalar variable pressure-like is deﬁned at the center

of each cell, while velocity components are deﬁned at the cell faces. For the spatial

derivatives, the well-know central diﬀerence scheme is employed.

(2) Corrector step:

The corrected velocity ﬁeld can be obtained by solving the following equation with

the force term f ,

∂u

∂t

= f . (15)

30 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

In the conventional IBM, the force density f is determined in advance and then the

corrected velocity u can be obtained from Eq. (15). However, there is no guarantee

that the velocity at the Lagrangian point interpolated from this corrected velocity u

always satisﬁes the non-slip boundary condition. This drawback can be overcome by

the present method. The basic idea is that the force density f is considered as unknown,

which is evaluated implicitly in such a way that the velocity u

(X (s) , t)

at the boundary

(Lagrangian) point interpolated from the corrected velocity u at the Eulerian grids equals

the given boundary velocity U

By applying forward Euler method to approximated

the temporal derivative in Eq. (15), it is clear that adding a forcing term f is equivalent

to making a correction δu in the velocity ﬁeld. By deﬁning

u = u

∗

+ δu . (16)

The Eq. (15) gives

δu

δt

= f , (17)

where

δt

is the time step size, u is the corrected velocity, u

∗

is the intermediate velocity,

is the velocity correction. Suppose that the immersed boundary is represented by a set of

Lagrangian points X

(i = 1, 2, · · · , M)

, and the ﬂuid ﬁeld is covered by a ﬁxed uniform

Cartesian mesh x

(j = 1, 2, · · · , N)

with mesh spacing ∆

= ∆

. Furthermore, let

δ (x − X (s, t)) be smoothly approximated by a continous kernel distribution

= D



− X



− X

− Y

. (18)

Here,

is the mesh size of Eulerian mesh, and

δ (r)

is proposed by Lai and Peskin [

]

δ (r) =













3 − 2 |r| +

1 + 4 |r| − 4r



, |r| ⩽ 1



5 − 2 |r| +

−7 + 12 |r| − 4r



, 1 < |r| ⩽ 2

0 , |r| > 2 .

(19)

Note from Eq. (7) that the force density f at the Eulerian point is distributed from

the boundary force F through the Dirac delta function

(δ (x − X (s, t)))

interpolation.

Eq. (7) can be written in the discrete form as

f (x

, t) =



, t





− X



∆s

, (j = 1, 2, · · · , N; i = 1, 2, · · · , M) .

(20)

Substituting Eq. (20) into Eq. (17), we have

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 31

δu (x

, t) =

F (X

, t) δt



− X



∆s

, (i = 1, 2, · · · , M) , (21)

where

, t)

is the velocity correction at the

Eulerian point, ∆

is the

boundary

segmente length. It is noted that the known quantities in Eq. (21) are the boundary

forces F

(abbreviation for F

, t)

). To satisfy the non-slip boundary condition, the

velocity at the boundary point interpolated from the corrected ﬂuid velocity ﬁeld via

the smooth delta function

must be equal to the boundary velocity U

, t)

at the

same position



, t



u (x

, t) D



− X



, (j = 1, 2, · · · , N) . (22)

Here, u

, t)

is the corrected velocity at the

Eulerian grid point, expressed as

the sum of predicted velocity u

∗

, t) and velocity correction δu (x

, t)

u (x

, t) = u

∗

, t) + δu (x

, t) . (23)

Substituting Eqs. (23) and (21) into Eq. (22), gives



, t



∗

, t) D



− X



F (X

, t) δt



− X



∆s



− X



, (24)

for i = 1, 2, · · · , M; j = 1, · · · , N; k = 1, · · · , M .

The system of equations for the boundary force can be written in the following matrix

form

A X = B , (25)

where

A =

δt







∆s

· · · D

∆s

· · · D

∆s

· · · D

∆s













· · · D







(26)

32 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

B =













− h







· · · D













∗







(27)

X =













(28)

and U

(i = 1, · · · , M)

, u

∗

(j = 1, · · · , N)

are the abbreviations for U

, t)

and

∗

, t)

, respectively. In this way, the boundary force F

(i = 1, · · · , M)

is evaluated

directly and implicity so that the satisfying of no-slip condiction is accurately achieved. It

can be observed that the entries of matrix A depend only on the coordinate information

of the Lagrangian boundary points and their adjacent Eulerian points. By solving the

equation system (26), using a direct method or iteractive method, the unknow boundary

force F

(i = 1, · · · , M)

at all Lagrangian boundary points are obtained simultaneously,

which are then substituted into Eqs. (21) and (23) to calculate the velocity correction

δu

and the corrected velocity u

(j = 1, · · · , N).

A remarkable advantage of the present method is that the force on the Lagrangian

points can be directly derived, since it is a part of the solution. We know that

(i = 1, · · · , M)

is the boundary force density acting on the ﬂuid, and due to the fact

that action and reaction are equal and opposite, it is obvious that the force exerceted

by ﬂuid on the boundary should be −F

(i = 1, · · · , M) at each Lagrangian point.

2.3 Temperature correction procedure

Similarly, using the fractional step approach, the solution of Eq. (3) can be also

obtained by predictor-corrector steps. The predictor step solves the following energy

equation without the heat source term q ,

∂T

∂t

+ (u • ∇) T

= k∇

T . (29)

The resultant solution from Eq. (29) is noted as predicted temperature

∗

(x, t) ,

which can be obtained by solving the following equation obtained by the Crank-Nicolson

scheme

∗

− T

δt

= −

h

n+1

· ∇



∗



n+1

· ∇





∇

∗

+ ∇



. (30)

Similar to the velocity correction procedure, the corrected temperature

T (x, t)

obtained by solving the following equation,

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 33

∂T

∂t

= q . (31)

It can be clearly observed from Eq. (31) that adding a heat source/sink is equivalente

to making a temperature correction. Applying Euler explicit scheme to Eq. (31) gives

δT

δt

= q , (32)

with

δT (x, t) = T (x, t) − T

∗

(x, t) , (33)

where

δT (x, t)

is the temperature correction and

T (x, t)

is the corrected temperature.

Note from Eq. (8) that the heat source/sink,

at the Eulerian grid point is distributed

from the boundary heat ﬂux,

through Dirac delta function interpolation, which can

be expressed in the discret form as

q (x

, t) =



, t





− X



∆s

, (i = 1, · · · , M; j = 1, · · · , N) . (34)

Substituting Eq. (33) into Eq. (32) leads to

δT (x

, t) =

Q (X

, t) δt

ρc



− X



∆s

, (i = 1, · · · , M; j = 1, · · · , N) .

(35)

It is noted that the unknows in Eq. (35) are the boundary heat ﬂuxes

Q (X

, t)

. To

satisfy the physical boundary condiction, we have to make sure that the temperature

at the boundary point interpolated from the corrected temperature ﬁeld by the delta

function D

is equal to the speciﬁed temperature T

, t), that is,



, t



T (x

, t) D



− X



(i = 1, · · · , M; j = 1, · · · , N) . (36)

Substituting Eqs. (33) and (35) into Eq. (36) gives



, t



∗

, t) D



− X





, t



δt

ρc



− X



∆s



− X



, (37)

with

= 1

, · · · , M

;

= 1

, · · · , N .

The Eq. (37) can be put in the following matrix form

34 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

C Y = D , (38)

where

C =

δt

ρc







∆s

· · · D

∆s

· · · D

∆s

· · · D

∆s













· · · D







(39)

D =













− h







∆s

· · · D

∆s

· · · D

∆s

· · · D

∆s













∗







(40)

Y =













(41)

with

(i = 1, · · · , M)

∗

(j = 1, · · · , N)

and

(i = 1, · · · , M)

, representing

T (X

, t)

∗

, t) and Q (X

, t), respectively.

As shown above, the essence of temperature correction procedure is that the heat

source/sink is determined implicitly in such a way that the temperature

T (X (s) , t)

the boundary point interpolated from the corrected temperature ﬁeld

T (x, t)

equals to

the speciﬁed boundary temperature

. In other words, the physical boundary condition

is enforced. The system of linear equations (38) can be solved by a direct method or

an interactive method. After

at all Lagrangian points are obtained, they are then

substituted into Eq. (35) to obtain the temperature correction

δT

, which are further

substituted into Eq. (33) to get ﬁnally the corrected temperature T

(j = 1, · · · , N).

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 35

3 Evaluating the average number of Nusselt

Nusselt number is a key parameter in the heat-transfer problem. The local Nusselt

number Nu is deﬁned as

Nu (X (s) , t) =

h (X (s) , t) L

, (42)

where

is a typical characteristic length. According to Newton’s cooling law and

Fourier’s law, the heat conducted away from the wall by the ﬂuid is equal to the heat

convection from the wall, that is

−k

∂T

∂n

(X (s) , t) = h (X (s) , t) (T

− T

∞

) . (43)

Substituting Eq. (43) into Eq. (42) one obtains

Nu (X (s) , t) =

− T

∞

−k

∂T

∂n

(X (s) , t) L

= −

− T

∞

∂T

∂n

(X (s) , t) . (44)

The average Nusselt number

is an important parameter to examine the rate of

heat-transfer from the heated surface and is deﬁned by

Nu =

Nu · ds =

−

− T

∞

∂T

∂n

(X (s) , t) ds , (45)

where

is the total length of the immersed boundary. As shown in Eq. (44) and

(45), the evaluation of the local and average Nusselt numbers involves the calculation of

the temperature gradient at the boundary points.

3.1 Estimation of the average Nusselt at Eulerian points

As shown in Eq. (8), the heat source at one Eulerian points is from the heat ﬂux at

a few Lagrangian points. If we consider all the Eulerian points which receive the heat

from the immersed boundary, then from energy conservation law, we have

Ω

q (x, t) dx =

Q (X (s) , t) ds , (46)

in which the l.h.s. (left hand side) is the volume integral in the whole ﬂuid domain, and

the r.h.s. (right hand side) is the surface integral along the cylinder wall. According to

Fourier’s law, Q (X (s) , t) can be written as

Q (X (s) , t) = −k

∂T

∂n

(X (s) , t) . (47)

36 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

Substituting Eqs. (32) and (47) into Eq. (46) gives

Ω

ρc

δT

δt

dx =

− k

∂T

∂n

(X (s) , t) ds . (48)

The r.h.s. of Eq. (48) can be also be written in terms of local Nusselt numbers

Nu (X (s) , t) as

− k

∂T

∂n

(X (s) , t) ds =

k · Nu (X (s) t) (T

− T

∞

)

ds =

Nu (T

− T

∞

) kL

(49)

Eq. (48) can then be simpliﬁed to

Nu (T

− T

∞

)

kL =

Ω

ρc

δT

δt

dx . (50)

Thanks to the approximation of the volume integral oﬀered by the mean theorem, the

equality (50) can be ﬁnally written as

Nu =

ρc

δT

δt

∆x

∆y

kL (T

− T

∞

)

, (j = 1, · · · , N) . (51)

In Figs. (2) and (3), we present the behaviour of the ﬂow around the neighborhood

of the ishotermal square cylinder (bluﬀ body) at diﬀerent values of

(

= 1

2 and

Re = 10) and Ri (Ri = 0 and Ri = 0, 5).

(a) (b)

Fig. 2. Streamlines for ﬂow around ishotermal square cylinder at Re = 1 and Re = 2 for different

Richardson numbers (Ri).

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 37

Fig. 3. Streamlines for ﬂow around ishotermal square cylinder at

= 10 for different Richardson

numbers (Ri).

The ﬂuid begins the process of separation of bluﬀ body when

Re ≥

2, for diﬀerent

values. A detachment process from the circulation bubble is also observed in Fig. (2-(c))

for Re = 10 with Ri = 0.

Other works using diﬀerent numerical methodologies have presented a similar scenario;

for example, see Chatterjee [

], where it is conducted a study for a pair of square

cylinders in a tandem vertical channel conﬁguration, with low Reynolds numbers. Here,

the author concluded that the beginning of separation of the ﬂow occurs for

= 1

−

= 2

−

3 and

= 3

−

4, considering the same values for the Richardson number

presented in our work. Therefore, in this work, it was observed that, as the Reynolds

number is (gradually) increased, the separation of the ﬂow occurs at the trailing edge of

the immersed ishotermal body and the two (symmetrical) vortices of the recirculation

bubble begin to be formed downstream of the obstacle, even under the permanent

outﬂow regime. This can be attributed to the increase in buoyancy, which leads to a

higher velocity gradient on the surface of the cylinder. Consequently, the pressure on

the immersed body’s surface decreases. Other signiﬁcant studies, including Singh [

]

and Moulay [

], in addition to presenting similar scenarios, there is a critical value to

be studied for the adimensionless number

, i.e. the authors take into account the

onset of the vortex detachment and analyze the inﬂuence of this parameter with the

increase of Re.

Table 1. Comparison of drag coefﬁciet (

) and average Nusselt (

N u

) numbers for isothermal

boundary condition, for different Reynolds (Re) numbers.

Re Cd Nu

Present 1 4,8500 0,6917

Anjaiah et. al. [17] 1 4,8600 0,6916

Dhiman et. al. [18] 1 4,8714 0,6916

Sharma et. al. [19] 1 4,8714 0,6928

Present 10 0,9300 1,5532

Anjaiah et. al. [17] 10 0,9437 1,5624

Continued on next page

38 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

Table 1 – continued from previous page

Dhiman et. al. [18] 10 0,9437 1,5624

Sharma et. al. [19] 10 0,9343 1,5573

Present 40 0,2497 2,6802

Anjaiah et. al. [17] 40 0,2538 2,6969

Dhiman et. al. [18] 40 0,2538 2,6969

Sharma et. al. [19] 40 0,2493 2,6012

Note that the r.h.s. of Eq. (51) are dimensionless. The actual form of Eq. (51) for

dimensionless variables depend on the reference characteristics taken for the problem

considered. In the numerical studies reported below, the speciﬁc forms of Eq. (51) will

be given.

The heat ﬂux at Lagrangian points obtained from equation system (38) can be directly

used to compute the average Nusselt number, that is,

Nu =

k (T

− T

∞

)

Q (X (s) , t) ds . (52)

The the discretized form of Eq. (52) is

Nu =

∆s

kL (T

− T

∞

)

, (i = 1, · · · , M) . (53)

Note that Eq. (53) is also based on dimensionless, variables. Its speciﬁc form for

dimensionless variables will be provided in Section 5.

4 Onset of turbulence

In addition to laminar regime ﬂow, we have also simulated (onset of turbulence) ﬂows

for Reynolds numbers going up to 5

. Thus, it is necessary to use turbulence

models to close the Navier-Stokes equations. A model has been implemented, namely

the Spalart-Allmaras model, with URANS model (Unsteady Reynolds Averaged Navier-

Stokes).

The URANS model is used to refer to RANS (Reynolds Averaged Navier-Stokes), where

in this model, the dependent variables of the Navier-Stokes equation are decomposed

into ﬁltered components and ﬂoating components, and then all terms are ﬁltered.

4.1 Spalart-Allmaras model

The Spalart-Allmaras model is a one-equation model that solves the transport equation

for the kinematic eddy turbulent viscosity [

], without involving turbulent energy,

dissipation or vorticity calculations, available in other models. Thus, the Spalart-

Allmaras model concentrates into a single parameter the behavior of turbulence, being

classiﬁed as a closed model.

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 39

The equation for the turbulent viscosity is constructed using mainly ﬂow empirical

considerations, dimensional analysis and Galileo’s principle of relativity. The model uses

a working variable,

ν, given by the following transport equation

∂

∂t

∂

∂x

ν) = c

(1 − f

)

ν +

∂

∂x

(ν +

ν)

∂

∂x

+ c

∂

∂x

∂

∂x

−



−





+ f

∆U

(54)

where the terms on the right side are, respectively: the production of turbulent viscosity,

the molecular and turbulent diﬀusion of

, the dissipation of

, the destruction of

which reduces the turbulent viscosity near the wall, and, ﬁnally, the terms that model

the transition eﬀects, indicated by sub-index

. The turbulence viscosity,

, is calculated

from the auxiliary variable of the Spalart-Allmaras model and damped by the function

along the wall, and is given by:

νf

, (55)

where

+ C

, (56)

with

χ =

. (57)

5 Results

5.1 Mixed convection around a isothermal cylinder

Forced convection over an isothermal cylinder has been extensively studied and used

as a benchmark to examine the capability of numerical methods. The experimental [

]

and numerical [

]-[

] results are available for square cylinder. The ﬂuid and heat ﬂow

are characterized by Reynolds number

ρU

∞

/µ

and Prandtl number

P r

µc

where

is the ﬂuid density,

∞

is the free stream velocity,

is the cylinder diameter,

is the dynamic viscosity,

is heat at constant pressure and

the thermal diﬀusivity.

In this work, numerical simulations are conducted for diﬀerent low Reynolds numbers

= 1

−

500, while keeping the Prandtl number ﬁxed at

P r

= 0

7. Both heat and

ﬂuid ﬂow characteristics like the drag coeﬃcient

, recirculation behind the cylinder,

streamline and isotherm patterns, average Nusselt number on the cylinder surface are

presented and compared with previous results in the literature.

A schematic view of two-dimensional heat and ﬂuid ﬂow across a heated square

cylinder is shown in Fig. (1). In the simulation, governing equations in dimensionless

form are used. The free-stream velocity

∞

and the cylinder diameter

are chosen as

the reference velocity and reference length, respectively. The temperature is normalized

∗

T − T

∞

− T

∞

, (58)

40 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

where

is the temperature on the cylinder surface and

∞

is the free-stream temperature.

The speciﬁc dimensionless forms to compute the average Nusselt number by the methods

1 and 2 abovementioned, are

Nu =

Re · P r

δT

δt

∆x

∆y

, (j = 1, · · · , N) , (59)

and

Nu =

∆s

, (i = 1, · · · , M) , (60)

respectively. Note that in Eq. (59), the dimensionless form of

is actually

−

∂T

∂n

, t)

which can be directly obtained from the solution of equation system (38). The drag

coeﬃcient is deﬁned as

/2 ρU

∞

, (61)

where

is the drag force that, in the present work, it is computed from the solution of

equation system (25), which makes the drag force calculation very easy. In the order

to validade our code, simulations with a single heated square cylinder were compared

with the results found in the literature. To mimic free boundary conditions, it was

considered

= 55

and width

= 30

. The heated isothermal square cylinder was

placed at x = 16.5d and y = 15d, to minimize boundary inﬂuence. The grid resolution

for these simulations was 318

164, 840

600 and 1120

800 points, a non-uniform

grid was used to better capture the eﬀects with a total of 202 Lagrangian grid points

inside the immersed body. The grid is uniform in the region of the cylinder, maintaining

a minimum number of 30 grid inside. The time step used in the calculation process

is in the range 1

−6

to 1

−4

, which is dynamically calculated by the

Courant-Friedrichs-Lewy (CFL) condition:

CF L =

max (|u| + c)

−0.5, x

+0.5]

× ∆t

∆x

. (62)

To avoid numerical problems, the maximum ratio of 3% grid expansion was employed

in this region in both directions. In the present work, the ratio between Lagrangian

and Eulerian grid size was maintained at 0

98. The boundary condition were prescribed

as follow: 1) uniform and nonuniform ﬂow with velocity

∞

, pressure

= 0, and

temperature

input

= 0 at the entrance of the domain (left side); 2) ﬂow fully developed

in the output

(∂u/∂x = 0, ∂T/∂x = 0)

; 3) condition of symmetries at the upper and

lower boundaries

(∂u/∂y = v = 0)

(∂T/∂y = 0)

. Concerning the initial condition, we

have considered

= 0, at

= 0 (time), for the entire computational domain. The

square cylinder is maintained at a constant dimensionless temperature equal to 1, i. e.,

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 41

= 1, while the ﬂuid has an initial temperature equal to zero (

= 0). Although the

boundary condition at the output is not a reﬂexive condition, the simulation results

successfully corroborate the formation of large (vortex) structure outside the domain

without any reﬂection. The current pressure boundary conditions at the inlet and outlet

are imposed to ensure consistency with the equations for velocities, which is achieved

through the arrangement of the staggered grid.

5.2 Mixed convection around isothermal cylinders in tandem

In this section, the ﬂow around a pair of heated circular cylinders with diﬀerent

conﬁgurations is studied, where the heat transfer process around obstacles has its

importance and relevance in engineering problems. Here, two cylinder conﬁgurations

will be considered, where in both cases the two cylinders have equal diameters and the

same center-to-center distance (

). In the ﬁrst case, the angle formed by the segment

joining the centers of the two cylinders and the axis of the abscissa is zero.

5.2.1 Description of the problem and boundary conditions

In the Fig. (4) the two cylinders are identical and ﬁxed with the same diameter,

maintained in tandem conﬁguration, the cylinder B, downstream of the cylinder A.

Fig. 4. Ilustration of the computational domain with two cylinders in tandem conﬁguration.

In this case, the cylinders is conﬁned to a channel with free ﬂow, with uniform velocity

(

∞

) and constant temperature (

(> T

∞

)

). The horizontal and vertical spacing between

the cylinders are ﬁxed in

= 16

and

= 19

, repectively. These values are

chosen to reduce the eﬀect of boundary condictions on the inlet and outlet relative to

the ﬂow pattern at the cylinder boundary. Moreover, these choices are also consistent

with other contemporary works available in the literature, such as, for example, Mahir

and Altaç [

], Sohankar and Etminan [

], Santos [

], Laidoudi and Bouzit [

] and

Santos et. al. [

]. The drag (

) and lift (

) coeﬃcients for the calculation of the each

cylinder are performed as follows:

42 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

= C

+ C

ρU

∞

, (63)

= C

+ C

ρU

∞

, (64)

where

and

represent the lift coeﬁcients due to pressure and viscous forces,

respectively.

In a similar way,

and

represent the drag coeﬁcients due to the pressure and

viscous forces.

and

are forces of drag and lift, respectively, acting on the surface of

the cylinder. Thus, the drag and lift coeﬃcients can be obtained from the expressions:











= 2

− p

) dy ,





(

∂u

∂y

∂u

∂y

)

dx +







∂u

∂x

∂u

∂x











(65)











= 2

− p

) dy ,











∂v

∂x

∂v

∂x







dx +

(

∂v

∂y

∂v

∂y

)





(66)

The subscripts

and

refer to "front", "back", "bottom" and "top" surfaces of the

cylinders, respectively.

The main results for the simulations can be summarized as follows:

•

A wake is formed upstream of the second cylinder, and it is necessary to investigate

whether it can be reduced or eliminated by increasing the distance between the

cylinders;

•

The isothermal lines reﬂect the same behavior of the pattern of the streamlines

(current lines);

•

The average Nusselt number increase for

= 500 for diﬀerent values of

, even

keeping the distance between the cylinders;

•

The thermal buoyancy is suppressed in the recirculation zones of the tandem

cylinders, even with a mounting angle;

•

The thermal buoyancy tends to increase the coeﬃcient of drag and the average

Nusselt number of the cylinder more than the second;

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 43

5.2.2 Variations of the Nusselt number

One of the main objectives of heat-transfer calculations involving cylinders is to

determine the local and overall heat-transfer around isothermal cylinders. The eﬀect of

the ﬂow, especially regarding the heat-transfer, can be better observed by analyzing the

local heat-transfer coeﬃent, also known as the Nusselt local number.

In the Fig. (5), for diﬀerent Richardson numbers, the distributions of Nusselt numbers

along the perimeter of the upstream and downstream cylinders are provided. For

= 3,

= 100,

= 200 and

= 500, for diﬀerent Richardson numbers, the

local distributions of the Nusselt number along the perimeter of the upstream and

downstream cylinders are provided. For

= 3, the heat transfer characteristics

between two closely spaced cylinders diﬀer signiﬁcantly. While the local Nusselt number

proﬁle of the upstream cylinder resembles that of an isolated cylinder, the downstream

cylinder exhibits distinct behavior. This discrepancy arises due to the strong connection

between heat transfer and ﬂuid ﬂow. Speciﬁcally, we observed localized regions of

minimal heat transfer at the front and back stagnation points of the downstream

cylinder. These points experience relatively low ﬂow velocities, which impact the heat

transfer process.

(a) (b)

Fig. 5. Local variation of the Nusselt number to the same dimensionless instant. (a) -

= 100,

200 and

= 500 for

= 0 (forced convection) and (b) -

= 500 for

= 1

= 2

and Ri = 5.0 (natural convection).

Thus, in Fig. (5-(a)), the maximum heat-transfer from the downstream cylinder is

exhibited with a double protuberance occuring in

θ ≈

◦

and

θ ≈

265

◦

from the

cylinder wall, where thermal layers (also known as thermal plumes) and hydrodynamics

becomes thinner. The formation of vortices in the downstream region of the cylinder

coincides with the oscillation of the average Nusselt number from large amplitude to

low amplitude during a vortex release period for

= 3 and

= 500 for diﬀerent

values of Ri as see in Fig. (5-(b)).

It is important to note that although the Nusselt’s local distribution of the downstream

44 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

cylinder resembles that of the upstream cylinder, typiﬁed as a large protuberance, its

magnitude is smaller than of the upstream cylinder, indicating smaller heat-to-cylinder

convection to downstream.

6 Conclusions

In this work, a boundary condition-enforced immersed boundary method is developed

for simulations of heat and mass transfer problems. The eﬀect of thermal boundaries

in the ﬂow and temperature ﬁelds is considered through the velocity and temperature

corrections. The temperature correction is evaluated implicit in such a way that the

temperature at the immersed boundary, interpolated from the corrected temperature

ﬁeld, satisﬁes the physical boundary conditions.

For the momentum transfer between the immersed body and the surrounding ﬂuid,

the additional momentum forcing obtained by using the forcing term is added to the

ﬂuid-body equations. To model the onset of turbulence, the Spalart-Allmaras model is

used. This model uses URANS concept, with only one transport equation for turbulence

viscosity, being calibrated in pressure gradient layers. A computational code was

developed to implement the mentioned methodology and analyze the combined heat-

transfer phenomena and the onset of turbulence for thermoﬂuid dynamics interactions

around complex isothermal geometries. The agreement of the results with the available

data in the literature validates the numerical method.

Abbreviations

IBM = Immersed Boundary Method

l.h.s. = left hand side

r.h.s. = right hand side

RANS = Reynolds Averaged Navier-Stokes

Disclosure statement. The authors declare no conﬂict of interest in the writing of the

manuscript, or in the decision to publish the results.

ORCID

Rômulo D. C. Santos https://orcid.org/0000-0002-9482-1998

Quétila G. Silva https://orcid.org/0009-0000-5928-507X

Samia R. L. Tananta https://orcid.org/0009-0004-6842-5221

References

Peskin, C. S. (1972). Flow Patterns Around Heart Valves: A Numerical Method. J. Com-

put. Phys., 10(2), 252–271. https://doi.org/10.1016/0021-9991(72)90065-4

Peskin, C. S. (1977). Numerical Analysis of Blood Flow in the Heart. J. Comput. Phys.,

25, 220–252.

Park, S. G., Chang, C. B., Kim, B. and Sung, H. J. (2017). Simulation of Fluid-

Flexible Body Interaction with Heat Transfer. Int. J. Heat Mass Transfer, 110, 20–33.

https://doi.org/10.1016/j.ijheatmasstransfer.2017.03.012

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 45

Ashraﬁzadeh, A. and Hosseinjani, A. A. (2017). A Phenomenological Study

on the Convection Heat Transfer around Two Enclosed Rotating Cylinders via

an Immersed Boundary Method. Int. J. Heat Mass Transfer, 107, 667–685.

https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.078

Mittal, R. and Iaccarino, G. (2005). Immersed Boundary Method. Annu. Rev. Fluid Mech.,

37, 239–261. https://doi.org/10.1146/annurev.ﬂuid.37.061903.175743

Lai, M. C., & Peskin, C. S. (2000). An Immersed Boundary Method with Formal Second-

Order Accuracy and Reduced Numerical Viscosity. J. Comput. Phys., 160(2), 705–719.

https://doi.org/10.1006/jcph.2000.6483

Goldstein, D., Handler, R., & Sirovich, L. (1993). Modeling a No-slip Flow

Boundary with an External Force Field. J. Comput. Phys., 105(2), 354–366.

https://doi.org/10.1006/jcph.1993.1081

Yang, Q. and Cao, S. (2013). Numerical Simulation of Flow around Bluﬀ Bodies Based on

Virtual Boundary Method. The 8th Asia–Paciﬁc Conference on Wind Engineering, Chennai,

pp. 10-14, 582-591. https://doi.org/10.3850/978-981-07-8012-8_154

Mohd-Yusof, J. (1997). For simulations of ﬂow in complex geometries. Annual Research

Briefs, 317.

10.

Ye, T., Mittal, R., Udaykumar, H. S., & Shyy, W. (1999). An Accurate Cartesian Grid

Method for Viscous Incompressible Flows with Complex Immersed Boundaries. J. Com-

put. Phys., 156(2), 209–240. https://doi.org/10.1006/jcph.1999.6356

11.

Santos, R. D., Gama, S. M., & Camacho, R. G. (2018). Two-Dimensional Simu-

lation of the Navier-Stokes Equations for Laminar and Turbulent Flow around a

Heated Square Cylinder with Forced Convection. Applied Mathematics, 9(03), 291–312.

https://doi.org/10.4236/am.2018.93023

12.

Silva, A. L. E., Silveira-Neto, A. and Damasceno, J. J. R. (2003). Numerical Simulation of

Two-Dimensional Flows over a Circular Cylinder using the Immersed Boundary Method.

J. Comput. Phys., 189, 351–370. https://doi.org/10.1016/s0021-9991(03)00214-6

13.

Chorin A. J. (1968). Numerical Solution of the Navier-Stokes Equations. Mathematics of

computation;22(104):745–762. https://doi.org/10.1090/s0025-5718-1968-0242392-2

14.

Chatterjee, D. (2010). Mixed Convection Heat Transfer from Tandem Square Cylinder in

Vertical Channel at Low Reynolds Numbers. Numerical Heat Transfer, Part A: Applications,

58(9), 740–755. https://doi.org/10.1080/10407782.2010.516703

15.

Singh, S., Panigrahi, P. & Muralidhar, K. (2007). Eﬀect of Bouyancy on the Wakes of

Circular and Square Cylinders: A Schlieren-Interferometric Study. Experiments in ﬂuids,

43(1):101–123. https://doi.org/10.1007/s00348-007-0329-8

16.

Moulay, M. A., Belkady, M., Aounallah, A. et.al. (2017). Eﬀect of Opposing Bouyancy on

Two-dimensional Laminar Flow and Heat Transfer Across a Conﬁned Circular Cylinder.

Mechanics, 23(6):859–865. https://doi.org/10.5755/j01.mech.23.6.17291

17.

Anjaiah, N., Dhiman, A. & Chhabra, R. (2006). Mixed Convection Heat Transfer from

a Square Cylinder to Power-Law Fluids in Cross-Flow. In AMSE 2006 2nd Joint US-

Europen Fluids Engineering Summer Meeting Collocated with the 14

International Con-

ference on Nuclear Engineering, pp. 1435–1444. American Society of Mechanical Engineers.

https://doi.org/10.1115/fedsm2006-98072

18.

Dhiman, A., Anjaiah, N., Chhabra, R. & Eswaran, V. (2007). Mixed Convection from

Heated Square Cylinder to Newtonian and Power-Fluids. Journal of Fluids Engineering,

129(4):506-513. https://doi.org/10.1115/1.2436586

19.

Sharma, N., Dhiman, A. K., & Kumar, S. (2012). Mixed Convection Flow and

Heat Transfer Across a Square Cylinder Under the Inﬂuence of Aiding Buoy-

ancy at Low Reynolds Numbers. Int. J. Heat Mass Transfer, 55(9-10), 2601–2614.

https://doi.org/10.1016/j.ijheatmasstransfer.2011.12.034

20.

Breuer, M., Bernsdorf, J., Zeiser, T. & Durst. F. (2000). Accurate Computations of

46 | https://doi.org/10.22481/intermaths.v4i1.12683 R D. C. Santos; Q. G. Silva; S. R. L. Tananta

the Laminar Flow Past a Square Cylinder Based on Two Diﬀerent Methods: Lattice-

Boltzmann and Finite Volume. International Journal of Heat and Fluid Flow, 21(2):186–196.

https://doi.org/10.1016/s0142-727x(99)00081-8

21.

Kim, J., Kim, D. & Cho, H. (2001). An Immersed-Boundary Finite-Volume Method

for Simulations of Flow in Complex Geometries. J. Comput. Phys., 171(1):132–150.

https://doi.org/10.1006/jcph.2001.6778

22.

Kim, J. & Cho, H. (2004). An Immersed-Boundary Finite-Volume Method for Simulation

of Heat Transfer in Complex Geometries. KSME International Journal, 18(6):1026–1035.

https://doi.org/10.1007/bf02990875

23.

Perumal, D. A., Kumar, G. V., & Dass, A. K. (2012). Numerical Simulation of Viscous Flow

over a Square Cylinder Using Lattice Boltzmann Method. ISRN Mathematical Physics,2012.

https://doi.org/10.5402/2012/630801

24.

Sharma, A. (2003). Numerical Investigation of Unconﬁned and Channel-Conﬁned Flow

Across a Square Cylinder with Forced and Mixed Convection Heat Transfer. Ph.D. thesis,

Indian Institute of Technology Kanpur, India.

25.

Sharma, A. & Eswaran, V. (2004). Eﬀect of Aiding and Opposing Buoyancy on the Heat

and Fluid Flow Across a Square Cylinder at

= 100. Numerical Heat Transfer, Part A:

Applications, 45(6):601–624. https://doi.org/10.1080/10407780490277798

26.

Sohankar, A. & Etminan, A. (2009). Forced-Convection Heat Transfer from Tandem Square

Cylinders in Cross Flow at low Reynolds Numbers. International Journal for Numerical

Methods in Fluids, 60(7):733–751. https://doi.org/10.1002/ﬂd.1909

27.

Yen, S., San, K., & Chuang, T. (2008). Interactions of Tandem Square Cylinders

at Low Reynolds Numbers. Experimental Thermal and Fluid Science, 32(4):927–938.

https://doi.org/10.1016/j.expthermﬂusci.2007.07.001

28. Smagorinsky, J. (1963) General Circulation Experiments with the Primitive Equations: I.

The Basic Experiment. Monthly Weather Review, 91, 99–164. https://doi.org/10.1175/1520-

0493(1963)091%3C0099:gcewtp%3E2.3.co;2

29.

Spalart, P. & Allmaras, S. (1992). A One-Equation Turbulence Model for Aerodynamics

Flows. Recherce Aeropatiale, No. 1, 5–21.

30.

Mahir, N. & Altaç, Z.(2008). Numerical Investigation of Convective Heat Transfer in

Unteady Flow Past Two Cylinders in Tandem Arragements. International Journal of Heat

and Fluid Flow, 29(5):1309–1318. https://doi.org/10.1016/j.ijheatﬂuidﬂow.2008.05.001

31.

Sohankar, A., Norberg, C. & Davidson, L. (1999). Simulation of Three-Dimensional Flow

Around a Square Cylinder at Moderate Reynolds Numbers. Physics of Fluids, 11(2):

288–306. https://doi.org/10.1063/1.869879

32.

Santos, R. D. C. d. (2014). Análise Bidimensional Termo-Fluido Dinâmica de Cilindros

Rotativos com o Método da Fronteira Imersa/Modelo Físico Virtual. Master’s thesis.

Universidade Federal de Itajubá - Minas Gerais, Brazil.

33.

Laidoudi, H. & Bouzit, M. (2016). Mixed Convection in Poiseuille Fluid from an Asymmet-

rically Conﬁned Heated Circular Cylinder. Thermal Science, (00): 172–172.

34.

Zdravkovich, M. (1985). Flow Induced Oscillations of Two Interfering Circular Cylin-

ders. Journal of Sound and Vibrations, 101(4): 511–521. https://doi.org/10.1016/s0022-

460x(85)80068-7

CC BY-NC 4.0

R D. C. Santos; Q. G. Silva; S. R. L. Tananta INTERMATHS, 4(1), 25–47, June 2023 | 47