Two-Dimensional Isothermal and Newtonian Flow in Complex Geometries

Rômulo D. C. Santos¹, Sílvio M. A. Gama²

¹Faculty of Civil Engineering, Campinas State University, UNICAMP, São Paulo, Brazil.
²Department of Mathematics and CMUP, Faculty of Sciences, University of Porto, Porto, Portugal.

DOI: 10.4236/am.2021.122008                

Abstract

In this paper, we investigate the thermal and turbulent behaviour of incompressible Newtonian flow, by numerical simulation, combining two physical phenomena, namely, the heat-transfer by mixed convection and the onset of turbulence, around different isothermal complex geometries, using the immersed boundary method coupled with the virtual physical model, in order to model the presence of the isothermal body. Boundary conditions of the Dirichlet and Neumann type are implemented. For turbulence modelling, the Smagorinsky and Spalart-Allmaras models are used, for Reynolds and Richardson numbers ranging up to 5000 and 5, respectively. This work confirms that, downstream of the immersed body, the recirculation: 1) increases with the increase in the number of Reynolds, keeping the number of Richardson constant, and 2) decreases with the increase in the number of Richardson, preserving the number of Reynolds constant. It also confirms the generation of thermal plumes moving upwards. For Reynolds numbers in the order of a few hundred and Richardson numbers around 5, it is observed, for tandem cylinders, the vortex wake being established in the downstream region. Interactions within the vortex wake, with the shear layer separated from the downstream cylinder, create two vortices near the downstream cylinder. The shear layer separating from the upstream cylinder creates a vortex behind the downstream cylinder.

Keywords

Immersed Boundary Method, Mixed Convection, Onset Turbulence

1. Introduction

Many physical phenomena in fluid mechanics can be described by mathematical modelling, using a set of partial differential equations, often non-linear, known as the laws of conservation of fluid mechanics. They are: conservation of momentum, conservation of mass (continuity) and conservation of energy. These laws model the dynamics of the forces acting on the fluid with the so-called Navier-Stokes equations, as well as the energy exchanges that occur in the different regions of the flow.

Traditional methods of domain discretization present some difficulties in terms of computational implementation and CPU resources, requiring successive remeshing processes (construction of new meshes) for each new iterative process, in addition to the need to introduce a new generalized coordinate system. Thus, considering these physical and mathematical problems, a computational code was developed for the immersed boundary methodology for the thermofluid dynamics interaction. Next, we review some of the most significant works in the area of immersed boundary method [1], considering the heat-transfer by different convection processes. These references aim the preparation of the theoretical foundation of the methodology used in this work.

Among the predecessor works involving the thermal part, we mention Badr and Dennis [2], and Badr et al. [3]. These works do not present the immersed boundary method, but they are a precursor concerning the flow around a rotating circular cylinder involving forced convection, being the theoretical basis used in other works for the implementation of the thermal part ( [4] [5] [6] ). In summary, the experimental and numerical works of Badr and Dennis consider the problem of laminar flow with heat-transfer by forced convection with a rotating circular cylinder around its own axis, initially with the Reynolds number equal to 100, at a specific rotation rate $\alpha \in \left[\mathrm{0.5,3}\right]$. Meanwhile, the second work uses the Reynolds number range varying between ${10}^3\le Re\le {10}^4$, with specific rotation rate equal to $\alpha =4$. In these both works, the cylinder is located in a uniform flow. Their findings report that the temperature fields are strongly influenced by the rotation of the cylinder and that the heat-transfer coefficient tends to decrease, as the rotation of the cylinder increases. An attribution of the presence of a layer of rotating fluid around the cylinder is also verified, being separated from the main flow of the flow. Also, their results demonstrate convergence and numerical stability compared to others available works (see, Sharma et al. [6] and Lai and Peskin [7] ).

Santos [8] presents the immersed boundary method with a virtual physical model, based on the equations of conservation of momentum, mass and energy, to model the physics of incompressible laminar flow. The motivation was to continue the work developed by Silva et al. [9], seeking to contribute to the resolution of engineering problems of flow around obstacles, implementing a computational code destined to present the heating in a rotating cylinder subject to forced convection; in other words, for the simulation of a flow around non-stationary cylinders, rotating with different values for the specific rotation varying in the interval $\left[\mathrm{0,4}\right]$, and the Reynolds number in the interval $\left[\mathrm{47,250}\right]$. The vorticity, temperature and pressure fields are obtained, as well as the drag, lift, pressure coefficients, Strouhal (which calculates the vortex shedding frequency) and Nusselt numbers. The results are compared based on the experimental work of Carvalho [10] and other numerical works.

Ashrafizadeh and Hosseinjani [11] present an extension of the immersed boundary method for the analysis of the natural and mixed convection process, for the case of two cylinders side by side with rotation velocity alternated with heating, both in a cold and square cavity containing air, with Prandtl number (Pr) equal to 0.7. The variables in this study involve the rotation rates, $\alpha$, for different numbers of Rayleigh number ( $Ra\in \left\{{10}^4{\mathrm{,10}}^5{\mathrm{,10}}^6\right\}$ ), and Richardson number ( $Ri\in \left\{\mathrm{1,2,}\cdots \mathrm{,100}\right\}$ ). Cases of pure natural convection are also considered. Isothermal lines (or isotherms), streamlines, Nusselt numbers (Nu), local and average around the solid boundary, is studied and presented in each case. An extension of the immersed boundary method is used for the boundary conditions and in the resolution of the governing equations, for a two-dimensional flow, using an orthogonal Cartesian mesh.

There are several flows inside the square cavity, called ascending, descending and circulatory. It is verified the location of ascending and descending plumes that are affected by the rotation of the cylinders. The authors find that the average Nusselt numbers on the side walls of the immersed body are higher than on the upper and lower surfaces. These extreme points depend on both the number of Richardson (Ri), and the number of Rayleigh (Ra), as well as, the specific rotation. With respect to thermal convection, the rotation shows a smaller effect compared to the effects of Rayleigh and Richardson number. In the analysis for $Re={10}^6$, the flow patterns became unstable and oscillating, resulting from the different rotation rates and the number of Richardson. It was also verified the existence of a mixed convection regime, where for Richardson numbers greater than 1, there is a decrease in the heat transfer rate in comparison with the natural convection regime for the relevant Rayleigh number. The results showed a good numerical approximation with others available in the literature. For the case of forced convection, considering an isothermal body with flow around a square cylinder, the work developed by Santos et al. [12] presents more details.

The contributions presented in this work, in contrast to some studies, such as Silva et al. [9], and Ashrafizadeh and Hosseinjani [11], are the two-dimensional simulations of flows, with the Reynolds number between 1 and 5000, associated with the heat transfer process by mixed convection. Heat transfer is carried out around bodies, immersed in an incompressible Newtonian fluid, with constant temperature on its surfaces. These contributions involved the development of a computational code that brings together both physical phenomena that, in different works, including those mentioned above, are done in a segregated way. The methodology, which we propose here, performs the calculation of the force field over a sequence of Lagrangian points, which represent the interface, using the Navier-Stokes and heat equations. Two turbulence models were implemented, namely, the Smagorinsky sub-mesh model and the Spalart-Allmaras model with heat transfer, both in the immersed boundary methodology coupled with the virtual physical model. The term “virtual” refers to the non-slip conditions that must be modelled without the direct imposition of velocity on the interface (for more details, see [12] ). Flow simulations were carried out around an isothermal square cylinder and around tandem and angulated isothermal circular cylinders. It was found that the simulations presented numerical convergence compared with published results.

This paper is organized as follows. In Section 2, we introduce the mathematical methodology. The Section 3 is devoted to the URANS (unsteady Reynolds-averaged Navier-Stokes), Smagorinsky and Spalart-Allmaras turbulence models. Section 4 details the numerical schemes necessary for the spatial discretization of the various differential operators, and the numerical algorithms for the temporal evolution of the different fields involved. In Section 5, we report the results, followed by the conclusions.

2. Mathematical Methodology

2.1. Mathematical Formulation for the Fluid Motion and Temperature

Consider an incompressible and two-dimensional flow for a Newtonian fluid with a domain $\Omega$ and border $\partial {\Omega }_b$, with the surface of the immersed body being heated with isothermal temperature (or constant temperature on the surface of the immersed body), modelled through discretized (Lagrangian) points. Since the effect of the frontier is taken into account through the introduction of a forcing term in the momentum and heat equations, the description of heat transfer by mixed convection in the immersed boundary methodology is written by (Pope [13] ):

$\nabla \cdot u=0\mathrm{,}$ (1)

${\rho }_0\left(\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u\right)=-\nabla p+\mu {\nabla }^{2}u+f\mathrm{,}$ (2)

${\rho }_0\left(\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u\right)=-\nabla p+\mu {\nabla }^{2}u+{\rho }_0\left(1-\beta \left(T-T_{\infty }\right)\right)g+f\mathrm{,}$ (3)

${\rho }_0{c}_p\left(\frac{\partial T}{\partial t}+\left(u\cdot \nabla \right)T\right)=k{\nabla }^{2}T+q\mathrm{,}$ (4)

where Equations (1), (2) and (4) are used in forced convection, while Equations (1), (3) and (4) are for natural convection; in Equation (3), the Boussinesq approximation is used. The fields $u$, p, T and the scalar $T_{\infty }$ denote the incompressible velocity vector, the pressure, the temperature, and the reference temperature, respectively. In (2)-(4), the quantities ${\rho }_0$, $\mu$, $\beta$, k and $c_p$ are, respectively, the fluid density at temperature $T=T_{\infty }$, the viscosity, the thermal diffusivity, the thermal expansion coefficient, and the specific heat at constant pressure; $g$ is the downward gravitational acceleration; the term ${\rho }_0\beta \left(T-T_{\infty }\right)g$ accounts for the effects of the fluid temperature on the fluid flow.

The force $f$ appearing in the both Equations (2) and (3) is the Euler force field, being intended to model the interface immersed in the flow

$f\left(x\mathrm{,}t\right)=\underset{\partial {\Omega }_b}{{\displaystyle \int }}F\left(X_k\mathrm{,}t\right)\delta \left(x-X_k\right)\text{d}{X}_k\mathrm{,}$

where $F\left(X_k\mathrm{,}t\right)$ is the Lagrangian force density, calculated on the interface points $x$ and $X_k$ which are the position vectors of the Eulerian and Lagrangian points, respectively. The presence of the Dirac delta function, $\delta \left(x-X_k\right)$, represents the interaction between the fluid and the immersed boundary. Similarly, the thermal source, represented by q, is added to Equation (4), being responsible for making the flow feel the presence of the heated solid interface; in other words, it is the heating source at the Lagrangian point on the immersed border, being expressed by

$q\left(x\mathrm{,}t\right)=\underset{\partial {\Omega }_b}{{\displaystyle \int }}Q\left(X_k\mathrm{,}t\right)\delta \left(x-X_k\right)\text{d}{X}_k\mathrm{.}$

Here, $Q\left(X_k\mathrm{,}t\right)$ is the heat flux on the Lagrangean point $X_k$.

2.2. Dimensionless Equations for Mixed Convection

Equations (1)-(4) in dimensionaless form assume the writing (Pope [13] )

$\nabla \cdot U=0\mathrm{,}$ (5)

$D_{t}U=-\nabla p+\left(\frac{1}{Re}\right){\nabla }^{2}U+\underset{\partial {\Omega }_b}{{\displaystyle \int }}F\left(X_k\mathrm{,}t\right)\delta \left(x-X_k\right)\text{d}{X}_k\mathrm{,}$ (6)

$D_{t}U=-\nabla p+\left(\frac{1}{Re}\right){\nabla }^{2}U+\left(\frac{Gr}{R{e}^2}\right)\Theta j+\underset{\partial {\Omega }_b}{{\displaystyle \int }}F\left(X_k\mathrm{,}t\right)\delta \left(x-X_k\right)\text{d}{X}_k\mathrm{,}$ (7)

$D_t\Theta =\left(\frac{1}{Pe}\right){\nabla }^2\Theta +\underset{\partial {\Omega }_b}{{\displaystyle \int }}Q\left(X_k\mathrm{,}t\right)\delta \left(x-X_k\right)\text{d}{X}_k\mathrm{,}$ (8)

where $D_t•={\partial }_t•+\left(U\cdot \nabla \right)•$ is the material derivative, $U$, p and $\Theta$ are the dimensionless velocity, pressure and temperature fields, respectively, and $Re=\left(u_{\infty }D\right)/\nu$ is the Reynolds number ( $u_{\infty }$ is a reference flow velocity and D is a characteristic linear dimension of the fluid).

In this work, we will also refer to the number of Richardson

$Ri=\frac{g\beta \left(T_c-T_{\infty }\right)D}{u_{\infty }^2}\mathrm{.}$

which expresses the ratio between natural convection and forced convection (g is the gravitational acceleration, $\beta$ is the thermal expansion coefficient, $T_c$ is the boundary/wall temperature, and $T_{\infty }$ is a reference temperature).

2.3. The Indicator Function

The indicator function, represented by $\mathcal{I}\left(x\mathrm{,}t\right)$, aims to signal and identify the position of the immersed body. This function was proposed by Unverdi and Tryggvason [14], being an interface monitoring method, where the function is calculated in the whole domain or part of it. In this work, the indicator function is defined by

${\nabla }^2\mathcal{I}\left(x\mathrm{,}t\right)=\nabla \cdot G\left(x\mathrm{,}t\right)\mathrm{,}$ (9)

being $G\left(x\mathrm{,}t\right)$ given by

$G\left(x\mathrm{,}t\right)={\displaystyle \underset{k=1}{\overset{N}{\sum }}}{D}_{ij}\left(x-X_k\right)n\left(X_k\right)\Delta s\left(X_k\right)\mathrm{,}$

where N is the number of Lagrangian points, $D_{ij}\left(x-X_k\right)$ is an interpolation/distribution function (to be defined for convenience), $n\left(X_k\right)$ is the vector normal to the interface at the Lagrangian point $X_k$, and $\Delta s\left(X_k\right)$ is the distance between two consecutive Lagrangian points. After solving the Poisson Equation (9), we obtained the indicator function, $\mathcal{I}\left(x\mathrm{,}t\right)$, for the whole calculation domain. To solve the linear system obtained by the discretization of Equation (9), we used the modified strongly implicit procedure (MSIP), developed by Schneider and Zedan [15].

2.4. The Virtual Physical Model

The virtual physical model, developed by Silva et al. [9], is a methodology for calculating the forces that act on the discrete points of a given border, also called the interfacial force or Lagrangian force. The characterization of the Lagrangian force represents the difference between the various immersed boundary methodologies. In this work, only rigid boundaries were treated (no elasticity), but the model can be used, or extended, to other types of interface (for example, elastic boundaries, boundaries between different fluids, etc). The virtual physical model uses the diffusion of interfacial forces on the interior of the flow. Thus, the Eulerian force field is applied in the vicinity of the immersed boundary, and its value is minimized as the distance to the interface increases. This model dynamically assesses not only the force that the fluid exerts on the solid surface immersed in the flow, but takes into account the thermal exchange between them.

The Lagrangian force, $F\left(X_k\mathrm{,}t\right)$, and the thermal source, $Q\left(X_k\mathrm{,}t\right)$, are evaluated separately. In other words, for the Lagrangian force, a balance of momentum was carried out on a fluid particle that is close to the fluid-solid interface, while, for the thermal part, the dimensionless heat Equation (8) was applied, which shows the interaction between the particle-fluid and the interface, as shown in Figure 1. Thus, the density of the interfacial force can be evaluated using the principle of conservation of momentum and energy. Therefore, taking the particle-fluid, as illustrated in Figure 1, crossing an arbitrary immersed interface, we obtain

$\begin{array}{l}F\left(X_k\mathrm{,}t\right)=\underset{F_a}{\underset{︸}{\rho \frac{\partial U\left(X_k\mathrm{,}t\right)}{\partial t}}}+\underset{F_i}{\underset{︸}{\rho \nabla \left[U\left(X_k\mathrm{,}t\right)U\left(X_k\mathrm{,}t\right)\right]}}\\ +\underset{F_p}{\underset{︸}{\nabla p\left(X_k\mathrm{,}t\right)}}-\underset{F_v}{\underset{︸}{\mu {\nabla }^{2}U\left(X_k\mathrm{,}t\right)}}\mathrm{.}\end{array}$ (10)

Here, from left to right, we have the acceleration force $F_a$, the inertial force

Figure 1. Control volume illustration located in a Lagrangian fluid particle.

$F_i$, the pressure force $F_p$, and the viscous force $F_v$, respectively. In Equation (10), the terms $F_a$, $F_i$, $F_p$ and $F_v$ are analysed at the Lagrangian mesh interface itself, considering the pressure and velocity fields, where, simultaneously, the calculation of these forces is performed in Eulerian mesh, through Lagrange interpolation.

Similarly, the calculation of the thermal source in the particle-fluid in contact with the interface, an energy balance is performed as follows

$Q\left(X_k\mathrm{,}t\right)=\frac{\partial \Theta \left(X_k\mathrm{,}t\right)}{\partial t}+\nabla \left[U\left(X_k\mathrm{,}t\right)\Theta \left(X_k\mathrm{,}t\right)\right]-\left(\frac{1}{Pe}\right){\nabla }^2\Theta \left(X_k\mathrm{,}t\right)\mathrm{,}$ (11)

where, from left to right, we have the local temperature variation rate, the thermal dissipation rate due to convection, and the diffusive thermal energy transport rate.

2.5. Calculation of Velocity, Pressure and Temperature

2.5.1. Auxiliary Point Allocation Process

The first step is to arbitrate an initial Lagrangian point for calculating the interfacial force, $F\left(X_k\mathrm{,}t\right)$. Then, two mutually orthogonal auxiliary lines are drawn on this point, one of which is parallel to one of the Eulerian axes. Two auxiliary points are marked on each of the lines, on the outside of the solid body, at a distance $\Delta x$ and $2\Delta x$ of the Lagrangian point considered. This distance is necessary in order to prevent two auxiliary points from being allocated within the same Eulerian cell. The meshes that are more than $2\Delta x$ distance away from the Lagrangian point, do not contribute to the interpolation. The internal and external regions of the solid body were identified with the aid of the normal unitary vector on the surface, which has its positive direction facing outside the immersed body. The auxiliary points are always located in the region of interest of the flow. Thus, the velocity, the pressure and the temperature at these points, in general, are not known, but they can be obtained from the neighbouring cells by interpolation.

Therefore, the velocity at the Lagrangian and auxiliary points is written

$U\left(X_k\right)={\displaystyle \underset{i}{\sum }}D\left(x_i-X_k\right)U\left(x_i\right)\mathrm{,}$ (12)

where $U\left(X_k\right)$ are the Lagrangian velocities, calculated at the auxiliary points and at the point $X_k$ by the interpolation of the Eulerian velocities.

The general equations for obtaining the pressure and temperature at the auxiliary points, or at the interface, and at the Lagrangian point, in x and y directions, are given respectively by

$\{\begin{array}{l}p\left(X_k\right)={\displaystyle \underset{i}{\sum }}D\left(x_i-X_k\right)p\left(x_i\right)\\ \Theta \left(X_k\right)={\displaystyle \underset{i}{\sum }}D\left(x_i-X_k\right)\Theta (\; x\; i\; )\end{array}$

and

$\{\begin{array}{l}p\left(Y_k\right)={\displaystyle \underset{i}{\sum }}D\left(y_i-Y_k\right)p\left(y_i\right)\\ \Theta \left(Y_k\right)={\displaystyle \underset{i}{\sum }}D\left(x_i-Y_k\right)\Theta (\; y\; i\; )\end{array}$

where $p\left(X_k\right)$ and $p\left(Y_k\right)$ are the pressure values at the interface, or at auxiliary points, and $p\left(x_i\right)$ and $p\left(y_i\right)$ are the pressure values at the nearest Eulerian mesh points, in x and y directions, respectively. Similarly, $\Theta \left(X_k\right)$ are the temperatures at the auxiliary points, and $\Theta \left(x_i\right)$ are the temperatures at the nearest Eulerian points. Since this work deals with two-dimensional flows, the transverse direction was disregarded.

The distribution/interpolation function, $D(\cdot )$, adopted in this work, is used for the interpolation of variables in the Eulerian mesh. Equation (13) shows the specific mass calculation process in the $X_k$ position

$\rho \left(X_k\mathrm{,}t\right)={\displaystyle \underset{i=1}{\overset{n_x}{\sum }}}{\displaystyle \underset{j=1}{\overset{n_y}{\sum }}}D\left(x_{i\mathrm{,}j}-X_k\right)\rho \left(x_{i\mathrm{,}j}\mathrm{,}t\right)\text{d}x\text{d}y\mathrm{,}$ (13)

where $n_x$ and $n_y$ represent the totality of cells in the x and y directions, respectively. Concerning the computational cost involved, it was reduced when considering non-zero $D(\cdot )$ for distances less than $2\Delta x$ from the interpolation point, which is also valid for the $F\left(X_k\mathrm{,}t\right)$ distribution. Therefore, this work has considered (Peskin [16] [17], Juric [18] ):

$D\left(x-X_k\right)={\displaystyle \underset{m=1}{\overset{N}{\prod }}}\frac{g\left(r_x\right)g\left(r_y\right)}{h^2}\mathrm{,}$

with

$g\left(r\right)=\{\begin{array}{ll}{g}_1\left(r\right)\hfill & 0\le r<1\hfill \\ 0.5-g_1\left(2-r\right)\hfill & 1\le r<2\hfill \\ 0\hfill & r\ge 2\hfill \end{array}$

where $g_1\text{}\left(r\right)=\frac{1}{8}\left(3-2r\sqrt{1+4r-4r^2}\right)$, being r the radius of influence of the

distribution function, given here by $h^{-1}\left(x-x_k\right)$ or $h^{-1}\left(y-y_k\right)$. The quantity $h=\Delta x=\Delta y$, is the size of the Eulerian mesh and $\left(x\mathrm{,}y\right)$ the coordinates of an Eulerian point in the domain. According to the definition of $g(\cdot )$, the interpolation covers a distance of, at most, two meshes from the point considered, since the value of the function $D(\cdot )$ is cancelled for greater distances. This can be explained because the region of coverage of the function $D(\cdot )$ must be finite in order to keep the computational cost reasonable. Without this observation, a bunch of points would contribute, causing a high computational cost. This distribution function is of Gaussian type. It is important to note that $D(\cdot )$ also allows the transfer of quantities between the Eulerian and Lagrangian meshes. This function acts as an interpolation function, distributing the Lagrangian force outside the interface. A convenient procedure to tune the interpolation region is to allocate a square centred on each auxiliary point (green dotted box), represented by dot lines in Figure 2, to obtain the velocity, pressure and temperature.

Thus, only the points inside this “dotted box” in Figure 2 are important in the calculation process, thereby reducing the computational cost, because for Eulerian points very distant from the analysed point $X_k$, the function $D(\cdot )$ takes the value zero. The use of internal points of the interface during the interpolation procedure (see Figure 3) is physically coherent, since the internal flow is also solved by the Navier-Stokes equations.

Figure 2. Illustration of the velocity interpolation procedure at auxiliary point 3: (a) in x-direction, and (b) in y-direction.

Figure 3. Illustration of auxiliary points in the interpolation scheme for the calculation of Lagrangian forces.

After the interpolation of velocity, pressure and temperature at the interface at auxiliary points, the derivatives at the Lagrangian source terms, in (10) and (11), are determined using Lagrange polynomials of first and second order.

2.5.2. Calculation of Lagrangian Force Distribution and Thermal Source

After calculating the terms of Equation (10) and Equation (11), and obtaining the values for $F\left(X_k\mathrm{,}t\right)$ and $Q\left(X_k\mathrm{,}t\right)$, then the Eulerian terms are calculated for $f$ and q. Figure 4 illustrates the distribution of interfacial force and thermal field.

The calculation of $f$ and q is done by

$\{\begin{array}{l}f\left(x_i\right)={\displaystyle \underset{i}{\sum }}D\left(x_i-X_k\right)F\left(X_k\right)\Delta s\\ q\left(x_i\right)\mathrm{<mo>\; =\; </mo>}{\displaystyle \underset{i}{\sum }}D\left(x_i-X_k\right)Q\left(X_k\right)\Delta s\end{array}$ (14)

$\{\begin{array}{l}f\left(y_i\right)={\displaystyle \underset{i}{\sum }}D\left(y_i-Y_k\right)F\left(Y_k\right)\Delta s\\ q\left(y_i\right)={\displaystyle \underset{i}{\sum }}D\left(y_i-Y_k\right)Q\left(Y_k\right)\Delta s\end{array}$ (15)

where $f\left(x_i\right)$ and $f\left(y_i\right)$ are the forces at each Eulerian node, while $F\left(X_k\right)$ and $Q\left(X_k\right)$ are the force in each Lagrangian node. In (14)-(15), (i) $q\left(x_i\right)$ and $f\left(y_i\right)$ are the thermal sources for each Eulerian node, and (ii) $Q\left(X_k\right)$ and $Q\left(Y_k\right)$ are the thermal source for each Eulerian node. Thus, there is a thermal field of Eulerian forces acting on the Lagrangian grid, both acting together with the fluid particles in the vicinity of the border. The remaining terms are the average distances between two adjacent Lagrangian points, represented by $\Delta s$, and D is the distribution function.

In this work, the boundary conditions, for the Navier-Stokes equation, are Dirichlet and Neumann conditions (Dirichlet condition defines the values at the

Figure 4. Distribution of the interfacial force, $f\left(x_i\right)$ and the thermal force, $q\left(x_i\right)$ at each Lagrangian node in the body immersed in the flow, resulting in a heated interfacial force field, allowing the particle-fluid to recognize the heated body.

border; Neumann condition fixes the values for the derivatives at the border). Dirichlet and Neumann conditions are used to specify the pressure, temperature and components of the velocity vector at the boundaries of the calculation domain.

3. Mathematical Modelling of Turbulence

Among the methodologies available for numerical resolution of turbulence models, the most common is the direct numerical simulation (DNS), which provides full resolution of the equations (all scales of the flow are solved), requiring, in general, high temporal and spatial resolutions. As a consequence, extremely refined meshes must be used, so that all scales present in the flow can be properly calculated. For this reason, studies involving DNS, in general, are restricted to very simple generic configurations and to flows involving low or moderate Reynolds numbers.

3.1. The URANS Formulation for Turbulence

In the URANS method, the turbulent vortices are not directly resolved and the entire turbulence behaviour is then modelled. In fact, with the temporal average technique, only the relatively low frequency variations of the flow field are captured by the URANS approach, in contrast to the turbulent high-frequency fluctuations, being fully modelled by other methods. In this respect, the velocity field is decomposed into

$u_i\left(t\right)={\bar{u}}_i\left(t\right)+{u^{\prime }}_i\left(t\right)\mathrm{,}$

where ${\bar{u}}_i\left(t\right)$ represents the coherent part of the velocity in the turbulent flow, that involves all organized low-frequency movement, and ${u^{\prime }}_i\left(t\right)$ indicates the high frequency floating part of the turbulent velocity. This decomposition is also applied to the pressure. The URANS equations for incompressible flows, in Einstein notation for indices, are

$\frac{\partial {\bar{u}}_i}{\partial t}+\frac{\partial \left({\bar{u}}_i{\bar{u}}_j\right)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x_i}+\nu \frac{{\partial }^2{\bar{u}}_i}{\partial x_j\partial x_j}-\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}\mathrm{,}$ (16)

where ${\tau }_{ij}=\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}$ is the Reynolds stress tensor (Pope [13] ). This tensor acts as a tension on the fluid element and, physically, is interpreted as a correlation of velocity fluctuation. To make Equation (16) “solvable”, ${\tau }_{ij}$ must be expressed based on coherent flow variables. The turbulent viscosity models relate the Reynolds stress tensor with a new variable, defined as turbulent viscosity, ${\nu }_t$ (Lesieur [19] ).

3.2. The Smagorinsky Turbulence Model

The Smagorinsky model [13] [19] [20] is based on the local equilibrium hypothesis for small scales, so that the injected energy in the spectrum equals the dissipated energy by the viscous effects. This model offers

${\nu }_t={\left(C_s\mathcal{l}\right)}^2\sqrt{2{\bar{S}}_{ij}{\bar{S}}_{ij}}\mathrm{.}$ (17)

Here, ${\bar{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i}\right)$ is the strain rate tensor, $\mathcal{l}$ is the length of the

sub-mesh, and $C_s$ is Smagorinsky constant.

3.3. The Spalart-Allmaras Turbulence Model

According to this model, the turbulent viscosity, ${\nu }_t$, is calculated thanks to the auxiliary work variable $\tilde{\nu }$ [13] [21]. Then

${\nu }_t=\tilde{\nu }{f}_{v_1}\mathrm{,}$

where $f_{v_1}={\chi }^3/\left({\chi }^3+C_{v_1}^3\right)$, with $\chi =\tilde{\nu }/\nu$, and $\tilde{\nu }$ obeying the transport equation:

$\begin{array}{c}\frac{\partial \tilde{\nu }}{\partial t}+\frac{\partial }{\partial x_j}\left(u_j\tilde{\nu }\right)=c_{b_1}\left(1-f_{t_2}\right)\tilde{S}\tilde{\nu }+\frac{1}{\sigma }\left[\frac{\partial }{\partial x_j}\left(\left(\nu +\tilde{\nu }\right)\frac{\partial \tilde{\nu }}{\partial x_j}\right)+c_{b_2}\frac{\partial \tilde{\nu }}{\partial x_j}\frac{\partial \tilde{\nu }}{\partial x_j}\right]\\ -\left(c_w{f}_w-\frac{c_{b_1}}{k^2}{f}_{t_2}\right){\left(\frac{\tilde{\nu }}{d_w}\right)}^2+f_{t_1}\Delta U^2.\end{array}$ (18)

Here, the installments on the right-hand side of Equation (18) represent the production of turbulent viscosity, the turbulent and molecular diffusions of turbulent viscosity, the dissipation of turbulent viscosity, the destruction of turbulent viscosity, responsible for its decrease, and terms that model transition effects to turbulence. The various coefficients (and corresponding numerical values) that appear in this turbulence model can be found in the 2015 Kostic review article [22].

4. Numerical Method

Exact solutions of differential equations are only known in restricted situations, since, in general, these equations involve nonlinear terms. This implies the need to use numerical algorithms and demanding computational facilities to perform the calculations. In our case, the spatial (computational) domain was discretized, that is, divided into control volumes and, for each volume, a nodal point was associated. The set of theses discrete points (nodal points) is called the mesh, where the solution was obtained for each of these nodal points. The conservation equations were discretized by the finite difference method. Among the various existing numerical methods, we opted for the fractional-step method (FSM), developed by Chorin [23], and later improved by Kim and Moin (1985) [24], for the temporal discretization.

The numerical solution of Equations (5)-(7) was performed in a segregated manner, requiring the treatment of the coupling between pressure and velocity. The method FSM is used in a non-iterative way, initially solving the momentum equations, using the velocity fields of the previous time, which leads to an estimate for the velocity field, which, in general, does not satisfy the continuity equation. The difference of this estimate is used as a source term for the solution of a Poisson equation for the correction of pressure, being used to update the velocity field and, finally, guaranteeing the conservation of mass to the corrected velocity. The importance of the Poisson equation for the correction of pressure is due to the fact that it makes the connection between the equations of momentum and continuity.

The Navier-Stokes equation in forced convection is discretized, and for natural convection the process is analogous. Thus, Equation (6) can be rewritten in the following way

$\frac{u_i^{n+1}-u_i^n}{\Delta t}+\left[\frac{\partial }{\partial x_j}\left(u_i^n{u}_j^n\right)\right]=-\frac{1}{\rho }\frac{\partial p^{n+1}}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\nu +{\nu }_t\right)\left(\frac{\partial u_i^n}{\partial x_j}+\frac{\partial u_j^n}{\partial x_i}\right)\right]+f_i^n\mathrm{.}$ (19)

In the FSM method, the velocity, the pressure and the forcing term of the predictive instant, n, are used to calculate, in the predictive step, an estimate for the velocity in the current time ${\bar{u}}_i^{n+1}$, offered by the equation

$\frac{{\bar{u}}_i^{n+1}-u_i^n}{\Delta t}+\left[\frac{\partial }{\partial x_j}\left(u_i^n{u}_j^n\right)\right]=-\frac{1}{\rho }\frac{\partial p^n}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\nu +{\nu }_t\right)\left(\frac{\partial u_i^n}{\partial x_j}+\frac{\partial u_j^n}{\partial x_i}\right)\right]+f_i^n\mathrm{.}$ (20)

The next step, in the FSM method, is to subtract Equation (20) from Equation (19), resulting in

$\frac{{\bar{u}}_i^{n+1}-u_i^n}{\Delta t}=\frac{1}{\rho }\frac{\partial {p^{\prime }}^{n+1}}{\partial x_i}\mathrm{,}$ (21)

where

${p^{\prime }}^{n+1}=p^{n+1}-p^n\mathrm{.}$ (22)

Knowing that the velocity field must satisfy the continuity equation, one obtains

$\frac{1}{\Delta t}\frac{\partial {\bar{u}}_i^{n+1}}{\partial x_i}=\frac{1}{\rho }\frac{{\partial }^2{p^{\prime }}^{n+1}}{\partial x_i\partial x_i}\mathrm{.}$ (23)

The Equation (23) is a Poisson equation for the correction of pressure, $p^{\prime }$. Note that the source term is a function of the estimated velocity, which is known from the predictor step, Equation (21). Finally, with the pressure field calculated, we return to Equation (21), obtaining, then, the equation corrected for the velocity

$u_i^{n+1}={\bar{u}}_i^{n+1}-\frac{\Delta t}{\rho }\frac{\partial {p^{\prime }}^{n+1}}{\partial x_i}\mathrm{.}$ (24)

The reader can find more details about this numerical scheme in [8] [12].

4.1. The Fractional Method Algorithm

This algorithm consists in carrying out the following steps:

1) Estimate the velocity field, Equation (19);

2) With the estimated velocity field, solve the linear system for the pressure correction, Equation (23);

3) Correct the velocity field, Equation (24), and pressure, Equation (22);

4) Check the conservation of mass within the specified tolerance;

5) Advance to the next time step.

4.2. Temporal Discretization

The temporal discretization techniques used in the present work are presented in this section.

4.2.1. Second-Order Runge-Kutta Method

In two-dimensional simulations around a single immersed isothermal body, with or without border movement, the second-order Runge-Kutta method was used for the temporal discretization of the Navier-Stokes equations. The method consists of a predictor step that evaluates the field in an intermediate time, called the corrector, and then the update to the current instant is performed. In other words,

· Predictor step:

$\{\begin{array}{l}{\tilde{u}}_{I,J}^{n+\frac{1}{2}}=u_{i,j}^n+\frac{\Delta t}{2}\left[-A_{xi,j}^n+D_{xi,j}^n-\frac{1}{\rho }{P}_{xi,j}^n+f_{xi,j}^n\right]\\ {\tilde{v}}_{I,J}^{n+\frac{1}{2}}=v_{i,j}^n+\frac{\Delta t}{2}\left[-A_{yi,j}^n+D_{yi,j}^n-\frac{1}{\rho }{P}_{yi,j}^n+f_{yi,j}^n\right]\end{array}$

· Corrector step:

$\{\begin{array}{l}{\tilde{u}}_{I,J}^{n+\frac{1}{2}}=u_{i,j}^n+\Delta t\left[-A_{xi,j}^{n+\frac{1}{2}}+D_{xi,j}^{n+\frac{1}{2}}-\frac{1}{\rho }{P}_{xi,j}^{n+\frac{1}{2}}+f_{xi,j}^{n+\frac{1}{2}}\right]\\ {\tilde{v}}_{I,J}^{n+\frac{1}{2}}=v_{i,j}^n+\Delta t\left[-A_{yi,j}^{n+\frac{1}{2}}+D_{yi,j}^{n+\frac{1}{2}}-\frac{1}{\rho }{P}_{yi,j}^{n+\frac{1}{2}}+f_{yi,j}^{n+\frac{1}{2}}\right]\end{array}$

where A represents the advective source, D is for the diffusive term, P the pressure correction gradient, and f the force term.

4.2.2. Second-Order Adams-Bashforth Method

This method was used in the simulations for the cases of a bluff isothermal body (square cylinder) and two tandem cylinders (in a row, with or without frontier movement), for the time advance of the momentum equations. For the first iteration, the first order Euler method was used:

$\{\begin{array}{l}{\tilde{u}}_{I,J}^n=u_{i,j}^n+\Delta t\left[1.5\left(-A_{xi,j}^n+D_{xi,j}^n\right)-0.5\left(-A_{xi,j}^{n-1}+D_{xi,j}^{n-1}\right)\right]-\Delta t\left(\frac{1}{\rho }{P}_{xi,j}^n+f_{xi,j}^n\right)\\ {\tilde{v}}_{I,J}^n=v_{i,j}^n+\Delta t\left[1.5\left(-A_{yi,j}^n+D_{yi,j}^n\right)-0.5\left(-A_{yi,j}^{n-1}+D_{yi,j}^{n-1}\right)\right]-\Delta t\left(\frac{1}{\rho }{P}_{yi,j}^n+f_{yi,j}^n\right)\end{array}$

4.3. The Spatial Discretization of Equations

4.3.1. Discretization of Navier-Stokes Equations

The discretization of the Navier-Stokes and heat equations used finite differences centred of second order for the spatial discretization, and the explicit Euler method of first order for the temporal discretization. The velocities are located on the control volume. In Figure 5, a control volume of the Eulerian point $\left(i\mathrm{,}j\right)$ is illustrated. This scheme was adopted to provide more stability in the pressure-velocity coupling. Therefore, considering a cell with coordinate $\left(i\mathrm{,}j\right)$, we have the distribution: (i) the pressure ( $p_{i\mathrm{,}j}$ ), and temperature ( ${\theta }_{i\mathrm{,}j}$ ) in the centre of the Eulerian cell; (ii) thex components of Eulerian velocity, and force on the left side face ( $\Delta x/2$ away from the centre of the cell), and (iii) the y components Eulerian velocity, and force on the left bottom face ( $\Delta y/2$ away from the centre of the cell). All simulations have been carried in non-uniform meshes. For non-uniform mesh (Figure 6), we have:

$\frac{\partial p^{\prime }}{\partial x}\approx \frac{{p^{\prime }}_{i\mathrm{,}j}-{p^{\prime }}_{i-\mathrm{1,}j}}{\Delta x_i}\mathrm{,}$

and

$\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)\approx \frac{1}{\Delta x_j}\left[\left({\bar{u}}_p{\bar{u}}_p\right)-\left({\bar{u}}_w{\bar{u}}_w\right)\right]+\frac{1}{\Delta y_i}\left[\left({\bar{u}}_{nw}{\bar{v}}_{nw}\right)-\left({\bar{u}}_{sw}{\bar{v}}_{sw}\right)\right]\mathrm{.}$

The velocities defined by the average operator must be interpolated, being functions of the velocity on the faces of the mesh (using the mesh displacement scheme) given by

Figure 5. Control volume of an Eulerian point.

Figure 6. Illustrative scheme of two non-uniform Cartesian cells.

$\begin{array}{lll}{\bar{u}}_p=f\left(u_{i,j+1},u_{i,j}\right),\hfill & {\bar{u}}_w=f\left(u_{i,j},u_{i,j-1}\right),\hfill & {\bar{u}}_{nw}=f\left(u_{i,j},u_{i+1,j}\right),\hfill \\ {\bar{u}}_{sw}=f\left(u_{i,j},u_{i-1,j}\right),\hfill & {\bar{v}}_w=f\left(v_{i+1},v_{i+1,j-1}\right),\hfill & {\bar{v}}_{sw}=f\left(v_{i,j},v_{i,j-1}\right).\hfill \end{array}$

The approximation (25) discretizes the diffusive term,i.e.

$\begin{array}{c}\frac{\partial }{\partial x_j}\left[{\nu }_{ef}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]\approx \frac{1}{\Delta x_j}\left(2{\nu }_{e{f}_{i\mathrm{,}j}}\frac{u_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}}{\Delta x_j}-2\nu f_{i\mathrm{,}j-1}\frac{u_{i\mathrm{,}j}-u_{i\mathrm{,}j-1}}{\Delta x_j}\right)\\ +\frac{1}{\Delta y_i}[{\tilde{\nu }}_{e{f}_N}\left(\frac{u_{i+1}\mathrm{,}{u}_{i\mathrm{,}j}}{\Delta y_{i+1}}+\frac{v_{i+1}-v_{i+\mathrm{1,}j-1}}{\Delta x_j}\right)\\ -{\tilde{\nu }}_{e{f}_S}\left(\frac{u_{i\mathrm{,}j}-u_{i-\mathrm{1,}j}}{\Delta y_i}+\frac{v_{i\mathrm{,}j}-v_{i\mathrm{,}j}-v_{i\mathrm{,}j-1}}{\Delta x_j}\right)]\mathrm{.}\end{array}$ (25)

The effective viscosity, which must be interpolated on the faces, is a function of the viscosity in the neighbouring cells

${\bar{\nu }}_{e{f}_N}=f\left({\nu }_{e{f}_{i\mathrm{,}j}}\mathrm{,}{\nu }_{e{f}_{i\mathrm{.}j-1}}\mathrm{,}{\nu }_{e{f}_{I+\mathrm{1,}j}}\mathrm{,}{\nu }_{e{f}_{I+\mathrm{1,}j-1}}\right)\mathrm{,}$

${\bar{\nu }}_{e{f}_S}=f\left({\nu }_{e{f}_{i\mathrm{,}j}}\mathrm{,}{\nu }_{e{f}_{i\mathrm{.}j-1}}\mathrm{,}{\nu }_{e{f}_{I-\mathrm{1,}j}}\mathrm{,}{\nu }_{e{f}_{I-\mathrm{1,}j-1}}\right)\mathrm{.}$

4.3.2. Discretization of the Pressure Correction Equation

The Equation (23) for the correction of pressure is written as follows

$\frac{{\partial }^2{p}^{\prime }}{\partial x^2}+\frac{{\partial }^2{p}^{\prime }}{\partial y^2}=\frac{\rho }{\Delta t}{\left[\frac{\partial \bar{u}}{\partial x}+\frac{\partial \bar{v}}{\partial y}\right]}^{n+1}\mathrm{.}$

Using the finite difference centred scheme, it reads

$\begin{array}{l}\frac{1}{\Delta x}\left[\frac{{p^{\prime }}_{i\mathrm{,}j+1}^{n+1}-{p^{\prime }}_{i\mathrm{,}j}^{n+1}}{\Delta x}-\frac{{p^{\prime }}_{i\mathrm{,}j}^{n+1}-{p^{\prime }}_{i\mathrm{,}j-1}^{n+1}}{\Delta x}\right]+\frac{1}{\Delta y}\left[\frac{{p^{\prime }}_{i+\mathrm{1,}j}^{n+1}-{p^{\prime }}_{i\mathrm{,}j}^{n+1}}{\Delta y}-\frac{{p^{\prime }}_{i\mathrm{,}j}^{n+1}-{p^{\prime }}_{i-\mathrm{1,}j}^{n+1}}{\Delta y}\right]\\ =\frac{\rho }{\Delta t}\left[\frac{1}{\Delta x}\left(\frac{{\bar{u}}_{i\mathrm{,}j+1}^{n+1}+{\bar{u}}_{i\mathrm{,}j}^{n+1}}{2}-\frac{{\bar{u}}_{i\mathrm{,}j}^{n+1}+{\bar{u}}_{i\mathrm{,}j-1}^{n+1}}{2}\right)+\frac{1}{\Delta y}\left(\frac{{\bar{v}}_{i+\mathrm{1,}j}^{n+1}+{\bar{v}}_{i\mathrm{,}j}^{n+1}}{2}-\frac{{\bar{v}}_{i\mathrm{,}j}^{n+1}+{\bar{v}}_{i-\mathrm{1,}j}^{n+1}}{2}\right)\right]\mathrm{.}\end{array}$ (26)

The Equation (26) presents the terms for the pressure correction at iteration time $n+1$. This way, all equations are coupled, generating a linear system of equations, which can be represented in the compact form

$a_p{p^{\prime }}_p+a_e{p^{\prime }}_e+a_w{p^{\prime }}_w+a_n{p^{\prime }}_n+a_s{p^{\prime }}_s=b_p\mathrm{,}$ (27)

where ${p^{\prime }}_p={p^{\prime }}_{i,j}$, ${p^{\prime }}_e={p^{\prime }}_{i+1,j}$, ${p^{\prime }}_w={p^{\prime }}_{i-1,j}$, ${p^{\prime }}_n={p^{\prime }}_{i,j+1}$ and ${p^{\prime }}_s={p^{\prime }}_{i,j-1}$. The coefficients $a_p\mathrm{,}{a}_e\mathrm{,}{a}_w\mathrm{,}{a}_n$ and $a_s$ are constant, for the entire uniform mesh. For the non-uniform mesh, they change along the domain, depending on the dimensions of each cell. The sub-indices p, e, w, n and s are illustrated in Figure 7, and the coefficients take the values:

$\begin{array}{l}{a}_e=a_w=\frac{1}{\Delta x^2},\\ a_n=a_s=\frac{1}{\Delta y^2},\\ a_p=-2\left(a_e+a_n\right),\end{array}$

Figure 7. Collocation points p, e, w, n and s.

and

$b_p=\frac{\rho }{\Delta t}\left(\frac{{\bar{u}}_{i+\mathrm{1,}j}^{n+1}-{\bar{u}}_{i\mathrm{,}j}^{n+1}}{\Delta x}-\frac{{\bar{v}}_{i\mathrm{,}j+1}^{n+1}-{\bar{v}}_{i\mathrm{,}j}^{n+1}}{\Delta y}\right)\mathrm{.}$

To solve the linear system (27), it was used the modified strongly implicit procedure—MSIP—proposed by Schneider and Zedan [15].

4.3.3. Discretization of the Smagorinsky Model

The spatial discretization of the algebraic equation for the turbulent viscosity of the Smagorinsky sub-mesh model, Equation (17), is

$S\equiv \sqrt{2S_{ij}{S}_{ij}}\mathrm{,}$

where

$S=\sqrt{2\underset{S_{11}}{\underset{︸}{{\left(\frac{\partial u}{\partial x}\right)}^2}}+\underset{S_{12}}{\underset{︸}{{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2}}+2\underset{S_{22}}{\underset{︸}{{\left(\frac{\partial v}{\partial y}\right)}^2}}}$

with

$\{\begin{array}{l}{S}_{11}=\left(\frac{u_{i,j+1}-u_{i,j}}{\Delta x}\right)\\ S_{12}=\left(\frac{u_{mn}-u_{ms}}{\Delta y}+\frac{v_{me}-v_{mw}}{\Delta x}\right)\\ S_{22}=\left(\frac{v_{i+1}-v_{i,j}}{\Delta y}\right)\end{array}$

where

$\begin{array}{c}\{\begin{array}{l}{u}_{ms}=\frac{1}{4}\left(u_{i,j}+u_{i,j+1}+u_{i-1,j}+u_{i-1,j+1}\right)\\ v_{me}=\frac{1}{4}\left(v_{i+1,j-1}+v_{i+1,j}+v_{i,j}+v_{i+j+1}\right)\\ v_{mw}=\frac{1}{4}\left(v_{i+1,j-1}+v_{i+1,j}+v_{i,j-1}+v_{i,j}\right)\end{array}\end{array}$

is for uniform meshes. With regard to non-uniform meshes, the next section presents the interpolation process in a concise way.

4.3.4. Interpolation of properties for non-uniform meshes

With the use of non-uniform meshes, the quantities defined on the faces must be properly interpolated. The work of Patankar [25] proposes the use of an interpolation scheme. For the interpolation of velocities and density on the faces, a linear approximation between the points is used as follows

${\varphi }_e=f_e{\varphi }_P+\left(1-f_e\right){\varphi }_E\mathrm{,}$

where the interpolation factor, $f_e$, is the ratio between the distances of the non-uniform mesh, given by

$f_e=\frac{{\left(\delta x\right)}_e^{+}}{{\left(\delta x\right)}_e}\mathrm{,}$

and for the viscosity, it is used

${\varphi }_e={\left(\frac{1-f_e}{{\varphi }_p}+\frac{f_e}{{\varphi }_E}\right)}^{-1}\mathrm{.}$

4.3.5. Discretization of the Spalart-Allmaras Model

The spatial discretization of the transport equation for the turbulent viscosity of the Spalart-Allmaras model, Equation (18), is presented below. The same spatial discretization method was used to discretize the Navier-Stokes equation:

· Advective term:

$\frac{\partial }{\partial x_j}\left(\tilde{\nu }{u}_l\right)\approx \frac{1}{\Delta x_j}\left(u_{i\mathrm{,}j+1}{\overline{\stackrel{˜}{\nu }}}_e-u_{i\mathrm{,}j}{\overline{\stackrel{˜}{\nu }}}_w\right)+\frac{1}{\Delta y_i}\left(v_{i+\mathrm{1,}j}{\overline{\stackrel{˜}{\nu }}}_n-u_{i\mathrm{,}j}{\overline{\stackrel{˜}{\nu }}}_s\right)\mathrm{.}$

The auxiliary variable, $\tilde{\nu }$, is located at the centre of the mesh cell. Thus, to calculate the values on the faces, one must interpolate with values of the neighbourhood points, represented by

$\begin{array}{l}{\overline{\stackrel{˜}{\nu }}}_e=f\left({\tilde{\nu }}_{i,j+1},{\tilde{\nu }}_{i,j}\right),{\overline{\stackrel{˜}{\nu }}}_w=f\left({\tilde{\nu }}_{i,j},{\tilde{\nu }}_{i,j-1}\right),\\ {\overline{\stackrel{˜}{\nu }}}_n=f\left({\tilde{\nu }}_{i+1,j},{\tilde{\nu }}_{i,j}\right),{\overline{\stackrel{˜}{\nu }}}_s=f\left({\tilde{\nu }}_{i,j},{\tilde{\nu }}_{i-1,j}\right).\end{array}$

The turbulent viscosity production term of the Spalart-Allmaras model is proportional to the norm of the strain rate tensor, given by Equation (17) and is discretized as follows

$S_{ij}{S}_{ij}=S_{11}^2+2S_{12}^2+S_{22}^2\mathrm{,}$

where

$\{\begin{array}{l}{S}_{11}=\frac{\partial u}{\partial x}\approx \frac{u_{i,j+1}-u_{i,j}}{\Delta x_j}\\ S_{22}=\frac{\partial v}{\partial y}\approx \frac{v_{i+1,j}-v_{i,j}}{\Delta y_j}\\ S_{12}=\frac{1}{2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\approx \frac{1}{2}\left[\left(\frac{{\bar{u}}_N-{\bar{u}}_S}{\Delta y_i}\right)+\left(\frac{{\bar{v}}_E-{\bar{v}}_W}{\Delta x_j}\right)\right]\end{array}$

The interpolated velocities are functions of the following variables:

$\begin{array}{l}{\bar{u}}_N=f\left(u_{i,j},u_{i,j+1},u_{i+1,j},u_{i+1,j+1}\right),{\bar{u}}_S=f\left(u_{i,j},u_{i,j+1},u_{i-1,j},u_{i-1,j+1}\right)\\ {\bar{v}}_E=f\left(v_{i,j},v_{i,j+1},v_{i,j+1},v_{i+1,j+1}\right),{\bar{v}}_W=f\left(v_{i,j},v_{i+1,j},v_{i,j-1},v_{i+1,j-1}\right).\end{array}$

· Diffusive term:

The conservative diffusive term for the two-dimensional incompressible flow can be rewritten as follows

$\begin{array}{l}\left(\nu +\tilde{\nu }\right)\frac{{\partial }^2\tilde{\nu }}{\partial x_j\partial x_j}+{\left(\frac{\partial \tilde{\nu }}{\partial x_j}\right)}^2\\ \approx \left({\nu }_{i\mathrm{,}j}+{\tilde{\nu }}_{i\mathrm{,}j}\right)[\frac{{\tilde{\nu }}_{i\mathrm{,}j}-2{\tilde{\nu }}_{i\mathrm{,}j}+{\tilde{\nu }}_{i\mathrm{,}j-1}}{\Delta x_j^2}+\frac{{\tilde{v}}_{i+\mathrm{1,}j}-2{\tilde{\nu }}_{i\mathrm{,}j}+{\tilde{\nu }}_{i-\mathrm{1,}j}}{\Delta y_i^2}]\\ +\left(\frac{{\tilde{\nu }}_{i\mathrm{,}j+1}-{\tilde{\nu }}_{i\mathrm{,}j-1}}{\Delta x_j+\Delta x_{j+1}}\right)\left(\frac{{\tilde{\nu }}_{i\mathrm{,}j+1}-{\tilde{\nu }}_{i\mathrm{,}j-1}}{\Delta x_j+\Delta x_{j+1}}\right)+\left(\frac{{\tilde{\nu }}_{i+\mathrm{1,}j}-{\tilde{\nu }}_{i-\mathrm{1,}j}}{\Delta y_i+\Delta y_{i+1}}\right)\left(\frac{{\tilde{\nu }}_{i+\mathrm{1,}j}-{\tilde{\nu }}_{i-\mathrm{1,}j}}{\Delta y_i+\Delta y_{i+1}}\right)\mathrm{.}\end{array}$

· Diffusive term (non-conservative):

$\frac{\partial \tilde{\nu }}{\partial x_j}\frac{\partial \tilde{\nu }}{\partial x_j}\approx \left(\frac{{\tilde{\nu }}_{i\mathrm{,}j+1}-{\tilde{\nu }}_{i\mathrm{,}j}}{\Delta x_j}\right)\left(\frac{{\tilde{\nu }}_{i\mathrm{,}j+1}-{\tilde{\nu }}_{i\mathrm{,}j}}{\Delta x_j}\right)+\left(\frac{{\tilde{\nu }}_{i+\mathrm{1,}j}-{\tilde{\nu }}_{i\mathrm{,}j}}{\Delta y_j}\right)\left(\frac{{\tilde{\nu }}_{i+\mathrm{1,}j}-{\tilde{\nu }}_{i\mathrm{,}j}}{\Delta y_j}\right)\mathrm{.}$

4.3.6. Discretization of the Indicator Function

As for pressure, the indicator function and variables $G_x$ and $G_y$ located in the centre of the Eulerian mesh, Equation (9), for the two-dimensional case, can be rewritten as

$\frac{{\partial }^{2}I}{\partial x^2}+\frac{{\partial }^{2}I}{\partial y^2}=\frac{\partial G_x}{\partial x}+\frac{\partial G_y}{\partial y}\mathrm{.}$

For its discretization, the following approximations were used

$\frac{{\partial }^{2}I}{\partial x^2}\approx \frac{I_{i+1,j}-2I_{i,j}+I_{i-1,j}}{\Delta x_j^2}=\frac{I_e-2I_p+I_w}{\Delta x_j^2},$

$\frac{{\partial }^{2}I}{\partial y^2}\approx \frac{I_{i,j+1}-2I_{i,j}+I_{i-1,j}}{\Delta y_j^2}=\frac{I_n-2I_p+I_s}{\Delta y_j^2},$

$\frac{\partial G_x}{\partial x}\approx \frac{\left(G_{x_{i,j+1}}+G_{x_{i,j}}\right)-\left(G_{x_{i,j-1}}+G_{x_{i,j}}\right)}{2\Delta x_j},$

$\frac{\partial G_y}{\partial y}\approx \frac{\left(G_{y_{i\mathrm{,}j}}+G_{y_{i+\mathrm{1,}j}}\right)-\left(G_{y_{i-\mathrm{1,}j}}+G_{y_{i\mathrm{,}j}}\right)}{2\Delta x_j}\mathrm{.}$

4.3.7. Discretization of the Heat Equation

The heat equation, Equation (8), can be rewritten as

$\frac{\partial \theta }{\partial t}+\frac{\partial u\theta }{\partial x}+\frac{\partial v\theta }{\partial y}=\frac{1}{RePr}\left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}\right)+q\mathrm{.}$ (28)

Performing the discretization of Equation (28), results in

$\begin{array}{l}\frac{\partial \theta }{\partial t}\approx \frac{{\theta }_{i\mathrm{,}j}^{n+1}-{\theta }_{i\mathrm{,}j}^n}{\Delta t}\\ \frac{\partial u\theta }{\partial x}\approx \frac{1}{\Delta x}\left[u_{i+\mathrm{1,}j}\frac{{\theta }_{i+\mathrm{1,}j}+{\theta }_{i\mathrm{,}j}}{2}-u_{i\mathrm{,}j}\frac{{\theta }_{i\mathrm{,}j}+{\theta }_{i-\mathrm{1,}j}}{2}\right]\\ \frac{\partial v\theta }{\partial y}\approx \frac{1}{\Delta y}\left[v_{i+\mathrm{1,}j}\frac{{\theta }_{i+\mathrm{1,}j}+{\theta }_{i\mathrm{,}j}}{2}-v_{i\mathrm{,}j}\frac{{\theta }_{i\mathrm{,}j}+{\theta }_{i-\mathrm{1,}j}}{2}\right]\\ \frac{{\partial }^2\theta }{\partial x^2}\approx \frac{1}{\Delta x}\left[\frac{{\theta }_{i+\mathrm{1,}j}-{\theta }_{i\mathrm{,}j}}{\Delta x}-\frac{{\theta }_{i\mathrm{,}j}-{\theta }_{i-\mathrm{1,}j}}{\Delta x}\right]\\ \frac{{\partial }^2\theta }{\partial y^2}\approx \frac{1}{\Delta y}\left[\frac{{\theta }_{i+\mathrm{1,}j}-{\theta }_{i\mathrm{,}j}}{\Delta y}-\frac{{\theta }_{i\mathrm{,}j}-{\theta }_{i-\mathrm{1,}j}}{\Delta y}\right]\\ q\approx q_{i\mathrm{,}j}\end{array}$

5. Results

This section is dedicated to the simulations and analyses of three different cases, namely, isothermal square cylinder, isothermal tandem cylinders and angled isothermal cylinders, involving the modelling of the onset of turbulence associated with heat-transfer by mixed convection. As methodology, the immersed boundary method was used, using unsteady Reynolds averaged Navier-Stokes equation - URANS together with DNS. To calculate the turbulent viscosity, the Smagorinsky sub-mesh model was used, implemented in the context of the DNS methodology, and the Spalart-Allmaras turbulence model, in the URANS context, both adopted for the two-dimensional case. For the validation of the methodology, flow simulations were carried out around isolated bluff and circular bodies, and in tandem (in line), without boundary movement. The study of the main parameters related to the flow dynamics was carried out, thus demonstrating the reliability, relevance and potential of the methodology of this work for problems involving thermo-fluid dynamics around bodies immersed in incompressible Newtonian fluids. The results obtained in the simulations considered low and moderate Reynolds numbers, that is, ( $1\le Re\le 5000$ ), considering the heat-transfer process by mixed convection (natural and forced convections).

The numerical results obtained were compared and validated with other numerical results from different authors, available in the literature. In the following sections, the results are presented, where all the simulations were performed until the flow field regime was established.

5.1. Simulation and Analysis for a Bluff Isothermal Body

In this section, the results of numerical simulations are presented for low and moderate Reynolds numbers in $\left[\mathrm{1,100}\right]$. The interval $\left[\mathrm{1,40}\right]$ was used to test the code, where we have obtained the results presented in Anjaiah [4], Dhiman [5] and Sharma [6]. For these simulations, a non-uniform mesh with 318 × 164 points was used. The most refined mesh was used in regions with higher gradients, close to the bluff immersed body, where better captures the hydrodynamic fields, while in regions with low gradient, a less refined mesh was used.

The use of a uniform mesh in all computational domains, would considerably increase the computational cost. The simulation of flow develops from left to right. The velocity profile in the input field is uniform, with unit value. In all simulations, the domain dimensions were $L_x=55d$ and $L_y=30d$. The bluff body was positioned at $x=16.5d$ and $y=15.0d$, where d represents the length of the side of the bluff immersed body. The Cartesian mesh used is shown in Figure 8.

The number of points in the mesh, Figure 8, aims to capture thermal and rotation effects (for cases related to border movement, without considering elasticity). For all simulations, a total of 201 points was used for the Lagrangian mesh, with the non-uniform mesh in the region of the isothermal square cylinder, maintaining a minimum amount of 30 meshes inside. For all simulations, the time step used in the calculation process was ${10}^{-a}$, with $a=\mathrm{4,5,6}$ being dynamically calculated by the Courant-Friedrichs-Lewy stability criterion, also known as the CFL condition.

The Prandtl number (Pr) was kept constant, namely, 0, 7 (in the air), in all simulations. For the lateral boundaries of the domain, free boundary conditions

Figure 8. Illustration of a non-uniform mesh.

were used, i.e. zero velocity at the border. At the entrance of the domain, a velocity profile ( $U_{\infty }$ ) was imposed and, at the exit, the derivative is zero for the velocity. In summary, the boundary conditions are as follows

· Domain entry:

$u=U_{\infty },v=0.$

· Domain exit:

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial x}=0.$

· Upper and lower boundaries of the domain:

$\frac{\partial u}{\partial y}=\frac{\partial v}{\partial y}=0.$

For the pressure, the boundary conditions used were Neumann at the entrance and Dirichlet at the exit and on the sides of the domain. The boundary conditions at the entrance of the domain, at the exit, and at the lower border were defined by

$p=0,$

$\frac{\partial p}{\partial x}=0,\frac{\partial p}{\partial y}=0.$

For the temperature, the conditions are analogous, given by

$\theta =0,$

$\frac{\partial \theta }{\partial x}=0,\frac{\partial \theta }{\partial y}=0.$

5.2. Two-Dimensional Flow with Re = 1 Up to Re = 40, for Different Ri Numbers

Figure 9 shows the temperature, pressure, isothermal temperature lines, and the aerodynamic coefficients (drag and lift) around the blunt isothermal body, for $Re=1$ and $Ri=0$ (here, $x\in \left[\mathrm{0,33}\right]$ ). It was found that there is no buoyancy or flow separation around the immersed geometry. Figure 10 and Figure 11 correspond to $Re=1$ with $Ri=0.5$ and $Ri=1.0$, respectively. The goal was to verify the influence of heat around an immersed square cylinder.

Throughout these simulations, it was found that the heating can result in a small “delay” with respect to the separation of the flow around the isothermal cylinder. In other words, the flow is separated more quickly (after numerical observations) even for low Reynolds numbers, when heat transfer is neglected in the simulations which, in turn, may require more heating for the immersed body, as the Reynolds number is increased, as observed, for example, in $Re\ge 2$, for $Ri=0.5$ and $Ri=1.0$. It was also observed that with an increase in the Richardson number (Ri), for any fixed Reynolds number, the region of the so-called von Kármán wake decreases. Therefore, from these simulations, even fixing the Reynolds number, the length of the wake for $Ri=0$ is greater than

Figure 9. (a) Temperature, (b) pressure, (c) isothermal lines, and (d) drag and lift coefficients, for Re = 1 and Ri = 0.

Figure 10. (e) Temperature, (f) pressure, (g) isothermal lines, and (h) drag and lift coefficients, for $Re=1$ e $Ri=0.5$.

Figure 11. (i) Temperature, (j) pressure, (k) isothermal lines, and (l) drag and lift coefficients, for $Re=1$ e $Ri=1$.

the length for $Ri\in \left[\mathrm{0.5,1.0}\right]$. Although there is no vortex shedding, which would be something common in simulations involving moderate and high Reynolds numbers, the damping function used in this work also acts for low Reynolds numbers, for different of Richardson.

Regarding the pressure fields, they were obtained for the same dimensionless time in all the simulations performed, and upstream of the isothermal body, at the point of stagnation the velocity is zero and the static pressure is maximum. Regarding the aerodynamic coefficients (drag coefficients ( $C_d$ ) and lift ( $C_l$ )), observed in Figure 9(d), Figure 10(h) and Figure 11(l), those refer to changes with respect to Re and Ri numbers. Concerning the simulated interval under certain conditions, it appears that the value of the drag coefficient ( $C_d$ ) is quite insensitive to the value of Ri values, while the lift coefficient ( $C_l$ ) is sensitive to Ri. These cases coexist for low values of Prandtl number (Pr), presumably due to the boundary layer being thicker in this case.

What is found, is that the value of the lift coefficient ( $C_l$ ) in the different simulations is zero, even considering the permanent flow. This occurs because the Richardson number is different from zero, which results in the increase of the so-called shear and pressure forces. About the pattern of the isothermal lines presented for $Ri=0$ up to $Ri=1.0$ for different Reynolds values, see Figure 9(c) up to Figure 11(l), we observe that the temperature field not only becomes asymmetrical with the introduction of buoyancy (Ri), but its asymmetry increases with the increase of the Richardson number (Ri).

It was also observed that buoyancy is more present downstream of the isothermal bluff body (isothermal square cylinder) than in relation to the other parts of its surface. In addition, it was verified that the Nusselt number presents a higher value. This is due to the accumulation of the so-called isothermal lines, when compared to the other parts of the immersed body surface. Indeed, in Figure 12, the difference is subtle, since the number of Richardson goes from 0 to 0.5, but it is possible to notice the narrowing of the current lines in the upper and lower parts of the cylinder, i.e. they approach the isothermal body, this is due to the variation of both dimensionless numbers, that is, Reynolds and Richardson numbers.

Other studies using different numerical methodologies have presented a similar situation; e.g. Chatterjee [26], who carried out a study for a pair of square cylinders in a vertical channel configuration in tandem (cylinders in line position) for low Reynolds numbers. The author concluded that the beginning of the flow separation occurs for $Re\in \left[\mathrm{1,2}\right]$, $Re\in \left[\mathrm{2,3}\right]$ and $Re\in \left[\mathrm{3,4}\right]$, for $Ri=0$, $Ri=0.5$ and $Ri=1.0$, respectively. Then, he observes that, if Reynolds number is gradually increased, the flow separation occurs at the trailing edge of the isothermal blunt body and the two (symmetrical) vortices of the recirculation bubble begin to be formed downstream of the obstacle, even in permanent flow regime. However, the opposite trend is observed with the increase of Richardson number at a fixed Reynolds number. This is probably due to the fact that, as the fluctuation (Ri) increases, the velocity gradient on the cylinder surface increases, resulting in reduced pressure on the surface of the immersed body. This behaviour is similar to the flow patterns, with respect to the Reynolds number and the Richardson number, around a circular cylinder under the influence of buoyancy, which is also observed in the works of Badr [2] and Patnaik et al. [27].

Figure 12. Streamlines around a square cylinder for (a) $Re=\mathrm{1,}Ri=0$, and (b) $Re=\mathrm{2,}Ri=0.5$.

In other works, Singh [28] and Moulay [29], in addition to presenting similar issues, there is a critical value for Ri, that is, the authors take into account the beginning of the detachment of vortices, and begin to analyse the influence of this number with the increase of Re, considering the influence of buoyancy in an unstable flow regime. Next, the contour lines with respect to the vorticity are presented, which provide useful information about the flow behaviour, especially close to the square cylinder for the interval considered in Figure 13.

It is observed that the recirculation increases downstream of the immersed body, as the Reynolds number increases for the Richardson number fixed, and decreases when Richardson number increases, keeping the Reynolds number constant. The separation points are located with the help of the vorticity lines. Therefore, it is concluded that as the number of Reynolds or the number of Richardson increases, there is an elongation of the so-called vorticity lines. The standard isotherm lines involving the parameters Re and Ri near the immersed obstacle in the flow are shown in Figure 14, following to the Prandtl number equal to 0.7.

Figure 13. Vorticity fields for (a) $Re=\mathrm{1,}Ri=0$, and (b) $Re=\mathrm{1,}Ri=0.5$.

Figure 14. Isothermal fields for (a) $Re=\mathrm{40,}Ri=0.5$, and (b) $Re=\mathrm{40,}Ri=1$.

5.2.1. Validation of Average Drag Coefficient and Nusselt Number

In Table 1, it is presented a comparison of our mean values (time averages) of the drag coefficients, for $Ri=0$ and $Pr=0.7$, versus different Reynolds numbers, with the numerical results of other authors, for the case of one isothermal bluff body ( [4] [5] [6] ).

In Anjaiah et al. [4], the authors use a numerical scheme, called SMAC-Implicit, implemented in a staggered mesh to solve the Navier-Stokes and heat equations. Their scheme is implicit in the diffusion terms and explicit in the convective terms.

In Dhiman et al. [5] the researchers use the numerical method of finite volumes for a non-scaled mesh. In summary, a semi-explicit method is used to solve the Navier-Stokes equations, in which the momentum equation is explicitly discretized, excepted for the pressure gradient term that is treated implicitly. Consequently, the pressure-velocity coupling is reduced to a Poisson equation for the pressure correction (to prevent oscillations due to the pressure-velocity decoupling in the mesh, they use the interpolation scheme of Rhie and Chow [30] ).

Sharma et al. [6], based on the work of Dhiman [5], use a non-uniform computational mesh, developed in Matlab^Ò. The semi-explicit finite volume method is implemented in the staggered mesh arrangement, used to solve the equations laws, where the momentum equation is explicitly discretized, while the pressure gradient term is treated implicitly.

Considering these three previous works, the algorithm we propose, and corresponding numerical convergence, is valid in view of the corresponding numerical values we found (they agree within 2% - 3% margin).

Table 1. Comparison of the average values of the drag coefficients and the Nusselt number, for $Ri=0$ and $Pr=0.7$, with different values of Re.

5.2.2. Two-Dimensional Flow Fields for Re = 5000 with Ri= 0

In this section, Figure 15 and Figure 16 are discussed. Figure 15(c) and Figure 16(i) present the effective viscosity. It was necessary to adjust the original Smagorinsky model, because in this model there is no damping of the turbulent viscosity close to the walls, which is not physically consistent, because in the regions of the boundary layer, the turbulence models must have relatively low velocities and, consequently, must not have high values for effective viscosity, which can result in premature displacement of the boundary layer.

Figure 15. Smagorinsky model for (a) temperature, (b) pressure, (c) effective viscosity, (d) vorticity, (e) isothermal lines, and (f) drag and lift coefficients ( $Re=\mathrm{5000,}Ri=0$ ).

Figure 16. Spalart-Allmaras model for (g) temperature, (h) pressure, (i) effective viscosity, (j) vorticity, (k) isothermal lines, and (l) drag and lift coefficients ( $Re=\mathrm{5000,}Ri=0$ ).

The Figure 15(c), Figure 15(d) and Figure 16(i), Figure 16(j) show the effective viscosity and the vorticity fields for $Re=5000$. The effective viscosity close to the wall is calculated by the models DNS and URANS, respectively, and has the same approximate magnitude as the effective viscosity for the regions downstream of the immersed body. In the wake regions, the effective viscosity can assume high values, inhibiting the release of vortices as can be observed in Figure 16(j). Regarding the aerodynamic coefficients, presented by Figure 15(f) and Figure 16(l), sharper oscillations are observed in the drag coefficient ( $C_d$ ). In other words, this means that the process of generation and release of vortices is totally swirling in the flow. The same happens with the support coefficient ( $C_l$ ).

5.2.3. Frequency of Vortex Shedding for the Isothermal Square Cylinder

The frequency of vortex shedding is measured by the Strouhal number, St. In Table 2, some numerical values for St are presented, as well the value obtained by Sohankar [31] for the isothermal square cylinder. Note that the results presented in Table 2, corresponding to $Re=500$, have a good numerical agreement. The URANS model is more suitable for the study involving the mean flow behaviour. Therefore, this model captures less information related to transient physical phenomena. In contrast, the LES models have greater resolution in capturing the phenomena associated with large-scales. Therefore, this class of turbulence models is characterized by calculation of lower levels of turbulent viscosity. For this reason, the detachment frequency is slightly higher in LES model, when compared with URANS model.

5.3. Forced and Natural Convection around Isothermal Cylinders in Tandem

In this section, the flow around a pair of heated circular cylinders with different configurations is studied, where the heat transfer process around obstacles has its importance and relevance in engineering problems. Here, two cylinder configurations will be considered, where, in both cases, the two cylinders have equal diameters and the same centre-to-centre distance ( $L_{cc}$ ). In the first case, the angle formed by the segment joining the centres of the two cylinders and the x-axis is zero.

5.3.1. Description of the Problem and Boundary Conditions

In Figures 17-19, the two cylinders are identical and with the same diameter, maintained in tandem configuration, and with angulation (or angular inclination) for the cylinder B. In both cases, the cylinders are confined to a channel with free flow, with uniform velocity ( $U_{\infty }$ ) and constant temperature $T_c\left(>T_{\infty }\right)$. The horizontal and vertical spacing between the cylinders are fixed to $L_u=16.5d$ and $L_d=19.5d$, respectively. These values are chosen to reduce the effect of boundary conditions on the inlet and outlet due to the flow pattern at the cylinder boundary. Moreover, these choices are also consistent with existing works available in the literature, e.g. Santos [8], Santos et al. [12], Mahir and Altaç [32], Sohankar and Etminan [33], Laidoudi and Bouzit [34]. The drag ( $C_d$ ) and lift ( $C_l$ ) coefficients, for the calculation of each cylinder, are evaluated

Table 2. Strouhal number for different Re and Ri values.

Figure 17. Illustration of a two-cylinder configuration: (a) cylinders in tandem ( ${\phi }_m=0^{\circ }$ ), and (b) cylinders in angulation ( $0^{\circ }<{\phi }_m<{90}^{\circ }$ ).

Figure 18. Computational domain with two cylinders in tandem configuration.

Figure 19. Computational domain with two cylinders in angular configuration.

by:

$C_d=C_{dp}+C_{dv}=\frac{2F_d}{\rho U_{\infty }^{2}d},$ (29)

$C_l=C_{lp}+C_{lv}=\frac{2F_l}{\rho U_{\infty }^{2}D},$ (30)

where $C_{lp}$ and $C_{lv}$ represent the lift coefficients due to pressure and viscous forces, respectively. Similarly, $C_{dp}$ and $C_{dv}$ represent the drag coefficients due to the pressure and viscous forces. $F_d$ and $F_l$ are forces of drag and lift acting on the surface of the cylinder, respectively. Thus, the drag and lift coefficients can be obtained from the expressions:

$\{\begin{array}{l}{C}_{dp}=2{\displaystyle {\int }_0^1}\left(p_f-p_r\right)\text{d}y\\ C_{dv}=\frac{2}{Re}{\displaystyle {\int }_0^1}\left(\left({\frac{\partial u}{\partial y}|}_s+{\frac{\partial u}{\partial y}|}_i\right)\text{d}x+\left({\frac{\partial u}{\partial x}|}_f+{\frac{\partial u}{\partial x}|}_r\right)\text{d}y\right)\end{array}$ (31)

$\{\begin{array}{l}{C}_{lp}=2{\displaystyle {\int }_0^1}\left(p_i-p_s\right)\text{d}y\\ C_{lv}=\frac{2}{Re}{\displaystyle {\int }_0^1}\left(\left({\frac{\partial v}{\partial x}|}_f+{\frac{\partial v}{\partial x}|}_b\right)\text{d}x+\left({\frac{\partial v}{\partial y}|}_s+{\frac{\partial v}{\partial y}|}_i\right)\text{d}y\right)\end{array}$ (32)

The subscripts f, b, i and s refer to “front”, “back”, “bottom” and “top” surfaces of the cylinders, respectively.

5.3.2. Flow Fields for Re = 500 and Ri = 0 in Cylinders in Tandem with ${\phi }_m=0$ by Forced Convection

The Figure 20 and Figure 21 present the temperature, pressure, effective viscosity, vorticity, isothermal lines, and aerodynamical drag ( $C_d$ ) and lift ( $C_l$ ) coefficients, for the flow around tandem cylinders, for $Re=500$ and $Ri=0$. The Table 3 represents the distribution of the Nusselt and Strouhal numbers, setting the value of the spacing ratio ( $L_{cc}/d=3$ ) between the cylinders with mounting angles equal to 65˚. In general, it is observed that the average Nusselt number on cylinder A increases more rapidly compared to cylinder B. This is due to the thinning of the boundary-layers around the surface of the cylinder.

5.3.3. Flow Fields for Re = 500 and Ri = 5.0 in Cylinder in Tandem with ${\phi }_m=65^{\circ }$ by Mixed Convection

Concerning Figure 21 and Figure 22, the simulations are based on the work

Figure 20. Smagorinsky model for (a) temperature, and (b) pressure (cylinders A and B, $Re=\mathrm{500,}Ri=0$ ).

Figure 21. Smagorinsky model for (c) effective viscosity, (d) vorticity, (e) isotherms lines, and (f) drag and lift coefficients (cylinders A and B, $Re=\mathrm{500,}Ri=0$ ).

of Zdravkovich [37], where the author relates three flow regimes for tandem cylinders [38]. The author endorses several values for the spacing between the cylinders. The vortex structures and isothermal lines are described for $Re=500$, $Ri=0$ and $L_{cc}/d=3$. For $Re=500$ and $Ri=0$, the formation of vortices between the two cylinders is not observed; however, the oscillation propagates further downstream of the second cylinder, forming a vortex wake in the downstream region. As both vorticity and thermal energy are transported through the flow field, the vorticity and isothermal contours exhibit similar characteristics. The spacing between the centres is still below the critical value. The shear layer separated from the upstream cylinder forms a vortex behind the downstream cylinder. The heat is similarly carried by the flow to the downstream region. Then, what is observed is the vortex wake being established; however, interactions within the vortex wake with the shear layer separated from the downstream cylinder create two separated vortices near the downstream cylinder. The magnitude of the vortices increases due to the faster movement of the fluid and the vortex detachment. The main results for the simulations can be summarized in 1) a wake forms upstream of the second cylinder, but it needs to be checked whether it can be decreased or suppressed with the increase of the distance between the cylinders; 2) the isothermal lines reflect the same behaviour of the pattern of the streamlines (current lines); 3) the average Nusselt number increases with $Re=500$ and for different values of Ri, even keeping the distance

Figure 22. Spalart-Allmaras model for (a) temperature, (b) pressure, (c) effective viscosity, (d) vorticity, (e) isotherms lines, and (f) drag and lift coefficients (cylinders A and B, $Re=\mathrm{500,}Ri=5$ ).

Table 3. Flow parameters: in this work, $Re\in \left\{\mathrm{100,200,500}\right\}$ (ca: cylinder A, cb: cylinder B); ^aMittal and Laccarino [1], ^bMeneghini [35], ^cSlaouti and Sansby [36].

between the cylinders; 4) the thermal buoyancy is suppressed in the recirculation zones of the tandem cylinders, even with a mounting angle; and 5) the thermal buoyancy tends to increase the drag coefficient and the average Nusselt number of the first cylinder more than the second one.

5.3.4. Variations of the Nusselt Number

One of the main purposes of the heat transfer calculations involving cylinders is to determine the local and total transfer around isothermal cylinders. The effect of the flow, especially with respect to the heat transfer, can be better observed by analysing the local heat transfer coefficient, the Nusselt local number. In the Figure 23, for different Richardson numbers, the distribution of the Nusselt number along the perimeter of the upstream and downstream cylinders are provided. For $L_{cc}/d=3$, $Re=100$, $Re=200$ and $Re=500$, for different Richardson numbers, the local distributions of the Nusselt number along the perimeter of the upstream and downstream cylinders is provided. For $L_{cc}/d=3$, although the local profile of the Nusselt number of the upstream cylinder is similar to that of an isolated cylinder, the downstream cylinder has completely different characteristics. As the heat transfer rate is closely related to the flow, the local minimum rates of heat transfer appear at the front and back stagnation points of the downstream cylinder, where the magnitudes of velocities are relatively small.

Thus, in Figure 23(a), the maximum heat transfer from the downstream cylinder is exhibited with a double protuberance occurring at $\theta \approx {57}^{\circ }$ and $\theta \approx {265}^{\circ }$ from the cylinder wall, where thermal layers (also known as thermal plumes) become 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 $L_{cc}/d=3$ and $Re=500$ for different values of Ri, as observed in Figure 23(b). It is important to note that, although the Nusselt’s local distribution of the downstream cylinder resembles that of the upstream cylinder, typified as a large protuberance,

Figure 23. Local variation of Nusselt number: (a) $Re\in \left\{\mathrm{100,200,500}\right\}$ and $Ri=0$ (forced convection), and (b) $Re=500$ and $Ri\in \left\{\mathrm{1,2,5}\right\}$ (natural convection).

its magnitude is smaller than of the upstream cylinder, indicating smaller heat-to-cylinder transfer to downstream areas.

6. Conclusions

This work was devoted to the study of two-dimensional incompressible Newtonian fluids subject to heat transfer processes by mixed convection (natural and forced) at moderate Reynolds ( $Re\in \left[\mathrm{1,5000}\right]$ ) and Richardson numbers ( $Ri\in \left\{\mathrm{0,0.5,1,5}\right\}$ ), with the Smagorinsky and Spalart-Allmaras models of turbulence. Our goal was to continue exploring the potential of the immersed boundary method coupled with the virtual physical model for the analysis of turbulent flows, around cylinders, combined with heat transfer processes by mixed convection in stationary/mobile boundary problems without elasticity. In this context, different results were obtained throughout this work. In general, the methodology used here is applied to the Navier-Stokes equations with a force term that performs the modelling of the immersed interface.

The discretization was performed by finite differences and the boundary conditions were properly imposed on the fluid domain and on the immersed body. The updating of the different numerical fields (velocity, temperature, etc.) to the current time is carried out through the information of previous times. For the first iteration, the first order Euler method was used. For time advancing, the second-order Adams-Bashforth scheme was used, with a second-order spatial centred method, intended to contribute to the application of this methodology to flow problems, around complex bodies, with heat transfer by mixed convection. Without this numerical procedure, the kinetic energy of the so-called physical instabilities would accumulate, resulting in a numerical divergence.

Flow simulations were also carried out around an isothermal square cylinder (with constant temperature on its surface), and around tandem isothermal circular cylinders. It was found that the simulations showed numerical convergence, when compared with other works available in the literature. This paper presents an algorithm whose computational implementation has low processing time (even in mixed convection and turbulence).

Flow parameters, such as Strouhal numbers, lift and drag coefficients, are obtained and compared with the available literature, thus verifying the robustness of the results obtained by our code. The temperature, pressure, turbulence, vorticity and isothermal lines were obtained to understand and interpret the flow and heat transport. In addition to the average Nusselt numbers, the local Nusselt numbers along the perimeter of the cylinders, the upstream and downstream were obtained. It was observed that, for moderate Reynolds numbers, the local Nusselt numbers of the downstream of the cylinder exhibit a generation of thermal plumes that move upwards, which is related to the increase in the number of Richardson.

Acknowledgements

This work was supported by Project STRIDE NORTE-01-0145-FEDER-000033, funded by ERDF NORTE 2020; by project MAGIC POCI-01-0145-FEDER-032485, funded by FEDER via COMPETE 2020 - POCI and by FCT/MCTES via PIDDAC; and by CMUP, which is financed by national funds through FCT - Fundação para a Ciência e a Tecnologia, I.P., under the project with reference UIDB/00144/2020.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Mittal, R. and Iaccarino, G. (2005) Immersed Boundary Method. Annual Review of Fluid Mechanics, 37, 239-261. https://doi.org/10.1146/annurev.fluid.37.061903.175743
[2]	Badr, H. and Dennis, S. (1985) Laminar Forced Convection from a Rotating Cylinder. International Journal of Heat and Mass Transfer, 28, 253-264. https://doi.org/10.1016/0017-9310(85)90027-4
[3]	Badr, H., Countanceau, M., Dennis, S. and Menard, C. (1990) Unsteady Flow past a Rotating Circular Cylinder at Reynolds Numbers 103 and 104. Journal of Fluid Mechanics, 220, 490-484. https://doi.org/10.1017/S0022112090003342
[4]	Anjaiah, N., Dhiman, A.K. and Chhabra, R.P. (2006) Mixed Convection Heat Transfer from a Square Cylinder to Power-Law Fluids in Cross-Flow. Fluids Engineering Division Summer Meeting, Miami, 17-20 July 2006, 1435-1444. https://doi.org/10.1115/FEDSM2006-98072
[5]	Dhiman, A.K., Anjaiah, N., Chhabra, R.P. and Eswaran, V. (2007) Mixed Convection from a Heated Square Cylinder to Newtonian and Power-Law Fluids. Journal of Fluids Engineering, 129, 506-513. https://doi.org/10.1115/1.2436586
[6]	Sharma, N., Dhiman, A.K. and Kumar, S. (2012) Mixed Convection Flow and Heat Transfer across a Square Cylinder under the Influence of Aiding Buoyancy at Low Reynolds Numbers. International Journal of Heat and Mass Transfer, 55, 2601-2614. https://doi.org/10.1016/j.ijheatmasstransfer.2011.12.034
[7]	Lai, M.C. and Peskin, C.S. (2000) An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical Viscosity. Journal of Computational Physics, 160, 705-719. https://doi.org/10.1006/jcph.2000.6483
[8]	Santos, R.D. (2014) Análise Bidimensional Termo-Fluido Dinamica de Cilindros Rotativos com o Método da Fronteira Imersa/Modelo Fsico Virtual. Master’s Thesis, Universidade Federal de Itajubá-Minas Gerais, Brazil.
[9]	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. Journal of Computational Physics, 189, 351-370. https://doi.org/10.1016/S0021-9991(03)00214-6
[10]	Carvalho, G.B. (2003) Estudo Experimental do Escoamento em Torno de Cilindros Circulares em Movimento de Rotao. Master’s Thesis, Universidade Estadual Paulista-UNESP, So Paulo.
[11]	Ashrafizadeh, A. and Hosseinjani, A.A. (2017) A Phenomenological Study on the Convection Heat Transfer around Two Enclosed Rotating Cylinders via an Immersed Boundary Method. International Journal of Heat and Mass Transfer, 107, 667-685. https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.078
[12]	Santos, R.D., Gama, S.M. and Camacho, R.G. (2018) Two-Dimensional Simulation of the Navier-Stokes Equations for Laminar and Turbulent Flow around a Heated Square Cylinder with Forced Convection. Applied Mathematics, 9, 291-312. https://doi.org/10.4236/am.2018.93023
[13]	Pope, B.P. (2011) Turbulent Flows. Cambridge University Press, Cambridge.
[14]	Unverdi, S.O. and Tryggvagson, G. (1994) A Front-Tracking Method for Viscous, Incompressible, Multi-Fluid Flows. Journal of Computational Physics, 100, 25-37. https://doi.org/10.1016/0021-9991(92)90307-K
[15]	Schneider, G. and Zedan, M. (1981) A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems. Numerical Heat Transfer, 4, 1-19. https://doi.org/10.1080/01495728108961775
[16]	Peskin, C.S. (1972) Flow Patterns around Heart Valves: A Numerical Method. Journal of Computational Physics, 10, 252-271. https://doi.org/10.1016/0021-9991(72)90065-4
[17]	Peskin, C.S. (1977) Numerical Analysis of Blood Flow in the Heart. Journal of Computational Physics, 25, 220-252. https://doi.org/10.1016/0021-9991(77)90100-0
[18]	Juric, D. (1996) Computations of Phase Change. PhD Thesis, University of Michigan, Ann Arbor.
[19]	Lesieur, M. (2012) Turbulence in Fluids (Vol. 40). Springer Science & Business Media, Berlin.
[20]	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<0099:GCEWTP>2.3.CO;2
[21]	Spalart, P. and Allmaras, S. (1992) A One-Equation Turbulence Model for Aerodynamics Flows. Recherce Aeropatiale, No. 1, 5-21. https://doi.org/10.2514/6.1992-439
[22]	Kostić, C. (2015) Review of the Spalart-Allmaras Turbulence Model and Its Modifications to Three-Dimensional Supersonic Configurations. Scientific Technical Review, 65, 43-49. https://doi.org/10.5937/STR1501043K
[23]	Chorin A.J. (1968) Numerical Solution of the Navier-Stokes Equations. Mathematics of Computation, 22, 745-762. https://doi.org/10.1090/S0025-5718-1968-0242392-2
[24]	Kim, J. and Moin, P. (1985) Application of Fractional-Step Method to Incompressible Navier-Stokes Equation. Journal of Computational Physics, 59, 308-323. https://doi.org/10.1016/0021-9991(85)90148-2
[25]	Patankar, S. (1980) Numerical Heat Transfer and Flow. CRC Press, Boca Raton.
[26]	Chatterjee, D. (2000) Mixed Convection Heat Transfer from Tandem Square Cylinder in Vertical Channel at Low Reynolds Numbers. Numerical Heat Transfer, Part A: Applications, 58, 740-755. https://doi.org/10.1080/10407782.2010.516703
[27]	Patnaik, B.V., Narayana, P.A. and Seetharamu, K. (1999) Numerical Simulation of Vortex Shedding Past a Circular Cylinder under the Influence of Buoyancy. International Journal of Heat and Mass Transfer, 42, 3495-3507. https://doi.org/10.1016/S0017-9310(98)00373-1
[28]	Singh, S., Panigrahi, P. and Muralidhar, K. (2007) Effect of Buoyancy on the Wakes of Circular and Square Cylinders: A Schlieren-Interferometric Study. Experiments in Fluids, 43, 101-123. https://doi.org/10.1007/s00348-007-0329-8
[29]	Moulay, M.A., Belkady, M. and Aounallah, A. (2017) Effect of Opposing Buoyancy on Two-Dimensional Laminar Flow and Heat Transfer across a Confined Circular Cylinder. Mechanics, 23, 859-865. https://doi.org/10.5755/j01.mech.23.6.17291
[30]	Rhie, C.M. and Chow, W.L. (1983) Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation. AIAA Journal, 21, 1525-1532. https://doi.org/10.2514/3.8284
[31]	Sohankar, A., Norberg, C. and Davidson, L. (1999) Simulation of Three-Dimensional Flow around a Square Cylinder at Moderate Reynolds Numbers. Physics of Fluids, 11, 288-306. https://doi.org/10.1063/1.869879
[32]	Mahir, N. and Alta, Z. (2008) Numerical Investigation of Convective Heat Transfer in Unsteady Flow past Two Cylinders in Tandem Arrangements. International Journal of Heat and Fluid Flow, 29, 1309-1318. https://doi.org/10.1016/j.ijheatfluidflow.2008.05.001
[33]	Sohankar, A. and 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, 733-751. https://doi.org/10.1002/fld.1909
[34]	Laidoudi, H. and Bouzit, M. (2018) Mixed Convection in Poiseuille Fluid from an Asymmetrically Confined Heated Circular Cylinder. Thermal Science, 22, 821-934. https://doi.org/10.2298/TSCI160424172L
[35]	Meneghini, J.R. (1993) Numerical Simulation of Bluff Body Flow Control Using a Discrete Vortex Method. PhD Thesis, University of London, London. https://doi.org/10.1007/978-3-662-00414-2_57
[36]	Slaouti, A. and Stansby, P. (1992) Flow around Two Circular Cylinders by the Random-Vortex Method. Journal of Fluids and Structures, 6, 641-670. https://doi.org/10.1016/0889-9746(92)90001-J
[37]	Zdravkovich, M. (1985) Flow Induced Oscillations of Two Interfering Circular Cylinders. Journal of Sound and Vibrations, 101, 511-521. https://doi.org/10.1016/S0022-460X(85)80068-7
[38]	Yen, S., San, K. and Chuang, T. (2008) Interactions of Tandem Square Cylinders at Low Reynolds Numbers. Experimental Thermal and Fluid Science, 32, 927-938. https://doi.org/10.1016/j.expthermflusci.2007.07.001