Koren B., Vuik K. (Editors) Advanced Computational Methods in Science and Engineering

Подождите немного. Документ загружается.

Yunus Hassen and Barry Koren

In Figure 18, we show the error E, which is computed as the difference between

the exact and numerical solutions, for the solutions given in Figure 17. As can be

seen, there is relatively more discrepancy near the EBs. This near-EB discrepancy is

of the same order for both test cases. However, as mentioned earlier, the discrepancy

is signiﬁcantly larger, about ﬁve times more, in the cosine-with-cavity case at and

near the extrema, due to the limiters. In the former case, the results obtained with

the RK3b scheme are superior to those obtained with the Forward Euler scheme, for

the obvious reason. In the latter case, both schemes, limited and/or unlimited, yield

almost the same accuracy.

0 0.2 0.4 0.6 0.8 1

−1

−0.5

0.5

x 10

−3

(a) On a 20-cell grid, for (5a)

0 0.2 0.4 0.6 0.8 1

−0.05

−0.03

−0.01

0.01

0.02

(b) On a 20-cell grid, for (5b)

0 0.2 0.4 0.6 0.8 1

−1

−0.5

0.5

x 10

−3

0 0.2 0.4 0.6 0.8 1

−0.02

−0.01

0.005

0.01

(d) On a 40-cell grid, for (5b)

Fig. 18 Errors after one full-period, for the initial solutions (5).

-blue: unlimited higher-order

upwind-biased with Forward Euler,

∗-red: limited higher-order upwind-biased with Forward Euler,

⋄-green: unlimited higher-order upwind-biased with RK3b, ×-black: limited higher-order upwind-

biased with RK3b.

6 Extension to more general cases

First extensions to the method presented so far, which must and will be made, are:

(i) to higher dimensions and (ii) to higher-order accuracy in time. We already per-

264

Finite-Volume Discretizations and Immersed Boundaries

formed some work into this direction. Here follows a brief account of the ideas we

are currently pursuing.

6.1 Extension to higher dimensions

Extension to 2D and 3D of the current 1D space discretization is done by dimen-

sional splitting. For this purpose, a multi-D embedded boundary is to be projected

ﬁrst on each of its separate coordinate directions (Figure 19). The details of the pro-

jection step are crucial. Dimensional splitting has been applied with success to the

standard

-scheme and, as such, is widely spread in CFD. We presume that it will

be successful here as well.

(a) control volume with EB (in red)

(b) 1D projections of the EB

Fig. 19 Example of 1D projections of a 2D embedded boundary.

6.2 Higher-order accuracy in time

Forward Euler can readily be made second-order accurate by following the Modiﬁed

Euler approach. Given (43), Modiﬁed Euler reads:

ˆc

n+1

= c

F(c

), (68a)

n+1

= c



F( ˆc

n+1

) + F(c

)



. (68b)

Forward Euler step (68a) is the predictor and (68b) the corrector step. Modiﬁed

Euler is still explicit. Extension of the local adaptivity in time, introduced in § 4.2,

from Forward Euler to Modiﬁed Euler is rather straightforward, given the close

similarity of the two schemes. Details about our local time adaptivity method and

Modiﬁed Euler will be given in future work.

265

Yunus Hassen and Barry Koren

7 Conclusion

A novel immersed-boundary approach, for solving advection problems, has been

introduced. The essence of the approach is that moving bodies are embedded in a

regular ﬁxed grid and speciﬁc ﬂuxes in the vicinity of the embedded boundary (EB)

are computed in such a way that they accurately and monotonously accommodate

the boundary conditions valid on the moving body. To suppress the wiggles that

exhibit near discontinuities, tailor-made limiters are introduced for the ﬂuxes that

are especially modiﬁed. Then, over the majority of the domain, where we do not

have inﬂuence of the embedded boundaries, we can readily use standard methods

on the underlying regular ﬁxed grid. Excellent results are achieved, without much

computational overhead.

In summary:

• A generalized

-scheme that uses EB information (EB location and EB solu-

tion values), and which is an optimally accurate upwind-biased ﬁnite-volume

discretization, has been proposed. This near-EB spatial discretization is a gen-

eralization of a well-proven ﬁnite-volume discretization, which allows us to ac-

commodate EBs.

• Generalized limiters that use EB information have been proposed. These limiters

satisfy the spatial monotonicity requirement and Harten’s TVD requirement. To

be consistent with standard limiters, the generalized limiters are made indepen-

dent of the CFL number.

• Locally adaptive splitting of the time step, near EBs, has been proposed; a two-

stage approach which requires the least computational time and memory for a

given gain in accuracy.

We foresee that the numerical methods, developed so far and still to be developed,

can readily be extended to realistic ﬂow problems.

Acknowledgements The ﬁrst author’s research is funded by the Delft Centre for Computational

Science and Engineering (DCSE).

Appendix

Referring to (25) and (29), the local truncation error terms e in cells i + 1 and i + 2,

when the net ﬂuxes in cells i + 1 and i + 2 are considered to optimize

, respec-

tively, are:

i+1

3 −2

∂

, (69a)

and

i+2

−1

∂

. (69b)

266

Finite-Volume Discretizations and Immersed Boundaries

For a given grid size h, (69) can be rewritten, as a function of

∈ [0,1], as:

i+1

= 9 −6

, (70a)

and

i+2

= 2

−1, (70b)

where

b =

∂

. (70c)

The scaled error terms (70a) and (70b) are plotted in the (e/b,

)-diagram given in

Figure 20. We see that |e

i+1

| ≥ 3 | e

i+2

|, ∀

. Also notice that e

i+2

= 0 for

;

i.e., we get a third-order accurate net ﬂux in cell i+ 2 when the EB is situated at the

center of cell i.

Fig. 20 Variation of the

scaled, leading local trunca-

tion error terms in cells i + 1

and i + 2 with

i+2

i+1

−1

References

1. Calhoun, D.: A Cartesian grid method for solving the two-dimensional streamfunction-

vorticity equations in irregular regions. Journal of Computational Physics 176(2), 231–275

(2002)

2. Fadlun, E.A., Verzicco, R., Orlandi, P., Mohd-Yusof, J.: Combined immersed-boundary meth-

ods for three-dimensional complex ﬂow simulations. Journal of Computational Physics 161,

35–60 (2000)

3. Godunov, S.K.: Finite difference method for numerical computation of discontinuous so-

lutions of the equations of ﬂuid dynamics. Matematicheskii Sbornik 44, 271–306 (1959).

Translated from Russian at the Cornell Aeron. Lab.

4. Goldstein, D., Handler, R., Sirovich, L.: Modeling a no-slip ﬂow boundary with an external

force ﬁeld. Journal of Computational Physics 105, 354–366 (1993)

267

Yunus Hassen and Barry Koren

5. Harten, A.: On a class of high resolution total-variation-stable ﬁnite-difference schemes.

SIAM Journal on Numerical Analysis 21, 1–23 (1984)

6. Hundsdorfer, W., Koren, B., van Loon, M., Verwer, J.G.: A positive ﬁnite-difference advection

scheme. Journal of Computational Physics 117, 35–46 (1995)

7. Kim, J., Kim, D., Choi, H.: An immersed-boundary ﬁnite-volume method for simulations of

ﬂow in complex geometries. Journal of Computational Physics 171, 132–150 (2001)

8. Koren, B.: A robust upwind ﬁnite-volume method for advection, diffusion and source terms.

In: Vreugdenhil, C.B., Koren, B. (eds.) Notes on Numerical Fluid Mechanics, 45, pp. 117-

138. Vieweg, Braunschweig (1993)

9. Leer, B. van: Upwind-difference methods for aerodynamic problems governed by the Euler

equations. In: Lectures in Applied Mathematics, 22 - 2, pp. 327-336. American Mathematical

Society, Providence, RI (1985)

10. Mittal, R., Iaccarino, G.: Immersed boundary methods. Annual Review of Fluid Mechanics

37, 239–261 (2005)

11. Mohd-Yusof, J.: Combined immersed-boundary/B-spline methods for simulations of ﬂow in

complex geometries. In: CTR Annual Research Briefs, pp. 317-327. Center for Turbulence

Research, NASA Ames/Stanford University (1997)

12. Mohd-Yusof, J.: Development of immersed boundary methods for complex geometries.

In: CTR Annual Research Briefs, pp. 325-336. Center for Turbulence Research, NASA

Ames/Stanford University (1998)

13. Peskin, C.S.: Flow patterns around heart valves: a numerical method. Journal of Computa-

tional Physics 10, 252–271 (1972)

14. Peskin, C.S.: Numerical analysis of blood ﬂow in the heart. Journal of Computational Physics

25, 220–252 (1977)

15. Peskin, C.S.: The ﬂuid dynamics of heart valves: experimental, theoretical and computational

methods. Annual Review of Fluid Mechanics 14, 235–259 (1982)

16. Saiki, E.M., Biringen, S.: Numerical simulation of a cylinder in uniform ﬂow: application of

virtual boundary method. Journal of Computational Physics 123, 450–465 (1996)

17. Su, S.-W., Lai, M.-C., Lin, C.-A.: An immersed boundary technique for simulating complex

ﬂows with rigid boundary. Computers & Fluids 36, 313–324 (2007)

18. Sweby, P.: High resolution schemes using ﬂux limiters for hyperbolic conservation laws.

SIAM Journal on Numerical Analysis 21, 995–1011 (1984)

19. Tseng, Y.H., Ferziger, J.H.: A ghost-cell immersed boundary method for ﬂow in complex

geometry. Journal of Computational Physics 192, 593–623 (2003)

20. Zhang, N, Zheng, Z.C.: An improved direct-forcing immersed-boundary method for ﬁnite

difference applications. Journal of Computational Physics 221, 250–268 (2007)

268

Large Eddy Simulation of Turbulent

Non-Premixed Jet Flames with a High Order

Numerical Method

S. van der Hoeven, B.J. Boersma, and D.J.E.M. Roekaerts

Abstract In this chapter we will report on the Large Eddy Simulation (LES) of a

turbulent non premixed jet ﬂame. The numerical model used for the LES is based on

a discretization with high order compact schemes. These schemes have a negligible

amount of numerical dissipation. The subgrid terms in the compressible Navier-

Stokes equations are modelled with simple eddy viscosity models. The LES model

is used to simulate the Sandia Flame D.

1 Introduction

Large-scale industrial burners typically have a cross section of the order of 0.1 meter

or larger. The smallest length scale in the reacting ﬂow generated using such a burner

is in general associated with the ﬂame thickness. The laminar ﬂame thickness of a

premixed ﬂame is to a good approximation independent of the size of the geometry

and depends on material properties such as viscosity, thermal diffusivity, activation

energy, etc. The ﬂame thickness is in general less than a millimeter. A full three

dimensional numerical simulation (or direct numerical simulation) in which both

the smallest and largest scales of the ﬂow are well resolved is not possible with

present-day computers. A possible solution to this problem is the use of Large Eddy

Simulation (LES), in which the small scales are removed from the spectrum by a

ﬁltering procedure. This ﬁltering splits the energy spectrum, see Figure 1, into two

parts, a resolved part (large scales) and an unresolved part (small scales). The large

scales are resolved on a computational grid. The effect of the small scales which

are removed by the ﬁltering is compensated for by a subgrid model. Depending on

S. van der Hoeven & D.J.E.M. Roekaerts

Department of Multi-Scale Physics , Delft University of Technology, the Netherlands

B.J. Boersma

Laboratory for Aero and Hydrodynamics, Mekelweg 2, 2628 CA, the Netherlands,

e-mail: b.j.boersma@tudelft.nl

269

Lecture Notes in Computational Science and Engineering 71,

DOI 10.1007/978-3-642-03344-5_9,

B. Koren and C. Vuik (eds.), Advanced Computational Methods in Science and Engineering,

S. van der Hoeven, B.J. Boersma, and D.J.E.M. Roekaerts

E(k)

LES DNS

large scales small scales

Fig. 1 A sketch of a 1D energy spectrum in a turbulent ﬂow. In Direct Numerical Simulation

(DNS) the full energy spectrum is resolved. In Large Eddy Simulation (LES) only the large scale

part of the energy spectrum is resolved, and the small scale part is modelled with a subgrid model.

the cut-off length of the ﬁlter we can have a well resolved LES or a poorly resolved

LES. The fact that ﬂame structure in general belongs to the smallest scales, makes

LES modeling of combustion more difﬁcult than modeling of inert ﬂow.

In general dissipation is associated with the small scales in a ﬂow. These small

scales have been removed by the LES ﬁlter and the subgrid scale model has to

provide the required dissipation. Most numerical models also introduce a certain

amount of dissipation. The amount of numerical dissipation strongly depends on the

type of scheme. Low order central schemes have in general a higher numerical dis-

sipation than high order schemes. Central (symmetric) schemes have in general less

dissipation than upwind (asymmetric) schemes. In an actual Large Eddy Simulation

the dissipation is always a combination of dissipation given by the sub-grid model

and by the numerics of the LES model. This makes the validation of a sub-grid

model complicated because dissipative effects can be contributed to the numerics

or the sub-grid model. For the proper validation of the sub-grid model numerical

dissipation should be minimized.

In the literature many low dissipative methods are reported, using high order

ﬁnite differences, compact ﬁnite differences and orthogonal polynomials, see for

instance Lele (1992), Chu & Fan (1998). The disadvantage of these schemes is that

they are not very stable in general. Recently, Boersma (2005) proposed a staggered

variant of the compact ﬁnite difference scheme as developed by Lele (1992). This

new scheme is stable because it has certain symmetry properties, see for instance

Verstappen & Veldman (2003). In this chapter we will explore the possibility of this

scheme for combustion simulations.The model is used to simulate the Sandia Flame

270

Large Eddy Simulation of Jet Flames

D, which is one of the test cases of the International Workshop on Measurement and

Computation of Turbulent Non-Premixed Flames [6].

2 Governing equations

The derivation of the governing equations for combustion is discussed in many text-

books, such as those of Peters [16] and Kuo [17]. Here we will present the Large-

Eddy ﬁltered version of these equations. A LES ﬁltered quantity is given by an

overbar and is obtained using the following ﬁlter function:

f (x) =

G(x −x

′

) f (x

′

)dx

′

. (1)

In this equation f is a ﬂow variable, V is the volume of the LES ﬁlter and G is a

ﬁlter kernel. The exact form of G is at this stage not important. In a compressible

ﬂow it is customary to use Favre-averaged (or density weighted) variables. This is

also the approach we will follow in the present paper. The Favre average, denoted

by a tilde, is deﬁned as follows

f =

. (2)

Variations with respect to the Favre average are denoted with a double prime, i.e.

′′

= f −

With help of the LES and Favre ﬁlters we can write for the equation for conservation

of mass:

∂

˜u

= 0, (3)

where

is the LES averaged density, and ˜u

the Favre-averaged velocity vector. The

LES ﬁltered momentum equations read:

∂

˜u

∂

˜u

∂

= −

∂

¯p

∂

i j

∂

(

˜u

−

). (4)

In this equation ¯p is the averaged pressure,

˜u

−

the so-called subgrid

stress, which will be discussed in more detail later, and

i j

the averaged stress tensor

with components

i j



∂

˜u

∂

˜u

∂



−

µδ

i j

∂

˜u

∂

, (5)

in which

is the averaged dynamic viscosity of the ﬂuid and

i j

the Kronecker delta

function. The dynamic viscosity is a function of the temperature. The temperature

dependence of

will be speciﬁed later. The equation for the mean of the total non-

chemical energy E = C

T + u

/2, with T the temperature and C

the speciﬁc heat

271

S. van der Hoeven, B.J. Boersma, and D.J.E.M. Roekaerts

at constant volume, reads [18]:

∂

˜u

H =

∂



˜u

i j

(

∂

˜x



H − ˜u





. (6)

Here H is the speciﬁc enthalpy,

the thermal conductivity which is a function of

the temperature and

is the heat release rate of the combustion process. In this

equation again a subgrid term appears, i.e.

H − ˜u

H, which has to be modelled.

In combustion modelling of non-premixed ﬂames it is customary to use an equation

for the mixture fraction

Z. This equation reads:

∂

˜u

∂



Z − ˜u



, (7)

where D is the diffusion coefﬁcient of the mixture, and again a subgrid term appears.

The mean pressure is obtained from the mean thermal equation of state (8):

¯p =

∑

i=1

, (8)

where at the right hand side unclosed terms containing correlations between tem-

perature and composition have been left out. In this equation R is the universal gas

constant, N

is the number of species and

∑

i=1

is the inverse of the mean molar

mass of the mixture.

The subgrid ﬂuxes in the LES ﬁltered Navier Stokes equation, equation (4) are

closed with a simple model based on an eddy viscosity assumption:

[gu

− ˜u

˜u

] = −2

ρν

i j

, (9)

where

i j

is the strain-rate tensor given by

i j



∂

˜u

∂

˜u

∂



. (10)

The turbulent viscosity

is given by a model proposed by Vreman [19]:

= C

i j

. (11)

The model constant is given by C

= 2.5C

, where C

is the Smagorinsky constant.

The functional dependence of the eddy viscosity on local conditions is given by:

Z =

−

, (12)

where

i j

∆

m j

, (13)

272

Large Eddy Simulation of Jet Flames

i j

≡ (∇

i j

∂

. (14)

In the present work we use an implicit LES ﬁlter. The ﬁlter width

∆

is the same in

all three Cartesian directions and proportional to the cubic root of the grid volume,

i.e.

∆

= (

∆

1/3

. (15)

The Vreman model was found to perform better than a standard Smagorinsky model,

because its dissipation is relatively small in the transitional region in the lower part

of the jet.

The subgrid term in the energy equation

H − ˜u

H is modelled via the temper-

ature gradient, i.e.

H − ˜u

H = −

∂

, (16)

/Pr

, (17)

where d

is the turbulent thermal diffusivity,

the eddy viscosity given by equation

(12), and Pr

the turbulent Prandtl number. Similarly in the resolved mixture fraction

equation

Z − ˜u

Z = −

ρΓ

∂

, (18)

/Sc

. (19)

In this study the Smagorinsky constant, C

was set to 0.08. Following [7] the tur-

bulent Prandtl number was set to 0.7 and the turbulent Schmidt number was set to

0.4.

The dynamic viscosity

and the thermal conductivity

are a function of the

temperature. We use the following simple relations





0.76





0.76

. (20)

In regions with large temperatures the increase in dynamic viscosity or thermal dif-

fusivity can be substantial. For the stability of the energy equations it is essential

that the increased thermal conductivity is taken into account.

3 Chemistry model

For the chemistry we use a steady ﬂamelet model [16], with a library consisting

of only one ﬂamelet with strain rate a = 100s

−1

. The low-Mach number resolved

273