# A Three-Dimensional, Immersed Boundary, Finite Volume Method for the Simulation of Incompressible Heat Transfer Flows around Complex Geometries

Authors: Hassan Badreddine; Yohei Sato; Matthias Berger; Bojan Ničeno
Journal: International Journal of Chemical Engineering
 (2017)

---

## Abstract

The current work focuses on the development and application of a new finite volume immersed boundary method (IBM) to simulate three-dimensional fluid flows and heat transfer around complex geometries. First, the discretization of the governing equations based on the second-order finite volume method on Cartesian, structured, staggered grid is outlined, followed by the description of modifications which have to be applied to the discretized system once a body is immersed into the grid. To validate the new approach, the heat conduction equation with a source term is solved inside a cavity with an immersed body. The approach is then tested for a natural convection flow in a square cavity with and without circular cylinder for different Rayleigh numbers. The results computed with the present approach compare very well with the benchmark solutions. As a next step in the validation procedure, the method is tested for Direct Numerical Simulation (DNS) of a turbulent flow around a surface-mounted matrix of cubes. The results computed with the present method compare very well with Laser Doppler Anemometry (LDA) measurements of the same case, showing that the method can be used for scale-resolving simulations of turbulence as well.


---

## Body

## 1. Introduction

Computational fluid dynamics (CFD) has reached a development level at which it is routinely applied in the industrial environment, at least for single-phase flows. Grid generation techniques and methods for discretization of governing equations all play important roles in successful application of CFD. Clearly, flow simulation is facilitated by the ever-increasing computational power on the one hand and the development of more efficient numerical methods on the other. Prominent examples of these algorithmic advances are the multigrid methods [1] and direct solvers based on Fast Fourier Transforms (FFTs). The most efficient numerical algorithms are usually defined for structured orthogonal grids, suitable only for problems occurring in simple, canonical geometries. Compared to unstructured grids, structured grids take less computational time and data storage.Flow simulations around complex geometries can be achieved in two ways: either by using unstructured codes already applicable for complex geometries or by tailoring the existing highly efficient and/or highly accurate structured methods to nontrivial geometries. The former approach has already been adopted by commercial CFD vendors and, as such, represents the state of the art. Indeed, leading commercial CFD solvers are based on unstructured grid methods and can be applied to complex geometries. The latter approach can go in two directions. One possibility is to define the underlying highly efficient/accurate method to boundary fitted grids [2, 3]. The other possibility is to leave the basic numerical method in its original form, that is, to discretize the system of governing equations using the original method on orthogonal grids and then to make the necessary changes to include the effects of a complex surface immersed in the computational domain. Such methods, usually referred to as immersed boundary methods (IBMs), can conveniently be divided into two broad classes (forcing andcut-cells methods), depending on the way the boundary conditions are imposed on the immersed body.In theforcing methods, one can find the continuous forcing approaches [4–7] where a continuous force is added to the governing equations before discretization, making the formulation of such approaches independent of the spatial discretization. On the other hand, in the discrete forcing approaches [8–12], the governing equations are first discretized on the grid neglecting the presence of the immersed body. Then interpolation schemes are used for the grid points at the vicinity of the immersed body in such a way that the required boundary condition at the interface is satisfied. Compared to continuous forcing approaches, discrete forcing approaches have the advantage of allowing direct control over stability and numerical accuracy.In the second class of IBM, thecut-cells methods, once the body is immersed in the domain, and the intersections between computational cells and immersed body are evaluated and the fraction of the cell volumes and cell face surfaces are computed. This procedure is required in order to modify the discretized system of equations for the truncated cells for the purpose of imposing the boundary conditions on the immersed body. Developments in this area include [13–16]. A comparative overview on IBM approaches is given in [17, 18].Several works exist in the literature regarding extending IBM approaches to solve the energy equation around complex geometries by imposing Dirichlet and Neumann boundary conditions [19–27]. Conjugate heat transfer simulations, where the energy equation is also solved inside the immersed body, using IBM were also reported in [28–31]. However, to the best of the authors’ knowledge, most of these works are based on the extension of the IBMforcing method to heat transfer. In this work, the implementation and development of a new IBM based on acut-cell approach to simulate heat transfer flow will be presented.


## 2. Numerical Method

A new cut-cell approach is described in this section and implemented in the in-house code PSI-BOIL [32]. In this work, we assume incompressibility of fluid and constant physical properties. The governing equations are the Navier-Stokes equations and the energy balance equation written here in integral formulation:(1)∮SudS=0,(2)∫V∂u∂tdV+∮SuudS=-1ρ∮SpdS+ν∮S∇udS+FBoussinesq,(3)∫V∂T∂tdV+∮SuTdS=λρcp∮S∇TdS.In the above equations, ρ is the density,u is the velocity vector, ν is the kinematic viscosity, T is the temperature, λ is the thermal conductivity, and cp is the specific heat capacity. FBoussinesq represents the Boussinesq approximation. V is the cell volume and S is the area vector of cell surface.The governing equations are integrated in time using a semi-implicit projection method [33]. For the spatial discretization, the orthogonal staggered finite volume method is used, with the second-order central difference being used for the diffusive terms and advection terms. Momentum equations are first solved for a tentative velocity field, u⋆, which is not divergence-free (i.e., ∇·u⋆≠0):(4)∫Vu⋆-unΔtdV+Hn-1/2=Pn+Dn+1/2+Fn.Here, H, D, and P represent the advective, viscous, and pressure terms in the momentum equations (see (4)), respectively. The superscript n represents current time step and n-1/2 indicates Adams-Bashforth time discretization scheme, whereas n+1/2 indicates an implicit time level. Once the tentative velocity field u⋆ is available, the algorithm proceeds by the solution of the pseudo-pressure p~ elliptic equation:(5)∮S∇p~dS=ρΔt∮Su⋆dS.The final step of the algorithm is projection of the velocity onto a divergence-free field un+1:(6)un+1=u⋆-1ρ∇p~Δt,(7)pn+1=pn+p~.As argued in [34], the pressure in this method converged only at a first-order rate. To converge at a second-order term, the pressure should be updated as follows:(8)pn+1=pn+p~-νΔt2∇2p~.Diagonally preconditioned conjugate gradient method is used to solve (4) and a variant of the additive correction multigrid method [35] is used to solve (5). However, it should be noted that multigrid methods, in certain circumstances, could be inefficient and require large computational time [22]. Other solvers such as FFT [36] could be also used.


### 2.1. Modification for an Immersed Body

#### 2.1.1. Cutting the Cells

In the first step of our immersed body approach, the immersed body is imported from a file in stereolithographic (STL) format. STL files define triangulated surfaces by the coordinates of the triangle’s vertices and normal surface vectors. The next step is to construct the intersection between the edges of the Cartesian grid cell with the triangles of the immersed body. The edge-triangle intersections are then used to define the cutting plane for cellsi,j,k (Figure 1).Figure 1
Procedure for immersing a body into a Cartesian grid: (a) intersecting the edges of celli,j,k with triangulated surface t1,t2; (b) creation of the cutting plane πi,j,k; and (c) definition of the final computational cell.
(a)
(b)
(c)Since we are using a staggered grid approach, we have four finite volume grids: one for the scalar variables and one each for the three velocity components. To cope with immersed boundaries, we cut all three grids separately, using the procedure outlined in Figure1, as illustrated in Figure 2.Figure 2
Each grid is cut separately: (a) scalar cell and (b), (c), and (d) momentum cells inx-, y-, and z-directions, respectively.
(a)
(b)
(c)
(d)


#### 2.1.2. Cell Classification

Once the scalar cells (pressure and temperature cells) are cut, they are classified into three sets:(i)
Cut cells, belonging to set C, are the cells which are intersected by the immersed boundary. The fluid volume fraction fV for these cells is 0<fV<1.(ii)
Fluid cells, belonging to set F, are the cells whose center remains in the fluid part of the domain, no matter if they are cut or not. This means that fV>0.5 for these cells.(iii)
Solid cells, belonging to set S, are all those which are not in set F, that is the cells whose center lies in the solid part of the domain.The classification of momentum (velocity) cells, whether they are ON (in the fluid) or OFF (in the immersed body), is related to the scalar cells. For instance, theu-component velocity cell is considered to be in the fluid if the two surrounding scalar cells in I-direction (i.e., cell (i) and cell (i-1) are both in fluid). A similar procedure is applied to the other two velocity components.Figure3(a) shows a 2D computational grid, with each scalar cell center being denoted with a black dot. Immersed body is represented by a shaded area. The cells whose cell center lies in fluid belong to set F, while the rest of the cells belong to set S as illustrated in Figure 3(b). The cells which are cut with the immersed body also belong to set C.Figure 3
A 2D computational domain with classified cells: (a) scalar cell centers and (b) sets. The purple shaded area represents immersed body.
(a)
(b)


#### 2.1.3. Inertial Terms

Cuts performed on the computational cells change their geometrical properties, and this has to be incorporated in the discretized system of equations. Modification of the inertial terms for the presence of an immersed body is straightforward. For a cell containing an immersed body, the inertial term, in terms of a generalized variableϕ, is multiplied by the volume fluid fraction of cell fVi,j,k:(9)∫Vϕ⋆-ϕnΔtdV≈I^i,j,k⋆-I^i,j,kn=ϕi,j,k⋆-ϕi,j,knVi,j,kfVi,j,kΔt.


#### 2.1.4. Advective Terms

Advective terms for the immersed body are affected by the changes in cell face areas. For each cell face which is cut, the fraction remaining in the fluid part of the cell is first estimated. As an example, the cell face ati+δ,j,k before and after immersing a body is illustrated in Figure 4.Figure 4
A cell face,i+δ,j,k before (a) and after (b) cutting. Cell face area before the cutting was Si+δ,j,k and after S^i+δ,j,k.
(a)
(b)The fraction of the cell face area remaining in the fluid is(10)fi+δ,j,kS=S^Si+δ,j,k.Using equations like (10) for all faces in the computational domain which are cut, the advective term for the general scalar variable ϕ for the immersed body can be evaluated as follows:(11)∮SuϕdS≈∑δ=-1/21/2H^i+δ,j,k+H^i,j+δ,k+H^i,j,k+δ,where H^ represents those advective terms modified for the presence of the immersed boundary, computed from(12)∀i,j,k∈F:H^i+δ,j,k=Hi+δ,j,kfi+δ,j,kSifi+2δ,j,k∈C,where Hi+δ,j,k=uϕSxi+δ,j,k represents the unmodified advective term. Analogous expressions could be easily formulated for the y- and z-directions.


#### 2.1.5. Diffusion Terms

The diffusion terms for the case of an immersed body can be written as follows:(13)∮S∇ϕdS≈∑δ=-1/21/2D^i+δ,j,k+D^i,j+δ,k+D^i,j,k+δ.For the evaluation of the diffusion terms, modified for the presence of the immersed body (D^), two different cases have to be considered. The first occurs if both cells surrounding the face for which the diffusive flux is estimated belong to the fluid set (Figure 5(a)). In such a case, the distance between the cells remains unchanged, but the cell face area between them changes, affecting the diffusive flux through the cell face, and is accounted for with the factor introduced in (10). The second case occurs if one cell, say i+2δ,j,k, is in the solid (Figure 5(b)). In such a case, the cell in the solid becomes effectively a boundary cell, and distance between the cells decreases. That is accounted for by defining the ratio between the new d^ and the old distances between the cells:(14)fi+δ,j,kd=d^di+δ,j,k.To summarize, corrections which need to be made for the diffusive fluxes are as follows:(15)∀i,j,k∈F:D^i+δ,j,k=Di+δ,j,kfi+δ,j,kSfi+δ,j,kdifi+2δ,j,k∈C,where Di+δ,j,k=(ϕi+2δ,j,k-ϕi,j,k)Sx/di+δ,j,k represents the unmodified diffusive flux. It should be noted that same procedure is also used to discretize the elliptic equation for the pseudo-pressure (see (5)).Figure 5
(a) Both cells are in the fluid. (b) Celli+2δ,j,k in the solid. The effective diffusive flux through the connection diminishes in proportion to the reduced face area and reduced distance between cells.
(a)
(b)The above treatment for the diffusion terms is valid for the momentum equations and for the heat equation when fixing the velocity or the temperature to a predefined value is required. Another approach for the diffusion term in the heat equation should be taken in case conduction is allowed inside the solid. For instance, when both cellsi,j,k and i+2δ,j,k are cut but their centers are in the fluid (similar to Figure 5(a)), the diffusive flux at the cell face i+δ,j,k is evaluated as follows:(16)D^i+δ,j,k=λfS^i+δ,j,k+λsSx-S^i+δ,j,kdi+δ,j,kϕi+2δ,j,k-ϕi,j,k,where λf and λs are the thermal conductivity of the fluid and solid, respectively. The second case is when the center of cell i,j,k is in the fluid and the center of i+2δ,j,k is in the solid (similar to Figure 5(b)). In this case, the diffusive flux at cell face i+δ,j,k is multiplied by(17)D^i+δ,j,k=Sxλfλsλfdi+δ,j,k-d^i+δ,j,k+λsd^i+δ,j,kϕi+2δ,j,k-ϕi,j,k.The above equation can be proved by considering the continuity of the diffusion flux at the interface between the solid and the fluid:(18)λfTI-Tijkd^i+δ,j,k=λsTi+2δ,j,k-TIdi+δ,j,k-d^i+δ,j,k,where TI is the temperature at the interface. The diffusion flux at cell face i+δ,j,k is(19)λ∂T∂xi+δ,j,k≈λfTI-Tijkd^i+δ,j,k.From (18), TI could be formulated in terms of Tijk and Ti+2δ,j,k and then inserted in (19) to get (17).The corrections for the presence of the immersed body, defined for the inertial term by (9), the advection term by (12), and the diffusive term by (15), are performed before compiling the final discretized system of equations. Once this is done, additional care has to be taken for the cells that belong to the solid set. Their coefficients in the discretized system have to be modified in a way which ensure that the values taken by the variables (temperature and velocity) in these cells remain equal to the predefined value (in case no conduction is solved in the solid).


#### 2.1.6. Velocity Projection

The discrete form of the velocity projection equation (6) is(20)ui+δ,j,kn+1=ui+δ,j,k⋆+Δtρp~i,j,k-p~i+2δ,j,kdi+δ,j,k.It should be noted that the velocity is only updated for the momentum cells which are inside the fluid (i.e., when the pressure scalar cellsi,j,k and i+2δ,j,k are both in the fluid).


## 2.1. Modification for an Immersed Body

### 2.1.1. Cutting the Cells

In the first step of our immersed body approach, the immersed body is imported from a file in stereolithographic (STL) format. STL files define triangulated surfaces by the coordinates of the triangle’s vertices and normal surface vectors. The next step is to construct the intersection between the edges of the Cartesian grid cell with the triangles of the immersed body. The edge-triangle intersections are then used to define the cutting plane for cellsi,j,k (Figure 1).Figure 1
Procedure for immersing a body into a Cartesian grid: (a) intersecting the edges of celli,j,k with triangulated surface t1,t2; (b) creation of the cutting plane πi,j,k; and (c) definition of the final computational cell.
(a)
(b)
(c)Since we are using a staggered grid approach, we have four finite volume grids: one for the scalar variables and one each for the three velocity components. To cope with immersed boundaries, we cut all three grids separately, using the procedure outlined in Figure1, as illustrated in Figure 2.Figure 2
Each grid is cut separately: (a) scalar cell and (b), (c), and (d) momentum cells inx-, y-, and z-directions, respectively.
(a)
(b)
(c)
(d)


### 2.1.2. Cell Classification

Once the scalar cells (pressure and temperature cells) are cut, they are classified into three sets:(i)
Cut cells, belonging to set C, are the cells which are intersected by the immersed boundary. The fluid volume fraction fV for these cells is 0<fV<1.(ii)
Fluid cells, belonging to set F, are the cells whose center remains in the fluid part of the domain, no matter if they are cut or not. This means that fV>0.5 for these cells.(iii)
Solid cells, belonging to set S, are all those which are not in set F, that is the cells whose center lies in the solid part of the domain.The classification of momentum (velocity) cells, whether they are ON (in the fluid) or OFF (in the immersed body), is related to the scalar cells. For instance, theu-component velocity cell is considered to be in the fluid if the two surrounding scalar cells in I-direction (i.e., cell (i) and cell (i-1) are both in fluid). A similar procedure is applied to the other two velocity components.Figure3(a) shows a 2D computational grid, with each scalar cell center being denoted with a black dot. Immersed body is represented by a shaded area. The cells whose cell center lies in fluid belong to set F, while the rest of the cells belong to set S as illustrated in Figure 3(b). The cells which are cut with the immersed body also belong to set C.Figure 3
A 2D computational domain with classified cells: (a) scalar cell centers and (b) sets. The purple shaded area represents immersed body.
(a)
(b)


### 2.1.3. Inertial Terms

Cuts performed on the computational cells change their geometrical properties, and this has to be incorporated in the discretized system of equations. Modification of the inertial terms for the presence of an immersed body is straightforward. For a cell containing an immersed body, the inertial term, in terms of a generalized variableϕ, is multiplied by the volume fluid fraction of cell fVi,j,k:(9)∫Vϕ⋆-ϕnΔtdV≈I^i,j,k⋆-I^i,j,kn=ϕi,j,k⋆-ϕi,j,knVi,j,kfVi,j,kΔt.


### 2.1.4. Advective Terms

Advective terms for the immersed body are affected by the changes in cell face areas. For each cell face which is cut, the fraction remaining in the fluid part of the cell is first estimated. As an example, the cell face ati+δ,j,k before and after immersing a body is illustrated in Figure 4.Figure 4
A cell face,i+δ,j,k before (a) and after (b) cutting. Cell face area before the cutting was Si+δ,j,k and after S^i+δ,j,k.
(a)
(b)The fraction of the cell face area remaining in the fluid is(10)fi+δ,j,kS=S^Si+δ,j,k.Using equations like (10) for all faces in the computational domain which are cut, the advective term for the general scalar variable ϕ for the immersed body can be evaluated as follows:(11)∮SuϕdS≈∑δ=-1/21/2H^i+δ,j,k+H^i,j+δ,k+H^i,j,k+δ,where H^ represents those advective terms modified for the presence of the immersed boundary, computed from(12)∀i,j,k∈F:H^i+δ,j,k=Hi+δ,j,kfi+δ,j,kSifi+2δ,j,k∈C,where Hi+δ,j,k=uϕSxi+δ,j,k represents the unmodified advective term. Analogous expressions could be easily formulated for the y- and z-directions.


### 2.1.5. Diffusion Terms

The diffusion terms for the case of an immersed body can be written as follows:(13)∮S∇ϕdS≈∑δ=-1/21/2D^i+δ,j,k+D^i,j+δ,k+D^i,j,k+δ.For the evaluation of the diffusion terms, modified for the presence of the immersed body (D^), two different cases have to be considered. The first occurs if both cells surrounding the face for which the diffusive flux is estimated belong to the fluid set (Figure 5(a)). In such a case, the distance between the cells remains unchanged, but the cell face area between them changes, affecting the diffusive flux through the cell face, and is accounted for with the factor introduced in (10). The second case occurs if one cell, say i+2δ,j,k, is in the solid (Figure 5(b)). In such a case, the cell in the solid becomes effectively a boundary cell, and distance between the cells decreases. That is accounted for by defining the ratio between the new d^ and the old distances between the cells:(14)fi+δ,j,kd=d^di+δ,j,k.To summarize, corrections which need to be made for the diffusive fluxes are as follows:(15)∀i,j,k∈F:D^i+δ,j,k=Di+δ,j,kfi+δ,j,kSfi+δ,j,kdifi+2δ,j,k∈C,where Di+δ,j,k=(ϕi+2δ,j,k-ϕi,j,k)Sx/di+δ,j,k represents the unmodified diffusive flux. It should be noted that same procedure is also used to discretize the elliptic equation for the pseudo-pressure (see (5)).Figure 5
(a) Both cells are in the fluid. (b) Celli+2δ,j,k in the solid. The effective diffusive flux through the connection diminishes in proportion to the reduced face area and reduced distance between cells.
(a)
(b)The above treatment for the diffusion terms is valid for the momentum equations and for the heat equation when fixing the velocity or the temperature to a predefined value is required. Another approach for the diffusion term in the heat equation should be taken in case conduction is allowed inside the solid. For instance, when both cellsi,j,k and i+2δ,j,k are cut but their centers are in the fluid (similar to Figure 5(a)), the diffusive flux at the cell face i+δ,j,k is evaluated as follows:(16)D^i+δ,j,k=λfS^i+δ,j,k+λsSx-S^i+δ,j,kdi+δ,j,kϕi+2δ,j,k-ϕi,j,k,where λf and λs are the thermal conductivity of the fluid and solid, respectively. The second case is when the center of cell i,j,k is in the fluid and the center of i+2δ,j,k is in the solid (similar to Figure 5(b)). In this case, the diffusive flux at cell face i+δ,j,k is multiplied by(17)D^i+δ,j,k=Sxλfλsλfdi+δ,j,k-d^i+δ,j,k+λsd^i+δ,j,kϕi+2δ,j,k-ϕi,j,k.The above equation can be proved by considering the continuity of the diffusion flux at the interface between the solid and the fluid:(18)λfTI-Tijkd^i+δ,j,k=λsTi+2δ,j,k-TIdi+δ,j,k-d^i+δ,j,k,where TI is the temperature at the interface. The diffusion flux at cell face i+δ,j,k is(19)λ∂T∂xi+δ,j,k≈λfTI-Tijkd^i+δ,j,k.From (18), TI could be formulated in terms of Tijk and Ti+2δ,j,k and then inserted in (19) to get (17).The corrections for the presence of the immersed body, defined for the inertial term by (9), the advection term by (12), and the diffusive term by (15), are performed before compiling the final discretized system of equations. Once this is done, additional care has to be taken for the cells that belong to the solid set. Their coefficients in the discretized system have to be modified in a way which ensure that the values taken by the variables (temperature and velocity) in these cells remain equal to the predefined value (in case no conduction is solved in the solid).


### 2.1.6. Velocity Projection

The discrete form of the velocity projection equation (6) is(20)ui+δ,j,kn+1=ui+δ,j,k⋆+Δtρp~i,j,k-p~i+2δ,j,kdi+δ,j,k.It should be noted that the velocity is only updated for the momentum cells which are inside the fluid (i.e., when the pressure scalar cellsi,j,k and i+2δ,j,k are both in the fluid).


## 2.1.1. Cutting the Cells

In the first step of our immersed body approach, the immersed body is imported from a file in stereolithographic (STL) format. STL files define triangulated surfaces by the coordinates of the triangle’s vertices and normal surface vectors. The next step is to construct the intersection between the edges of the Cartesian grid cell with the triangles of the immersed body. The edge-triangle intersections are then used to define the cutting plane for cellsi,j,k (Figure 1).Figure 1
Procedure for immersing a body into a Cartesian grid: (a) intersecting the edges of celli,j,k with triangulated surface t1,t2; (b) creation of the cutting plane πi,j,k; and (c) definition of the final computational cell.
(a)
(b)
(c)Since we are using a staggered grid approach, we have four finite volume grids: one for the scalar variables and one each for the three velocity components. To cope with immersed boundaries, we cut all three grids separately, using the procedure outlined in Figure1, as illustrated in Figure 2.Figure 2
Each grid is cut separately: (a) scalar cell and (b), (c), and (d) momentum cells inx-, y-, and z-directions, respectively.
(a)
(b)
(c)
(d)


## 2.1.2. Cell Classification

Once the scalar cells (pressure and temperature cells) are cut, they are classified into three sets:(i)
Cut cells, belonging to set C, are the cells which are intersected by the immersed boundary. The fluid volume fraction fV for these cells is 0<fV<1.(ii)
Fluid cells, belonging to set F, are the cells whose center remains in the fluid part of the domain, no matter if they are cut or not. This means that fV>0.5 for these cells.(iii)
Solid cells, belonging to set S, are all those which are not in set F, that is the cells whose center lies in the solid part of the domain.The classification of momentum (velocity) cells, whether they are ON (in the fluid) or OFF (in the immersed body), is related to the scalar cells. For instance, theu-component velocity cell is considered to be in the fluid if the two surrounding scalar cells in I-direction (i.e., cell (i) and cell (i-1) are both in fluid). A similar procedure is applied to the other two velocity components.Figure3(a) shows a 2D computational grid, with each scalar cell center being denoted with a black dot. Immersed body is represented by a shaded area. The cells whose cell center lies in fluid belong to set F, while the rest of the cells belong to set S as illustrated in Figure 3(b). The cells which are cut with the immersed body also belong to set C.Figure 3
A 2D computational domain with classified cells: (a) scalar cell centers and (b) sets. The purple shaded area represents immersed body.
(a)
(b)


## 2.1.3. Inertial Terms

Cuts performed on the computational cells change their geometrical properties, and this has to be incorporated in the discretized system of equations. Modification of the inertial terms for the presence of an immersed body is straightforward. For a cell containing an immersed body, the inertial term, in terms of a generalized variableϕ, is multiplied by the volume fluid fraction of cell fVi,j,k:(9)∫Vϕ⋆-ϕnΔtdV≈I^i,j,k⋆-I^i,j,kn=ϕi,j,k⋆-ϕi,j,knVi,j,kfVi,j,kΔt.


## 2.1.4. Advective Terms

Advective terms for the immersed body are affected by the changes in cell face areas. For each cell face which is cut, the fraction remaining in the fluid part of the cell is first estimated. As an example, the cell face ati+δ,j,k before and after immersing a body is illustrated in Figure 4.Figure 4
A cell face,i+δ,j,k before (a) and after (b) cutting. Cell face area before the cutting was Si+δ,j,k and after S^i+δ,j,k.
(a)
(b)The fraction of the cell face area remaining in the fluid is(10)fi+δ,j,kS=S^Si+δ,j,k.Using equations like (10) for all faces in the computational domain which are cut, the advective term for the general scalar variable ϕ for the immersed body can be evaluated as follows:(11)∮SuϕdS≈∑δ=-1/21/2H^i+δ,j,k+H^i,j+δ,k+H^i,j,k+δ,where H^ represents those advective terms modified for the presence of the immersed boundary, computed from(12)∀i,j,k∈F:H^i+δ,j,k=Hi+δ,j,kfi+δ,j,kSifi+2δ,j,k∈C,where Hi+δ,j,k=uϕSxi+δ,j,k represents the unmodified advective term. Analogous expressions could be easily formulated for the y- and z-directions.


## 2.1.5. Diffusion Terms

The diffusion terms for the case of an immersed body can be written as follows:(13)∮S∇ϕdS≈∑δ=-1/21/2D^i+δ,j,k+D^i,j+δ,k+D^i,j,k+δ.For the evaluation of the diffusion terms, modified for the presence of the immersed body (D^), two different cases have to be considered. The first occurs if both cells surrounding the face for which the diffusive flux is estimated belong to the fluid set (Figure 5(a)). In such a case, the distance between the cells remains unchanged, but the cell face area between them changes, affecting the diffusive flux through the cell face, and is accounted for with the factor introduced in (10). The second case occurs if one cell, say i+2δ,j,k, is in the solid (Figure 5(b)). In such a case, the cell in the solid becomes effectively a boundary cell, and distance between the cells decreases. That is accounted for by defining the ratio between the new d^ and the old distances between the cells:(14)fi+δ,j,kd=d^di+δ,j,k.To summarize, corrections which need to be made for the diffusive fluxes are as follows:(15)∀i,j,k∈F:D^i+δ,j,k=Di+δ,j,kfi+δ,j,kSfi+δ,j,kdifi+2δ,j,k∈C,where Di+δ,j,k=(ϕi+2δ,j,k-ϕi,j,k)Sx/di+δ,j,k represents the unmodified diffusive flux. It should be noted that same procedure is also used to discretize the elliptic equation for the pseudo-pressure (see (5)).Figure 5
(a) Both cells are in the fluid. (b) Celli+2δ,j,k in the solid. The effective diffusive flux through the connection diminishes in proportion to the reduced face area and reduced distance between cells.
(a)
(b)The above treatment for the diffusion terms is valid for the momentum equations and for the heat equation when fixing the velocity or the temperature to a predefined value is required. Another approach for the diffusion term in the heat equation should be taken in case conduction is allowed inside the solid. For instance, when both cellsi,j,k and i+2δ,j,k are cut but their centers are in the fluid (similar to Figure 5(a)), the diffusive flux at the cell face i+δ,j,k is evaluated as follows:(16)D^i+δ,j,k=λfS^i+δ,j,k+λsSx-S^i+δ,j,kdi+δ,j,kϕi+2δ,j,k-ϕi,j,k,where λf and λs are the thermal conductivity of the fluid and solid, respectively. The second case is when the center of cell i,j,k is in the fluid and the center of i+2δ,j,k is in the solid (similar to Figure 5(b)). In this case, the diffusive flux at cell face i+δ,j,k is multiplied by(17)D^i+δ,j,k=Sxλfλsλfdi+δ,j,k-d^i+δ,j,k+λsd^i+δ,j,kϕi+2δ,j,k-ϕi,j,k.The above equation can be proved by considering the continuity of the diffusion flux at the interface between the solid and the fluid:(18)λfTI-Tijkd^i+δ,j,k=λsTi+2δ,j,k-TIdi+δ,j,k-d^i+δ,j,k,where TI is the temperature at the interface. The diffusion flux at cell face i+δ,j,k is(19)λ∂T∂xi+δ,j,k≈λfTI-Tijkd^i+δ,j,k.From (18), TI could be formulated in terms of Tijk and Ti+2δ,j,k and then inserted in (19) to get (17).The corrections for the presence of the immersed body, defined for the inertial term by (9), the advection term by (12), and the diffusive term by (15), are performed before compiling the final discretized system of equations. Once this is done, additional care has to be taken for the cells that belong to the solid set. Their coefficients in the discretized system have to be modified in a way which ensure that the values taken by the variables (temperature and velocity) in these cells remain equal to the predefined value (in case no conduction is solved in the solid).


## 2.1.6. Velocity Projection

The discrete form of the velocity projection equation (6) is(20)ui+δ,j,kn+1=ui+δ,j,k⋆+Δtρp~i,j,k-p~i+2δ,j,kdi+δ,j,k.It should be noted that the velocity is only updated for the momentum cells which are inside the fluid (i.e., when the pressure scalar cellsi,j,k and i+2δ,j,k are both in the fluid).


## 3. Results

### 3.1. Steady-State Heat Conduction with a Source Term

As a first step to assess the performance of our immersed boundary method, we solve the heat conduction equation with a source term:(21)∫V∂T∂tdV=∮Sλ∇TdS+∫Vq˙dVW.The computational domain (Figure 6) is a rectangular channel of dimension 1.22 × 0.125 × 1 in the x-, y-, and z-directions, respectively. A rectangular immersed body of width 0.22 divides the domain into two parts: the x-coordinates of the “solid" domain range from -0.72 to -0.5, while the x-coordinates of the “solid-source" domain range form -0.5 to 0.5. The full domain is meshed with 64 × 2 × 64 grid cells. Given the fact that the solution of this case is two-dimensional, only two grid cells are used for the y-direction. The heat equation (see (21)) is solved with the source term set to 100W/m3 inside the “solid-source” domain and to zero inside the “solid” domain. The following boundary conditions are used for temperature T: (i)
T=289 [K] at x=-0.72 and T=300 at x=0.5(ii)
Periodic boundary condition in they-direction(iii)
∂T/∂n=0 in the z-direction.Figure 6
The computational domain used to solve the heat conduction equation. The green cells belong to the “solid-source” domain and the purple cells to the immersed body (“solid” domain).The thermal conductivityλ for the “solid-source” domain was set to 0.5. Three different cases are considered: the first one without the immersed body (i.e., the domain size in the x-direction ranges from -0.5 to 0.5 with T being set to 300 for the west and east boundaries). In this case, (21) has an analytical solution:(22)Tx=q˙2λx-0.5-x-0.52+300.For the other cases, the immersed body is included with λ(body) taking two values: one and two. For λ(body) = 1, the temperature at the interface between the fluid and body domains would be equal to 300. This will allow us to compare the results with the analytical data.The diffusion term was discretized in time using Crank-Nicolson scheme and in space using the second-order central difference scheme. Equation (21) is solved until steady-state solution is reached. Temperature profiles for all the cases are shown in Figure 7. Excellent agreement with the analytical solution is obtained for the case simulated without IBM and for λ(body) = 1. Inside the body zone, a linear profile for the temperature is observed as expected since the heat equation is solved without a source term there. This profile is also steeper for λ=1 compared to λ=2 as expected.Figure 7
Temperature profile inside the domain. The vertical black line corresponds to the interface between the “solid-source” and body domains.


### 3.2. Natural Convection in a Square Cavity

In this section, we validate our immersed boundary method for natural convection flow in square domain named “fluid.” The computational domain is similar to the one used in Section3.1. The only difference is that the width of the immersed body is 0.1 instead of 0.22. The fluid governing equations (see (1) and (2)) and the energy balance equation (see (3)) are solved inside the fluid domain and the temperature is fixed to Thot=301 K inside the body zone. The following boundary conditions are used:(i)
Tcold=300 K at x=0.5(ii)
∂T/∂n=0 for the top and bottom faces (z-direction)(iii)
Periodic boundary condition forT and u at the south and north faces (y-direction)(iv)
No-slip boundary condition foru at x- and z-directions,(v)
∂p~/∂n=0,where n is the normal vector at corresponding boundary. The whole domain is meshed with 128 × 2 × 128 grid cells. Only two cells are used for the y-direction because the solution is two-dimensional for this case. This problem can be characterized with two nondimensional numbers: the Rayleigh (Ra=gβΔTL3/να) and the Prandtl (Pr=ν/α) number. Four cases are simulated with different Rayleigh numbers (103 to 106). Prandtl number is set to 0.71 for all the cases. The governing equations are solved until steady state is reached.Results of averaged Nusselt number and maximum vertical dimensionless velocity on the horizontal mid-plane and its location and maximum horizontal dimensionless velocity on the vertical mid-plane and its location are compared against the benchmark data published in [37]. The dimensionless velocity is obtained by scaling the velocity with α/L, where α is the thermal diffusivity and L is the length of square cavity. The average Nusselt number Nu¯ is defined as follows:(23)Nu¯=∫hot  wall∂T∂xdz.Results are presented in Table1 and very good agreement is achieved compared to the benchmark data. The temperature profile for Ra=105 case is shown in Figure 8.Table 1
Comparison between the IBM and the benchmark (data) solutions for the natural convection simulation in a square cavity.
Ra
1
0
3
1
0
4
105
106
IBM
data
IBM
data
IBM
data
IBM
data
N
u
-
1.119
1.118
2.249
2.243
4.528
4.519
8.859
8.800
u
m
a
x
3.645
3.649
16.220
16.178
34.755
34.730
64.854
64.630
z
m
a
x
0.816
0.813
0.824
0.823
0.855
0.855
0.847
0.850
w
m
a
x
3.732
3.697
19.766
19.617
68.709
68.590
219.696
219.360
x
m
a
x
0.179
0.178
0.119
0.119
0.067
0.066
0.0418
0.0379Figure 8
Temperature profile for natural convection flow in square cavity with Ra =105 simulated using the IBM approach. The vertical black line corresponds to the interface between the immersed body and the fluid domain.


### 3.3. Natural Convection in a Square Enclosure with a Circular Cylinder

In this section, the simulation of natural convection in a square enclosure of length one with an immersed circular cylinder of radius0.2 centered at the center of the enclosure is presented. The temperature of the walls of the square enclosure is fixed to Tcold=300 K, while Thot=301 K is imposed on the cells that are inside the circular cylinder. Prandtl number was set to 0.71 and four different cases were simulated depending on Rayleigh number that ranges from 103 to 106.The domain (Figure9) is meshed with 256 cells in the x- and z-directions, while only two cells were used in y-directions, where periodic boundary condition was used, since the solution is two-dimensional (x-z plane).Figure 9
A section of the mesh used to simulate natural convection of heated cylinder. The green cells belong to the fluid domain and the purple cells to the immersed cylinder.The governing equations are solved until steady state is reached. The local Nusselt number (Nu) distribution along the surface of the square walls (Figure 10), isothermals (Figure 11), and streamlines (Figure 12) for all the cases were compared against results published in [25]. The local Nusselt number is defined as follows:(24)Nu=∂T∂nwall,where n is the wall-normal vector. As can be seen, the results obtained using our IBM approach are in very good agreement with the results of Kim et al. [25].Figure 10
Local Nusselt number distribution along the surface of the square enclosure for Ra =103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)Figure 11
Comparison of isothermals contours simulated by our IBM approach (coloured figures) and by Kim et al. [25] for Ra = 103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)Figure 12
Comparison of velocity streamlines coloured by temperature simulated by our IBM approach (cylinder shown in red) and by Kim et al. [25] for Ra = 103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)


### 3.4. Turbulent Flow around a Matrix of Cubes

In this section, we validate our immersed boundary method for isothermal turbulent flow around a surface-mounted matrix of cubes. As reference data, LDA measurements from [38] have been used. This case was used for 6th and 8th ERCOFTAC/IAHR/COST Workshops on refined flow and turbulence modeling [39, 40] and has been analyzed by a number of groups. The physics of the flow is well understood and will not be explained here in great detail. Our focus in this work is to compare computed mean and fluctuating velocity profiles against the measured data.The representation of the problem domain including all relevant dimensions is sketched in Figure13. An equidistantly spaced matrix of cubes is mounted on the lower surface a rectangular channel. The channel height is H=51 mm, height of the cubes is h=15 mm, and the pitch between the cubes is S=60 mm. The experimental setup consists of 25×10 cubes in the streamwise and spanwise directions, respectively. Measurements were taken around the 18th row of cubes, where the flow was expected to be spatially fully developed. This aspect makes this test case particularly attractive, as it allows analysis to be carried out in a periodic box, avoiding the specification of inlet conditions. The dashed lines in Figure 13 depict the computational domain. Dimensions are S×S×H in the streamwise (x), spanwise (y), and wall-normal directions (z), respectively. The origin of the coordinate system is located on the lower surface, centrally with respect to the computational domain and the cube under consideration, as shown in Figure 13.Figure 13
Schematic representation of the matrix of cubes test case with all relevant dimensions. The computational box is represented by the dashed lines, and coordinate system origin is indicated.The working fluid is air, with a bulk velocity ofUb=3.86 m/s, resulting in a Re based on cube height of Reh=3854. This is a modest Re number, allowing simulations to be performed without the use of turbulence model, that is, as a DNS. Grid resolution is set to 128×128×108 cells in x-, y-, and z-directions, respectively. The grid is uniform in all directions; that is, no stretching towards the wall was used in any coordinate direction. Cell dimensions in x- and y-directions are Δx=Δy=0.46875 mm, and the number of cells in normal direction (z) was chosen in such a way as to make Δz as close as possible to Δx and Δy; in the present case, Δz=0.472˙. The total number of cells was close to 1.8 million. For stability reason, the time step was set to 2×10-5, resulting in a maximum Courant-Friedrich-Levy (CFL) number around the value of 0.3. The total number of time steps performed in this simulation was 105, but only the last 8×104 were used for computation of turbulent statistics.The computed mean and fluctuating streamwise and spanwise velocity profiles are compared with measurements from [38] in Figures 14–16. Figure 14 shows profiles in the vertical (top row) and horizontal (bottom row) mid-plane of the cube at x=-12 mm. Figures 15 and 16 show the same, but for x=3 mm and x=12 mm, respectively. Computed mean velocity profiles match measured values well at all three measuring locations. The same can also be claimed for the streamwise velocity fluctuations. The most striking agreement is probably the one in the vertical plane at x=-3 mm, shown in Figure 15(b) at the top. The sharp peak in the streamwise velocity fluctuations is captured exactly and is very close for the spanwise component as well. Overall, the agreement is good, and we can conclude that the immersed boundary approach we describe can be used for large-scale simulations of turbulence, at least for the modest Re number considered in this work. With increasing Re, a subgrid scale (SGS) model would need to be used, as well as wall functions. However, implementation of state-of-the-art SGS models on a Cartesian grid is a straightforward task, and this is an added benefit of the current approach.Figure 14
Comparison of computed and measured velocity atx=-12 mm: (a) position of the measuring plane and line for LDA, (b) mean streamwise, (c) fluctuating streamwise, and (d) fluctuating spanwise. The top row shows comparison in vertical (x-z) plane and the bottom row in horizontal (x-y) plane. Code predictions are given by the solid line (DNS) and the measurements by open circles (EXP).
(a)
(b)
(c)
(d)Figure 15
Same as in Figure14 but for x=-3 mm.
(a)
(b)
(c)
(d)Figure 16
Same as in Figure14 but for x=12 mm.
(a)
(b)
(c)
(d)


## 3.1. Steady-State Heat Conduction with a Source Term

As a first step to assess the performance of our immersed boundary method, we solve the heat conduction equation with a source term:(21)∫V∂T∂tdV=∮Sλ∇TdS+∫Vq˙dVW.The computational domain (Figure 6) is a rectangular channel of dimension 1.22 × 0.125 × 1 in the x-, y-, and z-directions, respectively. A rectangular immersed body of width 0.22 divides the domain into two parts: the x-coordinates of the “solid" domain range from -0.72 to -0.5, while the x-coordinates of the “solid-source" domain range form -0.5 to 0.5. The full domain is meshed with 64 × 2 × 64 grid cells. Given the fact that the solution of this case is two-dimensional, only two grid cells are used for the y-direction. The heat equation (see (21)) is solved with the source term set to 100W/m3 inside the “solid-source” domain and to zero inside the “solid” domain. The following boundary conditions are used for temperature T: (i)
T=289 [K] at x=-0.72 and T=300 at x=0.5(ii)
Periodic boundary condition in they-direction(iii)
∂T/∂n=0 in the z-direction.Figure 6
The computational domain used to solve the heat conduction equation. The green cells belong to the “solid-source” domain and the purple cells to the immersed body (“solid” domain).The thermal conductivityλ for the “solid-source” domain was set to 0.5. Three different cases are considered: the first one without the immersed body (i.e., the domain size in the x-direction ranges from -0.5 to 0.5 with T being set to 300 for the west and east boundaries). In this case, (21) has an analytical solution:(22)Tx=q˙2λx-0.5-x-0.52+300.For the other cases, the immersed body is included with λ(body) taking two values: one and two. For λ(body) = 1, the temperature at the interface between the fluid and body domains would be equal to 300. This will allow us to compare the results with the analytical data.The diffusion term was discretized in time using Crank-Nicolson scheme and in space using the second-order central difference scheme. Equation (21) is solved until steady-state solution is reached. Temperature profiles for all the cases are shown in Figure 7. Excellent agreement with the analytical solution is obtained for the case simulated without IBM and for λ(body) = 1. Inside the body zone, a linear profile for the temperature is observed as expected since the heat equation is solved without a source term there. This profile is also steeper for λ=1 compared to λ=2 as expected.Figure 7
Temperature profile inside the domain. The vertical black line corresponds to the interface between the “solid-source” and body domains.


## 3.2. Natural Convection in a Square Cavity

In this section, we validate our immersed boundary method for natural convection flow in square domain named “fluid.” The computational domain is similar to the one used in Section3.1. The only difference is that the width of the immersed body is 0.1 instead of 0.22. The fluid governing equations (see (1) and (2)) and the energy balance equation (see (3)) are solved inside the fluid domain and the temperature is fixed to Thot=301 K inside the body zone. The following boundary conditions are used:(i)
Tcold=300 K at x=0.5(ii)
∂T/∂n=0 for the top and bottom faces (z-direction)(iii)
Periodic boundary condition forT and u at the south and north faces (y-direction)(iv)
No-slip boundary condition foru at x- and z-directions,(v)
∂p~/∂n=0,where n is the normal vector at corresponding boundary. The whole domain is meshed with 128 × 2 × 128 grid cells. Only two cells are used for the y-direction because the solution is two-dimensional for this case. This problem can be characterized with two nondimensional numbers: the Rayleigh (Ra=gβΔTL3/να) and the Prandtl (Pr=ν/α) number. Four cases are simulated with different Rayleigh numbers (103 to 106). Prandtl number is set to 0.71 for all the cases. The governing equations are solved until steady state is reached.Results of averaged Nusselt number and maximum vertical dimensionless velocity on the horizontal mid-plane and its location and maximum horizontal dimensionless velocity on the vertical mid-plane and its location are compared against the benchmark data published in [37]. The dimensionless velocity is obtained by scaling the velocity with α/L, where α is the thermal diffusivity and L is the length of square cavity. The average Nusselt number Nu¯ is defined as follows:(23)Nu¯=∫hot  wall∂T∂xdz.Results are presented in Table1 and very good agreement is achieved compared to the benchmark data. The temperature profile for Ra=105 case is shown in Figure 8.Table 1
Comparison between the IBM and the benchmark (data) solutions for the natural convection simulation in a square cavity.
Ra
1
0
3
1
0
4
105
106
IBM
data
IBM
data
IBM
data
IBM
data
N
u
-
1.119
1.118
2.249
2.243
4.528
4.519
8.859
8.800
u
m
a
x
3.645
3.649
16.220
16.178
34.755
34.730
64.854
64.630
z
m
a
x
0.816
0.813
0.824
0.823
0.855
0.855
0.847
0.850
w
m
a
x
3.732
3.697
19.766
19.617
68.709
68.590
219.696
219.360
x
m
a
x
0.179
0.178
0.119
0.119
0.067
0.066
0.0418
0.0379Figure 8
Temperature profile for natural convection flow in square cavity with Ra =105 simulated using the IBM approach. The vertical black line corresponds to the interface between the immersed body and the fluid domain.


## 3.3. Natural Convection in a Square Enclosure with a Circular Cylinder

In this section, the simulation of natural convection in a square enclosure of length one with an immersed circular cylinder of radius0.2 centered at the center of the enclosure is presented. The temperature of the walls of the square enclosure is fixed to Tcold=300 K, while Thot=301 K is imposed on the cells that are inside the circular cylinder. Prandtl number was set to 0.71 and four different cases were simulated depending on Rayleigh number that ranges from 103 to 106.The domain (Figure9) is meshed with 256 cells in the x- and z-directions, while only two cells were used in y-directions, where periodic boundary condition was used, since the solution is two-dimensional (x-z plane).Figure 9
A section of the mesh used to simulate natural convection of heated cylinder. The green cells belong to the fluid domain and the purple cells to the immersed cylinder.The governing equations are solved until steady state is reached. The local Nusselt number (Nu) distribution along the surface of the square walls (Figure 10), isothermals (Figure 11), and streamlines (Figure 12) for all the cases were compared against results published in [25]. The local Nusselt number is defined as follows:(24)Nu=∂T∂nwall,where n is the wall-normal vector. As can be seen, the results obtained using our IBM approach are in very good agreement with the results of Kim et al. [25].Figure 10
Local Nusselt number distribution along the surface of the square enclosure for Ra =103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)Figure 11
Comparison of isothermals contours simulated by our IBM approach (coloured figures) and by Kim et al. [25] for Ra = 103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)Figure 12
Comparison of velocity streamlines coloured by temperature simulated by our IBM approach (cylinder shown in red) and by Kim et al. [25] for Ra = 103 (a), Ra = 104 (b), Ra = 105 (c), and Ra = 106 (d).
(a)
(b)
(c)
(d)


## 3.4. Turbulent Flow around a Matrix of Cubes

In this section, we validate our immersed boundary method for isothermal turbulent flow around a surface-mounted matrix of cubes. As reference data, LDA measurements from [38] have been used. This case was used for 6th and 8th ERCOFTAC/IAHR/COST Workshops on refined flow and turbulence modeling [39, 40] and has been analyzed by a number of groups. The physics of the flow is well understood and will not be explained here in great detail. Our focus in this work is to compare computed mean and fluctuating velocity profiles against the measured data.The representation of the problem domain including all relevant dimensions is sketched in Figure13. An equidistantly spaced matrix of cubes is mounted on the lower surface a rectangular channel. The channel height is H=51 mm, height of the cubes is h=15 mm, and the pitch between the cubes is S=60 mm. The experimental setup consists of 25×10 cubes in the streamwise and spanwise directions, respectively. Measurements were taken around the 18th row of cubes, where the flow was expected to be spatially fully developed. This aspect makes this test case particularly attractive, as it allows analysis to be carried out in a periodic box, avoiding the specification of inlet conditions. The dashed lines in Figure 13 depict the computational domain. Dimensions are S×S×H in the streamwise (x), spanwise (y), and wall-normal directions (z), respectively. The origin of the coordinate system is located on the lower surface, centrally with respect to the computational domain and the cube under consideration, as shown in Figure 13.Figure 13
Schematic representation of the matrix of cubes test case with all relevant dimensions. The computational box is represented by the dashed lines, and coordinate system origin is indicated.The working fluid is air, with a bulk velocity ofUb=3.86 m/s, resulting in a Re based on cube height of Reh=3854. This is a modest Re number, allowing simulations to be performed without the use of turbulence model, that is, as a DNS. Grid resolution is set to 128×128×108 cells in x-, y-, and z-directions, respectively. The grid is uniform in all directions; that is, no stretching towards the wall was used in any coordinate direction. Cell dimensions in x- and y-directions are Δx=Δy=0.46875 mm, and the number of cells in normal direction (z) was chosen in such a way as to make Δz as close as possible to Δx and Δy; in the present case, Δz=0.472˙. The total number of cells was close to 1.8 million. For stability reason, the time step was set to 2×10-5, resulting in a maximum Courant-Friedrich-Levy (CFL) number around the value of 0.3. The total number of time steps performed in this simulation was 105, but only the last 8×104 were used for computation of turbulent statistics.The computed mean and fluctuating streamwise and spanwise velocity profiles are compared with measurements from [38] in Figures 14–16. Figure 14 shows profiles in the vertical (top row) and horizontal (bottom row) mid-plane of the cube at x=-12 mm. Figures 15 and 16 show the same, but for x=3 mm and x=12 mm, respectively. Computed mean velocity profiles match measured values well at all three measuring locations. The same can also be claimed for the streamwise velocity fluctuations. The most striking agreement is probably the one in the vertical plane at x=-3 mm, shown in Figure 15(b) at the top. The sharp peak in the streamwise velocity fluctuations is captured exactly and is very close for the spanwise component as well. Overall, the agreement is good, and we can conclude that the immersed boundary approach we describe can be used for large-scale simulations of turbulence, at least for the modest Re number considered in this work. With increasing Re, a subgrid scale (SGS) model would need to be used, as well as wall functions. However, implementation of state-of-the-art SGS models on a Cartesian grid is a straightforward task, and this is an added benefit of the current approach.Figure 14
Comparison of computed and measured velocity atx=-12 mm: (a) position of the measuring plane and line for LDA, (b) mean streamwise, (c) fluctuating streamwise, and (d) fluctuating spanwise. The top row shows comparison in vertical (x-z) plane and the bottom row in horizontal (x-y) plane. Code predictions are given by the solid line (DNS) and the measurements by open circles (EXP).
(a)
(b)
(c)
(d)Figure 15
Same as in Figure14 but for x=-3 mm.
(a)
(b)
(c)
(d)Figure 16
Same as in Figure14 but for x=12 mm.
(a)
(b)
(c)
(d)


## 4. Conclusions

In this paper, a new immersed boundary method (IBM) based on thecut-cell approach to simulate three-dimensional, incompressible flows with heat transfer has been described. The governing equations are first discretized using the Cartesian, staggered, finite volume method. The resulting systems of linear equations are then corrected for the presence of the immersed body by recomputing individual entries where needed. The simple structure of the Cartesian grid facilitated usage of an efficient algebraic multigrid solver, but special care had to been taken with the modifications performed on the discretized system of equations not to hinder its convergence properties.The verification tests comprise steady-state heat conduction equation with a source term solved inside a cavity coupled with the solution of the heat conduction equation without source term inside the immersed body. The parabolic temperature profile obtained inside the cavity compares very well to the analytical solution, and temperature profile inside the immersed body was as expected linear. Next simulations of natural convection in cavity with and without immersed cylinder for different Rayleigh numbers were presented. Results in terms of Nusselt number distribution, isothermal contours, and velocity streamlines were in excellent agreement with results published in [25, 37]. Finally, the new IBM approach has been also validated in the turbulent flow regime, using Direct Numerical Simulation (DNS) of the flow around a surface-mounted matrix of cubes at Re=3854. The computed results for velocity and velocity fluctuations are in good agreement with experimental Laser Doppler Anemometry (LDA) measurements, showing that the proposed method can be used for large-scale simulations of turbulence.The general benefits of the IBM approach over body-fitted grid methods are also shared by the method proposed in this work. The structure of the Cartesian grid allows the usage of efficient linear solvers (multigrid) and easier optimization of the computer code for various computational platforms. Furthermore, the implementation of state-of-the-art physical models is a straightforward task on a Cartesian grid. Finally, using an IBM rids the analyst from the burden of a body-fitted mesh generation procedure.



---

*Source: 1726519-2017-06-19.xml*