Chapter 7
Turbulent Combustion Modeling
7.1 Regimes of turbulent combustion
Numerous physical phenomena are involved in turbulent (premixed and nonpremixed)
combustion and different have been developed to capture them in the modeling depending on
the combustion regimes.The premixed combustion regimes are identi?ed according to key pa
rameter diagram proposed in [32].In the case of high intensity turbulence,the dispersion of
Figure 7.1:Diagramof turbulent combustion regimes.
?uid elements becomes suf?ciently large to cause different portions of wrinkled?ame sheets to
approach each other.Laminar?ame propagation may cause the reaction fronts to intersect and
cut off pockets of unburnt gas.The?ame front reconnection can be interpreted as a turbulent
?ame propagation,and formation of an archipelago of unburned gas pockets surrounded by
combustion products could be a main consequences of chemistryturbulence interaction.The
transition from one regime of turbulent combustion to another is responsible to?ame accel
eration phenomenon because of the increase in?ame surface area of the?amelets inside the
turbulent structure.For suf?ciently high level of turbulence,the?amelet structure maybe
destroyed and replaced by a?distributed?reaction zone.Thus,under these conditions,turbu
lence exerts the main in?uence on the motion of the reaction zone,but chemical kinetics and
molecular transport effects are never negligible.
38
7.2 Multicomponent gas phase chemistry
The averaged gas phase equations for a chemically reacting mixture of ideal gases with
embedded condensed liquid droplets can be written as
@Y
k
@t
+r[Y
k
~u
t
Sc
k
rY
k
] = _
c
k
+ _
s
k
;k = 1;:::;N
s
(7.1)
where _
c
k
is the chemical source termde?ned by combustion mechanism.
The N
r
elementary chemical reactions involving N
s
species X
s
are represented in the form
N
s
X
s=1
0
sr
X
s
k
r
f
*
)
k
r
b
N
s
X
s=1
00
sr
X
s
;r = 1;N
r
:(7.2)
Here X
s
stands for one mole of sspecies,and
0
sr
and
00
sr
are stoichiometric coef?cients for the
rth reaction.
Formal kinetics requires
0
sr
00
sr
= 0,i.e.species?s?is formed only fromthe other species.
The chemical mass production terms in the mass balance eqs (4) are
_
s
= M
s
N
r
X
r=1
(
00
sr
0
sr
) _!
r
;(7.3)
and the chemical heat release termin the internal energy equation is given by
_
Q
c
=
N
r
X
r=1
q
r
_!
r
;(7.4)
where q
r
is the heat of rreaction at T = T
ref
,
q
r
=
X
s
(
0
sr
00
sr
)(h
s
)
o
f
:(7.5)
The variable _!
r
is the rate of progress of the rth reaction:
_!
r
= k
r
f
(T)
N
s
Y
k=1
(
k
M
k
)
0
kr
k
r
b
(T)
N
s
Y
k=1
(
k
M
k
)
00
kr
;(7.6)
where T is the temperature,k
r
f
and k
r
b
are the rate coef?cients for forward and backward stages
of the rreaction.The rate multipliers k
r
f
and k
r
b
are used,as a rule,in a generalized Arrhenius
form
k
r
f
= a
r
f
T
r
f
exp(T
r
af
=T)M
r
;
k
r
b
= a
r
b
T
r
b
exp(T
r
ab
=T)M
r
;(7.7)
with the equilibriumconstraint
k
r
f
= k
r
b
K
r
e
(T);(7.8)
where K
r
e
is the equilibrium constant in concentration units,and the third body factor,M
r
,
de?ned as
M
r
=
(
P
N
s
m
m
=W
m
r;m
for third body reactions
1:0 otherwise
(7.9)
where
r;m
are the third body ef?ciencies for mspecies in rreaction.
39
7.2.1 Mechanismof hydrogen combustion in oxygen (air)
Hydrogen/oxygen combustion,which powers e.g,the space shuttle,is nearly the most pow
erful chemical reaction known and the combustion mechanism has been studied for decades.
The?niterate H
2
=O
2
combustion mechanism,(8 species,26 reactions),and reaction rate pa
rameters in the Arrhenius form are listed in Table 1 and based mainly on the data in [?].
The chaperon ef?ciencies for the third body concentrations M
r
are de?ned by following expres
sions:
1;H
2
=
2;H
2
=2.5,
1;H
2
O
=12.0,
2;H
2
O
=15.0.Rate coef?cients for the backward reactions
are computed with the help of linear regression by satisfying of the constraint (7.8).Stages
2731 included in the Table 1 represent chemically equilibriummechanism used in the rocket
propulsion modeling.
The mechanism was tested by comparing computed results with the experimental data [?]
for ignition delay times of stoichiometric H
2
=Air mixtures at 2 atm,(see Fig.7.2a).Although
experimental data for rocket chamber conditions are not available in literature,relevant data
for highly diluted (by Ar) stoichiometric mixtures obtained using a new highpressure shock
Table 1
Reaction Mechanismfor H
2
=O
2
Combustion
Reaction a
r
f
r
f
E
r
af
1:H
2
+O
2
*
)
OH +OH 1:7 10
13
0:0 47780:
2:H
2
+OH
*
)H
2
O +H 1:2 10
09
1:30 3626:
3:O+OH
*
)O
2
+H 4:0 10
14
0:50 0:
4:O+H
2
*
)OH +H 5:0 10
04
2:67 6290
5:H +HO
2
*
)O +H
2
O 3:1 10
10
0:0 3590:
6:O+OH +M
*
)HO
2
+M 1:0 10
16
0:0 0:
7:H +O
2
+M
2
*
)HO
2
+M
2
3:6 10
17
0:72 0:
8:OH +HO
2
*
)H
2
O +O
2
7:5 10
12
0:0 0:
9:H +HO
2
*
)OH +OH 1:4 10
14
0:0 1073:
10:O+HO
2
*
)O
2
+OH 1:4 10
13
0:0 1073:
11:OH +OH
*
)O +H
2
O 6:0 10
08
1:3 0:
12:H +H +M
1
*
)H
2
+M
1
1:0 10
18
1:0 0:
13:H +H +H
2
*
)H
2
+H
2
9:2 10
16
0:6 0:
14:H +H +H
2
O
*
)H
2
+H
2
O 6:0 10
19
1:25 0:
15:H +OH +M
*
)H
2
O +M 1:6 10
22
2:0 0:
16:H +O +M
*
)
OH +M 6:2 10
16
0:6 0:
17:O+O +M
*
)
O
2
+M 1:9 10
13
0:0 1788:
18:H +HO
2
*
)
H
2
+O
2
1:2 10
13
0:0 0:
19:HO
2
+HO
2
*
)
H
2
O
2
+O
2
2:0 10
12
0:0 0:
20:H
2
O
2
+M
*
)
OH +OH +M 1:3 10
17
0:0 45500:
21:H
2
O
2
+H
*
)
HO
2
+H
2
1:6 10
12
0:0 3800:
22:H
2
O
2
+H
*
)
H
2
O +OH 1:0 10
13
0:0 3590:
23:H
2
O
2
+OH
*
)
H
2
O +HO
2
1:0 10
13
0:0 1800:
23:H
2
O
2
+H
*
)H
2
O +OH 1:0 10
13
0:0 3590:
24:H
2
O
2
+O
*
)H
2
O +O
2
8:4 10
11
0:0 4260:
25:H
2
O
2
+O
*
)OH +HO
2
2:0 10
13
0:0 5900:
26:HO
2
+H
2
*
)H
2
O +OH 6:5 10
11
0:0 18800:
27:O
2
*
)
O +O 1:0
28:H
2
*
)
H +H 1:0
29:H
2
+O
2
*
)
OH +OH 1:0
30:H
2
O +OH
*
)
HO
2
+H
2
1:0
31:H
2
O +H
2
O
*
)
H
2
O
2
+H
2
1:0
The table consists of two blocks for elementary and equilibrium reactions,respectively.Rate
coef?cients,third body concentrations M
i
and chaperon ef?ciencies are taken in the formof [?].
40
0.70 0.80 0.90 1.00 1.10 1.20
1000/T, K
10
1
10
2
10
3
10
4
10
5
Ignition delay time, sec
Ignition Delays for H
2
/Air Mixtures
Experiment, 1 atm
Experiment, 2 atm
This paper, 1 atm
This paper, 2 atm
Balakrish.&Williams
Balakrish.&Williams
Miller et al.
Miller et al.
=1
(a)
1 10
Pressure, bars
10
1
10
2
10
3
10
4
10
5
Ignition delay, sec
Ignition of H
2
−O
2
Mixtures
T
o
= 950 K
T
0
=1000 K
T
o
=1100 K
=1
(b)
Figure 7.2:a) Calculated ignition delays for the stoichiometric H
2
/O
2
mixture vs shocktube
experimental data;b) nonmonotonous dependence of ignition delays on pressure for different
initial temperatures.
5 6 7 8 9 10 11
10000/T, K
10
0
10
1
10
2
10
3
10
4
10
5
Ignition delays, sec
Auto−ignition H
2
−O
2
−Ar Mixtures
(Comparison with shock−tube data)
33 atm, calc.
33 atm, exp.
5 atm, calc.
5 atm, exp.
Figure 7.3:Predicted and measured (in shock tube experiments) ignition delay times for stoi
chiometric H
2
/O
2
/Ar mixtures at elevated pressures.
tube [49] (stoichiometric mixtures for pressures between 33 and 87 atm
1
and temperatures
from1175 to 1880 K) were compared with predictions,and the results of such a comparison are
presented in Fig.7.3.The measured ignition delays are insensitive to inertgas dilution for Ar
levels (95.099.85 %) used in the experiments.At high temperatures,the ignition delays are
practically independent on pressure,but they turn out to be very sensitive to pressure levels at
high temperatures.Ignition delays are less sensitive to pressure level at low temperatures.
A sensitivity analysis (based on the SENKIN code of the Chemkin2 library) applied to this
reaction mechanismstudy shows that at lowpressures the ignition process is well controlled by
the classic H
2
=O
2
combustion mechanism,(reaction 3 in Table 1 dominates),while at high pres
sures,the role of the trimolecular reaction 7 becomes more pronounced during the induction
period.The competition between reactions 3 and 7 explains the nonmonotonous dependence
1
The 5 atmdata of Skinner and Ringrose are added for completeness
41
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
1000
0
1000
2000
3000
4000
reaction number, j
sensitivity coefficient, T/Aj
T
o
=1100 K, P
o
=1 bar, =1.0
(a)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
100
0
100
200
300
400
500
sensitivity coefficients, T/Aj
reaction numbers, j
T
o
=1100 K, P
o
=100 bar, =1.0
(b)
Figure 7.4:Absolute sensitivity coef?cients for the H
2
=O
2
reaction mechanismof Table 1 at the
end of ignition period (t
id
=3.010
5
sec):a) atmospheric pressure case,b) high pressure case.
of ignition delays on the pressure presented in Fig.7.2b and the other peculiarities of ignition
at high pressures.
7.2.2 Motor fuel surrogates models
The word fuels,see [17],refers to a range of products re?ned from crude oil,which are
usually liquids (occasionally gasses),that burn in the presence of air and enable the operation
of heat engines,either rocket,piston (mainly gasoline and diesel) or turbine (aircraft engines).
Practical diesel fuels consist of a great number of aliphatic and aromatic compounds,and their
combustion is too complex to be modeled using a comprehensive chemical mechanism.Instead
of this,the model or surrogate fuels are used in numerical simulations.Aliphatic components
can be represented by long chain hydrocarbons such as nheptane or ndodecane to their cetane
number of approximately 56,which is similar to the cetane number of conventional diesel
fuel.Aromatic components signi?cantly contribute to soot formation.The diesel fuel surrogate,
0.7 0.9 1.1 1.3 1.5
1000 K/T
10
−2
10
−1
10
0
10
1
10
2
Ignition delays, ms
Autoignition of Diesel Fuel Surrogate
(Simulation of shock−tube experiments)
diesel oil surrogate
n−heptane, C
7
H
16
toluene, C
7
H
8
=1.0
Figure 7.5:Ignition delays for fuel surrogate and its components simulating shock tube exper
iments;p
o
=41.0 bar, = 1 for all compounds.
42
which can suf?ciently represent the cetane number as well as the other properties of real fuel
is assumed to be a 70/30 % mixture of nheptane,C
7
H
15
,and toluene,C
7
H
8
.The alternative
model consisting of the mixture of ndodecane and methylnaphthalene was also proposed.If
the physical properties of the model fuel are represented by properties of real diesel oil,the
difference in combustion models seems not critical.
At?rst,detailed reaction mechanisms for nheptane and toluene oxidation have been devel
oped and validated using shocktube autoignition experimental data.Such a detailed mecha
nismintegrating the nheptane and toluene oxidation chemistry with the kinetics of aromatics
(up to four aromatic rings) formation for rich acetylene?ames by [15] was,then,reduced on the
bases of sensitivity analysis (using the software of the Chemkin2 [4] package) to the mecha
nismconsisting of 68 species participating in 280 reactions which implemented into the KIVA
3V code [42][43] to simulate diesel spray combustion.This mechanismwas assumed to predict
the fuel autoignition,combustion development and formation of polycyclic aromatic hydrocar
bons (PAHs) till gaseous soot precursors based on acenaphthylene,A2R5,which are nonplanar
aromatic rings.The mechanismfor both aliphatic and aromatic components includes principal
fragments of C
1
C
7
and O/Hchemistry taken fromdifferent sources,but mainly fromLLL [44]
and MIT [45] databases.More detailed discussion of the fuel surrogate models can be found in
the Appendix III.
7.2.3 Soot formation model based on PAHconcept
Soot formation mechanism consists of series of elementary reaction steps leading from
acetylene and hydrogen to the formation of the?rst aromatic ring,A
1
.Reaction steps lead
ing to the formation of the phenyl,A
1
,radical and the?rst aromatic ring are followed by
the successive stages of H abstraction,C
2
H
2
 addition (HACA mechanism),thus,yielding a
chain of aromatic rings.Toluene can formA
1
without the presence of acetylene.A special reac
tion sequence was included leading via benzene dimerization to acenaphthylene.The incipient
soot was formed fromlongchain acetylene,C
6
H
2
,and acenaphthylene via?graphitization?pro
cesses
C
6
H
2
*H
2
+6C(s) (7.10)
A
2
R5 *4H
2
+12C(s)
The notation C(s) + H
2
can be formally attributed to incipient soot as the latter contains a
signi?cant amount of hydrogen.Thermal properties of C(s) have been taken as properties of
graphite.The complete reaction mechanismand species thermodynamic property are provided
by [46].The physical properties of the fuel surrogate were taken as properties of real fuel
compiled in the DI model of the KIVA3V fuel library with the original fuel chemical formula
replaced by C
14
H
28
.
Soot formation mechanism consists of series of elementary reaction steps leading from
acetylene and hydrogen to the formation of the?rst aromatic ring,A
1
.Reaction steps lead
ing to the formation of the phenyl,A
1
,radical and the?rst aromatic ring are followed by
the successive stages of H abstraction,C
2
H
2
 addition (HACA mechanism),thus,yielding a
chain of aromatic rings.Toluene can formA
1
without the presence of acetylene.A special reac
tion sequence was included leading via benzene dimerization to acenaphthylene.The incipient
soot was formed fromlongchain acetylene,C
6
H
2
,and acenaphthylene via?graphitization?pro
cesses
C
6
H
2
*H
2
+6C(s) (7.11)
A
2
R5 *4H
2
+12C(s)
43
The notation C(s) + H
2
can be formally attributed to incipient soot as the latter contains a
signi?cant amount of hydrogen.Thermal properties of C(s) have been taken as properties of
graphite.The complete reaction mechanismand species thermodynamic property are provided
by [46].The physical properties of the fuel surrogate were taken as properties of real fuel
compiled in the DI model of the KIVA3V fuel library with the original fuel chemical formula
replaced by C
14
H
28
.
7.3 Multicomponent molecular diffusion
Diffusion transport in the Ncomponent gas mixture is governed by the systemof equations:
N
X
j6=i
(x
i
x
j
)=D
ij
(u
j
u
i
) = rx
i
;i = 1;:::;N (7.12)
where u
i
is the ith species speci?c velocity,x
i
is the mole fraction of ith species,D
ij
is the
binary diffusivity of the pair (i,j).
If speci?c velocities are de?ned fromthe system(7.12),the diffusion?uxes J
i
relative to the
mass averaged velocity u needed for closure of the mass and energy conservation laws of?uid
mechanics can be calculated as
J
i
=
i
(u
i
u);i = 1;:::;N (7.13)
where
i
is the partial mass density of ith species.
In the system(7.12),only N1 equations are linearly independent,as their sumover i yields
0=0.Thus,the determinant of the system (7.12) is equal to 0,and a unique solution for this
system does not exist.The lack of the solution uniqueness is due to the fact that Eqs (7.12)
involve only velocity differences,and there is no information about absolute velocities.The
latter information is contained in the equation de?ning the massaveraged velocity.If this
velocity is de?ned as
u =
N
X
i
(
i
=)u
i
;(7.14)
the diffusion?uxes identically sumto zero.,i.e
N
X
i
J
i
= 0 (7.15)
7.4 Chemistry/turbulence interaction
The direct implementation of the approach described above for the particular chemical
mechanism in the CFD codes sometimes is restrictive due to the computer time and mem
ory constraints.Moreover,this approach known as the quasilaminar model does not included
the effect of turbulence on chemical reactions.Till recently,there was no explicit coupling of
chemical kinetics and the rate of combustion in entirely turbulent mixingcontrolled regime.
In particular,in [33] was proposed a method which requires only the evaluation of averaged
species mass fractions to describe chemical kinetics,albeit in a very simple form,in the course
of the precombustion period.Let us analyze this model (known as the Eddy Dissipation Con
cept,EDC,in more details.
44
The original EDC model of [33] assumes a onestep global reaction such as
(1kg)fuel +(kg)oxidizer!(1 +kg)products;
where is the stoichiometric ratio of mass of oxidizer to that of fuel.
Although an application of the original EDC model to multireaction mechanismis,in prin
ciple,possible,it meets some dif?culties.They are not only because of the fact that in the past,
this has been dictated only by the fact that till recently,the integrations of largescale chemical
mechanisms with 3DCFDmodels were far beyond the capabilities of practical simulations.The
Eddy Dissipation Concept (EDC) is the model which can facilitate the efforts making this in
tegration a reality.However,the EDC combustion model was not yet analyzed in detail,albeit
the model is said to be implemented in practically all commercial codes available.For example,
in the latest review by Veynante and Vervisch published recently in [31],the model has been
referred as belonging to the family of?in?nitely fast chemistry?models and summarized in only
15 halflines of the following text:?The eddy dissipation model (EDC) is a direct extension to
nonpremixed?ame of the eddy break up (EBU) closure initially devoted to turbulent premixed
combustion...The fuel burning rate is calculated according to:
_w
f;edc
=
k
min(
~
Y
f
;
~
Y
o
s
;
~
Y
p
1 +s
);(7.16)
where and are adjustable parameters of the closure?.In the above equation,?...the reaction
rate is limited by a de?cient species.To account for the existence of burnt gases bringing the
energy to ignite the fresh reactants,this species may be the reaction product (when is non
zero)...Actually,for large ,the model is dif?cult to justify.?
It is clear that the above review attributes the EDC model only to its?rst,Magnussen and
Hjertager,MH,formulation [33],but even this version of the model was not amply covered in
the literature.The reaction product term (when is nonzero) is introduced into the rate ex
pression to account for?niterate chemistry effects (ignition/quenching) that makes it incorrect
to place the EDCapproach in the family of?in?nitely fast chemistry?models.In fact,the actual
reaction rate is taken to be in a very controversial form
_!
j
= min( _!
c
;_!
edc
) (7.17)
argued by?heuristic?considerations to extend the approach to multistep reaction mechanisms.
When the reaction rates are determined,the source terms in the species conservation equa
tions are de?ned as
_
c
k
= M
k
Nr
X
j=1
(
00
k;j
0
k;j
) _!
j;k
;(7.18)
and the chemical heat release termin the energy equations becomes
_
Q
c
=
Nr
X
j=1
Q
j
_!
j
;(7.19)
where Q
j
is the negative of the jreaction heat at the reference temperature de?ned as
Q
j
=
X
k
(
00
k;j
0
k;j
)(h
o
f
)
k
;(7.20)
where (h
o
f
)
k
is the heat of formation of kspecies at the reference temperature.
45
Some other versions of so called EDC model,e.g.,[35] were formulated without a proper
connection with the original EDC approach.Also,Fluent [36] claimed the implementation of
a new version of the EDC model?...for the modeling of turbulent?niterate chemistry.It is
completely different fromthe eddy dissipation model that has long been available.?
The section below is aimed to?ll a gap in the EDC based models analysis starting with an
outline of the classic EDC model features.The basic assumption of the EDC is that chemical
reactions occur only in the turbulence?ne structure,i.e.,Kolmogorov size eddies and that the
reactions are quenched,if the characteristic chemical times for limiting species are longer than
the Kolmogorov time.Thus,the classic EDC model does not allowthe reaction zone broadening
beyond the Kolmogorov limit.
7.4.1 Formulation of the newEDC Model
The new model formulation is based on the operatorsplitting procedure (see Appendix II)
applied to the mass conservation equations for species participating in any multistep reaction
mechanism.In terms of this approach,the time differencing was performed in three steps:the
?rst step was assumed to be the convection contribution,the second one  diffusion effect with
out a contribution of micromixing,and the third step was the chemical kinetics effect coupled
with micromixing.The last step of such a mass balance can be interpreted as representing
combustion in a constant volume partially stirred reactor (PaSR) of a computational cell size,
where reactions occur in a fraction of its volume.Since the reaction zone parameters,as a rule,
can not be resolved on a computational grid,the subgrid diffusion term due to micromixing
was approximated with the help of introduction of a micromixing time as suggested in the
?interaction with the mean?approach [37].Then,assuming that chemical processes proceed in
such a way that the shortest chemical time associated with particular (reference) species partic
ipating in the reaction was constrained by the micromixing time,a system of PaSR equations
consisting of two sets (?rst  differential,second  algebraic) can obtained
dc
1
dt
=
c
1
c
o
=
c
c
;
c c
1
mix
=
c
c
;(7.21)
where is the time step,
c
is the chemical reaction time,and
mix
is the micromixing time.
The model distinguishes (see Fig.2) between the concentration (in mean molar density) at
the reactor exit,c
1
,the concentrations in the reaction zone,c,and in the feed,c
0
.When time
proceeds,c
1
trades place for c
0
.
After algebraic manipulations,one can yield the analytical solutions of the above problem
CC
0 1
C
*
Figure 7.6:Partially stirred reactor,PaSR,consisting of nonpremixed,premixed,and pre
mixed and reactive regions.
46
which,?nally,can be represented as
c
1
c
o
= (
c
1
c
) =
1
2
H(
c
1
c
;
c
1
mix
);(7.22)
where =
c
=(
mix
+
c
),and His a harmonic mean.
The solution (7.22) shows that the turbulent combustion time is a sum of the mixing and
reaction times,if the process is expressed in terms of the reactor output parameters.This is the
main consequence of the PaSR model.
For the detailed chemical mechanism,the species production rate due to the particular r
chemical reaction (separated into creation and destruction rates) contrary to the?rst equation
of the system(7.21) can be rewritten as
dc
1
d
= f
r
(c) =
c c
1
mix
(7.23)
where
f
r
(c) = (
00
r
0
r
) _!
r
(c) = term
r
term
r
where
00
r
and
0
r
are stoichiometric coef?cients of the backward and forward stages,_!
r
is a
progress variable rate of the rreaction.Species indices are omitted for simpli?cation.
The reaction progress variable rate _!
r
(c) is calculated fromthe massaction law
_!
r
(c) = k
r
f
(T)
N
s
Y
s=1
(c
s
)
0
sr
k
r
b
(T)
N
s
Y
s=1
(c
s
)
00
sr
;
where k
r
f
and k
r
b
are the rate coef?cients for the forward and backward stages of the rreaction,
T is the temperature,and N
s
is the total number of species in the mixture.here is a substantial
difference between Eqs (7.21) and (7.23)  the original model has the analytical solution,while
the application of the second one requires numerical integration.Combination of the terms in
the model (7.23) leads to a de?nition of the reactive volume as introduced in [48]:
c
1
= c
+c
o
(1
);(7.24)
where the multiplier
= =( +
mix
),while is the integration step.
Rewriting the?rst relation of (7.23) in terms of the reactor exit parameters,one can get
c
1
c
o
= f
r
(c)'f
r
(c
1
) +(@f
r
=@c)
jc=c
1
(c c
1
)
= f
r
(c
1
)
c c
1
c
The rhs of the above equation is constructed using Taylors expansion at the state c
1
,assuming
that the reaction times can be estimated as reciprocal values of elements of the Jacobian matrix
evaluated at unknown values c=c
1
,i.e.,
c
[@f
r
=@c]
1
and that (@f
r
=@c)
jc=c
1 < 0.Elimination
of cvalues using c
1
values taken from Eq.(7.24)
c
1
c
o
= f
r
(c
1
) [
c
1
c
o
(1
)
c
1
]=
c
;(7.25)
if the reaction time is de?ned as
c
1
= f
r
o
=(c
s
o
+term
r
);
47
immediately puts the reaction rate into the form
c
1
c
o
= f
r
(c
1
)
c
c
+
mix
(7.26)
=
c
s
o
f
r
o
c
s
o
+term
r
+term
r
mix
;
containing no reaction zone parameters c which cannot be resolved on a grid,but replacing
their effect with the rate multiplier .The superscript
o
and the subscript
s
denotes the value
at the start of the integration step and the?reference?species index de?ned below.
This expression looks different from the original?niterate formulation of the EDC model
known as
_w
s
=
M
1
M
c
s
c
1
s
mix
;(7.27)
where
M
is the?ne structure volume,and
denotes the?nestructure quantities,but a generic
proximity of the approximations can be established using the PaSR formulation.
The chemical reaction times are formally de?ned as characteristic times of the destruction
rates,if the species chemical production rates are presented as a sumof creation and destruc
tion terms.In this way,for each reaction of the mechanism,chemical times can be attributed
to each species participating in the depletion stages,e.g.,k
r
c
i
c
k
c
i
=
c;i
c
k
=
c;k
,where
i;k are species indices.In terms of the EDC model,the shortest of these times is restricted
from below by the micromixing time that explains the usage of de?cient (limiting) species in
the rate expression (7.16).
Application of the rate expression (7.26) to a linear approximation of the Arrhenius reaction:
f
o
r
= c
o
s
=
c
leads to the formula
c
o
s
f
o
r
c
o
s
+term
r
+term
r
mix
=
c
o
s
+
c
+
mix
=
c
1
s
c
+
mix
=
1
2
H(
c
1
s
c
;
c
1
s
mix
)
which is perfectly consistent with the PaSR analytical solution (7.22).In a limiting case of
chemical equilibrium,the model automatically approaches the MHmodel.
7.4.2 Reference species denition
The reference species is similar to the limiting species used in the MHmodel [33] by de?n
ing the reaction rates.Let us consider the balance equation of sspecies participating in the
rreaction:
c
1
s
c
o
s
= (
00
sr
0
sr
)[k
r
f
(T)
1
k
r
b
(T)
2
] =
[
0
sr
k
r
b
2
+
00
sr
k
r
f
1
]

{z
}
term
r
[
0
sr
k
r
f
1
+
00
sr
k
r
b
2
]

{z
}
term
r
;
where
1
=
N
s
Y
n=1
(c
1
n
)
0
nr
;
2
=
N
s
Y
n=1
(c
1
n
)
00
nr
;
48
The above equation written in a formseparating the net chemical production rate into creation
and destruction terms is assumed to be linearized as
c
1
s
c
o
s
= f
r
(t;c
1
s
) = term
r
term
r
c
1
s
c
o
s
(7.28)
c
1
sj=0
= c
o
s
;
This linearized expression is mostly accurate for the species which concentration is less than
those of the other reaction partners.Thus,such a species (called reference species) is playing
a special role in the analysis.The factor before this species has a dimensional representation
of an inverse time.This time is the shortest time for the particular reaction.By algebraic
manipulation with Eq.(7.28),one can get the expressions for unknown terms in Eq.(7.26),
viz.,f
r
(c
1
) and
c
:
f
r
(c
1
s
) =
c
o
s
f
r
o
c
o
s
+term
r
;
c
1
= (@f
r
=@c)
j
c=c
1
f
r
o
=(c
s
o
+term
r
)
As shown in [48],introduction of reference species assures the equilibrium conditions imple
mentation at !1for each reaction in the mechanism.
7.4.3 Micromixing time denition
The correct de?nition of the micromixing time is a matter of importance for any EDC tur
bulent combustion model.The micromixing time de?nition taken in this study is applicable to
the systemwith?natural?initial mixture nonuniformities and based on the model
2
proposed
by Frisch [40,p.135140].In terms of this approach,at each stage of the Richardson cascade,
the number of eddies formed froma given?parent?eddy is chosen such that the fraction of vol
ume occupied by?active?eddies is decreased by a factor < 1.In this way,the intermittency
is introduced in a geometric form.The factor is an adjustable parameter of the model.If l
o
is
the initial eddy size,and l is the size of the eddy formed on the nth step of fragmentation,the
fraction p
l
of the active space can be calculated as
p
l
=
n
= (
l
l
o
)
3D
;(7.29)
where the exponent D can be interpreted as a fractal dimension of the turbulent?eld.The
expression (7.29) can be derived by considering the analytical representation of the Richardson
cascade in the form:
f(x)f(y=x) = f(y)f(1);(7.30)
where f(1) is a number of eddies at the initial stage,f(x) = f(l
o
=) is a number of fragments
at the particular stage of fragmentation,x = l
o
=,y = l
o
=, and , < are fragment sizes.
he above equation has a general solution the?fragmentation law?
N f(x) = Cx
D
;(7.31)
where C and D are the solution parameters,N is the number of fragments.The fraction of the
?active?volume can be calculated as
V ol
N
= f(x)
3
= C(l
o
=)
D
3
= Cl
D
o
3D
= Cl
3
o
(=l
o
)
3D
2
This has nothing in common with that mentioned in the review paper [31].
49
This relation is identical to the expression (7.29).Since the active eddies of size l?ll only a
fraction p
l
of a total volume,the speci?c energy associated with this scale is
E
l
= v
l
2
p
l
= v
l
2
(
l
l
o
)
3D
(7.32)
assuming the energy?ux fromscales l to smaller scales is E
l
=
l
,where
l
= l=v
l
is the eddy
turnover time.
In the inertial range of scales,the energy?ux is independent of l,i.e.,
= v
3
o
=l
o
;
where is a dissipation rate of turbulent kinetic energy.After some algebra,one can get the
expression for the eddy turnover time
l
=
l
o
v
o
(
l
l
o
)
2
3
+
3D
3
(7.33)
The viscous cutoff scale of the model is obtained by equating the eddy turnover time and
the viscous diffusion time
d
=
2
d
= that gives an expression for the dissipation scale
d
= l
o
Re
3
1+D
;(7.34)
where Re = l
o
v
o
=,and is molecular viscosity.When D = 3,the classic Kolmogorov expression
is recovered;indeed,this assumption neglects the intermittency.
f l
o
is de?ned in terms of k model of turbulence,from Eq.(7.34),the viscous dissipation
time can be derived as
d
=
mix
= (
k
)
[(
)
1=2
]
(1)
(7.35)
= (
k
) (c
=Re
t
)
1
2
;
where =
3(D3)
1+D
,and k and are the turbulent kinetic energy and energy dissipation rate,
respectively.Fromexperiments reported in [41] follows that the fractal dimension for turbulent
dissipation is about D=2.7,i.e
mix
= (
k
) (c
=Re
t
)
0:
621
In particular,the eddy breakup time k= corresponds to D=5,and Kolmogorov microscale
time
k
= (=)
1=2
 to D=3.
Another approach developed is based on the usage of Kolmogorov expression for the micro
mixing time with the molecular viscosity replaced by a fraction of the effective viscosity as
sumed to be responsible for micromixing.If,for example,the RNG k model is employed,
the turbulent viscosity related to k,the turbulent kinetic energy,and ,the dissipation rate of
k,is given by the general expression.
t
=
l
(1 +
s
c
k
2
=
l
)
2
(7.36)
= (
l
+
s
t
+2
p
l
s
t
);
where
s
t
is the?standard?k value of the kinematic viscosity.
Provided that the?rst two terms in the expression (7.36) contribute to conventional diffu
sion transport,the geometrical mean
3
term can be assumed to determine the characteristic
3
The micromixing time of such a formhas been proposed (without deviation) in [26].
50
time scale for micromixing in the turbulence/chemistry interaction model,if written in a form
similar to the Kolmogorov time de?nition,
k
= (
l
=)
1=2
replacing the molecular viscosity by a
corresponding expression for the turbulent viscosity,i.e.,
mix
= (2
p
l
s
t
=)
1=2
= (2
p
c
k
k
)
1=2
;(7.37)
where c
=0.09 is the constant of the k model of turbulence.
The above expression formally corresponds to D=3.8 in the general formula (7.35) giving
only a provisional value for
mix
that must be re?ned in comparisons with experiments.The
model is assumed to be valid across a full range of?ow conditions from low to high Reynolds
numbers,if k and are determined from the generalized transport equations of the RNG k
model of turbulence implemented in the code [43].
51
Bibliography
[1] Sirignano,W.A.,Fluid Dynamics and Transport of Droplet and Sprays,Cambridge Uni
versity Press (2000)
[2] Kuo,K.K.,ed.,Recent Advances in Spray Combustion:Spray Combustion Measurements
and Model Simulations,Vol.II,in Progress in Astronautics and Aeronautics,Vol.171
(1996)
[3] Williams,F.A.,Combustion Theory,2nd ed.,AddisonWisley Publishing Co.,Inc,London
(1985)
[4] Kawano,D.,Senda,J.,Wada,Y.,et al.,Numerical Simulation of Multicomponent Fuel
Spray,SAE Paper 2003011838 (2003)
[5] Torres,D.J.,Rourke,P.J.,and Amsden,A.A.,A Discrete Multicomponent Fuel Model for
GDI Engine Simulations,ILASS Americas,14th Annual Conference on Liquid Atomiza
tion and Spray Systems,Deaborn,MI (May 2001)
[6] Veynante,D.,Vervisch,L.,Turbulence Combustion Modeling,Progress in Energy and
Combustion Science,28:193266 (2002)
[7] Borghi,R.,Turbulent Combustion Modeling,Progr.Energy Combust.Sciences.,14:245292
(1988)
[8] Magnussen,B.F.,Hjertager,B.H.,On the Numerical Modeling of Turbulent Combustion
with Special Emphasis on Soot Formation and Combustion,Sixteenth Symposium(Inter
national) on Combustion,Pittsburgh:The Combustion Institute,p.719729 (1976)
[9] Kralj,C.,Numerical Simulation of Diesel Spray Processes,Ph.D.thesis,Imperial College,
London (October 1995)
[10] Kong,S.F.,Marriot,C.D.,Reitz,R.D.,and Christensen,M.,Modeling and Experiments of
HCCI Engine Combustion Using Detailed Chemical Kinetics with Multidimentional CFD,
SAE Paper 2001011026 (2001) 4
[11] Fluent6:Applications:Reacting?ows,http://www.?uent.com/software/?uent/application
/reacting.htm(2002)
[12] Aubry,C.,and Villermaux,J.,J.Chemical Engineering Science,30,p.457 (1975)
[13] Golovitchev,V.I.,Nordin,N.,Jarnicki,R.,and Chomiak,J.,3D Diesel Spray Simulations
Using a New Detailed Chemistry Turbulent Combustion Model,SAE Paper 2000011891
(2000)
[14] Liang,P.Y.,and Ungewitter,R.J.,MultiPhase Simulations of Coaxial Injector Combustion,
AIAA Paper 920345 (1992)
52
[15] Appel,J.,Bockhorn,H.,and Frenklach,M.,http://www.me.berkeley.edu/soot/mecha
nisms/abf.html (1999)
[16] Kee,R.J.,Rupley,F.M.,and Miller,J.A.,ChemkinII:A Fortran Chemical Kinetics Package
for the Analysis of Gas Phase Kinetics,Sandia Report SAND898009B,UC706 (April
1992)
[17] Frisch,U.,Turbulence,the Legacy of A.N.Kolmogorov,Cambridge University Press (1995)
[18] Sreenivasan,K.R.,and Meneveau,C.,The fractal facet of turbulence,J.Fluid Mech.:vol.
173,pp.357386 (1986)
[19] Amsden A.A.,KIVA3:A KIVA Program with BlockStructured Mesh for Complex Ge
ometries,LA12503MS (1993)
[20] Amsden A.A.,KIVA3v:A Blockstructured KIVA Program for Engines with Vertical or
Canted Valves,LA13313MS,July (1997)
[21] Curran,H.J.,Gaffuri,P.,Pitz,W.J.,and Westbrook,C.K.,A Comprehensive Modeling Study
of nHeptane Oxidation,Combustion and Flame,114:149 177 (1998)
[22] Howard,J.B.,http://web.mit.edu/anish/www/MITcomb.html(2001)
[23] Golovitchev,V.I.,http://www.tfd.chalmers.se/?valeri/MECH(1998)
[24] Wiartalla,A.,B¤aker,H.,Br ¤uggeman,D.and Koss,H.,In?uence of Spray Development and
Fuel Quality on Soot Formation and Oxidation,JOULE0008D(AM) Report (1992)
[25] Golovitchev,V.I.,Nordin,N.,Detailed Chemistry SubGrid Scale Model of Turbulent
Spray Combustion for the KIVA Code,Paper 99ICE237,ICEVol.333:p.1725 (1999)
[26] Karlsson,J.,Modeling AutoIgnition,Flame Propagation and Combustion in Non Station
ary Turbulent Sprays,Chalmers University of Technology,PhD Thesis,G¤oteborg (May
1995)
[27] Edwards,C.F.,Siebers,D.L.,and Hoskin,D.H.,A Study of the Autoignition Process of a
Diesel Spray via High Speed Visualization.SAE Paper 920108 (1992)
[28] Onishi,S.,Jo,S.H.,Shoda,K.,Jo,P.D.,and Kato,S.,Active ThermoAtmospheric Com
bustion (ATAC)a New Combustion Process for Internal Combustion Engines,SAE Paper
790501 (1979)
[29] Frenklach,M.,et al,http://www.me.berkeley.edu/gri
mech (2000)
[30] Marinov,N.M.,et al,http:wwwcms.llnl.gov/combustion/nbutane
mech.txt (2000)
[31] Golovitchev.V.I.,and Chomiak,J.,Evaluation of Ignition Improvers for Methane Au
toignition,Comb.Sci.and Tech.:135,31 (1996)
[32] Lifshitz,A.,Scheller,K.,and Burcat,A.,ShockTube Investigation of Ignition in Methane
OxygenArgon Mixtures.Combustion and Flame,16,311 (1971)
[33] Petersen,E.L.,Davidson,D.F.,R¤ohrig,M.,and Hanson,R.K.,ShockInduced Ignition
of HighPressure H
2
O
2
Ar and CH
4
O
2
Ar Mixtures,31st AIAA/ASMe/SAE/ASEE Joint
Propulsion Conference and Exhibition,July 1012,AIAA Paper 953113 (1995)
53
Comments 0
Log in to post a comment