1. Introduction

The idea of Immersed Boundary Method (IBM) was proposed by [1] and [2] to simulate blood flows in mitral valves, and was later introduced to solve problems in computational fluid dynamics (e.g., [3] [4] and [5] ). This method provides a new numerical tool to investigate the flow around the immersed body with complex geometry in Cartesian coordinates. The immersed boundary method was reviewed in [6] , presenting there an extensive list of bibliographic references. Problems related to flow around immersed bodies have always been studied in many areas, especially in engineering applications such as aerodynamics or even in more complex configuration such as in aircraft and automobiles. In such cases, the involvement of complex geometries is inevitable in many problems. Thus, the classical simulation methods used have some disadvantages, especially in the study of fluid-structure interaction problems, for example, in the simulation of mobile or deformable border simulation. In these cases, traditional methods produce good results, but in general, with high computational costs. Basically, two methodologies are used for this type of problems.

The first methodology uses unstructured grids to describe complex geometries, with the remeshing process for the case of deformable bodies. The second methodology used the immersed boundary method, as proposed by [1] . The latter has some advantages, for example, the possibility of simulating complex geometries in Cartesian grids without the need to use the remeshing process in each time step of the iterative scheme, during deformation or movement of the border, without increasing the computational cost. An important work, developed by [7] , brings a study of a simulation of the flow around rigid bodies. The authors use a forcing term to reach the rigid boundary required to simulate the specified boundary conditions: the velocity of the fluid is equivalent to the velocity of the surface.

The work [8] 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 to form the boundary. The selection of constant values in the forcing term is discussed based on the error of velocity at the boundary. Then, the VBM method is used to simulate the flow around bluff bodies, where the qualitative description about this simulation is firstly presented by providing the current simulation’s streamlines and pressure contours. The phenomena illustrated are the same as others available in the literature. The VBM feasibility is verified by the good agreement between this simulation’s results and other reference’s outcomes.

In this paper, unlike other works, it is presented the Immersed Boundary Method as a method of capturing interfaces with the Virtual Physical Model (IBM/VPM) being used to simulate the flow in two-dimensional space, with the presence of immersed geometries in incompressible fluid with forced convection around the heated square cylinder (also known as bluff bodies, with its heated surface at constant temperature, and where there is no elasticity in the immersed body). The VPM was developed by [9] , where the model is based on Navier-Stokes equations, the presence of the body and its heating are calculated and included as a forcing term in the equations of momentum and energy. This methodology, that does not require high computational cost or facilities, can also be adapted to other researcher areas, such as systems involving renewable energy, energy converters of sea wave mechanics, solar heating systems, refrigeration systems for electronic components, cooling of nuclear reactors, etc.

The VPM model has the ability to self-adjust to the flow as the force required to “stop” the fluid particles near the interface is automatically calculated, without the need to adjust variables, unlike the original method developed by [1] and [8] . The interfacial force is calculated at Lagrangian points and distributed to neighboring Eulerian points, with the help of a type of Gaussian function.

In the IBM method, the force $f$ is introduced into the Navier-Stokes equations to model the solid-fluid interface. Analogously, the heated body is modeled by a forcing term q added to the energy equation. Thus, the method is based on a mixed formulation with a grid for the fluid (Eulerian grid) and another for fluid-solid interface (Lagrangian grid). The bulk of this paper is to verify the applicability and potentiality of IBM/VPM to model and simulate a flow over a heated square cylinder for $40\le \text{Re}\le 5\times {10}^3$ , considering the aerodynamic coefficients, this is, the drag and lift coefficients, the temperature, the Strouhal and Nusselt numbers. For the highest number of Reynolds here considered $\left(\text{Re}=5\times {10}^3\right)$ , two turbulence models are used, namely, Smagorinsky and Spalart-Allmaras models.

This paper is organized as follows. In the Section 2, it is presented the IBM methodology. Section 3 is devoted to the Virtual Physical Model used to solve the Navier-Stokes equations. The turbulence models are in Section 4. In the next Section, the numerical method is discussed. The results are in Section 6, and the concluding remarks are presented in 7.

2. The Methodology of Immersed Boundary

Mathematical Formulation

The physical phenomena of fluid mechanics are described by a set of partial differential equations, describing the conservation of mass (continuity), momentum and energy. For the non-isothermal case, theses equations, in the context IBM/VPM, are written under the following hypotheses:

Ÿ The fluid is Newtonian and incompressible with constant properties. The buoyancy term, based on the Boussinesq approximation, does not appear since we are only considering the case of forced convection;

Ÿ The immersed body (heated square cylinder) is describe by the VPM model instead of by the direct imposition of the boundary conditions of velocity and temperature;

Ÿ The energy generation term is neglected, because neither the effects of internal heat are considered (for example, absorption or emission radiation) nor the humidity which could be responsible for the latent heat exchange;

Ÿ There is no Coriolis force nor rotation effects of the coordinate system.

The equations for Newtonian incompressible flow with forced convection can be written, in dimensionless form, as

$\{\begin{array}{l}\nabla \cdot U=0\mathrm{,}\\ \frac{\partial U}{\partial t}+\nabla \cdot \left(UU\right)=-\nabla p+{\text{Re}}^{-1}\Delta U+f\mathrm{,}\\ \frac{\partial \Theta }{\partial t}+\nabla \cdot \left(U\Theta \right)={\text{Pe}}^{-1}\Delta \Theta +q\mathrm{,}\end{array}$ (1)

where U is the incompressible velocity field, p is the pressure, and Re is the Reynolds number. The Eulerian body-force $f$ is used to mimic the effects of the immersed body, W, in the flow. The Reynolds number is defined by $\text{Re}=UL/\nu$ , where U, L and v are the reference length, the reference velocity, and the kinematic viscosity, respectively. Pe is the Péclet number, defined by $\text{Pe}=\text{Re}\times \text{Pr}$ , where Pr is the Prandtl number, given by $\text{Pr}=c_p\mu /k$ . Here, $c_p$ is the heat capacity and k is the coefficient of thermal conductivity. The dimensonless temperature is defined by $\Theta =\left(\theta -{\theta }_0\right)/\left({\theta }_c-{\theta }_i\right)$ , being ${\theta }_c$ and ${\theta }_i$ the square cylinder temperature and the domain inlet temperature, respectively. The immersed boundary is represented by Lagrangian points, where the Lagrangian force is defined. The Eulerian and Lagrangian forces are related to each other through a regularized delta-function. The Eulerian force $f$ , added to the momentum equation, is mathematically represented by

$f\left(x\mathrm{,}t\right)={\displaystyle {\int }_{\partial \Omega }}F\left(x_k\mathrm{,}t\right)\delta \left(x-x_k\right)\text{d}{x}_k\mathrm{,}$ (2)

where, $x$ and $x_k$ are the position vectors of the Eulerian and Lagrangian points, respectively; $\partial \Omega$ is the domain of the immersed boundary, $f$ and $F$ are the Eulerian and Lagrangian forces, respectively. Similarly, in the energy equation (3rd equation in system (1)), the dimensionless external heat source q, which is the mutual interaction energy between fluid and immersed boundary, is expressed as follows

$q\left(x\mathrm{,}t\right)={\displaystyle {\int }_{\partial \Omega }}Q\left(x_k\mathrm{,}t\right)\delta \left(x-x_k\right)\text{d}{x}_k\mathrm{,}$ (3)

where $Q$ is the heat source on the Lagrangian point $x_k$ at the immersed boundary. In this paper, the boundary conditions of velocity and temperature are imposed by integration of the equations of momentum and energy, of the system (1). These terms are intended to model the immersed geometry by

$f\left(x\mathrm{,}t\right)={\displaystyle {\int }_{\partial {\Omega }_b}}F\left(x_k\mathrm{,}t\right)\delta \left(x-x_k\right)\text{d}{x}_k\simeq {\displaystyle \underset{k=1}{\overset{N}{\sum }}}F\left(x_k\mathrm{,}t\right)D\left(x-x_k\right)\Delta s^2\mathrm{,}$ (4)

$q\left(x\mathrm{,}t\right)={\displaystyle {\int }_{\partial {\Omega }_b}}Q\left(x_k\mathrm{,}t\right)\delta \left(x-x_k\right)\text{d}{x}_k\simeq {\displaystyle \underset{k=1}{\overset{N}{\sum }}}Q\left(x_k\mathrm{,}t\right)D\left(x-x_k\right)\Delta s^2\mathrm{,}$ (5)

where N is the number of Lagrangian points and $\Delta s$ is the discrete Lagrangian distance between two points. The term $D\left(x-x_k\right)$ is the distribution function suggested by [1] [10] and [11] , being approximated by

$D\left(x-x_k\right)\simeq h^{-2}\left(g_1\left[\left(x_k-x_i\right)h^{-1}\right]g_1\left[\left(y_k-y_i\right)h^{-1}\right]\right)\mathrm{,}$ (6)

with

$g_1\left(r\right)=\{\begin{array}{l}{g}_2(r),\left|r\right|<1\\ 0.5-g_2\left(2-\left|r\right|\right),1<\left|r\right|<2\\ 0,\left|r\right|>2\end{array}$ (7)

where h is the Euclidian grid size. The function $g_2\left(r\right)$ is given by

$g_2\left(r\right)=\frac{1}{8}\left(3-2\left|r\right|+\sqrt{1+4\left|r\right|-4{\left|r\right|}^2}\right)$ (8)

where r denotes $\left[\left(x_k-x_i\right)h^{-1}\right]$ or $\left[\left(y_k-y_i\right)h^{-1}\right]$ , and $\left(x_i\mathrm{,}{y}_i\right)$ are the coordinates of the Eulerian points. In this paper, the domain is discretized by a Cartesian grid, as shown in Figure 1.

3. The Virtual Physical Model

The Virtual Physical Model (VPM) was proposed by [9] to calculate the Lagrangian force field, based only in the momentum equation. All the Navier-Stokes terms are calculated over the Lagrangian points, given by

$F\left(x_k\mathrm{,}t\right)=\frac{\partial U}{\partial t}\left(x_k\mathrm{,}t\right)+\nabla \left(U\left(x_k\mathrm{,}t\right)\cdot U\left(x_k\mathrm{,}t\right)\right)-{\text{Re}}^{-1}\Delta U\left(x_k\mathrm{,}t\right)+\nabla p\left(x_k\mathrm{,}t\right)\mathrm{,}$ (9)

where the r.h.s. parcels, from left to right, are referred to as acceleration force, inertial force, viscous force and pressure force, respectively. Similarly, for the thermal source of the fluid particle 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)-{\text{Pe}}^{-1}\Delta \Theta \left(x_k\mathrm{,}t\right)\mathrm{,}$ (10)

where the r.h.s. parcels, from left to right, are called local rate of change of temperature, heat dissipation rate due to convection and diffusive transport rate of the thermal energy. The terms $F$ and $Q$ are calculated in the interface points, through interpolation of the pressure, velocity and temperature fields, being calculated in the Lagrangian grid. After calculating $Q\left(x_k\mathrm{,}t\right)$ , the Eulerian variable $q\left(x\mathrm{,}t\right)$ , in Equation (5), is obtained to finally calculate the new temperature through the 3rd equation in (1).

3.1. Calculation of Velocity, Pressure and Temperature

For the calculation of pressure and temperature derivatives at each Lagrangian point, it is necessary to obtain the pressure and temperature value on the interface point $x_k$ . For the calculation of pressure and temperature, it is necessary to use an auxiliary point, here called P, according to Figure 2.

The derivatives for the calculation of pressure force is performed by the finite difference method

$\{\begin{array}{l}\frac{\partial p}{\partial x}\simeq \frac{p_2-p_1}{x_2-x_1}\\ \frac{\partial p}{\partial y}\simeq \frac{p_4-p_3}{y_4-y_3}\end{array}$ (11)

The pressure and the temperature at the auxiliary point P belong to an Eulerian cell. So, to obtain the velocity, pressure and temperature, we use the following approximations (for details, see [12] ):

$\{\begin{array}{l}u\left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(x_i-x_k\right)u\left(x_i\right)\text{d}x\text{d}y\\ v\left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(y_i-y_k\right)v\left(x_i\right)\text{d}x\text{d}y\\ p\left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(x_i-x_k\right)p\left(x_i\right)\text{d}x\text{d}y\\ p\left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(y_i-y_k\right)p\left(x_i\right)\text{d}x\text{d}y\\ \Theta \left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(x_i-x_k\right)\Theta \left(x_i\right)\text{d}x\text{d}y\\ \Theta \left(x_k\right)\simeq {\displaystyle \underset{i}{\sum }}D\left(y_i-y_k\right)\Theta \left(x_i\right)\text{d}x\text{d}y\end{array}$ (12)

where the quantities $u\left(x_k\right)$ and $v\left(x_k\right)$ , in (12), represent the velocities of the fluid on the interface, taking into account the internal and external velocities to the interface on the Eulerian grid. The terms $u\left(x_i\right)$ , $v\left(x_i\right)$ , $p\left(x_i\right)$ and

$\Theta \left(x_i\right)$ are, respectively, the velocities in the x- and y-direction, the pressure and temperature in the Eulerian grid which will be interpolated; the terms $u\left(x_k\right)$ , $v\left(x_k\right)$ , $p\left(x_k\right)$ and $\Theta \left(x_k\right)$ , are the velocities, pressure and temperature on the auxiliary points, respectively. The derivative that composes the different terms are calculated using the second-order Lagrange polynomial.

The systems of Equations (13) and (14), in the variable f, represent the velocity and the pressure in the x- and y-direction, respectively. After the interpolation of velocity, pressure and temperature at the interface and in the auxiliary points, one determines the derivatives that constitute the terms for the calculation of the Lagrangian polynomial of the second order. Denoting the velocity or temperature components by (f), the calculation of the first and second derivatives in x and y directions, respectively, is given by

$\{\begin{array}{l}\frac{\partial \varphi }{\partial x}=\frac{\left(x_i-x_k\right)+\left(x_i-x_2\right)}{\left(x_1-x_2\right)\left(x_1-x_k\right)}{\varphi }_1+\frac{\left(x_i-x_k\right)+\left(x_i-x_1\right)}{\left(x_2-x_1\right)\left(x_2-x_k\right)}{\varphi }_2+\frac{\left(x_i-x_1\right)+\left(x_i-x_2\right)}{\left(x_k-x_1\right)\left(x_k-x_2\right)}{\varphi }_k\\ \frac{{\partial }^2\varphi }{\partial x^2}=\frac{2}{\left(x_1-x_2\right)\left(x_1-x_k\right)}{\varphi }_1+\frac{2}{\left(x_2-x_1\right)\left(x_2-x_k\right)}{\varphi }_2+\frac{2}{\left(x_k-x_1\right)\left(x_k-x_2\right)}{\varphi }_k\end{array}$ (13)

$\{\begin{array}{l}\frac{\partial \varphi }{\partial y}=\frac{\left(y_i-y_k\right)+\left(y_i-y_4\right)}{\left(y_3-y_4\right)\left(y_3-y_k\right)}{\varphi }_3+\frac{\left(y_i-y_k\right)+\left(y_i-y_3\right)}{\left(y_4-y_3\right)\left(y_4-y_k\right)}{\varphi }_4+\frac{\left(y_i-y_3\right)+\left(y_i-y_4\right)}{\left(y_k-y_3\right)\left(y_k-y_4\right)}{\varphi }_k\\ \frac{{\partial }^2\varphi }{\partial y^2}=\frac{2}{\left(y_3-y_4\right)\left(y_3-y_k\right)}{\varphi }_3+\frac{2}{\left(y_4-y_3\right)\left(y_4-y_k\right)}{\varphi }_4+\frac{2}{\left(y_k-y_3\right)\left(y_k-y_4\right)}{\varphi }_k\end{array}$ (14)

where ${\varphi }_1$ , ${\varphi }_2$ , ${\varphi }_3$ and ${\varphi }_4$ are obtained by interpolation of the nearest Eulerian variables. The coordinates of the auxiliary points (1, 2, 3 and 4) and the coordinates of the Lagrangian point $x_k$ are, respectively, the pairs $\left(x_k\mathrm{,}{y}_k\right)$ , $\left(x_1\mathrm{,}{y}_1\right)$ , $\left(x_2\mathrm{,}{y}_2\right)$ , $\left(x_3\mathrm{,}{y}_3\right)$ and $\left(x_4\mathrm{,}{y}_4\right)$ . The points 1, 2, 3 and 4 must be located outside the interface, in order to calculate the force to be independent of the flow properties in the interior. The distance between the points 1, 2, 3 and 4 is $\Delta x$ . Thus, the calculation of the inertial, viscous and pressure forces are independent of the flow inside the interface. The same process holds for the energy equation. The acceleration force, which is one of the terms of the total Lagrangian (or interface) force, is obtained by an approximation using the finite difference method, according to the expression,

$\frac{\partial U}{\partial t}\left(x_k\mathrm{,}t\right)\simeq \frac{U_k-U_{fk}}{\Delta t}\mathrm{,}$ (15)

where $U_k$ represents the interface velocity vector and $U_{fk}$ represents the velocity vector of the fluid at the same interface position. This acceleration force represents the most influential part in the calculation of the total Lagrangian force, and can be interpreted as the portion that guarantees the non-slip condition. Similarly, the time derivative for the temperature is approximated by

$\frac{\partial {\theta }_k}{\partial t}\simeq \frac{{\theta }_k-{\theta }_{fk}}{\Delta t}\mathrm{.}$ (16)

3.2. Calculation of Parameter L₂

Another important calculation is the determination of the parameter $L_2$ , which provides the difference between the dimensionless velocity of the fluid at the interface and the velocity at the interface itself and its temperature. The numerical value is around 10⁻³, being a coherent and acceptable value, as proposed by [13] . The calculation of $L_2$ is performed by

$L_2\simeq \sqrt{{\displaystyle \underset{k=1}{\overset{n_p}{\sum }}}\frac{{\left(u_k\left(x_k\right)-u_{fk}\left(x_k\right)\right)}^2+\left(v_k\left(x_k\right)-v_{fk}\left(x_k\right)\right)}{n_p{\left|U_{\infty }\right|}^2}}\mathrm{,}$ (17)

$L_2\simeq \sqrt{{\displaystyle \underset{k=1}{\overset{n_p}{\sum }}}\frac{{\left({\theta }_k\left(x_k\right)-{\theta }_{fk}\left(x_k\right)\right)}^2}{n_p{\left|U_{\infty }\right|}^2}}\mathrm{,}$ (18)

where $\left(u_k\mathrm{,}{u}_{fk}\right)$ , $\left(v_k\mathrm{,}{v}_{fk}\right)$ , $\left({\theta }_k\right)$ and $\left({\theta }_{fk}\right)$ are the components of velocities and temperatures at the interface, respectively; $n_p$ is the number of Lagrangian points in the immersed interface.

3.3. The Lagrangian Calculation of the Distribution of Force and Heat Source

After calculating the terms of Equations (9) and (10), and obtaining the values for $F$ and $Q$ , the Eulerian terms for $f$ and q are evaluated in order to obtain the distribution of the interfacial force and the thermal field. The calculation process is performed as follows:

$f\left(x_i\right)\simeq {\displaystyle \underset{k}{\sum }}{D}_{i\mathrm{,}j}\left(x_i-x_k\right)F\left(x_k\right)\Delta s\left(x_k\right)\mathrm{,}$ (19)

$q\left(x_i\right)\simeq {\displaystyle \underset{k}{\sum }}{D}_{i\mathrm{,}j}\left(x_i-x_k\right)Q\left(x_k\right)\Delta s\left(x_k\right)\mathrm{,}$ (20)

where $f\left(x_i\right)$ is the force on each Eulerian node, while $F\left(x_k\right)$ is the force of each Lagrangian node being distributed to the Eulerian nodes. Figure 3 illustrates the distribution of the interfacial force and the thermal source at each Lagrangian node in the body immersed in the flow, resulting in a heated interfacial force field, allowing the fluid particle to recognize the solid boundary and the heated body.

The average distance between two adjacent points is represented by $\Delta s\left(x_k\right)$ and $D_{i\mathrm{,}j}$ is a distribution function. In Equation (20), it is shown the thermal source term, $q\left(x_i\right)$ , for each Eulerian node, due to the presence of the heated immersed body, and $Q\left(x_k\right)$ , the thermal source in each Lagrangian node, being distributed to the Eulerian nodes, thus forming an Eulerian thermal field, that acts on the particles near the border.

4. Turbulence Modeling

In addition to laminar regime flow, we have also simulated (turbulent) flow for Reynolds numbers going up to $5\times {10}^3$ . Thus, it is necessary to use turbulence models to close the Navier-Stokes equations. Two models were implemented, namely the Smagorinsky model, with LES (Large Eddy Simulation) methodology, and the Spalart-Allmaras model (S-A), with URANS model (Unsteady Reynold Averange Navier-Stokes) (for further details, see [14] ).

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 filtered components and floating components, and then all terms are filtered. In the LES methodology the largest turbulence structures (eddy structures) are solved by filtered equations and only the smallest structures are modeled, since they are (statistically) homogeneous and isotropic. The scale of the small structures is calculated by the size of the grid used for the solution of the Navier-Stokes filtered equations. In other words, the width of the filter becomes a function of the grid. In this way, turbulent structures that are smaller than the grid solution are modeled by the so-called sub-grid models.

4.1. Smagorinky Model

The Smagorinsky sub-grid model [15] is based on the balance between the production of turbulent kinetic energy and isotropic dissipation of turbulent energy. The turbulent viscosity, ${\nu }_t$ , is calculated as a function of the shear rate, $S_{ij}$ , and the length scale, $\mathcal{l}$ , given by

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

where $C_S=0.18$ is the Smagorinsky constant and $\mathcal{l}=\sqrt{\Delta x\Delta y}$ is the length of the sub-grid scale.

4.2. Spalart-Allmaras Model

The Spalart-Allmaras model is a one-equation model that solves the transport equation for the kinematic eddy turbulent viscosity [16] , without involving turbulent energy, dissipation or vorticity calculations, available in other models. Thus, the model Spalart-Allmaras concentrates into a single parameter the behavior of the turbulence, being classified as a closed model.

The equation for the turbulent viscosity is constructed using mainly flow empirical considerations, dimensional analysis and Galileo’s principle of relativity. The model uses a working variable, $\tilde{\nu }$ , given by the following transport equation:

$\begin{array}{c}\frac{\partial \tilde{\nu }}{\partial t}+\frac{\partial }{\partial x_j}\left(u_j\tilde{\nu }\right)=c_{b1}\left(1-f_{t2}\right)\tilde{S}\tilde{\nu }+\frac{1}{\sigma }[\frac{\partial }{\partial x_j}\left(\left(\nu +\tilde{\nu }\right)\frac{\partial \tilde{\nu }}{\partial x_j}\right)+c_{b2}\frac{\partial \tilde{\nu }}{\partial x_j}\frac{\partial \tilde{\nu }}{\partial x_j}]\\ -\left[c_w{f}_w-\frac{c_{b1}}{k^2}{f}_{t2}\right]\left[\frac{\tilde{\nu }}{d_w}\right]-f_{t1}\Delta U^2\mathrm{,}\end{array}$ (22)

where the terms on the right side, respectively, represent: the production of turbulent viscosity, the molecular and turbulent diffusion of $\tilde{\nu }$ , the dissipation of $\tilde{\nu }$ , the destruction of $\tilde{\nu }$ , which reduces the turbulent viscosity near the wall and, finally, the terms that model the transition effects, indicated by sub-index t. The turbulence viscosity, ${\nu }_t$ , is calculated from the auxiliary variable of the Spalart-Allmaras model and damped by the function $f_{v1}$ along the wall, and is given by

${\nu }_t=\tilde{\nu }{f}_{v1},$ (23)

where

$f_{v1}=\frac{{\chi }^3}{{\chi }^3+C_{v1}^3},$ (24)

with

$\chi \equiv \frac{\tilde{\nu }}{\nu }\mathrm{.}$ (25)

The production term in the transport Equation (22) also needs a correction, which is performed by replacing the parameter S with a modified variable $\tilde{S}$ , which also depends on the influence of a damping function $f_{v2}$ . Thus, $\tilde{S}$ is defined by

$\tilde{S}\equiv S+\frac{\tilde{\nu }}{{\left(k{d}_w\right)}^2}{f}_{v2}\mathrm{,}$ (26)

where $d_w$ is the distance to the nearest wall, and S is the modulus of the strain rate given by

$S=\sqrt{2{\bar{S}}_{ij}{\bar{S}}_{ij}}.$ (27)

The term of destruction originally formulated by [16] presents problems, and this term decreases very slowly in certain regions of the boundary-layer. To correct this deficiency it is defined a dimensionless function $f_w$ , which adjusts the term of destruction. The function $f_w$ is adjusted to a unit value of the logarithmic region of the boundary-layer, increasing destruction of the term as it approaches the wall and going to zero for more distant regions of the wall, and is defined by the formulas

$f_w=g{\left(\frac{1+c_{w3}^6}{g^6+c_{w3}^6}\right)}^{1/6}\mathrm{,}$ (28)

$g=r+c_{w2}\left(r^6-r\right),$ (29)

with

$r\equiv \frac{\tilde{\nu }}{\tilde{S}{k}^2{d}_w^2}.$ (30)

(Here, k is an empirical constant of the model.) The influence of these terms allows the control of the transition in two different aspects: in the maintenance of the laminar flow in the desired regions or beginning the transition region. The control is performed with the addition of a source term, being controlled by the function $f_{t1}$ and a reduction in the production term of the turbulent viscosity controlled by the function $f_{t2}$ , being defined as follows

$f_{t1}=c_{t1}{g}_t\exp \left(-c_{t2}\frac{{\omega }_t^2}{\Delta U^2}\left[d_w^2+g_t^2{d}_t^2\right]\right),$ (31)

$f_{t2}=c_{t3}\exp \left(-c_{t4}\exp {\chi }^2\right),$ (32)

where $d_t$ is the distance to the start point of the transition, ${\omega }_t$ is the vorticity at the point of transition of the boundary-layer, and $\Delta U$ is the norm of the velocity difference between the flow and the transition point. The function $g_t$ is defined by $min\left\{\mathrm{0.1\; <mo>\; ;\; </mo>}\Delta U/\left(\omega \Delta x_t\right)\right\}$ , where $\Delta x_t$ is the size of the grid along the wall in the transition region. These terms were neglected, as the simulations were only performed for fully developed turbulence regimes. The other constants of the model are $\sigma =2/3$ , $c_{b1}=0.1355$ , $c_{b2}=0.622$ , $k=0.41$ ,

$c_{w1}=\frac{c_{b1}}{k^2}+\frac{1+c_{b2}}{\sigma }$ , $c_{w2}=0.3$ , $c_{w3}=2$ and $c_{v1}=7.1$ . These constants were determined empirically.

The mean temperature field $\bar{\Theta }$ is governed by

$\frac{\partial \bar{\Theta }}{\partial t}+\frac{\partial \left({\bar{u}}_j\bar{\Theta }\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\alpha +{\alpha }_t\right)\left(\frac{\partial \bar{\Theta }}{\partial x_j}\right)\right],$ (33)

where ${\alpha }_t$ is the turbulent diffusion coefficient, being determined by

${\alpha }_t=\frac{{\nu }_t}{{\mathrm{Pr}}_t},$ (34)

where ${\text{Pr}}_t$ is the Prandtl number, which controls the magnitude of the mean turbulent diffusion. The turbulent thermal diffusivity value is 0.9 in all simulations.

5. Numerical Method

In the order to validate 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 $L_x=55d$ and width $L_y=30d$ . The heated square cylinder was placed at $x=16.5d$ and $y=15d$ , in order to minimize the influence of the boundaries. The grid resolution for these simulations was $318\times 164$ points, a non-uniform grid was used to better capture the effects with a total of 201 Lagrangian grid points inside the immersed body. The 2nd and 3rd equations in (1) were discretized 1) in space, using the method of centered finite difference, and 2) in time, by the second order Adams-Bashforth scheme, and solved in the two-dimensional domain together with the fractional step method for calculating the iterative time step and the coupling between pressure, velocity and temperature, respectively. The discretized expression for the pressure correction results in a linear system that was solved using the MSIP (Modified Strongly Implicit Procedure Method), developed in [17] , which consists of solving the following system of equations:

$\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+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$ (35)

$\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_j\partial x_j},$ (36)

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

${p^{\prime }}^{n+1}=p^{n+1}-p^n,$ (38)

$\begin{array}{l}\frac{{\theta }_{i\mathrm{,}j}^{n+1}-{\theta }_{i\mathrm{,}j}^n}{\Delta t}+\left[u_{i+\mathrm{1,}j}\frac{{\theta }_{i+\mathrm{1,}j}+{\theta }_{i\mathrm{,}j}}{2\Delta x}-u_{i\mathrm{,,}j}\frac{{\theta }_{i\mathrm{,}j}+{\theta }_{i-\mathrm{1,}j}}{2\Delta x}\right]\\ +\left[v_{i+\mathrm{1,}j}\frac{{\theta }_{i+\mathrm{1,}j}+{\theta }_{i\mathrm{,}j}}{2\Delta y}-v_{i\mathrm{,,}j}\frac{{\theta }_{i\mathrm{,}j}+{\theta }_{i-\mathrm{1,}j}}{2\Delta y}\right]\\ ={\text{Pe}}^{-1}(\left[\frac{{\theta }_{i+\mathrm{1,}j}-{\theta }_{i\mathrm{,}j}}{\Delta x}-\frac{{\theta }_{i\mathrm{,}j}-{\theta }_{i-j}}{\Delta x}\right]\frac{1}{\Delta x}\\ +\left[\frac{{\theta }_{i+\mathrm{1,}j}-{\theta }_{i\mathrm{,}j}}{\Delta y}-\frac{{\theta }_{i\mathrm{,}j}-{\theta }_{i-j}}{\Delta y}\right]\frac{1}{\Delta y})+q_{i\mathrm{,}j}\mathrm{,}\end{array}$ (39)

where the calculation of the force field on the interface, and the solution of the momentum and energy equations are performed explicitly. The computational domain is shown in the schematic illustration of Figure 4, which has length 55d and width 30d.

The grid is uniform in the region of the heated square cylinder, maintaining a minimum number of 30 grid inside. The time step used in the calculation process is in the range $1.0\times {10}^{-6}\text{s}$ to $1.0\times {10}^{-4}\text{s}$ , which is dynamically calculated by the Courant-Friedrichs-Lewy condition (CFL condition):

${\text{CFL}}_j=\frac{\text{max}{\left(\left|u\right|+c\right)}_{\left[x_{j-0.5}\mathrm{,}{x}_{j+0.5}\right]}\times \Delta t}{\Delta x_j}\mathrm{,}$ (40)

which is necessary for the explicit solution over time. To avoid numerical problems, the maximum ratio of 3% grid expansion was employed in these regions in both directions. In the present paper, the ratio between Lagrangian and Eulerian grid size was maintained at 0.98.

6. Results and Discussion

In this section, the heat-transfer by forced convection and turbulence problems are simulated using the IBM/VPM method. The results of the simulations are compared with previous results. Isothermal IBM/VPM conditions are implemented and studied. The Reynolds numbers are 40, 50, 100, 500, and 5000. The Prandtl number is fixed at 0.7 in all cases. The staggered grid arrangement is applied in all cases. The square cylinder D is 1, with the center positioned at 16.5d from the domain entry and at 15d from top domain wall, and the free-stream velocity $U_{\infty }$ . The computational domain is $30d\times 55d$ , with the center of the cylinder located at $\left(16.5d\mathrm{,16.5}d\right)$ . These dimensions were determined numerically in order to minimize the influence of the domain on the flow around the square cylinder and, at the same time, to minimize the total number of nodes required. Newman boundary conditions were used.

The boundary condition were prescribed as follow: 1) uniform flow with velocity $U_{\infty }$ , pressure $p=0$ , and temperature ${\Theta }_{\text{input}}=0$ at the entrance of the domain (left side); 2) flow fully developed in the output $\left(\partial u/\partial x=0,\partial \Theta //\partial x=0\right)$ ; 3) condition of symmetries at the upper and lower boundaries $\left(\partial u/\partial y=v=0\right)$ , $\left(\partial \Theta /\partial y=0\right)$ . Concerning the initial condition, $u=v=0$ , at $t=0$ (time), for the entire computational domain. The square cylinder is maintained at a constant dimensionless temperature equal to 1, i.e. ${\Theta }_c=1$ , while the fluid has an initial temperature equal to zero $\left({\Theta }_f=0\right)$ . Although the boundary condition at the output is not a reflexive condition, the simulation results successfully corroborate the formation of large (vortex) structures outside the domain without any reflection. The current pressure limit conditions at the inlet and outlet limits are imposed to be consistent with the equations for velocities due to the arrangement of the staggered grid.

6.1. Recirculation Bubble for Re = 40 and Re = 50

Figure 5 shows the streamline patterns around the square cylinder (bluff body) for different Reynolds numbers (Re = 40 and Re = 50). At low Reynolds number (Re = 40), the flow is laminar, steady, and slightly separated from the cylinder depicted. With the increase of the Reynolds number, Re = 50, the flow separates and two vortices are arranged symmetrically (according to the Figure 5(a), in the center line of the channel). For Re = 40, the fluid is imprisoned in the first 20 s. For Re = 50, the fluid is released from the recirculation bubble in the first 10 s.

6.2. Vorticity Fields for Re = 40 and Re = 50

Figure 6 shows an extended view of the vorticity fields for these simulations. Statistically established regimes are displayed. The formation and detachment of

the vortices for the higher Reynolds numbers is clearly observed, whereas, for Re = 40, the flow is stable.

6.3. Lift and Drag Coefficients for Re = 40 and Re = 50

The time evolution of the drag, $C_D$ , and lift, $C_l$ , coefficients are shown in Figure 7 or Re = 40 and Re = 50, respectively. It is verified that, in $C_D$ , does not occur amplitude variations, which is a consequence of the periodic detachment of vortices of the upper and lower surfaces. The lift coefficient is zero, because there is no buoyancy in the flow. In the same figure, $C_l$ , despite presenting variations in time, does not present significant oscillation; besides the flow be laminar, C_l is not influenced by the distribution of pressure on the back side of the cylinder, which would be strongly influenced by the formation of vortices on the upper and lower sides of the square cylinder.

6.4. Vorticity Fields for Re = 100 and Re = 500

In Figure 8(a), Figure 8(b) the vorticity fields are shown around a heated square cylinder, at different times, for Re = 100 and Re = 500, respectively. The square cylinder was maintained at a constant dimensionless temperature, i.e. with ${\theta }_b=1$ . (The sub-index b refers to the immersed body.)

In these figures, the vortex detachment is observed in both cases (time » 70). After 1.7 million iterations (»4 hours of CPU time), the flow is fully developed (at low computational cost). In Figure 8(b), for Re = 500, it can be observed the total detachment of vortices in the wake (here, the Smagorinsky model was used). It is also observed that, with the increase of the Reynolds number, occurs an increase of the drag coefficient. This happens because the vortices collide with more intensity behind the cylinder.

6.5. Lift and Drag Coefficients for Re = 100 and Re = 500

Contrary to what was observed in Figure 7(a) and Figure 7(b), for Reynolds numbers equal to 100 and 500, the time evolution of the drag coefficients, $C_D$ ,

and the lift coefficient, $C_l$ , in the square cylinder, varies in time, i.e. occur amplitude variations. This is a consequence of the periodic detachment of vortices on the upper and lower surfaces of the square cylinder. The lift coefficient, $C_l$ , also presents oscillation, that is, there is buoyancy in the flow. In Figure 9(b), the lift coefficient $C_l$ presents great oscillations, and this is not a numerical problem,

but a physical consequence caused by the turbulence. This type of flow is unstable and contains fluctuations that are dependent on time and position in space. These disturbances, shown in Figure 9(b), are strong. A complete description of the transitions requires the analysis of nonlinear perturbation amplification process. This is a difficult theoretical issue, since we are facing with nonlinear problems. These perturbations are indeed physical instabilities, and the transition from laminar to turbulent flow is verified. Therefore, it was necessary to implement a turbulence model. In this case Re = 500, the Smagorinsky model was employed. In addition, Figure 9 shows that, from the laminar (Re = 100) to turbulent (Re = 500) flow, occurs the influence of the pressure distribution downstream of the cylinder, which is strongly influenced by the formation and detachment of vortices on the upper and lower sides of the square cylinder. The fluctuation value of the supporting force is now directly connected to the formation and detachment of vortices, and, therefore, their value vary between maxima and minima of equal magnitudes.

6.6. Vorticity Fields for Re = 5000―Spalart-Allmaras and Smagorinsky Models

In Figure 10(a) and Figure 10(b), the Spalart-Allmaras and Smagorinsky turbulence models were used for the energy transfer process between the larger and the smaller turbulence scales. The kinetic energy of the physical instabilities accumulates over the cut-off frequency (given by the discretization loop).

6.7. Lift and Drag Coefficients for Re = 5000

In Figure 11, the Reynolds number is 5000. It is observed more pronounced oscillations in the drag coefficient, $C_D$ , and, therefore, the vortex generation and detachment process is totally swirling in the flow. The same happens with the lift coefficient, $C_l$ .

6.8. Mean Drag Coefficients

In Table 1 it is presented the results obtained for the mean value of the drag coefficients, for different Reynolds numbers, and compare them with the drag coefficients presented in [18] and [19] , which were obtained by the Boltzmann method and finite volume discretizations. We observe here that, when the Reynolds number increases, the drag coefficient also increases, particularly for vortex spill flows.

6.9. Mean Nusselt

The mean Nusselt numbers are presented in Table 2 and they are compared with the numerical values found in [20] . Several correlation can be obtained for the mean Nusselt numbers. Here, we used the Hilpert correlation that considers the global mean conditions, given by

$\overline{Nu}=c\times {\text{Re}}_D^m\times {\text{Pr}}^{1/3}\mathrm{.}$ (41)

7. Conclusions

The immersed boundary method coupled with the virtual physical model was used to simulate incompressible two-dimensional flows around a heated square cylinder at constant temperature on its surface. In order to validate the code, the works [18] [19] and [20] were used. A good numerical convergence was obtained, being the margin of error, with respect to these works, less than 3%. The time evolution of the drag and lift coefficients, as well as the Nusselt number were obtained with this methodology, being the parameters obtained from the Eulerian fields, since the geometry used in this work has singularities, which were taken into account in the construction of the algorithm/code. The implementation process for the calculation of the drag and lift coefficients is simple. This fact is important, because it allows its applicability to other (less simple) geometries.

In all simulations, the results show that the influence of the surface of the heated body immersed in the flow increases as the Reynolds number increases.

For the temporal discretization, the second order Adams-Bashforth scheme was used together with the spatial centered scheme. A turbulence model was also used for the energy transfer process between the largest and smallest turbulent scales.

For future work, it is planned to extend this work to more complex geometries, more precisely, to simulate tandem cylinder configurations in stationary and non-stationary cases. This will allow investigating mixed convection (natural and forced).

Acknowledgements

SG and RS acknowledge the partial support by CMUP (UID/MAT/ 00144/2013), which is funded by FCT (Portugal) with national (MCTES) and European structural funds (FEDER), under the partnership agreement PT2020-ext. to 2018. RS acknowledges the financial support by CAPES (Brazil). SG acknowledges the Project STRIDE-NORTE-01-0145-FEDER-000033, funded by ERDF NORTE 2020. The authors thank Virtual Hydrodynamics Laboratory of the Institute of Mechanical Engineering of the Federal University of Itajubá (Brazil) for the computer facilities.