Simulation of turbulent combustion in aerobic chambers:

extension of Eulerian Monte Carlo methods

Nadezda PETROVA*,Vladimir SABELNIKOV

*Nadezda.Petrova@onera.fr

Nomenclature

Abbreviations

CEDRE:Code d'Ecoulements Diphasiques Reactifs pour l'

Energetique

(CFD package for diphasic reactive ows)

CFL:Courant - Friedrich - Levy condition

EMC:Eulerian Monte Carlo

ICT:Interaction"chemistry - turbulence"

O-U:Ornstein - Uhlenbeck

PDF:Probability Density Function

RANS:Reynolds Averaged Navier - Stokes Equations

SPDEs:Stochastic Partial Dierential Equations

GLM:Generalized Langevin Model

Abstract

In this work,we investigate an Eulerian Monte Carlo (EMC) method in order to solve the

joint probability density function (PDF) of turbulent reactive scalars and velocity eld.This

method has been recently proposed by V.Sabelnikov and O.Soulard [2].We use dierent

numerical schemes for solving the stochastic partial dierential equations (SPDEs) on which

the EMC method is based.To validate our approach we consider a one dimensional turbulent

reactive ow.The system of SPDEs is weakly hyperbolic and so the Godunov - like methods

(9) are not adapted for the SPDEs.Another problem which arises in SPDEs is the presence

of a multi - valued solution.We propose to nelgect the multi - valued solution of the SPDEs

and calculate the statistics using kinetic schemes in space [8].

Introduction

The modern energy and environmental context motivates the improvement of combustion

systems in aerobic chambers.The structure of turbulent ames in industrial combustion cham-

bers is complex:it is governed by turbulence,two - phase injection,chemical kinetics,acoustics

and radiation.Once validated by comparison with experimental data,numerical simulations

allow to improve the understanding of the dierent physical mechanisms,in particular,analyze

the combustion stability,explain the formation polluting emissions.The description of the in-

teraction"chemistry - turbulence"(ICT) is a key point in the development of such numerical

methods.Approaches based on the one-point joint velocity - scalar probability density functions

(PDF) is a natural and promising tool for the description of the ICT [1].But it has a severe nu-

merical constraint:the PDF possesses a potentially high number of dimensions,which induces

heavy computational costs.On the other hand,Monte Carlo methods yield a linearly growing

eort and so are well suited to solve PDF equations.For this reason the Monte Carlo methods

have been traditionally used to solve the PDF equations.

1

EMCmethods are based on stochastic Eulerian elds,which evolve according to the SPDEs

statistically equivalent to the PDF equation.EMC methods have already been proposed in [3]

and in [4],in order to compute the one - point PDF of turbulent reactive scalars.But these meth-

ods are numerically expensive,because the CFL criterion is similar to an advection/diusion

stability criterion.Afterwards,V.Sabelnikov and O.Soulard [2] have proposed SPDEs to solve

the Favre joint PDF of turbulent reactive scalars and velocity eld.The Ornstein - Uhlenbeck

process for uctuating velocity allows to use the time step which is proportional to a grid size

x divided by the total stochastic velocity j

e

U +u

00

j.The time step t for the RANS equations

with the SPDEs is expressed as:

t = min

"

CFL

x

j

e

U +u

00

j

;CFL

x

e

U +a

#

;(1)

where CFL is the CFL constant,

e

U is the Favre averaged velocity and u

00

is the Favre uctuations

of velocity and a is the speed of sound.

An EMC solver is implemented into the industrial code CEDRE of the ONERA.It is

tested for one and two dimensional turbulent ows.It can be shown that the SPDEs are weakly

hyperbolic,so the numerical schemes of Godunov type are not applicable to SPDEs.These

numerical spatial schemes cannot also capture the multi - valued solution of SPDEs.

The paper consists of four parts.Firstly,we recall the joint PDF of uctuating velocity

and scalar.Secondly,we write the SPDEs statistically equivalent to the PDF equation.Then

we describe numerical aspects and numerical validation tests in section 3.Finally,we adress the

following issues:a) the numerical schemes which are adapted to solve the weakly hyperbolic

equations;b) multi - valued solution of the SPDEs.

1 Favre joint uctuating velocity - scalar PDF

Let us consider the joint one point velocity - scalars PDF in turbulent reactive gas ows.

For simplicity writing we use only one reactive scalar c.For variable density ows,it is usual to

work with so - called Favre statistics,weighted by the density.The Favre and Reynolds PDFs

of velocity and scalar are related as

e

f

Uc

(U;c) =

< jU;c >

f

Uc

(U;c);(2)

where

e

f

Uc

and f

Uc

are the Favre and Reynolds PDFs;and < jU;c > represents the conditional

density for xed velocity and scalar.

Below,we work with the uctuating velocity u

00

and not the total stochastic velocity U:

u

00

= U

e

U;(3)

where

e

U is the Favre averaged velocity.eg is now the joint uctuating velocity - scalar PDF.It

is linked with the Favre joint velocity - scalar PDF by the relation:

eg(u

00

;c) =

e

f(

e

U+u

00

;c):(4)

It is supposed that the low Mach number assumption is veried,so the density can be expressed

as a function of the scalar concentration.In consequence,the model transport evolution equation

for the PDF eg(u

00

;c) is dened as [7]:

@

@t

(< > eg) +

@

@x

j

< >

f

U

j

+u

00

j

eg

=

2

@

@u

00

j

@

@x

k

< >

]

u

00

j

u

00

k

eg

+

@

@u

00

j

< >

C

1

!u

00

j

+

@

f

U

k

@x

k

!

u

00

j

eg

!

+

1

2

< > C

0

< >

@

2

eg

@u

00

j

@u

00

j

+

@

@c

(< > C

c

!(c ec) eg)

@

@c

(< > S(c)eg);j;k = 1;::;3;(5)

where < > is the mean density,< > is the mean turbulent energy dissipation rate,!is

the mean turbulent frequency,C

0

,C

1

and C

c

are empirical constants.The left - hand side of

(5) represents the temporal change and the advection of eg by the total velocity

e

U+u

00

.The

rst term on the right - hand side is the Reynolds stress tensor,the second is the generalized

Langevin model (GLM) which is used to model the in uence of molecular velocity diusion and

uctuating pressure.The third term on the right - hand side has its origin in the stochastic

force,the fourth term describes the eects of molecular mixing of the scalar c which is modeled

using the IEM model with scalar frequency!,and the last - the eects of chemical reaction.

Chemical reactions are treated exactly,and no closer is needed.This is a main advantage of

the PDF approach.The equation (5) need be completed by the mean continuity equation and

momentum equation.

2 SPDEs

We consider the conservative form of the SPDEs statistically equivalent to the equation

of the PDF (5).The uctuating velocity equation is given by

@ru

00

i

@t

+

@r(

f

U

k

+u

00

k

)u

00

i

@x

k

=

r

< >

@

ru

00

k

u

00

i

@x

k

1

T

L

ru

00

i

@

e

U

i

@x

k

ru

00

k

+r

p

C

0

< >

_

i

(t);k;i = 1;::;3:(6)

The equation for the mass fractions is

@rY

j

@t

+

@r

f

U

k

+u

00

k

Y

j

@x

k

= rC

Y

1

T

L

Y

j

e

Y

j

+S

j

(Y

l

;T);j;l = 1;::;N

esp

:(7)

The equation for the total enthalpy is formulated as

@rh

t

@t

+

@r

f

U

k

+u

00

k

h

t

@x

k

=

r

< >

@ < p >

@t

rC

h

1

T

L

h

t

e

h

t

:(8)

Here,r is a stochastic density (dierent from ),T

L

is a characteristic time of turbulence,

!=

1

T

L

.W

i

are standard independent Wiener noises,

_

i

=

dW

i

dt

are its time derivatives,i.e.

white noises.Y

j

are stochastic mass fractions,h

t

is a stochastic total enthalpy,c

h

and c

Y

are

model constants.Mean pressure is dened by the following formula:

< P >=< rR

0

T

N

esp

X

k=1

Y

k

M

k

>

s

;(9)

where <>

s

is the averaging operator related to

i

,R

0

is an absolute gas constant,M

k

is a

molecular weight of the species k.

It can be shown that the transport equation of the PDF eg

r

of stochastic elds (which

is dened from system of SPDEs (6),(7) and (8)) is identical to the equations (5),if the two

following conditions are veried:

< r >

s

=< > (10)

3

and

<

r

>

s

= 1:(11)

The rst condition (10) guarantees consistency at the level of Favre statistics,while the second

(11) guarantees consistency at the level of Reynolds statistics.

3 Numerical aspects and test cases

The EMC method for the joint mass - fraction,total enthalpy and uctuating velocity

PDF has been implemented into the CHARME solver of CEDRE.The CHARME allows to

solve hyperbolic equations.The SPDEs (6),(7) and (8)) are solved by schemes of Godunov

type in space and deterministic Runge - Kutta schemes in time.

It is known that the integration in time of SPDEs by deterministic explicit Runge - Kutta

methods reduces the formal weak temporal accuracy of time - space methods to rst order

globally,independently of the spatial operator [6].Only a generalization of the Euler method

for ordinary dierential equations to stochastic dierential equations,it is named as Euler {

Maruyama method,computes the SPDEs with theoretical temporal accuracy [12].

The issue of boundary conditions is a crucial one.The hyperbolic Euler equations require

temperature,mass - fractions,velocity at inlet of the domain and pressure at outlet.As can

be noted the SPDEs (6),(7) and (8) do not have a stochastic eld of pressure,and so all the

parameters of SPDEs need be specied at in ow boundaries (including the density).In this work

the uctuating velocity at in ow boundaries is xed at the time,and at out ow boundaries the

values of stochastic elds are computed from the interior of the domain.In general case when

the uctuating velocity changes the sign at the boundaries,the entrance of the domain can

become the exit,and inverse.

The solver EMC consists in resolving of the RANS equations with the turbulent model kl

and then the SPDEs which use the Favre averaged velocity,frequency and turbulent dissipation

obtained from the RANS equations.

In order to validate the solver EMC,we consider a one dimensional turbulent reactive

ow.C

1

= C

0

= C

h

= C

Y

= 1 is used.The turbulent frequency and the turbulent dissipation

are assumed to be constant.Within these assumptions,we perform several tests.

An auto - ignition problem.An inlet,a stoichiometric methane - air mixture is

injected,with temperature T

in

= 1500 K.The inlet velocity is random and follows a Gaussian

law with mean

e

U

in

= 10:4 m s

−1

and variance

f

u

002

in

= 1 m

2

s

−2

.At outlet,the pressure,which

is present in the momentum conversation equation,is xed at 1 bar.Instead of pressure,for

the SPDEs the density r is xed at inlet.The length of the 1D domain is L = 0:035 m.The

methane - air chemistry is dealt with a simple one step global reaction,taken from [10].

The calculation shows that in the case of the turbulent uctuations the auto - ignition

length increases in comparison of the case without uctuations (left gure 1).The convergence

rates uctuate around the N

0:5

r

theoretical value (right gure 1),where N

r

is a number of

realizations of stochastic elds.

A passive scalar Y transport in homogeneous stationary turbulence.At initial

time,the velocity proles are constant in space and randomly distributed.They follow a Gaus-

sian law with mean

e

U = 0 m s

−1

and variance

f

u

002

= 1 m

2

s

−2

.The turbulent frequency is set to

!= 20 s

−1

.The initial scalar proles are Heaviside functions all centered at the middle of the

domain.The stochastic density is constant in space,and is set to 0:22 kg.As a result,there is

no initial scalar variance or scalar ux.The domain has length L = 1 m.Neumann boundary

conditions are imposed at each end of the domain.

It can be shown that in the limits when the the time tends to innity,the diusion

coecient

4

Figure 1:The auto - ignition test in 1D(the CEDRE).Numerical schemes:the 2nd order Runge

- Kutta scheme in time and the 2nd order spatial\ODFI\scheme.The mesh si uniform with

300 internal cells;dt = 1 10

6

s.Mean temperature at T = 0:02 s (left),statistical convergence

of the mean temperature at T = 0:02 s (right).

t

=

]

u

00

Y

00

@

e

Y

@x

(12)

tends,as a rst approximation,to

d

=

f

u

002

2!

(13)

and,consequently,does not depend on x.The purpose of this test is to check that the asymptotic

behavior is obtained with the numerical scheme proposed in this work.

The solution obtained by the proposed numerical scheme conrms the hypothesis of the

asymptotic behavior.Figure 2 shows the temporal behavior of

t

at the center of the mixing

zone.It is seen that

t

tends to its theoretical limit value

d

after approximately!t

a

= 3.

Figure 2 compares the spatial prole of the variance of the passive scalar against the theoretical

approximated solution at time!t = 5.The approximated theoretical solution is given by:

g

Y

002

=

1

4!t

exp

(x x

0

)

2

2

d

t

;(14)

where x

0

= 0:5 m.The Favre mean and variance of the mass fraction converge to its analytical

solutions with rate around 0:5.

The tests proposed by V.Sabelnikov and O.Soulard in [2] for SPDEs (6),(7) and (8)

such as Riemann problem,return to Gaussiannity,passive scalar transport and auto-ignition

have been used to validate the solver EMC.But as it will be shown in the following section,in

spite of the fact that these tests yield theoretical results,the chosen numerical methods are not

suited to the SPDEs.

4 Diculties and methodology of its resolution for SPDEs

In the two - dimensional case,the validation tests have been conducted on the conguration

"backward facing step"for reactive and non - reactive ows.The physical domain is L = 1 m

5

Figure 2:The turbulent diusion test in 1D calculated by the EMC.Numerical schemes:the

2nd order stochastic Runge - Kutta in time and the 2nd order spatial (Godunov - Kolgan) scheme

(Matlab).The mesh is uniform with 100 internal cells.Evolution of the diusion coecient in

time (left),Favre variance of mass fraction at!T = 5 (right).

horizontally and H = 0:1 m vertically.The height of the step,placed at the lower wall,is

h = 0:035 m and its extremity is located at 0:1 m inside the computational domain.Positions

will hereafter be given using the step extremity as the origin.

At inlet,the prole of the mass fraction of methane is a gaussian function centered at

y = 0:0675 m with maximal value Y

CH

4

= 1.It is located in the interval y 2 [0:041 m;0:094 m],

outside this interval,the mass fraction of methane is zero.The mass fraction of nitrogen at inlet

is Y

N

2

= 1 Y

CH

4

.The mass fractions are injected at U = 58 m s

−1

and T = 525 K.The inlet

values of the turbulent quantities are K = 60 m

2

s

−2

and l = 0:01 m.At outlet,the pressure for

the RANS equations is xed P = 1 bar,for the SPDEs the density is xed at inlet.The upper

and lower walls are considered adiabatic.In the case of reactive ow the methane combustion in

air is modeled by a global single - step chemical reaction [10].At inlet,a methane - air mixture

is injected at U = 58 m s

−1

and T = 525 K,with a stoichiometrc equivalence ratio = 1.In the

interior of the domain the temperature is set to T = 1600 K.

In gure 3,the fraction of the methane is presented at time on the order of 10,i.e.at

L

U

0

= 0:15 s.Although the visual appearance of left and right gures of (3) is close enough,the

turbulent diusion calculated with the EMC are signicantly lower in comparison the RANS.

It persists also in the case of the reactive ow.We assume that the numerical schemes are not

adapted to solve the SPDEs when the turbulent intensity is large,e.g.in the recirculation

zone.These results led us to examine an another test of one dimensional ow.

Let us consider the case when the Favre averaged velocity is lower or has the same order

that the uctuating velocity.For illustration purposes,the Favre averaged velocity is set to zero

e

U = 0 m s

−1

.The transport equation for the density (a continuous eld is assumed) can be

written as

@r

@t

+r

@u

00

@x

+u

00

@r

@x

= 0;(15)

or

Dr

Dt

= u

00

@r

@x

:(16)

6

Figure 3:Calculation of methane in the injection chamber of conguration"backward faing

step"at the temperature of 525 K.Euler method with the spatial rst order scheme;2995

internal cells.Mass fraction of methane obtained by RANS at 0:15 s (left),Favre averaged eld

of mass fraction of methane obtained by EMC with 60 elds at 0:15 s (right).

where

D

Dt

is a total derivative.Along the trajectory,the stochastic density is governed by the

exponential function of the gradient of the uctuating velocity:

r exp

2

4

t

Z

0

@u

00

@x

d

3

5

:(17)

Calculations show that the gradient of the uctuating velocity

@u

00

@x

can be positive for a long time

and so the stochastic density tends to zero as represents (17).Let us illustrate this possibility

considering the following Riemann problem:

A vacuum test.This test demonstrates that r = 0 can be theoretically obtained for

the SPDEs.The Favre averaged velocity is zero at all domain.The initial uctuating velocity

prole is a Heaviside function H(x),i.e.

u

00

(x;t = 0) = u(x;t = 0) = H(x) =

(

0;if x < 0,

1;if x 0.

(18)

The density is constant in space and is taken as one.All sources are set to zeroes.The domain

is innite.

At the time t = T

> 0 the density is zero for x 2 [x

0

;x

0

+T

] (right gure 4).Numerically,

when we use hyperbolic schemes of the CHARME and when conservative variables are calculated,

the density r = 0 implies innite uctuating velocity:

ru

r

!1.The system of the SPDEs

becomes unstable,and the calculation fails.To overcome this diculty it is necessary to consider

schemes which are able to capture delta - shocks of density concentrations and treat data with

vacuum.One example of these schemes is given by F.Bouchut [8].

The Riemann problem with H(x) initial prole of velocity demonstrates a multi - valued

solution.Any entropic numerical scheme cannot ensure the correspondence between the SPDEs

(6),(7),(8) and the equation (5) of the Favre PDF eg,because the rst cannot capture the multi

- valued solution.

A test of multi - valued solution.In one dimension case we suppose a simplied

equation PDF for the velocity u.It is given by

7

Figure 4:The vacuum test.Initial prole of the uctuating velocity (left);prole of uctuating

velocity after time T

> 0 (right).

@ < >

e

f

@t

+

@ < > u

e

f

@x

=

@ < >

e

f

@t

+u

@ < >

e

f

@x

= 0;(19)

At initial time

e

f(u;x;t = 0) = H(x)(u 1) +H(x)(u);t = 0;(20)

< (x;t = 0) >= 1:(21)

The initial prole of velocity which corresponds to the PDF (19) is following

u(x;t = 0) = u

0

(x) = H(x):(22)

The solution of the Cauchy problem (19),(20) and (21) is

< >

e

f(u;x;t) = H(t x)(u 1) +H(x)(x);t > 0;(23)

where the density is dened by

< > (x;t) = H(t x) +H(x) =

(

1;if x < 0 and x > t,

2;if x < 0 t.

(24)

The equations (20),(21) and (23) show that at initial time t = 0 each point in space is associated

with one Dirac (u 1),if x < 0 and (u),if x 0 (left gure 5).With time,the information

associated with (u 1) propagates at the speed u = 1 and another associated with (u) does

not move.As a consequence,after the time T

> 0,the interval 0 < x T

contains two Dirac

(u 1) and (u) (gure 5,center).

Let us analyze the solution of the SPDEs which is chosen by the respecting the entropy

increase condition on the discontinuity shocks.The SPDEs statistically equivalent to (19),(20)

and (21) can be dened as

@r

@t

+

@ru

@x

= 0;(25)

@ru

@t

+

@ru

2

@x

= r

@u

@t

+ru

@u

@x

= 0;(26)

u(x;t = 0) = H(x);r(x;t = 0) = 1:(27)

The equation (26) can be rewritten as

8

Figure 5:The test of multi-valued solution.Initial prole of the velocity (left);solution of the

PDF (19) at t = T

(center);solution of the SPDEs by scheme of Godunov type at t = T

(right).

@u

@t

+u

@u

@x

= 0;r > 0:(28)

taking into account that r > 0.The solution of the SPDEs (28) and (27) is (right gure 5)

u(x;t) = H

t

2

+x

;t > 0:(29)

The Favre PDF corresponding to (29) is

e

f(u;x;t) = H

t

2

x

(u 1) +H

t

2

+x

(u):(30)

We can see that each point in space is associated with one Dirac:(u 1),if x <

t

2

and (u),

if x

t

2

.The solution (23) is not the same that the solution (29).So,the method used before

for the SPDEs cannot give a correct solution of the PDF.

A new path to solve the SPDEs (6),(7) and (8) consists in using a Runge - Kutta method

of the second order for numerical solution of stochastic dierential equations [11] and a second

- order kinetic scheme in space [8].We propose to neglect the multi - valued solution of the

SPDEs and nd the statistics of the joint velocity - turbulent scalars PDF using the methods

mentioned above.Nowadays the problem of numerical nding of multi - valued solution for

weakly hyperbolic equations is still open.

Conclusion

In this work,we show that the numerical schemes of Godunov type are not suited to solve

the SPDEs (6),(7) and (8).The numerical integration of the corresponding SPDEs is discussed

and a numerical spatial scheme is proposed.Further developments will include the validation of

this method for solving the SPDEs.

References

1.Pope S.B.,"PDF methods for turbulent reactive ows".Progress in Energy and Com-

bustion Science 27,119 - 192.

2.Sabel'nikov V.,Soulard O.,"Eulerian Monte Carlo methods for solving PDF transport

equations in turbulent reacting ows",Handbook of Combustion,Vol.3:Gases and

Liquids,p.75 - 119,2010.

9

3.Soulard O.,"Prise en compte d'un spectre d'echelles turbulentes dans la modelisation du

micro melange et elaboration d'une methode de MCE",PhD thesis,Universite de Rouen,

2005.

4.Ourliac M.,"Utilisation d'une methode eulerienne de resolution des equations PDF de

scalaires turbulents pour la simulation de l'allumage d'une chambre de combustion cry-

otechnique",PhD thesis,Universite de Lyon,2009.

5.Lemons D.S.,"Noise { induced instability in self { consistent Monte Carlo calculations",

Physical Rewiew,Vol.52,number 6,december 1995.

6.A.Donev,E.Vanden-Eijnden,A.Garcia,and John B.Bell,"On the Accuracy of Ex-

plicit Finite-Volume Schemes for Fluctuating Hydrodynamics",Communications in Ap-

plied Mathematics and Computational Science,5(2):p.149 - 197,2010.

7.S.B.Pope."Turbulent ows",Cambridge Univ.Press,2000.

8.Francoit Bouchut,Shi Jin,Xiantao Li"Numerical approximations of pressureless and

isothermal gas dynamics",SIAM J Numer.Anal.,2001.

9.Eleuterio F.Toro."Riemann solvers and numerical methods for uid dynamics",Springer,

2009.

10.C.K.Westbrook and F.L.Dryer."Chemical kinetic modeling of hydrocarbon combus-

tion",Prog.Energy Comb.Dci.,p.10:1 - 57,1984.

11.A.Tocino,R.Ardanuy."Runge - Kutta methods for numerical solution of stochastic

dierentail equtions",Journal of Computational and Applied Mathematics 138,p.219 -

241,2002.

12.P.E.Kloeden,E.Platen,and H.Schurz,"Numerical solution of SDE through computer

experiments",Springer,1994.

10

## Comments 0

Log in to post a comment