Simulation of turbulent combustion in aerobic chambers: extension of Eulerian Monte Carlo methods

monkeyresultMechanics

Feb 22, 2014 (3 years and 5 months ago)

73 views

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 Reactifs pour l'

Energetique
(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 Dierential 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 dierent
numerical schemes for solving the stochastic partial dierential 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 dierent 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
eort 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/diusion
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 veried,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 dened 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 diusion and
uctuating pressure.The third term on the right - hand side has its origin in the stochastic
force,the fourth term describes the eects of molecular mixing of the scalar c which is modeled
using the IEM model with scalar frequency!,and the last - the eects 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 (dierent 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 dened 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 dened from system of SPDEs (6),(7) and (8)) is identical to the equations (5),if the two
following conditions are veried:
< 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 dierential equations to stochastic dierential 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 specied 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 kl
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 proles 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 proles 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 innity,the diusion
coecient
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 conrms 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 prole 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 Diculties and methodology of its resolution for SPDEs
In the two - dimensional case,the validation tests have been conducted on the conguration
"backward facing step"for reactive and non - reactive ows.The physical domain is L = 1 m
5
Figure 2:The turbulent diusion 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 diusion coecient 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 prole 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 diusion calculated with the EMC are signicantly 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 conguration"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
prole 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 innite.
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 innite uctuating velocity:
ru
r
!1.The system of the SPDEs
becomes unstable,and the calculation fails.To overcome this diculty 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 prole 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 simplied
equation PDF for the velocity u.It is given by
7
Figure 4:The vacuum test.Initial prole of the uctuating velocity (left);prole 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 prole 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 dened 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 dened 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 prole 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 dierential 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 modelisation du
micro melange et elaboration d'une methode de MCE",PhD thesis,Universite de Rouen,
2005.
4.Ourliac M.,"Utilisation d'une methode eulerienne de resolution des equations PDF de
scalaires turbulents pour la simulation de l'allumage d'une chambre de combustion cry-
otechnique",PhD thesis,Universite 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.Francoit 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
dierentail 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