Flow Turbulence Combust (2008) 80:187–206 DOI 10.1007/s10494-007-9108-0
An Immersed Boundary Method for Complex Flow and Heat Transfer F. Paravento · M. J. Pourquie · B. J. Boersma
Received: 1 December 2006 / Accepted: 10 September 2007 / Published online: 12 October 2007 © Springer Science + Business Media B.V. 2007
Abstract The need to predict flow and heat transfer problems requires a flexible and fast tool able to simulate complex geometries without increasing the complexity of the flow solver architecture. Here we use a finite volume code that uses a direct solver with pressure correction. A new immersed boundary method (IBM) is used for a geometry consisting of a square body in a flow. The method is applied to flow cases with and without heat transfer. The obstacle simulated in the domain is implemented by local forcing of the flow with a procedure that adjusts locally the shear stress at the position of the object in conjunction with a non-penetration condition on the body walls. This approach has already been successfully applied by Breugem and Boersma (Phys. Fluids 17:15, 2005). We extend it for the case of heat transfer between body and flow. Comparison with other methods has been carried out as well. However, the proposed method can not be simply extended to immersed boundaries not aligned with the grid. Keywords Immersed boundary · Complex geometry · Low Mach number approximation
1 Introduction Most of the applications of computational fluid dynamics, heat transfer and combustion require modeling of complex geometries. Conventional numerical models generally use a complex (non-orthogonal) grid structure which requires a substantial computational effort. In order to retain the advantages of numerical accuracy
F. Paravento (B) · M. J. Pourquie · B. J. Boersma Laboratory for Aero and Hydrodynamics, Delft University of Technology, Leeghwaterstraat 21, 2628 CA Delft, The Netherlands e-mail:
[email protected],
[email protected]
188
Flow Turbulence Combust (2008) 80:187–206
and computational efficiency associated with simple orthogonal grids an immersed boundary method (IBM) can be used. The IBM can simulate the necessary geometry by adding extra forces to the momentum equations locally within the grid. In this way the grid does not have to follow the geometry. The term ‘immersed boundary method’ was first used by Peskin [2] to simulate cardiac mechanics and associated blood flow. However, this method presents severe limitations with respect to numerical stability. Since Peskin introduced his novel procedure, numerous modifications have been proposed and a number of variants of the IBM approach now exist [11]. The first successful improvements to Peskin’s approach, for the case of solid wall treatment, are ascribed to Mohd–Yusof [3], Verzicco et al. [4] and Fadlun et al. [9]. In these last papers, in order to apply the forcing to the flow equations, only local information is needed for the computation of the body force instead of the complete force distribution over the boundary as in the Peskin approach. In Fadlun and Verzicco’s methodology, [9], this is equivalent to the use of a local interpolation of the velocity directly on the walls and apply no-slip and nonpenetration conditions. Recently, a new method has been used by Breugem and Boersma [10] for simulation of porous media. In this geometry a grid of cubes mimics a permeable wall with a certain value of porosity. The cubes are aligned with the Cartesian mesh and this allows to apply exactly the forcing at their boundaries. In this method the shear stress on the boundary of the simulated obstacles is replaced in such a way that the no-slip velocity condition for the tangential component is applied at the wall. In conjunction, a non-penetration condition is also applied for the perpendicular velocity components at the boundary. The main contribution of the present paper is the extension of the Breugem and Boersma method for heat transfer with a low Mach number model and the study of the capability of this approach to keep the internal region of an obstacle well isolated under different conditions. The method is implemented on a simple equidistant grid where a square obstacle is simulated. A preliminary comparison with the data of Franke et al. [16] is presented and for the same geometry the L2 norm of the vertical velocity component is computed as well. Successively, a two dimensional comparison with the Fadlun and Verzicco approach is carried out, as well as with a conventional method (iterative solver). Afterwards, we test the method for a three dimensional laminar incompressible case with varying inflow velocity. We also perform a three dimensional test for a varying density case with varying heat flux from the cube. In this last case, although combustion is not considered here, we let the temperature vary between the common limits of a premixed combustion flame from the ambient temperature to the adiabatic one [6]. This study is intended as a preliminary step before the application of this new IBM to reactive flow cases including combustion. In fact, from preliminary tests on the interaction between a flame front and an obstacle we noted that it is important that the velocity components inside the body remain small to avoid non-physical oscillations of velocity, temperature and density in the surroundings of the body. Therefore, for all tests, we have also looked at the velocity profile inside the body as a measure for the quality of the method. The paper is organized in three main sections. In the first section the Low Mach number system of flow equations is introduced and the theoretical model for the IBM
Flow Turbulence Combust (2008) 80:187–206
189
is presented. The main characteristics of the numerical discretization are summarized as well. In the second section results are shown and discussed. They suggest that this method can be a good candidate for the simulation of flows with or without heat transfer and for combustion cases as well. The last section contains some considerations about the advantages and limitations of the proposed method and the conclusions.
2 Formulation and Numerical Method Many heat transfer phenomena and reacting flows, such as for instance the burning of natural gas, occur at low Mach numbers. In order to study these cases with direct numerical simulation (DNS), in which all flow scales are resolved, we need a system of equations that allows for large heat release, large temperature and density variations and substantial interaction of the reacting interface with the hydrodynamic flow field, including the effects of turbulence. A low Mach number approximation is suitable to characterize most of deflagration cases with appreciable advantages regarding the computational cost, because we do not have to resolve the acoustic oscillations and the set of equations is similar to the incompressible case, although the density may vary due to heat release [5]. After an asymptotic analysis, used for the first time by Rehm and Baum [1], the system of equations can be written in nondimensional form. In the following reference quantities are denoted by the subscript ‘0 ’, while the superscript ‘∗ ’ denotes dimensional quantities: ρ=
ρ∗ p∗ u∗ T ∗ − T0 ,p= ,u = ,T = , ρ0 ρ0 RT0 U0 Tf − T 0
x=
x∗ t∗ ∇∗ τ∗ μ∗ κ∗ , t = L 0 , ∇ = 1 , τ = μ0 U 0 , μ = ,κ = L0 μ0 κ0 L U L 0
0
(1)
0
T0 is the ambient temperature and Tf is the flame adiabatic temperature. The non-dimensional set of conservation equations, in conservative form, reads:
T −T
∂ρ + ∇ · (ρu) = 0 ∂t
(2)
1 ∂ρu + ∇ · (ρuu) = −∇ p + ∇ ·τ ∂t Re
(3)
1 ∂ρT + ∇ · (ρuT) = ∇ · (κ∇T) ∂t RePr
(4)
p0 = ρ (1 + τh T)
(5)
with τh = fT0 0 being the specific heat parameter. We have assumed an open system, thermodynamic pressure, p0 , flow properties and viscosity being constant and we have considered negligible influence of gravity, negligible influence of viscous dissipation in the energy equation and validity of the state equation for the ideal gas.
190
Flow Turbulence Combust (2008) 80:187–206
2.1 The pressure in the low Mach number context After the asymptotic analysis, once the acoustic effects are neglected, it can be shown that the static pressure can be written as 2 p(x, y, z, t) P = p0 (t) + M
(6)
We can interpret p0 as the thermodynamic pressure or as the mean value of the is static pressure and the quantity p is the deviation from the mean pressure [8]. M √ U0 √ the modified Mach number, M = γ M0 with M0 = γ RT . For the 1D case the total 0 pressure (static plus dynamic) in dimensional units is, 1 P∗Tot = P∗ (t) + ρ ∗ u∗2 2
(7)
making the last equation non-dimensional and using (6) we write the total pressure scaled with the ambient pressure of reference P∞ = ρ0 RT0 PTot =
2 P∗Tot 2 p + ρ0 U 0 1 ρu2 = p0 + M P∞ P∞ 2
(8)
2.2 Numerical method This section gives an outline of the main aspects of the numerical method we used for the time integration of the governing equations treated so far, while the details for the spatial discretization can be found in Appendix. The various quantities are defined on a staggered grid (Fig. 1).
Fig. 1 2D staggered grid. The velocities and the mass fluxes are defined at the edges of the cells while the scalars are at the center
Flow Turbulence Combust (2008) 80:187–206
191
2.2.1 Temporal discretization Following Treurniet [7], at first we advance the energy (4) (which has been put in non-conservative form by using the continuity equation in order to have T explicit): T n+1 − T n = t α (AT + DT )n − β (AT + DT )n−1 (9) where α = 1.55 and β = 0.55 are modified Adams–Bashforth second order (AB2) coefficients. One can show that the stability region for this scheme is larger than for the classical AB2 and comparable to a second order predictor–corrector scheme used by Treurniet. The operators AT and DT represent the advective and diffusive terms in (4). Subsequently, from the equation of state the density is calculated as: ρ n+1 =
p0 1 + τh T n+1
(10)
Then we integrate the momentum equations to an intermediate level indicated with: − (ρu)n (ρu) = α (Am + Dm )n − β (Am + Dm )n−1 − ∇ pn (11) t Then from the intermediate level ρ u we obtain a divergence free quantity (ρu)n+1 with the aid of a pressure correction defined by (ρu)n+1 − (ρu) = −∇ p∗ t
(12)
Applying the divergence on both sides the last equation becomes ∇ · (ρu)n+1 − ∇ · (ρu) = −∇ 2 p∗ t
(13)
Because of conservation of mass we have ∇ · (ρu)n+1 = −
∂ρ n+1 ∂t
(14)
The derivative of the density is calculated with the following second order backward discretization in time: ∂ρ n+1 1 n+1 − 4ρ n + ρ n−1 = 3ρ ∂t 2t
(15)
Now substituting (14) in (13) gives the Poisson equation, + ∇ · (ρu) t
∂ρ n+1 ∂t
= ∇ 2 p∗
(16)
Equation 16 must be solved. Once this is solved, (12) is used to calculate the mass flux at time level n + 1: − t∇ p∗ (ρu)n+1 = (ρu)
(17)
And the pressure is updated by adding its calculated correction value: pn+1 = pn + p∗
(18)
192
Flow Turbulence Combust (2008) 80:187–206
2.3 Immersed boundary method The main computational bottleneck in the numerical procedure outlined above is the solution of (16). For an arbitrary grid this can be done with an iterative solver, but it is a rather slow procedure, especially in 3D cases. For certain simple grid cases so-called Fast Poisson solvers exist, which are based on the separable nature of the Poisson equation. In this paper we have used a method which allows for complex geometries while still solving a simple Poisson problem. The main idea is to add a force f to the equation of motion, ∂ρu 1 + ∇ · (ρuu) = −∇ p + ∇ ·τ + f ∂t Re where f represents the body force. The force f can be prescribed on a Cartesian mesh so that the efficiency of the solution procedure on simple grids is maintained [9].
2.4 Fadlun and Verzicco’s method This was one of the first IBMs applied to a combustion problem but in the incompressible case. In [9] a 3D complex geometry is simulated (a IC piston) with moving boundary. Here we applied the method for a simpler configuration with a fixed square body in the flow. Moreover, in the framework of a staggered Cartesian grid, we put the body with one side aligned with the grid. In this way, the normal components of the velocity can be imposed exactly at the boundary of the object. The method computes the velocity value of each point closest to the boundary as linear interpolation between the zero velocity we want to simulate at the wall position and the velocity of a point further into the flow. The interpolation procedure is illustrated in Fig. 2a. The velocity vi−1 is known and the value vi is a linearization between vi−1 and zero. This procedure is equivalent to applying a body force to the momentum equations locally. The use of this method requires particular care for the treatment of the corner points because each of them receives the contribution from two faces (in 2D case) and this must therefore be taken into account in the interpolation procedure.
2.5 Breugem and Boersma’s stress method Compared to the method in Section 2.4, this method replaces the stress in the momentum equations (at first solved without body) to ensure that no-slip conditions exist at the boundary of the object. Here we show the main idea for a 2D case by applying the method on one of the walls of a body placed in the flow. In Fig. 2b a simple case is depicted with a body aligned along the mesh with its sides coinciding with the mesh points where the normal velocities are defined. We can imagine to apply a force ft (see Fig. 2b) to have no-slip at the position of the cross on the north boundary of the body. In this method three velocity values belonging to the body are involved: vi, j−1 , vi+1, j−1 , ui, j−1 . Looking in the spatial
Flow Turbulence Combust (2008) 80:187–206
193
a
b
c Fig. 2 a The Fadlun–Verzicco IB method: the velocity defined one grid point into the flow with respect to the body wall is a linear interpolation which ensures zero tangential velocity at the boundary. b The Breugem–Boersma IB method: at the position of the body the original tangential stress is removed and the new one is added to the momentum in order to impose no-slip at the wall. This is equivalent to applying a force ft which does not depend on a linear interpolation of velocity values outside the body but depends on the velocities at the boundary and inside it. c Temperature treatment for IB methods. The heat flux contained in the diffusive term of the energy equation is replaced with the correct flux we want to impose. This has also consequences for the convective term (par. 2.7)
discretization of the momentum (reported in Appendix) we see that these velocity values are involved in the following terms 1 vi+ 1 , j−1 (ρu)i, j− 1 vi, j−1 + vi+1, j−1 · 14 ρi, j + ρi+1, j ui, j + ρi, j−1 + ρi+1, j−1 ui, j−1 2 2 =2 yj yj
ui, j −ui, j−1 v j−1 −vi, j−1 μ τ yxi, j− 1 + i+1, xh Re yh j i 2 − =− yj yj contained in the advective and diffusive part of the momentum respectively.
194
Flow Turbulence Combust (2008) 80:187–206
These then form the terms which must be modified to make the flow ‘feel’ the body. In order to have non-penetration at the wall and no-slip condition at the position of the cross (Fig. 2b) vi, j−1 = vi+1, j−1 = 0 and ui, j−1 = −ui, j must hold and the above terms become, vi+ 1 , j−1 (ρu)i, j− 1 2
2
yj −
τ yxi, j− 1
2
yj
=−
μ Re
=0
2ui, j yh j
yj
The last equation expresses the fact that we update the stress term at the points half a grid cell away from the wall. This procedure is equivalent to subtracting the old flux term, Fold = vi+ 1 , j−1 (ρu)i, j− 1 2
yj
Fnew = −
2
μ Re
−
2u
i, j yh j
τ yx
i, j− 12
yj
, from the momentum (11) and to adding the new flux term
. Therefore a force ft is applied to the flow, ⎤ ⎡ u +u v
i, j i, j−1 i+1, j−1 −vi, j−1 − vi+ 1 , j−1 (ρu)i, j− 1 yh xh μ j i 2 2 ⎦ ⎣ − ft = −Fold + Fnew = − yj Re yj yj
(19) Furthermore, the non-penetration condition at the walls is enforced by imposing a zero value for the normal component of the mass flux at the intermediate time level, and for the related velocity. (ρu), This procedure for the stress is also applied to the inner side of the wall. 2.6 Penetration velocity treatment As it has already been seen in both the IBMs outlined above, the normal component of the velocity at the boundaries must be as close to zero as possible, because we do not want normal penetration at the boundary. Consider (17) again, − t∇ p∗ (ρu)n+1 = (ρu)
(20)
and this means that the What we enforce to be zero is the predicted mass flux (ρu) penetration at the wall is of O(−t∇ p∗ ). Therefore, it is better to keep both t and ∇ p∗ small. However, if the pressure correction does not remain small, for instance due to strong perturbations in the flow, then we have to reduce the time step. 2.7 Temperature treatment When temperature differences are introduced, (for example in the case of a hot body placed in the flow), the heat flux between boundaries and the flow can be well represented with a procedure similar to the stress replacement method for the momentum equation. In practice, we adjust the heat flux on the walls locally by modifying the diffusion and convective terms in the energy equation in such a
Flow Turbulence Combust (2008) 80:187–206
195
way to force the correct flux at the boundary. We write the energy equation in nonconservative form, 1 1 ∂T + u∇ · T = ∇ · (κ∇T) ∂t ρ RePr Let us consider the diffusion term at position i (Fig. 2c):
Ti+1 −Ti Ti −Ti−1 − κ κ xi xi−1 1 1 1 1 ∇ · (κ∇Ti ) = ρi RePr ρi RePr x
i −Ti−1 with the correct flux we want to impose at the We replace the term κ Tx i−1 boundary, Ti − Twall κ xi−1 /2 Similarly one has to replace the flux at (i − 1) as well. If we use a central differences scheme, we also have to consider the convective term ui + ui−1 Ti+1 + Ti Ti + Ti−1 − (u∇T)i = 2 2xi 2xi as at the wall we want T = Twall then we have to update the convective term in this way, ui + ui−1 Ti+1 + Ti Twall − (u∇T)i = 2 2xi xi the same applies to the convective term at the position (i − 1). This method allows full control of the temperature of a simulated body and it also allows to solve the heat equation inside the body.
3 Computational Results Three geometries have been used for the computation of a flow over a square obstacle. The first one (Fig. 3a) is the same two-dimensional configuration used by Franke [16] and it has been used for a comparison of the method with validated computational data. The second geometry (Fig. 3b) is also two-dimensional but is smaller than the first one. It considers a square cylinder placed symmetrically with respect to the height of the channel. The third geometry (Fig. 3c) is three-dimensional with a cube (representative for example of a motor vehicle in a tunnel) placed four grid cells away from the bottom of the channel. The grid resolutions of these geometries are given in Table 1. 3.1 Comparison with validated computational data In Franke et al. [16] the Strouhal number, St, and the drag coefficient Cd, are computed for different Reynolds numbers. These results are validated with experimental data. The geometry consisting of a square cylinder with diameter D placed in a
196 Fig. 3 a The Franke 2D geometry for the computation of a flow over a square cylinder. The cylinder is located along the central line of the channel. Free-slip boundary conditions are applied at the walls. b The 2D geometry used for the comparison between the Fadlun–Verzicco, Breugem–Boersma and iterative solver approaches. The cylinder is located along the central line of the channel. Periodic boundary conditions are applied at the walls. c The geometry for the 3D flow cases. The body is a cube of diameter D located 4 grid cells above the bottom of the channel. No-slip conditions are applied at the walls
Flow Turbulence Combust (2008) 80:187–206
a
b
c
channel flow is shown in Fig. 3a. The inflow velocity is unity. Free slip boundary conditions are applied at the walls. Firstly, for Re = 10, we compute the L2-norm (Fig. 4) of the vertical velocity component in one point along the center-line of the channel at a position 4D away from the back of the body. The values are calculated with respect to the number of points in the vertical direction by computing the flow for several grid resolutions, considering the solution on the finest grid as ‘exact’. The method is essentially accurate to the second order. In Table 2 some values of St and Cd for different Re numbers are compared with the data in [16]. The comparison is good. The Strouhal number was calculated by
Flow Turbulence Combust (2008) 80:187–206
197
Fig. 4 The L2-norm error of the vertical velocity versus the number of grid points in z−direction for the simulation of a flow past a square cylinder. Case for Recylinder = 10
counting the number of periods it takes for the z-velocity component to invert its sign at a point located 4D past the body along its center line. The drag was calculated by considering the contribution of all the horizontal forces acting on the cylinder. 3.2 Two dimensional flow: comparison between the IBMs and a standard method We have performed a 2D simulation of a cylinder in a flow, for the geometry illustrated in Fig. 3b, by using the methods described in paragraphs 2.4–2.5. The body is located symmetrically with respect to the z dimension, so that it is exactly aligned with the center line of the channel. In fact, due to the periodic boundary conditions for the flow in z−direction the body can have any position along z. For this case, the density is constant. The inflow mass flux in x-direction is given (with velocity uin = 1). Zero gradient condition for the pressure is used at the inlet. At the outlet the pressure is put to zero (or ambient pressure). The Reynolds number based on the side of the square cylinder is 1,400. Figure 5 shows a comparison between the new IBM, the Fadlun–Verzicco approach and a conventional iterative solver without immersed boundary. The averaged values of the velocity component u along the center line of the cylinder are plotted. The iterative solver was based on the SIP method [12] and it used 100 iterations for each time step to get a divergence of O(10−6 ) for the cell at the middle of the front of the body (Fig. 7). It was found (Fig. 5) that the results of new IBM and iterative solver coincide while for the Fadlun–Verzicco approach the recirculation length is smaller than for the other methods indicating a larger resistance to the flow. The velocity profiles have been averaged over a range of 15 shedding periods between the 25th and the 40th period.
Table 1 Grid resolution for 2D and 3D cases 2D case 3D case
x-Points
y-Points
z-Points
200 200
2 64
140 140
198
Flow Turbulence Combust (2008) 80:187–206
Table 2 Comparison of St and Cd for different Re numbers for square cylinder with validated data This work (grid 300 × 180 equidistant)
Franke et al. [16] (grid non-equidistant)
Re
St
Cd
Re
St
Cd
100 150 300
0.153 0.162 0.138
1.61 1.56 1.80
100 (grid: 88 × 76) 150 (grid: 88 × 76) 300 (grid: 186 × 156)
0.154 0.165 0.130
1.61 1.56 1.83
Furthermore, we have checked the new method described in paragraph 2.7, for the case with temperature and density coupled by the equation of state. The temperature is set with zero gradient on all boundaries. The body has a small temperature difference T = 0.1 with respect to the flow and the flux replacing method is used. The normal velocity at the wall of the body is of O(10−6 ÷ 10−5 ) with maximum penetration at the corner points and with minimum of O(10−7 ) inside the obstacle (Fig. 6). In these two dimensional simulations no enforcing of zero velocity has been applied inside the body. By enforcing it the magnitude of the velocity components inside the cylinder can vary between 10−7 ÷ 10−9 . In these simulations, with time step being T = O(10−4 ), the shedding appeared earlier for the Fadlun–Verzicco method (around 4,300 time steps) than for the stress method (around 40,000 time steps). This suggests that the numerical noise produced by the stress method is much smaller than for the case of the Fadlun–Verzicco approach. 3.3 Three dimensional flow: body in a laminar flow with periodically varying inflow conditions In the following all computations were conducted with the new stress IBM described in the paragraphs 2.5–7. This case refers to the geometry in Fig. 3c. The body is a cube of diameter D located symmetrically with respect to the y-direction and 0.2D
Fig. 5 2D square cylinder case, averaged x-velocity component versus body x-center line. Comparison between Fadlun–Verzicco, Breugem–Boersma and iterative solver approaches. The plots of the last two methods overlap while the Fadlun–Verzicco method gives a smaller recirculation area
Flow Turbulence Combust (2008) 80:187–206 Fig. 6 2D case, time-averaged x-velocity component versus x-center line inside the cylinder for the different IB methods and the standard (iterative) one
199
5e-07 Iterative solver New IBM Fadlun-Verzicco IBM
4e-07 3e-07 2e-07 1e-07 0 -1e-07 -2e-07 -3e-07 0.3
0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39
0.4
(4 grid cells) away from the bottom of the channel. The Reynolds number, based on the height l z of the box was Re = 2,100. For this test zero gradient condition for the pressure was used at the inlet, while at the outlet the pressure was put to zero. The flow was at constant temperature and a varying inflow velocity was applied defined by U inflow = U + 0.1U · sin (2π f · t) with U = 1. The inflow velocity was varied by 10% of its bulk value with frequency f . This means f cycles of sinusoidal perturbation occurring in one unit time (of O(D/U)). We performed tests with three frequencies: 100, 10 and 1. For these three cases, Figs. 8, 9 and 10 illustrate the x-mass flux component (ρu)x versus time at the body wall, in the cell adjacent to the center-point of the left body wall (Fig. 7), on
Fig. 7 τxz and mass flux (mx ) vectors at the cell P adjacent to the center point of a body wall
τxz P. mx
z x
200 Fig. 8 (ρu)x versus time at the body wall for the inlet mass flux perturbation of frequency f = 1
Flow Turbulence Combust (2008) 80:187–206 5e-07 0 -5e-07 -1e-06 -1.5e-06 -2e-06 -2.5e-06 -3e-06 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45 0.5
a range of half time unit. The results show that for the frequencies of 1 and 10 the mass flux at the body wall responded with the same frequency, on the contrary this was not the case for the highest frequency of 100. In this case a sort of ‘integral’ effect is noted and the body reacted like if affected by a global perturbation and the penetration was higher by an order of magnitude: from O(10−6 ) for the first two frequency cases to O(10−5 ) for the highest frequency case. However, for this case, the penetration was still acceptable. We can expect that for even higher frequencies the reaction time of the new IBM has a limit: for high frequencies it takes longer for the velocity field inside the body to adjust itself. In fact, as we have mentioned above, the penetration velocity depends on both time step and pressure correction (see (20)) and the higher the frequency of the perturbation becomes the higher the pressure correction. Therefore, in this case, we have to reduce the time step if we want small penetration and in particular, if f is the frequency of the perturbation applied to the system, we must have t · f 1.
Fig. 9 (ρu)x versus time at the body wall for the inlet mass flux perturbation of frequency f = 10
4e-06 3e-06 2e-06 1e-06 0 -1e-06 -2e-06 -3e-06 -4e-06 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Flow Turbulence Combust (2008) 80:187–206 Fig. 10 (ρu)x versus time at the body wall for the inlet mass flux perturbation of frequency f = 100
201
3.55e-05
3.5e-05
3.45e-05
3.4e-05
3.35e-05
3.3e-05
3.25e-05 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
3.4 Three dimensional flow: body with varying heat flux We wanted to prove that with the proposed extension for the heat transfer the Breugem–Boersma approach is also suitable for combustion problems. In particular, for premixed flames the interaction of the reacting front with an obstacle produces sharp temperature gradients at the body boundary. The method must be stable under these conditions and for an adiabatic obstacle no transport of energy or mass must take place at the walls to avoid non-physical oscillations of velocity, temperature and density in the surroundings of the body. For this case the Reynolds number, based on the height l z of the box, was Re = 2100 with inflow velocity on the left side of the domain uin = 1. The time step t was of O(10−4 ). The temperature of the entire body was perturbed with a positive where i represents the progressive number of time sine profile Tw = sin 2π it T steps during the calculation, t is the time and T is the period of the oscillations. Two frequencies have been applied ( f = 100 and f = 200) with periods of oscillations T = 0.01 and T = 0.005 respectively. The flow time scale is of O(D/uin ) = 0.1 therefore we apply an intense perturbation of the temperature with respect to the flow. The positive sine function has the property to go from 0 to 1 that in this model means from the ambient to the adiabatic flame temperature (here the dimensional values are 300 and 1,800 K respectively, comparable with a premixed air-methane flame). Moreover, after one period its change is sharp enough to check if the central difference scheme suffers of numerical instability for large temperature gradients in time. The influence of the different frequencies can be seen on several quantities, in particular, velocity, mass flux, pressure and boundary layers. We want to know if we are resolving the momentum and the thermal boundary layers. The ratio of the thickness of these two boundary layers is governed by the Prandtl number. Because in all our laminar simulations the Prandtl number was smaller than unity the thermal boundary layer was thicker than the momentum boundary layer. Therefore we only have to check if we are resolving the momentum boundary layer. We estimated
202 Fig. 11 Averaged z+ along the top x-center line of the body for two heat flux frequency cases ( f = 100 and f = 200)
Flow Turbulence Combust (2008) 80:187–206 0.65
lower frequency case higher frequency case
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15
0.3
0.31
0.32
0.33
0.34
0.35
0.36
0.37
0.38
0.39
0.4
the parameter z+ and the wall shear stress τw at the top face of the body with the formulae: 2 + τ2 τw = τzx zy z+ =
z p ν
τw ρ
where z p is the distance from the wall and ν is the kinematic viscosity. If z+ ≤ 5 we have sufficient grid points to resolve the viscous sublayer [13]. The data related to the first frequency case were averaged over 2,000 time steps while in the second case over 1,000 time steps with t = O(10−4 ). In both cases 20 periods were considered. In Fig. 11 plots of the averaged values of z+ are compared for the two heat flux frequencies cases. These pictures are plotted along the x-axis center line on the top of the body. We can see that in both cases we have z+ 5.
Fig. 12 Averaged x-velocity component, u, inside the body along its x-center line for two heat flux frequency cases ( f = 100 and f = 200)
0 lower frequency case higher frequency case
-5e-05 -1e-04 -0.00015 -0.0002 -0.00025 -0.0003 -0.00035 -0.0004 0.3
0.32
0.34
0.36
0.38
0.4
Flow Turbulence Combust (2008) 80:187–206 Fig. 13 x-Mass flux component (ρu)x versus time at the body wall for two heat flux frequencies cases ( f = 100 and f = 200)
203
1e-04 lower frequency case higher frequency case 5e-05
0
-5e-05
-0.0001
-0.00015
-0.0002 0.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
0.65
Figure 12 shows the profiles of the averaged x-velocity component, u, inside the body along its center line. The order of magnitude of u is between O(10−5 ) and O(10−4 ). Fig. 13 shows the x-mass flux component (ρu)x versus time at the wall of the body (in the grid cell depicted in Fig. 7). The mass flux penetration at the wall is of O(10−5 -10−4 ) as well. In Fig. 14 a time interval of an instantaneous plot of the total pressure (8) is shown as function of time. The values are taken at the yz plane which was 0.6l x away from the inlet and 0.2l x from the right side of the cube (l x being the x-size of the channel). The second frequency value for the perturbation of Tw was used in this case. We note periodic oscillations of the pressure with frequency of O(100) that is half of the frequency of the imposed heat flux. The picture shows the peaks of pressure oscillations corresponding to the minimum values of body temperature (when the sine profile starts a new cycle) but the oscillations do not create instabilities.
Fig. 14 Sample of total pressure versus time at a location two diameters past the body. The peaks correspond to the moment when the sinusoidal perturbation starts a new cycle
1.00006 1.00005 1.00004 1.00003 1.00002 1.00001 1 0.99999 0.99998 0.4
0.45
0.5
0.55
0.6
0.65
204
Flow Turbulence Combust (2008) 80:187–206
4 Final Considerations and Conclusions We summarize here the main advantages and limitations of the proposed method. Our code uses a pressure-correction method. The solution of the momentum and mass balance equations is satisfied in two steps. First the momentum balance is satisfied using the pressure at the old time level. Next the velocity is corrected to satisfy the mass balance. Correction is done with an update of the pressure for which we have to solve a Poisson equation. The boundary conditions are enforced by addition of extra momentum sources. This is done before the pressure-correction and is independent of it. As a result we can use so-called fast Poisson solvers which only work on a separable domain. Inclusion of obstacles makes the Poisson problem non-separable for a traditional approach and we have to resort to slower iterative solvers [14]. Therefore, as first advantage, the method we use is faster than traditional methods (like iterative solvers). The speedup with respect to Cartesian codes which do not use a direct solver is in general a factor of 2.5 or more (which of course depends on the divergence required and on the resolution) [14]. The code is cheaper by a factor of 10 or more per grid point when compared to curvilinear and unstructured codes [14]. Our IBM has been specifically designed for square geometries and to consider complicated geometries including many obstacles. In fact, by aligning a square body with the grid lines it is possible to apply almost exactly the non-penetration conditions for the vertical components of the velocity. Velocity interpolations are also avoided near the wall both for the non-penetration and the no-slip conditions. This is the second advantage. For non-aligned bodies or for curvilinear geometries our method can still be applied but its implementation becomes complicated because interpolations are required. In this case it loses the third advantage which is its simplicity to be implemented and to add as many square bodies as required which can be included independently of one another, or removed [10]. In conclusion, for non-aligned or curved obstacles the Fadlun–Verzicco method is certainly more suitable. However, it has higher wall normal leakage of the order of 10−3 of the bulk velocity [14]. In both methods the near wall quantities can not be well represented out of the subviscous layer, however in [15] has been shown that in this case the Fadlun–Verzicco approach presents larger errors due to its linear interpolation into the flow. In conclusion, we have performed calculations with an immersed boundary method for the case of a hot or insulated square body located in a box-domain open at the two faces perpendicular to the main stream direction. Two geometries have been used, a 2D geometry with few grid cells in the span direction and a fully 3D geometry. The object has been simulated with a new immersed boundary method. Several tests have been carried out by using different conditions for the inflow velocity and the heat flux from the body. Initially we have made a two dimensional simulation of a cold and hot body in a laminar channel flow with constant inflow conditions and comparison with a previous IBM and an iterative solver has been performed. The new IBM shows a very good agreement with respect to the results of an iterative solver. For the three dimensional geometry, at first we have simulated laminar cases with perturbation of the inflow mass flux for different frequencies. Secondly, we have studied the case of a laminar flow with a hot body whose heat flux and internal temperature can vary within the temperature limits of a common air-
Flow Turbulence Combust (2008) 80:187–206
205
methane premixed flame. For this last case two tests have been done with different frequencies of the temperature oscillations. The results show good behaviour of the stress IBM in all conditions considered with a velocity at the body walls and inside the object of an order of magnitude varying between 10−9 and 10−4 depending on the magnitude and the frequency of the disturbances applied. The method is able to adjust the velocity field inside the obstacle also for high frequency of mass or heat flux perturbations. We conclude that this new immersed boundary method is easy to implement, it requires less computational resources than standard methods, it has accuracy comparable with that one of an iterative solver approach and it is suited for heat transfer and combustion problems.
Appendix: Spatial Discretization The spatial discretization (second order finite volume) is outlined here, for pure central differences scheme of the terms A (advection) and D (diffusion) for the x-components of the energy and momentum equations. In the other directions the procedure is similar. For the energy equation (in non-conservative form) we have:
1 Ti+1, j,k − Ti−1, j,k 1 AT = − ui−1, j,k + ui, j,k 2 2 xi DT =
1 ρi, j,k
Ti+1, j,k + Ti−1, j,k − 2Ti, j,k 1 1 κ RePr xi xi
For the x-component of the momentum Considering u as the x-component of the velocity it is: Am = ∇ · (−u(ρu))i, j,k = −
−
−
Dm =
ui+ 1 , j,k (ρu)i+ 1 , j,k − ui− 1 , j,k (ρu)i− 1 , j,k 2
2
2
+
2
xhi
vi+ 1 , j,k (ρu)i, j+ 1 ,k − vi+ 1 , j−1,k (ρu)i, j− 1 ,k 2
2
2
2
yj wi+ 1 , j,k (ρu)i, j,k+ 1 − wi+ 1 , j,k−1 (ρu)i, j,k− 1 2
2
2
2
zk
τxxi+ 1 , j,k − τxxi− 1 , j,k 2
2
xhi
+
τ yxi, j+ 1 ,k − τ yxi, j− 1 ,k 2
2
yj
+
τzxi, j,k+ 1 − τzxi, j,k− 1 2
zk
2
206
Flow Turbulence Combust (2008) 80:187–206
with, τxxi+ 1 , j,k 2
τ yxi, j+ 1 ,k 2
τzxi, j,k+ 1 2
ui+1, j,k − ui, j,k μ 2 = − i+1, j,k 2 Re xi+1 3 vi+1, j,k − vi, j,k μ ui, j+1,k − ui, j,k = + Re yh j xhi wi+1, j,k − wi, j,k μ ui, j,k+1 − ui, j,k = + Re zhk xhi
and i, j,k =
ui+1, j,k − ui, j,k vi, j+1,k − vi, j,k wi, j,k+1 − wi, j,k + + xi yj zk
References 1. Rehm, R.G., Baum, H.R.: The equations of motion for thermally driven, buoyant flows. J. Res. Nat. Bur. Stds. 83, 297 (1978) 2. Peskin, C.S.: Flow patterns around heart valves: a numerical method. J. Comput. Phys. 10, 252–271 (1972) 3. Mohd-Yusof, J.: Combined Immersed-Boundary/B-Spline Methods for Simulations of Flow in Complex Geometries. CTR Annual Research Briefs, NASA Ames Research Center/Stanford Univ. Center for Turbulence Research. Stanford, CA (1997) 4. Verzicco, R., Mohd-Yusof, J., Orlandi, P., Haworth, D.: Large eddy simulation in complex geometric configurations using boundary body forces. AIAA J. 38(3), 12–40 (2000) 5. McMurtry, P.A., Jou, W.-H., Riley, J.J., Metcalfe, R.W.: Direct numerical simulation of a reactive mixing layer with chemical heat release. AIAA J. 24, 962–970 (1986) 6. Warnatz, J., Maas, U., Dibble, R.W.: Combustion, 3rd edn. Springer (2000) 7. Treurniet, T.: Direct numerical simulation of premixed turbulent combustion. Ph.D. thesis, TU Delft (2002) 8. Majada, A., Lamb, K.G.: Simplified equations for low Mach number combustion with strong heat release. Dynamical Issues in Combustion Theory. Springer (1991) 9. Fadlun, E.A., Verzicco, R., Orlandi, P., Mohd-Yusof, J.: Combined immersed-boundary finitedifference methods for three-dimensional complex flow simulations. J. Comput. Phys. 161, 35–60 (2000) 10. Breugem, W.P., Boersma, B.J.: Direct numerical simulations of turbulent flow over a permeable wall using a direct and a continuum approach. Phys. Fluids 17(025103), 15 (2005) 11. Mittal, R., Iaccarino, G.: Immersed boundary methods. Annu. Rev. Fluid Mech. 37, 239–61 (2005) 12. Ferziger, J.H., Peric, M.: Computational Methods for Fluid Dynamics. Springer-Verlag, Berlin (2002) 13. Versteeg, H.K., Malalasekera, W.: An introduction to computational fluid dynamics. The finite volume method. Longman Scientific & Technical, Longman Group Ltd (1995) 14. Pourquie, M.J., Nieuwstadt F.T.M.: Proceedings of McMAT2005, Joint ASME/ASCE/SES Conference on Mechanics and Materials Baton Rouge, Louisiana, USA, 1–3 June 2005 15. Pourquie, M., Paravento, F., Nieuwstadt, F.: Proceedings of EUROMECH Colloquium 471 on Turbulent Convection in Passenger Compartments, Göttingen, Germany, 13–14 October 2005 16. Franke, R., Rodi, W., Schönung, B.: Numerical calculation of laminar vortex-shedding flow past cylinders. J. Wind Eng. Ind. Aerodyn. 35, 237–257 (1990)