Center for Turbulence Research 91

Annual Research Briefs 2000

Stochastic modeling of scalar dissipation rate

fluctuations in non-premixed turbulent
combustion
By Heinz Pitsch AND Sergei Fedotov†

1. Motivation and objectives

In non-premixed combustion chemical reactions take place when fuel and oxidizer
mix on a molecular level. The rate of molecular mixing can be expressed by the scalar
dissipation rate, which is for the mixture fraction Z defined as
χ = 2DZ (∇Z)2 , (1.1)
where DZ is the diffusion coefficient of the mixture fraction. The scalar dissipation rate
appears in many models for turbulent non-premixed combustion as, for instance, the
flamelet model (Peters (1984), Peters (1987)), the Conditional Moment Closure (CMC)
model (Klimenko & Bilger (1999)), or the compositional pdf model (O’Brien (1980), Pope
(1985)). In common technical applications, it has been found that if the scalar dissipation
rates are much lower than the extinction limit, fluctuations of this quantity caused by
the turbulence do not influence the combustion process (Kuznetsov & Sabel’nikov (1990),
Pitsch & Steiner (2000)). However, it has been concluded from many experimental (Saitoh
& Otsuka (1976)) and theoretical studies (Haworth et al. (1988), Mauss et al. (1990),
Barlow & Chen (1992), Pitsch et al. (1995)) that there is a strong influence of these
fluctuations if conditions close to extinction or auto-ignition are considered. For instance,
in a system where the scalar dissipation rate is high enough to prohibit ignition, random
fluctuations might lead to rare events with scalar dissipation rates lower than the ignition
limit, which could cause the transition of the whole system to a burning state.
In this study, we investigate the influence of random scalar dissipation rate fluctua-
tions in non-premixed combustion problems using the unsteady flamelet equations. These
equations include the influence of the scalar dissipation rate and have also been shown
to provide very reasonable predictions for non-premixed turbulent combustion in a va-
riety of technical applications (Pitsch et al. (1996), Pitsch et al. (1998), Barths et al.
(1998)). However, it is clear that these equations are actually not capable of describing
all of the features which might occur in turbulent non-premixed flames. For instance,
in jet diffusion flames, local extinction events might occur close to the nozzle because
of high scalar dissipation rates. These extinguished spots might reignite downstream,
not by auto-ignition, but by heat conduction and diffusive mass exchange with the still
burning surroundings. It should be kept in mind that the motivation in this work is not
to predict actual turbulent reacting flows, but to study the dynamical system defined
by the equations described in the following section. The advantage of the present sim-
plified approach allows a study of the extinction process isolated from auto-ignition and
re-ignition events.
† UMIST, Manchester, England
92 Pitsch & Fedotov
The basic purpose of this paper is to analyze how random fluctuations of the scalar dis-
sipation rate can affect extinction of non-premixed combustible systems. The approach,
based on stochastic differential equations, allows us to take random extinction events
into account. In this case the critical conditions must be different from those involv-
ing deterministic situations. Here, we look at these phenomena in terms of noise-induced
transitions theory, where multiplicative noise of sufficient intensity can drastically change
the behavior of a system (Horsthemke & Lefever (1984)). In the present case, the proba-
bility density function for the temperature in the reaction zone may undergo qualitative
changes as the intensity of random fluctuations increases. It should be noted that a similar
analysis has been done in a series of works on the stochastic analysis of thermal ignition
of explosive systems in (Buyevich et al. (1993), Fedotov (1992), Fedotov (1993)).
Oberlack et al. (1999) have investigated the influence of Damköhler number fluctua-
tions in a well-stirred reactor. The fundamental difference compared to the present study
is that in a well stirred reactor the mixing process is assumed to be infinitely fast. The
Damköhler number, therefore, represents the residence time rather than the mixing time
and appears in the non-dimensional chemical source term. Hence, imposing stochastic
variations of the Damköhler number corresponds to a fluctuating chemical source term.
Here, however, the fluctuating quantity is the scalar dissipation rate, which appears as
a diffusion coefficient. The response of the mixing field to this fluctuating diffusion coef-
ficient and the interaction with the chemical source term are investigated. Moreover, in
the present formulation we allow for temporal changes of the fluctuating quantity and
also consider its pdf.
In this paper, we will first present the non-dimensional flamelet equations for a one-step
global reaction. With this assumption the system can be reduced to a single equation for
the temperature. We will then derive stochastic differential equations for the temperature
and the scalar dissipation rate. These equations lead to a partial differential equation for
the joint probability density function of the temperature and the scalar dissipation rate.
This equation will be discussed and numerical solutions will be presented.

2. Governing equations
2.1. Flamelet equations
Assuming an irreversible one-step reaction of the form νF F + νO O → P, where F, O, and
P denote fuel, oxidizer, and reaction product, respectively, the flamelet equations for the
mass fractions of fuel YF , oxidizer YO , reaction product YP , and the temperature T , can
be written as
∂Yi χ ∂ 2 Yi
− + νi Wi w = 0 , i = F, O, P (2.1)
∂t 2 ∂Z 2
∂T χ ∂2T Q
− 2
− w = 0. (2.2)
∂t 2 ∂Z cp
Here, νi are the stoichiometric coefficients, Wi the molecular weights, t is the time, ρ the
density, cp the specific
P heat capacity at constant pressure, and Q is the heat of reaction
defined as Q = − νi Wi hi , where hi denotes the enthalpy of species i. The mixture
i
fraction Z is defined as
ν̂YF − YO + YO,2 νO WO
Z= with ν̂ = , (2.3)
ν̂YF,1 + YO,2 νF WF
 Stochastic modeling of scalar dissipation rate fluctuations 93
where the subscripts 1 and 2 refer to the conditions in the fuel stream and the oxidizer
stream, respectively.
The parameter χ appearing in Eqs. (2.1) and (2.2) is the scalar dissipation rate, which
has already been defined by Eq. (1.1). The reaction rate per unit mass w is given by
 
YF YO E
w=ρ A exp − , (2.4)
WF WO RT
where A is the frequency factor and E the activation energy of the global reaction,
respectively. R is the universal gas constant.

2.2. Non-dimensionalization
In order to investigate the flamelet equations with respect to the relevant non-dimensional
parameters, it is convenient to introduce the non-dimensional temperature θ and mass
fractions of fuel yF , oxidizer yO , and reaction product yP as
T − Tst,u YF YO νF WF YP
θ= , yF = , yO = , yP = (ν + 1) , (2.5)
Tst,b − Tst,u YF,st,u YO,st,u νP WP YF,1
where the index st refers to stoichiometric conditions and the unburnt values of temper-
ature, fuel, and oxidizer at stoichiometric conditions are given by
Tst,u = T2 + (T1 − T2 ) Zst , Yi,st,u = Yi,2 + (Yi,1 − Yi,2 ) Zst , i = F, O . (2.6)
The adiabatic temperature for complete conversion of fuel to products is
L YF,1 Q YF,1
Tst,b = Tst,u + , L= ν = ν̂ . (2.7)
cp WF νF (ν + 1) YO,2
With these definitions and Eq. (2.3), the mixture fraction can be expressed as
νyF − yO + 1
Z= , (2.8)
ν +1
from which the stoichiometric mixture fraction
1
Zst = (2.9)
ν+1
follows .
The non-dimensional time τ is given by
χst,0 2∆Zν
τ= t with a = 2∆Z · Zst (1 − Zst ) = , (2.10)
a (1 + ν)2
where the reference value for the scalar dissipation rate χst,0 and the parameter ∆Z will
be introduced below.
The non-dimensional scalar dissipation rate x, the Damköhler number Da, and the
Zeldovich number Ze are defined as
χ ννF aρst,u YO,2 A
x= , Da = exp (−βref ) (2.11)
χst,0 (ν + 1) WO χst,0

Tst,b − Tu E
Ze = αβ , α=, β= . (2.12)
Tst,b RTst,b
With the assumption of constant molecular weight of the mixture, the density ρ can
94 Pitsch & Fedotov
be expressed in terms of the non-dimensional temperature θ as
(1 − α)
ρ= ρst,u . (2.13)
1 − α(1 − θ)
Introducing these definitions into the flamelet equations, Eqs. (2.1) and (2.2), yields
∂yF ax ∂ 2 yF 1
− 2
+ ω=0 (2.14)
∂τ 2 ∂Z ν +1
∂yO ax ∂ 2 yO ν
− 2
+ ω=0 (2.15)
∂τ 2 ∂Z ν +1
∂yP ax ∂ 2 yP
− −ω =0 (2.16)
∂τ 2 ∂Z 2
∂θ ax ∂ 2 θ
− −ω = 0, (2.17)
∂τ 2 ∂Z 2
where the non-dimensional chemical source term is given by
 
(ν + 1)2 (1 − α) exp (βref − β) 1−θ
ω = Da yF yO exp −Ze . (2.18)
ν 1 − α(1 − θ) 1 − α(1 − θ)
The boundary conditions for Eqs. (2.14) – (2.17) are
Z =0 : yF,2 = 0 , yO,2 = 1 , yP,2 = 0 , θ2 = 0 (2.19)

Z=1 : yF,1 = 1 , yO,1 = 0 , yP,1 = 0 , θ1 = 0 . (2.20)

2.3. Coupling functions
Adding Eqs. (2.14), (2.15), and (2.17) yields a conservation equation for yF + yO + θ as
∂ ax ∂ 2
(yF + yO + θ) − (yF + yO + θ) = 0 . (2.21)
∂τ 2 ∂Z 2
The boundary conditions for the conserved scalar can be determined from Eqs. (2.19)
and (2.20) to be unity at both sides. Then with the unburnt state as initial condition,
the solution of Eq. (2.21) is given by
yF + yO + θ = 1 . (2.22)
Note that this particular choice of the initial condition does not restrict the solution since
it is a requirement of every possible physical initial condition that it has to be realizable
from the unburnt state. Since the non-dimensional product mass fraction yP and tem-
perature θ are governed by a mathematically similar flamelet equation (Eqs. (2.16) and
(2.17)) and have the same boundary and initial conditions, it follows that
yP = θ , (2.23)
which shows that Eq. (2.22) represents the mass conservation condition.
With Eq. (2.22) and the definition of the mixture fraction, Eq. (2.8), the mass fractions
of fuel and oxidizer can be expressed in terms of mixture fraction and temperature as
ν
yO = 1 − Z − θ = 1 − Z − (1 − Zst ) θ (2.24)
ν +1
θ
yF = Z − = Z − Zst θ (2.25)
ν +1
 Stochastic modeling of scalar dissipation rate fluctuations 95
and the chemical reaction rate, defined in Eq. (2.18), as
    
(1 − α) exp (βref − β) Z 1−Z 1−θ
ω = Da −θ − θ exp −Ze .
1 − α(1 − θ) Zst 1 − Zst 1 − α(1 − θ)
(2.26)
With Eq. (2.26) the flamelet equation for the non-dimensional temperature given by
Eq. (2.17) depends only on the temperature itself and can be integrated without solving
the equations for the mass fractions of fuel, oxidizer, and product. If desired, these can
be computed from Eqs. (2.23), (2.24), and (2.25).
2.4. Stochastic differential equations
In this section we want to derive an equation for the joint pdf of the temperature and the
scalar dissipation rate. To complete Eq. (2.17) we need a stochastic differential equation
(SDE) that governs the evolution of the scalar dissipation rate.
We consider a Stratonovich SDE given by Horsthemke & Lefever (1984)
dχst = f (χst ) dt + σϕ (χst ) ◦ dW (t) , (2.27)
where W (t) denotes a Wiener process. In Eq. (2.27) the first term on the right-hand side
is a drift term, the second a random term. The stationary probability density function
corresponding to the Stratonovich SDE can be found to be
 χst 
Z
N 2f (z)
ps (χst ) = exp  dz  . (2.28)
σϕ (χst ) σ2 ϕ2 (z)
0

It is well known that a good approximation for a stationary pdf of χst (t) is a lognormal
distribution (Peters (1983)) given as
!
2
1 (ln χst − ln χst,0 )
ps (χst ) = √ exp − , (2.29)
χst 2πσ2 2σ2
from which it can easily be shown that the mean value of χst is
Z ∞  2
σ
χst = χst ps (χst ) dχst = χst,0 exp . (2.30)
0 2
To find f (χst ) and ϕ (χst ), one needs to equate Eqs. (2.29) and (2.28). From this we
obtain
s
χst 2 1
f (χst ) = − (ln χst − ln χst,0 ) , ϕ (χst ) = χst , N=√ . (2.31)
tχ tχ πtχ
For dimensional reasons a characteristic time scale tχ has been introduced, which appears
as a parameter of the problem. This time scale is associated with the time to reach a
statistically stationary state. Therefore, it does not appear in the stationary pdf given
by Eq. (2.29). The scalar dissipation rate χ (t) then obeys the following SDE
s
χst 2
dχst = − (ln χst − ln χst,0 ) dt + σ χst ◦ dW (t) . (2.32)
tχ tχ
In non-dimensional form, this equation can be rewritten as
r
xst 2
dxst = − ln xst dτ + σ xst ◦ dW (τ ) . (2.33)
δ δ
96 Pitsch & Fedotov
Here, δ = tχ χst,0 /a represents the ratio of the characteristic time scales of Eqs. (2.33)
and (2.2). In a turbulent flow, the time scale tχ would be modeled by the integral time
scale of the turbulence or the scalar (Pope (2000)). Hence, tχ can be expressed as

C0 Z 02
tχ = 2
(2.34)
χst,0 exp σ2
from which follows that
C0 Z 02
δ=
, (2.35)
2∆Z · Zst (1 − Zst )exp σ2
2

where C0 is a constant and Z 02 is the mixture fraction variance. This shows that δ is
independent of the mean scalar dissipation rate. Here, δ = 1 will be assumed, which for
C0 = 1, Zst = 0.5, and σ = 1 roughly corresponds to Z 0 = 0.2.
From a mathematical point of view, Eq. (2.17) with the source term Eq. (2.26) and
the random scalar dissipation rate is a nonlinear stochastic partial differential equation,
which can be solved but is very difficult to work with analytically. One way to analyze this
equation is to derive the corresponding equation for the probability density functional
for the temperature distribution θ(Z) (Fedotov (1992), Fedotov (1993)). However, since
the random parameter x(τ, Z) appears in Eq. (2.17) as a multiplicative noise, it would
be very difficult to obtain reasonable results. In order to simplify the problem, we will
derive ordinary stochastic differential equations (SDE) for these quantities by modeling
the diffusion term in Eq. (2.17).
It has been shown by Peters (1983) that these linear temperature profiles in the outer
non-reactive structure can be found as the first order solution of an asymptotic analysis
of the flamelet equations assuming one-step global chemistry. The assumption of linear
temperature profiles in the outer structure will now be used for an approximation of the
diffusion term appearing in Eq. (2.17).
The diffusion term evaluated at stoichiometric conditions can be written as a finite
difference approximation over the reaction zone of width ∆Z as
 + − !
∂ 2 T  1 ∂T  ∂T 
≈ − . (2.36)
∂Z 2  Zst ∆Z ∂Z  ∂Z 

If the temperature gradients appearing in this expression are evaluated with the as-
sumption of linear profiles in the non-reactive diffusion zones, the diffusion term can be
approximated as
  
∂ 2 T  1 Tst − T1 T2 − Tst Tst − Tst,u (Tst,b − Tst,u )
 ≈− − =− =− θst .
2
∂Z Zst ∆Z 1 − Zst Zst ∆ZZst (1 − Zst ) ∆ZZst (1 − Zst )
(2.37)
Here, it has to be assumed that the reaction zone thickness ∆Z is independent of the
scalar dissipation rate. Then, ∆Z is a constant which appears in the Damköhler number.
The actual choice of ∆Z then only changes the value of the Damköhler number and is
of no consequence for the conclusions of the paper. The validity of this assumption has
been numerically evaluated by Cha (2000).
Introducing Eqs. (2.37) and (2.26) into Eq. (2.17) formulated at Z = Zst yields
dθst
+ x (τ ) θst − ω (θst ) = 0 (2.38)
dτ
 Stochastic modeling of scalar dissipation rate fluctuations 97
with
 
(1 − α) exp (βref − β) 2 1 − θst
ω = Da (1 − θst ) exp −Ze . (2.39)
1 − α(1 − θst ) 1 − α(1 − θst )

2.5. Joint probability density function

Now we are in a position to analyze how random fluctuations of the scalar dissipation
rate can influence the non-premixed combustion process. It follows from Eqs. (2.38) and
(2.33) that the pair process (θst (τ ) , xst (τ )) is Markovian, and therefore their probability
density function p = p (τ, xst , θst ) is governed by the Fokker-Planck equation

∂p 1 ∂

∂   σ 2 ∂2

− ln xst − σ2 xst p + (−xθst + w (θst )) p = 2 x2st p (2.40)
∂τ δ ∂xst ∂θst δ ∂xst
with 0 < xst < ∞, 0 < θst < 1, and the boundary conditions
p (τ, 0, θst ) = p (τ, ∞, θst ) = p (τ, xst , 0) = p (τ, xst , 1) = 0. (2.41)

It is convenient to introduce the natural logarithm of the stoichiometric scalar dissi-

pation rate as a new independent variable
xln = ln xst . (2.42)
The pdf of xln can then be obtained by the normalization condition
plnx
p(xst ) = (2.43)
xst
and Eq. (2.40) can be written as
∂plnx 1 ∂ ∂   σ2 ∂ 2 p
lnx
− (xln plnx ) + (−exln θst + ω) plnx − = 0. (2.44)
∂τ δ ∂xln ∂θst δ ∂x2ln
The boundary conditions are given by
plnx (τ, −∞, θst ) = plnx (τ, ∞, θst ) = plnx (τ, xln , 0) = plnx (τ, xln , 1) = 0. (2.45)
Note that, as shown by Eq. (2.43), the distribution plnx is different from p and the
maximum will in general be at a different value of the scalar dissipation rate. However,
since both functions can easily be converted into each other, the conclusions do not
depend on the choice of the formulation used for the analysis.

3. Numerical solution
Equation (2.44) has been solved numerically using a 4th order Runge-Kutta scheme
with adaptive step-size control. The convection term in the xln -direction has been dis-
cretized using central differences, the convection term in the θst -direction by a robust,
globally second order upwind scheme as given by Koren (1996). The equations are solved
on a 300 × 300 equidistant grid. The numerical time-step is restricted by a CFL condi-
tion which is imposed by the high convection velocity in the θst -direction at high scalar
dissipation rate. This can be observed in Fig. 2, which will be described below. The ini-
tialization is performed with a numerical δ-function at some point in the xln –θst -space.
98 Pitsch & Fedotov
4. Results and discussion
In this section we will first provide a general discussion of Eq. (2.44) and the parameters
Da, Ze, and α appearing in this equation. Numerical solutions of Eq. (2.44) will then be
presented for a variation of the scalar dissipation rate variance σ, and the results will be
discussed.
Numerical solutions of Eq. (2.44) will then be presented. The results for different values
of the scalar dissipation rate variance σ will be discussed.

4.1. General discussion

Equation (2.44) is a two-dimensional unsteady partial differential equation depending
on xln and θst . In the direction of xln , the equation reveals a convective term and a
diffusion term. The convective term describes the relaxation to the mean. The mean
value is achieved when the convection velocity is zero. This implies that the mean value
of the non-dimensional scalar dissipation rate is xln = 0, which follows trivially from the
normalization of χst . However, it is interesting to note that only the scalar dissipation rate
itself determines the speed at which it relaxes to its mean. The diffusion term describes
the broadening of the pdf by fluctuations of the scalar dissipation rate with σ 2 appearing
as the diffusion coefficient.
In the direction of θst , Eq. (2.44) only reveals a convection term. Setting the convection
velocity Vθst = −exln θst + ω equal to zero yields the steady state relation between the
temperature and the scalar dissipation rate in the absence of scalar dissipation rate
fluctuations as
    
Ze 
2
Ze (1 − α) (1 − θst ) 1 − θst
x (θst ) = Da exp − exp −Ze . (4.1)
α ref α θst (1 − α(1 − θst )) 1 − α(1 − θst )
This relation describes the so called S-shaped curve for non-premixed combustion (Peters
(1984)), which depends on three parameters: the Damköhler number Da, the Zeldovich
number Ze, and the heat release parameter α, where Ze and α only depend on the
chemistry.
Figure 1 shows S-shaped curves from solutions of Eq. (4.1) for different values of these
parameters†. It is well known and will be shown in the following discussion that stable
solutions can only be achieved for the upper and the lower branch, whereas solutions
given by the middle branch are unstable. Considering the fact that the S-shaped curves
shown in Fig. 1 represent states with zero convection velocity in the direction of θst , it can
be seen from Eq. (2.44) that the convection velocity Vθst is positive for scalar dissipation
rates smaller than xln (θst ) as given from Eq. (4.1) and negative for larger values. The
consequence is that the convection velocity in the θ-direction is always directed away from
the intermediate branch, which shows that these solutions are unstable. It also shows that
starting from an unburnt solution, the scalar dissipation rate has to be decreased below
the value at the lower turning point of the curve to be able to auto-ignite. This value will,
therefore, be referred to as ignition scalar dissipation rate xig . Correspondingly, starting
from a burning solution, the flame can only be extinguished by increasing the scalar
dissipation rate over the value at the upper turning point. This value will, therefore, be
called extinction scalar dissipation rate xex .
For non-premixed methane flames, the activation energy of a one-step global reaction
† For constant scalar dissipation rate this relation would be plotted as function of the
Damköhler number, which would be proportional to x−1 ln . In the present paper these curves
are plotted over xln and therefore mirrored. However, we still use the phrase S-shaped curve.
 Stochastic modeling of scalar dissipation rate fluctuations 99
1

Da Ze α
0.8 100 4.91 0.679
10 4.91 0.679
100 6.95 0.679 Ze
0.6 100 6.95 0.866
θst

0.4 Vθst > 0

Da

0.2
α
Vθst < 0
0
-14 -12 -10 -8 -6 -4 -2 0 2
xln

Figure 1. S-shaped curve determined from Eq. (4.1) for different parameter variations.
Parameter changes are indicated by arrows.

can be assumed to be E = 150 kJ/kg (Seshadri (1999)). This implies a value of βref = 8.03
for a methane/air-system at ambient conditions. Then, the solid line in Fig. 1 corresponds
to a case with preheated air at T2 = 800 K and the dotted line to an air temperature
of T2 = 300 K. For both cases the fuel temperature is assumed to be T1 = 300 K and
the pressure to be 1 bar. It is clear from Eq. (4.1) and it can be seen in Fig. 1 that a
variation in the Damköhler number simply shifts the curve. In contrast, a variation of
the Zeldovich number leads to moderately lower scalar dissipation rate for extinction and
a strongly decreased ignition scalar dissipation rate.
The strongest influence however, can be seen by changing the heat release parameter.
Although by increasing α the extinction scalar dissipation rate is only slightly increased,
the ignition scalar dissipation rate is decreased very strongly to a value of approximately
xln,ig ≈ −40 corresponding to xig ≈ 10−17 for the example shown in Fig. 1. This merely
shows that auto-ignition of methane at ambient conditions is almost impossible.
Figure 2 shows a two-dimensional vector representation of the velocities of particles in
the θst -xln space, where the term particle is defined by a point and the associated velocity
in this space. This figure again clearly shows that the pdf tends to move to xln = 0 and
generally away from the unstable branch. However, at low temperature and low scalar
dissipation rate on the left side of the S-shaped curve, for instance, the driving force in
the direction of the mean scalar dissipation rate is so strong that particles might cross the
unstable branch. Even though these particles were initially in a regime which would for
constant xln lead to ignition, these particles will then be attracted by the lower branch.
In the present example this effect is not so obvious for particles originating from a
burning state with a scalar dissipation rate higher than the extinction limit, which would
be located in the upper right corner in Fig. 2. These particles can also change during the
extinction process to lower scalar dissipation rates and might cross the S-shaped curve.
This would lead to a recovery to the burning state. It has been discussed before and is
indicated in Fig. 1 that, in the absence of scalar dissipation rate changes, all particles
on the left side of the unstable branch of the S-shaped curve will change to the burning
100 Pitsch & Fedotov
1.00

0.75

row
θst

0.50

0.25

0.00
-4 -2 0 2 4
xcol
ln
Figure 2. Convection velocities in θst -xln space.
1.00
σ=0.5 σ=1 σ=2

0.75
θst

0.50

0.25

0.00
-4 -2 0 2 -4 -2 0 2 -4 -2 0 2
xln xln xln

0 12 25 38 0 1 2 4 5 6 8 9 0.0 0.2 0.4 0.6 0.8 1.0

Figure 3. Calculated joint θst -xln probability density function at τ = 5.

state, whereas particles on the right will change to the non-burning state. However, it has
clearly been demonstrated here that this is different in the case of a fluctuating scalar
dissipation rate, where the unstable branch does not uniquely separate these two regimes.

4.2. Numerical results

For the numerical solution of the equation for the joint pdf of θst and xln , Eq. (2.44),
the parameters appearing in this equation have been set to Da = 200, Ze = 4.91, and
α = 0.679. As mentioned above, this corresponds to a methane/air-system, where the
air is preheated to T2 = 800 K. Results of the numerical simulations at time τ = 5 are
presented exemplarily for σ = 1 in Fig. 3. All calculations have been started with a δ-
function at θst = 0.9 and xln = 0 as initial condition for the probability density function,
which is then given by p(τ = 0, xln , θst ) = δ(xln , θst − 0.9).
The joint pdfs of θst and xln are given in Fig. 3 for different values of the scalar
dissipation rate variance σ. It can be observed that even for the low variance case σ =
 Stochastic modeling of scalar dissipation rate fluctuations 101
0.5, the distribution of high probability density is already rather broad, extending from
approximately −1 < xln < 1 and mainly around the upper branch of the S-shaped curve.
Even though this cannot be seen in Fig. 3, the numerical results show that there is
already some probability to find the extinguished state around xln = 0.
It follows from the above discussion that extinguished particles originate from burning
particles, which, because of the fluctuations of the scalar dissipation rate, have experi-
enced a scalar dissipation rate high enough to completely extinguish the particle without
crossing the unstable branch of the S-shaped curve. This would result in re-ignition. The
low probability of finding these high scalar dissipation rates then forces the extinguished
particles to a state around xln = 0. In a real turbulent diffusion flame, these extinguished
areas could re-ignite by heat conduction from the surrounding, still burning gas. This
effect, however, is not included in the current analysis. Therefore, re-ignition can only
occur here if the scalar dissipation rate of an extinguished particle becomes smaller than
the ignition limit. This, however, is prohibited in the present simulations by choosing the
lower boundary for xln larger than the ignition scalar dissipation rate. This allows study
of the extinction process without the influence of auto-ignition.
It is important to recognize that because of this assumption the steady state solution
is always completely non-burning. This means that, for this dynamic system scalar dis-
sipation rate, fluctuations even of small amplitude lead to a phase transition from the
burning to the non-burning state. This dynamical character would not be observed in
the deterministic case.
For σ = 1 the distribution is even broader, revealing substantial probability for −2 <
xln < 2. Also, the probability of finding extinguished states is already of comparable
magnitude as for the burning states. As for σ = 0.5 the region of high probability is still
concentrated around the S-shaped curve, indicating that the chemistry is fast enough to
compensate scalar dissipation rate fluctuations. It is also interesting to note that similar
to the findings of Oberlack et al. (1999) there is only a very low probability of finding
states between burning and extinguished. This shows that the extinction process is fast
compared to other time scales of the system.
The solution for an even larger scalar dissipation rate variance of σ = 2 shows that
the probability distribution is further broadened and the fraction of extinguished states
is even higher. Most interesting here is the observation that, particularly at high scalar
dissipation rates close to extinction, the high probability region clearly departs from
the S-shaped curve. This can also be observed in Fig. 3 but to a smaller extent. The
departure from the S-shaped curve indicates that the chemistry is not fast enough to
relax the temperature in accordance with large scale scalar dissipation rate fluctuations
to the steady solution. At low scalar dissipation rate, the high probability region is still
very close to the S-shaped curve.

5. Conclusions and future work

In the present work the flamelet equations have been formulated for a one-step global
reaction and used for the investigation of the influence of scalar dissipation rate fluc-
tuations on non-premixed turbulent combustion. By modeling the diffusion term in the
flamelet equation, ordinary stochastic differential equations were derived for the temper-
ature and the scalar dissipation rate at stoichiometric mixture. From these, a Fokker-
Planck equation for the joint probability density function of temperature and the scalar
102 Pitsch & Fedotov
dissipation rate has been derived. The equation has been discussed and numerical solu-
tions for varying scalar dissipation rate variance provided.
The analysis shows that the S-shaped curve, which represents the steady-state solution
for a given scalar dissipation rate in the absence of scalar dissipation rate fluctuations,
separates the θst -xln space into two regimes, which will either lead to the burning or
the extinguished state. It is also shown that scalar dissipation rate fluctuations even of
small amplitude will under the present simplifications cause a phase transition from the
burning to the completely extinguished state.
Numerical solutions show an increasing fraction of extinguished states for increasing
scalar dissipation rate variance at a given time. It is also found that particles with a scalar
dissipation rate higher than the extinction limit can recover to a burning solution during
the extinction process. Therefore, for a fluctuating scalar dissipation rate, particles can
cross the S-shaped curve, which thereby no longer separates regimes that uniquely lead
to the extinguished or the burning state.
Furthermore, it is found that the low probability of finding high scalar dissipation rate
forces particles which have been extinguished at high scalar dissipation rate to rapidly
change to a state with lower scalar dissipation rate, where re-ignition could occur. For
higher scalar dissipation rate variance it is observed that the high probability region
clearly departs from the S-shaped curve. This indicates that the chemistry is not fast
enough to relax large scale scalar dissipation rate fluctuations to the steady state solution.
This has been shown to have an important implication in the application of flamelet type
models in non-premixed turbulent combustion.
The presented method has been shown to provide a useful tool to study the effect of
random scalar dissipation rate fluctuations. In future work, the model is to be corrobo-
rated with results from direct numerical simulations of turbulent reacting flows and the
re-ignition process is to be included. The investigation of the influence of scalar dissipa-
tion rate fluctuations on auto-ignition delay times and pollutant formation could also be
a worthwhile extension of the present work.

Acknowledgments
The authors gratefully acknowledge funding by the US Department of Energy in the
frame of the ASCI program and the Center for Turbulence Research. We are indebted
to Chong Cha for many inspiring discussions and for providing solutions of Monte Carlo
simulations for the investigated problem.

REFERENCES
Barlow, R. S. & Chen, J.-Y. 1992 On transient flamelets and their relationship to
turbulent methane-air jet flames. Proc. Combust. Inst. 24, 231-237.
Barths, H., Peters, N., Brehm, N., Mack, A., Pfitzner, M. & Smiljanovski, V.
1998 Simulation of pollutant formation in a gas turbine combustor using unsteady
flamelets. Proc. Combust. Inst. 27, 1841-1847.
Buyevich, Y. A., Fedotov, S. P. & Tret’yakov, M. V. 1993 Heterogeneous reac-
tion affected by external noise. Physica A. 198, 354-367.
Cha, C. M. 2000 Private communication.
Fedotov, S. P. 1992 Statistical model of the thermal ignition of a distributed system.
Comb. Flame. 91, 65-70.
 Stochastic modeling of scalar dissipation rate fluctuations 103
Fedotov, S. P. 1993 Stochastic analysis of the thermal ignition of a distributed ex-
plosive system. Phys. Lett. A. 176, 220-224.
Haworth, D., C., Drake, M., C., Pope, S., B. & Blint, R., J. 1988 The importance
of time-dependent flame structures in streched laminar flamelet models for turbulent
jet diffusion flames. Proc. Combust. Inst. 22, 589-597.
Horsthemke, W. & Lefever, R. 1984 Noise-Induced Transitions, Springer-Verlag.
Klimenko, A. Y. & Bilger, R. W. 1999 Conditional moment closure for turbulent
combustion. Prog. Energy Combust. Sci. 25, 595-687.
Koren, B. 1996 A robust upwind discretization method for advection, diffusion and
source terms. in C. B. Vreugdenhil & B. Koren (eds). Numerical methods for
advection-diffusion problems, Vol. 45 of Notes on Numerical Fluid Dynamics, Vieweg
Verlag, Braunschweig.
Kuznetsov, V. R. & Sabel’nikov, V. A. 1990 Turbulence and Combustion, Hemi-
sphere Publishing Corp.
Mauss, F., Keller, D. & Peters, N. 1990 A lagrangian simulation of flamelet
extinction and re-ignition in turbulent jet diffusion flames. Proc. Combust. Inst. 23,
693-698.
Oberlack, M., Arlitt, R. & Peter, N. 1999 On stochastic Damköhler number
variations in a homogeneous flow reactor. Comb. Theory Modeling, submitted.
O’Brien, E. E. 1980 The probability density function (pdf) approach to reacting turbu-
lent flows. in L. P. A. & F. A. Williams (eds). Turbulent Reacting Flows, Springer
Verlag, pp. 185-218.
Peters, N. 1983 Local quenching due to flame stretch and non-premixed turbulent
combustion. Combust. Sci. Technol. 30, 1.
Peters, N. 1984 Laminar diffusion flamelet models in non-premixed turbulent com-
bustion. Prog. Energy Combust. Sci. 10, 319-339.
Peters, N. 1987 Laminar flamelet concepts in turbulent combustion. Proc. Combust.
Inst. 21, 1231-1250.
Pitsch, H., Barths, H. & Peters, N. 1996 Three-dimensional modeling of nox and
soot formation in di-diesel engines using detailed chemistry based on the interactive
flamelet approach. SAE Paper 962057.
Pitsch, H., Chen, M. & Peters, N. 1998 Unsteady flamelet modeling of turbulent
hydrogen/air diffusion flames. Proc. Combust. Inst. 27, 1057-1064.
Pitsch, H. & Steiner, H. 2000 Large-eddy simulation of a turbulent piloted
methane/air diffusion flame (Sandia flame D). Phys. Fluids, submitted.
Pitsch, H., Wan, Y. P. & Peters, N. 1995 Numerical investigation of soot formation
and oxidation under diesel engine conditions. SAE Paper 952357.
Pope, S. B. 1985 Pdf methods for turbulent reactive flows. Prog. Energy Combust. Sci.
11, 119.
Pope, S. B. 2000 Turbulent Flows, Cambridge University Press, chapter 12.
Saitoh, T. & Otsuka, Y. 1976 Unsteady behavior of diffusion flames and premixed
flames for count flow geometry. Combust. Sci. Tech. 12, 135-146.
Seshadri, K. 1999 Rate-ratio asymptotics applied to flames. Paper presented at the
Western States Section of the Combustion Institute, Fall 1999 Meeting, University
of California at Irvine, Irvine, California.