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"
OU: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 onepoint 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 autoignition
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 multivalued 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.VandenEijnden,A.Garcia,and John B.Bell,"On the Accuracy of Ex
plicit FiniteVolume 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
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο