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.
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 and cut-cells methods), depending on the way the boundary conditions are imposed on the immersed body.
In the forcing 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, the cut-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 IBM forcing method to heat transfer. In this work, the implementation and development of a new IBM based on a cut-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:In the above equations, is the density, u is the velocity vector, is the kinematic viscosity, is the temperature, is the thermal conductivity, and is the specific heat capacity. represents the Boussinesq approximation. is the cell volume and 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, , which is not divergence-free (i.e., ):Here, , , and represent the advective, viscous, and pressure terms in the momentum equations (see (4)), respectively. The superscript represents current time step and indicates Adams-Bashforth time discretization scheme, whereas indicates an implicit time level. Once the tentative velocity field is available, the algorithm proceeds by the solution of the pseudo-pressure elliptic equation:The final step of the algorithm is projection of the velocity onto a divergence-free field :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:
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 cells (Figure 1).
(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 Figure 1, as illustrated in Figure 2.
(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 , are the cells which are intersected by the immersed boundary. The fluid volume fraction for these cells is .(ii)Fluid cells, belonging to set , are the cells whose center remains in the fluid part of the domain, no matter if they are cut or not. This means that for these cells.(iii)Solid cells, belonging to set , are all those which are not in set , 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, the -component velocity cell is considered to be in the fluid if the two surrounding scalar cells in -direction (i.e., cell (i) and cell (-1) are both in fluid). A similar procedure is applied to the other two velocity components.
Figure 3(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 , while the rest of the cells belong to set as illustrated in Figure 3(b). The cells which are cut with the immersed body also belong to set .
(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 :
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 at before and after immersing a body is illustrated in Figure 4.
(a)
(b)
The fraction of the cell face area remaining in the fluid isUsing 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:where represents those advective terms modified for the presence of the immersed boundary, computed fromwhere represents the unmodified advective term. Analogous expressions could be easily formulated for the - and -directions.
2.1.5. Diffusion Terms
The diffusion terms for the case of an immersed body can be written as follows:
For the evaluation of the diffusion terms, modified for the presence of the immersed body (), 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 , 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 and the old distances between the cells:To summarize, corrections which need to be made for the diffusive fluxes are as follows:where 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)).
(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 cells and are cut but their centers are in the fluid (similar to Figure 5(a)), the diffusive flux at the cell face is evaluated as follows:where and are the thermal conductivity of the fluid and solid, respectively. The second case is when the center of cell is in the fluid and the center of is in the solid (similar to Figure 5(b)). In this case, the diffusive flux at cell face is multiplied byThe above equation can be proved by considering the continuity of the diffusion flux at the interface between the solid and the fluid:where is the temperature at the interface. The diffusion flux at cell face isFrom (18), could be formulated in terms of and 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
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 cells and 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:The computational domain (Figure 6) is a rectangular channel of dimension 1.22 0.125 1 in the -, -, and -directions, respectively. A rectangular immersed body of width 0.22 divides the domain into two parts: the -coordinates of the “solid" domain range from to , while the -coordinates of the “solid-source" domain range form to . 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 -direction. The heat equation (see (21)) is solved with the source term set to inside the “solid-source” domain and to zero inside the “solid” domain. The following boundary conditions are used for temperature : (i) [K] at and at (ii)Periodic boundary condition in the -direction(iii) in the -direction.
The thermal conductivity for the “solid-source” domain was set to . Three different cases are considered: the first one without the immersed body (i.e., the domain size in the -direction ranges from to with being set to for the west and east boundaries). In this case, (21) has an analytical solution: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 compared to as expected.
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 Section 3.1. The only difference is that the width of the immersed body is instead of . 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 K inside the body zone. The following boundary conditions are used:(i) K at (ii) for the top and bottom faces (-direction)(iii)Periodic boundary condition for and at the south and north faces (-direction)(iv)No-slip boundary condition for at - and -directions,(v),where 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 -direction because the solution is two-dimensional for this case. This problem can be characterized with two nondimensional numbers: the Rayleigh () and the Prandtl () number. Four cases are simulated with different Rayleigh numbers ( to ). Prandtl number is set to 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 , where is the thermal diffusivity and is the length of square cavity. The average Nusselt number is defined as follows:
Results are presented in Table 1 and very good agreement is achieved compared to the benchmark data. The temperature profile for case is shown in Figure 8.
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 radius centered at the center of the enclosure is presented. The temperature of the walls of the square enclosure is fixed to K, while K is imposed on the cells that are inside the circular cylinder. Prandtl number was set to and four different cases were simulated depending on Rayleigh number that ranges from to .
The domain (Figure 9) is meshed with cells in the - and -directions, while only two cells were used in -directions, where periodic boundary condition was used, since the solution is two-dimensional (- plane).
The governing equations are solved until steady state is reached. The local Nusselt number () 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:where 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].
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(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 Figure 13. An equidistantly spaced matrix of cubes is mounted on the lower surface a rectangular channel. The channel height is mm, height of the cubes is mm, and the pitch between the cubes is mm. The experimental setup consists of 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 in the streamwise (), spanwise (), and wall-normal directions (), 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.
The working fluid is air, with a bulk velocity of m/s, resulting in a Re based on cube height of . 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 cells in -, -, and -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 - and -directions are mm, and the number of cells in normal direction () was chosen in such a way as to make as close as possible to and ; in the present case, . The total number of cells was close to 1.8 million. For stability reason, the time step was set to , resulting in a maximum Courant-Friedrich-Levy (CFL) number around the value of . The total number of time steps performed in this simulation was , but only the last 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 mm. Figures 15 and 16 show the same, but for mm and 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 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4. Conclusions
In this paper, a new immersed boundary method (IBM) based on the cut-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 . 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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.