Versteeg H., Malalasekra W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method
Подождите немного. Документ загружается.
12.17 CFD OF TURBULENT NON-PREMIXED COMBUSTION 377
The instantaneous velocity u is written as
u = U + u ″= +u″ (12.131)
In contrast to the Reynolds decomposition, where u′ represents a turbulent
velocity fluctuation, the quantity u″ also includes effects of density fluctu-
ations. If the flow is incompressible, the density is constant, so U = R and
u″≡u′.
In the convective term of the continuity equation we require the Favre
average of the product
ρ
u
i
. Multiplying equation (12.131) by
ρ
we obtain
ρ
u =
ρ
(U + u ″) =
ρ
U +
ρ
u″ (12.132)
Time averaging equation (12.132) we get
= 4U + (12.133)
By definition of Favre averaging (12.130) and equation (12.131) we have
= 0. Now the Favre-averaged continuity equation can be obtained:
+=0 (12.134)
Unlike Reynolds-averaged continuity equation (12.128), this equation has
the same form as the original continuity equation (12.127) and Reynolds-
averaged constant density equation (12.129), except that the mean velocity is
the density-weighted Favre-averaged velocity.
The Favre-averaging procedure considerably reduces the number of
additional terms in the other flow equations arising from density fluctua-
tions. For example, Reynolds averaging of the convective term
ρ
u
i
u
j
in the
momentum equation gives
=
= 4R
i
R
j
+ 4 + R
i
+ R
j
+ (12.135)
Favre averaging of the same term leads to
=
= 4U
i
U
j
+ (12.136)
This clearly highlights the reduced number of unknown correlations
(products of fluctuating quantities) which is the main advantage of Favre
averaging. However, care must be taken when comparing results obtained
by solving Favre-averaged equations with experimental data which are often
time averaged. Therefore, conversion of Favre-averaged quantities to time-
averaged quantities is necessary. In order to do this we need to know more
about the turbulent fluctuations first. In section 12.19 we look at descriptions
of turbulent fluctuations in terms of probability density functions and address
the conversion problem.
Favre averaging of other governing equations, i.e. momentum, energy,
scalar transport and species transport, yields the same form of equations as
those for turbulent constant density flows. Without presenting the deriva-
tions we give the set of Favre-averaged equations used to model turbulent
combusting flows as follows:
ρ
u″
i
u″
j
ρ
(U
i
+ u″
i
)(U
j
+ u″
j
)
ρ
u
i
u
j
ρ
′u′
i
u′
j
ρ
′u′
i
ρ
′u ′
j
u′v′
(
ρ
+
ρ
′)(R
i
+ u′
i
)(R
j
+ u ′
j
)
ρ
u
i
u
j
∂
(4U
i
)
∂
x
i
∂
4
∂
t
ρ
u″
ρ
u″
ρ
u
ρ
u
4
ANIN_C12.qxd 29/12/2006 04:44PM Page 377
378 CHAPTER 12 CFD MODELLING OF COMBUSTION
Continuity
+ 4U
i
= 0 (12.137)
Momentum
(4U
i
) + (4U
i
U
j
) =− + (8
ij
− ) (12.138)
Enthalpy
(4â) + (4U
j
â) =Γ
h
+ D
h
(12.139)
where Γ
h
= (
µ
/
σ
+
µ
t
/
σ
h
),
σ
h
= turbulent Prandtl number
Mixture fraction in conserved scalar transport model for combustion
(47) + (4U
j
7) =Γ
ξ
(12.140)
where Γ
ξ
= (
µ
/
σ
+
µ
t
/
σ
ξ
),
σ
ξ
= turbulent Schmidt number
In the two equations (12.139) and (12.140) the familiar gradient diffusion
treatment has been used to model averages of products of fluctuating quan-
tities, which we first encountered in section 3.7, where it was used to model
the diffusion of turbulence quantities k and
ε
.
Species conservation (used in more detailed combustion models) is
given by
(4J
k
) + (4U
j
J
k
) = 4D
k
−+4D
k
+ Ä
k
(12.141)
where Ä
k
is the Favre-averaged reaction rate.
Again applying gradient diffusion assumption gives
=
where
σ
k
is the turbulent Schmidt number for species k, and like other
equations we may write the transport equation for species as
(4J
k
) + (4U
j
J
k
) =Γ
k
+ 4D
k
+ Ä
k
(12.142)
where Γ
k
= (
µ
/
σ
+
µ
t
/
σ
k
). Note here that subscript k stands for species and
should not be confused with the similar effective diffusion coefficient used in
the turbulent kinetic energy equation.
D
E
F
∂
Y″
k
∂
x
j
A
B
C
∂
∂
x
j
D
E
F
∂
J
k
∂
x
j
A
B
C
∂
∂
x
j
∂
∂
x
j
∂
∂
t
∂
J
k
∂
x
j
µ
t
σ
k
ρ
Y″
k
u
i
″
D
E
F
∂
Y ″
k
∂
x
j
A
B
C
∂
∂
x
j
D
E
F
ρ
Y″
k
u
j
″
∂
J
k
∂
x
j
A
B
C
∂
∂
x
j
∂
∂
x
j
∂
∂
t
D
E
F
∂
7
∂
x
j
A
B
C
∂
∂
x
j
∂
∂
x
j
∂
∂
t
D
E
F
∂
â
∂
x
j
A
B
C
∂
∂
x
j
∂
∂
x
j
∂
∂
t
ρ
u″
i
u″
j
∂
∂
x
j
∂
Q
∂
x
i
∂
∂
x
j
∂
∂
t
∂
∂
x
i
∂
4
∂
t
ANIN_C12.qxd 29/12/2006 04:44PM Page 378
12.17 CFD OF TURBULENT NON-PREMIXED COMBUSTION 379
It should be noted that the Favre-averaged forms of the continuity,
enthalpy and mixture fraction equations do not include terms involving
density fluctuations. The Reynolds stress terms in the momentum equations
are modelled using the same turbulence modelling approach used for time-
averaged equations. For example, the k–
ε
model may be applied with Favre-
averaged transport equations for k and
ε
. The species equation has been
simplified as much as possible, but the involvement of species and density in
the reaction rate term still poses a problem.
The main problem in turbulent combusting flows arises from the averag-
ing of the species generation term N
k
. In section 12.11 it was shown that for
a simple single-step reaction of the form
Fuel + s Oxident → (1 + s)Products (12.143)
M
fu
= AT
α
C
fu
C
s
ox
exp (12.144)
For the simplest case where s = 1, and after converting concentrations into
mass fractions (see (12.73)), the generation or consumption of the species
fuel may be written as
M
fu
= A′
ρ
2
T
α
Y
fu
Y
ox
exp (12.145)
where ‘fuel’ and ‘oxidant’ are the reactants, and A′ is an appropriately
modified constant. This expression can be written in a slightly different
form as
M
fu
= A′
ρ
2
T
α
Y
fu
Y
ox
exp (12.146)
where T
A
= E
a
/R
u
is called the activation temperature and
ρ
is density.
As the reaction rate is highly non-linear, the average reaction rate N
fu
cannot
be easily expressed as a function of Favre-averaged mass fractions J
fu
, J
ox
,
the mean density 4 and mean temperature E. A Taylor series expansion may
be used for the exponential term as an attempt to obtain an averaged expres-
sion for N
fu
. Using such an approach, the averaged reaction rate may be
expressed as (Veynante and Vervisch, 2002)
N
fu
=−A4
2
H
α
J
fu
J
ox
exp
× 1 + + (P
1
+ Q
1
) +
+ (P
2
+ Q
2
+ P
1
Q
1
) ++. . . (12.147)
In this expression P
1
, Q
1
, P
2
, Q
2
etc. are terms of a series defined in Veynante
and Vervisch (2002). We do not intend to describe the details here; the
expression is shown just to highlight the complexity of the attempt to find an
J
K
L
D
E
F
Y ″
ox
T ″
2
J
ox
H
2
Y ″
fu
T ″
2
J
fu
H
2
A
B
C
D
E
F
Y″
ox
T ″
J
ox
H
Y ″
fu
T ″
J
fu
H
A
B
C
Y″
fu
Y ″
ox
J
fu
J
ox
G
H
I
D
E
F
−T
A
H
A
B
C
D
E
F
−T
A
T
A
B
C
D
E
F
−E
a
R
u
T
A
B
C
D
E
F
−E
a
R
u
T
A
B
C
ANIN_C12.qxd 29/12/2006 04:44PM Page 379
average. It can be seen that the averaged N
fu
term in this way involves corre-
lations such as Y″
fu
T ″, Y ″
ox
T ″, Y ″
fu
T ″
2
etc. and many others that are unknown
and have to be modelled. The above expression is for a single-step reaction.
When realistic chemical schemes are introduced involving many reactions
and species, it is not possible to model these correlations. Therefore, a
considerable amount of effort in turbulent combustion modelling has been
directed towards the development of models which avoid the solution of
averaged species mass fraction equations (12.142).
Given the difficulties associated with modelling of the source term in the
Favre-averaged species transport equation (12.142) it is useful to examine
the effect of turbulent fluctuations on the SCRS model. It was shown earlier
how the fast chemical reaction using a single-step reaction could be used to
develop a transport equation for the mixture fraction. In turbulent flows the
mixture fraction equation (12.102) applies in Favre-averaged form (12.140).
A Favre-averaged equation for enthalpy is also solved in cases where radi-
ation effects and other heat losses are significant. However, the calculation
of mean species and temperature using the field values of 7 and â is not
as straightforward as in the laminar case. The linear relationships between
species mass fractions and mixture fraction relate to instantaneous mixture
fraction
ξ
and not to 7. The same applies to the temperature–enthalpy
relationship. To calculate mean values of J
i
and H we need to know the
statistics of the variables (T, Y
i
,
ρ
) as a function of
ξ
. This is where an
approach known as the probability density function method is used in
turbulent combustion calculations.
We introduce a statistical approach to calculate mean quantities using the
probability density function (pdf) for fluctuating scalars (in this case the mix-
ture fraction). Mathematical descriptions of probability density functions
can be found in popular mathematics textbooks (see e.g. Evans et al., 2000)
and we do not attempt to describe the fundamentals here. The reader is also
referred to Kuo (1986) for further details and applications to the theory of
turbulent combustion.
For a random variable
φ
, the probability function F
φ
(
ψ
) is defined as
F
φ
(
ψ
) = Prob{
φ
<
ψ
} (12.148)
where Prob{
φ
<
ψ
} is the probability that
φ
is smaller than given value
ψ
.
Then the probability density function P
φ
(
ψ
) is defined as the deriva-
tive of the distribution function F
φ
(
ψ
):
P
φ
(
ψ
) = (12.149)
The properties of a probability density function or pdf are
(i) P
φ
(
ψ
) ≥ 0
(ii) P
φ
(
ψ
)d
ψ
= 1
∞
−∞
dF
φ
(
ψ
)
d
ψ
380 CHAPTER 12 CFD MODELLING OF COMBUSTION
SCRS model
for turbulent
combustion
12.18
Probability
density function
approach
12.19
ANIN_C12.qxd 29/12/2006 04:44PM Page 380
12.19 PROBABILITY DENSITY FUNCTION APPROACH 381
The probability density function P(
ψ
) can be used to calculate the mean of
any property q(
ψ
) which also depends on
ψ
in the following manner (see
Kuo, 1986):
á = q(
ψ
)P(
ψ
)d
ψ
(12.150)
In the context of combustion calculations we need to define the density-
weighted probability density function G(
ξ
) for the mixture fraction
ξ
at every
location. The density-weighted pdf G(
ξ
) and unweighted pdf P(
ξ
) are related
as follows:
G(
ξ
) = P(
ξ
) (12.151)
According to Jones and Whitelaw (1982), the density-weighted average of
any scalar quantity
φ
, which is itself a function of
ξ
, may now be obtained
from
É =
φ
(
ξ
)G(
ξ
)d
ξ
(12.152)
For example, the density-weighted average of the mass fraction of species i
is given by
J
k
= Y
k
(
ξ
)G(
ξ
)d
ξ
(12.153)
In the SCRC model the relationship between Y
k
and
ξ
, i.e. Y
k
= Y
k
(
ξ
), is
known from the fast chemistry assumptions (equations (12.103)–(12.106))
or in graphical form as shown in Figure 12.2. To obtain the density-weighted
average J
k
it is necessary to perform the integration (12.153), for which we
need to know the pdf G(
ξ
) at all locations. The form of the pdf varies for dif-
ferent types of flows. Figure 12.6 (see Bilger, 1980) shows observed pdfs at
1
0
1
0
ρ
(
ξ
)
4
∞
−∞
Figure 12.6 Schematic forms of
pdfs in a turbulent jet flame
ANIN_C12.qxd 29/12/2006 04:44PM Page 381
382 CHAPTER 12 CFD MODELLING OF COMBUSTION
various locations in a jet flame. Different analytical pdfs have been used
to approximate measured distributions, but the clipped Gaussian and beta
functions have been most successful. The interested reader is referred to
Bilger and Kent (1974), Lockwood and Naguib (1975), Bilger (1976), Pope
(1976), Lockwood and Monib (1980) and Pope (1985), among others, for
further details. The beta probability function for G(
ξ
) is currently the most
popular among modellers and has been incorporated in combustion model-
ling procedures of all the leading commercial CFD codes. Some details of
the beta function approach are given below.
A further point to note is that problems arise when CFD results from
Favre-averaged flow equations are compared with time-averaged experimen-
tal data. Suppose, however, that we know the pdf G(
ξ
) and the relationships
φ
(
ξ
) and
ρ
(
ξ
). Now we can replace u by
φ
/
ρ
in definition (12.130). This
enables us to compute the time-averaged mean of scalar 2 and time-averaged
density 4 as follows:
2 = 4G(
ξ
)d
ξ
(12.154)
and, using
φ
= 1 in (12.154),
4 = d
ξ
−1
(12.155)
The beta pdf (
β
-pdf ) is defined as
G(
ξ
) = (12.156)
Here Γ(z) denotes the gamma function (see Evans et al., 2000; Abramowitz
and Stegun, 1970). The shape of the
β
-pdf depends on the values of para-
meters
α
and
β
, which must both be positive (
α
> 0 and
β
> 0). This is
illustrated in Figure 12.7, which shows that it is possible to vary
α
and
β
to
give a reasonable match with measured pdfs such as those in Figure 12.6.
It can be shown that
α
and
β
can be determined from the mean and vari-
ance of
ξ
as follows:
α
= 7 − 1
β
= (1 − 7) (12.157)
where 7 is the Favre mean of
ξ
and 7″
2
is the Favre-averaged variance of
ξ
.
DNS calculations of pdfs in turbulent combusting flows also support the
use of a
β
-pdf: see for example Swaminathan and Bilger (1999). The reader
is also referred to Warnatz et al. (2001), Libby and Williams (1994), Lentini
α
7
J
K
L
7(1 − 7)
7″
2
G
H
I
ξ
α
−1
(1 −
ξ
)
β
−1
Γ(
α
+
β
)
Γ(
α
)Γ(
β
)
J
K
L
G(
ξ
)
ρ
(
ξ
)
1
0
G
H
I
φ
(
ξ
)
ρ
(
ξ
)
1
0
Beta pdf12.20
ANIN_C12.qxd 29/12/2006 04:44PM Page 382
12.20 BETA PDF 383
(1994) and Liu et al. (2002) for further details on pdf-based combustion
calculations.
To generate the
β
-pdf using (12.156) and (12.157) it is necessary to solve
the Favre-averaged equation for mean mixture fraction 7 (12.140) as well
as a further transport equation for mixture fraction variance 7″
2
. Without
derivation we write the modelled transport equation for 7″
2
in suffix notation
(see Lockwood, 1977):
(47″
2
) + (4U
j
7″
2
) =Γ
ξ
+ C
g1
µ
t
2
− C
g2
ρ
7″
2
(12.158)
On the left hand side of this equation we have the transient and convective
terms. All terms on the right hand side are modelled. The first term rep-
resents turbulent diffusion of mixture fraction variance along its gradients.
We note that this model takes the diffusion coefficient for mixture fraction
variance to be the same as the diffusion coefficient for mixture fraction. The
second and third terms are the source and sink terms, respectively, where
C
g1
and C
g2
are dimensionless model constants with values of 2.0 and 2.8
respectively.
Combining the
β
-pdf with the definition (12.152) of the density-weighted
mean É of any scalar flow variable
φ
we may now write
É =
φ
(
ξ
)G(
ξ
)d
ξ
=
φ
(
ξ
)
ξ
α
−1
(1 −
ξ
)
β
−1
d
ξ
(12.159)
If the relationship between variable
φ
and the instantaneous mixture frac-
tion
ξ
is somehow known (e.g. Figure 12.2 for species mass fractions),
1
0
Γ(
α
+
β
)
Γ(
α
)Γ(
β
)
1
0
6
ã
D
E
F
∂
7″
2
∂
x
j
A
B
C
D
E
F
∂
7″
2
∂
x
j
A
B
C
∂
∂
x
j
∂
∂
x
j
∂
∂
t
Figure 12.7 Behaviour of the
beta function for different values
of
α
and
β
ANIN_C12.qxd 29/12/2006 04:44PM Page 383
equation (12.159) can be numerically integrated for specified values of
α
and
β
. We normally use Romberg’s method with the midpoint approximation
(Press et al., 1993). At the end points of the integration interval (
ξ
= 0, 1) the
integrand becomes singular when the parameters
α
and
β
are less than 1.0.
This singularity can be eliminated analytically using the method suggested
by Bray et al. (1994) and Chen et al. (1994), which involves approximation of
the integration as follows:
É =
φ
(
ξ
)G(
ξ
)d
ξ
≅
φ
(0) +
φ
(
ξ
)
ξ
α
−1
(1 −
ξ
)
β
−1
d
ξ
+
φ
(1) (12.160)
where
η
is a very small number (say 10
−30
).
Another numerical difficulty is that computed values of
α
and
β
in equa-
tions (12.157) may approach magnitudes of several hundred thousand in the
iteration process (Chen et al., 1994). This problem leads to overflow in the
calculation of G(
ξ
). According to the characteristics of the beta function,
G(
ξ
) is close to a delta function when either values of exponents
α
and
β
are
sufficiently large. To avoid the overflow when
α
and
β
are large (say above
500) we approximate G(
ξ
) by a delta function: G(
ξ
) =
δ
(
ξ
− 7). The density-
weighted mean É of scalar flow variable
φ
is now given by
É =
φ
(
ξ
)G(
ξ
)d
ξ
=
φ
(
ξ
)
δ
(
ξ
− 7)d
ξ
=
φ
(7) (12.161)
Since the pdf is an assumed one – the
β
-pdf – this type of calculation is
sometimes called the presumed pdf approach. There are more elaborate
methods which solve transport equations to obtain the pdf. The interested
reader is referred to Pope (1990) and Dopazo (1994) for details. Although
these models are perceived to be more accurate they are currently too com-
putationally expensive to be used in engineering calculations. In contrast, the
simple fast chemistry model with presumed
β
-pdf provides reasonable esti-
mates of temperature and major species in turbulent combusting flows. The
main assumption in SCRS is that combustion takes place due to a one-step
irreversible reaction with fast chemistry and complete combustion. Fuel and
oxygen combine and burn completely on the stoichiometry surface: therefore
the model is also known as the flame-sheet model. The stoichiometric sur-
face is also the maximum temperature surface in the flame. The main disad-
vantage of the flame-sheet model is the fact that intermediate products such
as CO and H
2
are not predicted. The absence of endothermic dissociation
reactions generating the minor species may lead to a substantial overpredic-
tion of temperature and major species. A discussion of the modifications to
this model to include minor species is available in Peters (1984).
In section 12.8 we discussed how chemical equilibrium is used to calculate
species concentrations using simple chemical reactions. The methodology
1
0
1
0
η
β
α
1−
η
η
η
α
α
1
0
384 CHAPTER 12 CFD MODELLING OF COMBUSTION
The chemical
equilibrium model
12.21
ANIN_C12.qxd 29/12/2006 04:44PM Page 384
12.22 EDDY BREAK-UP MODEL OF COMBUSTION 385
can be extended to turbulent combustion calculations by considering inter-
mediate reactions. Dedicated chemical equilibrium calculation pro-
grams, such as CHEMKIN, can be employed to predict equilibrium species
concentrations including minor species. Species concentration profiles as
a function of mixture fraction can be generated with the aid of such equilib-
rium programs and used as an alternative to fast chemistry relationships
(Peters, 1984; Warnatz et al., 2001).
This method has been successfully used by Kent and Bilger (1973)
to predict hydrocarbon flames. Application of the equilibrium model to
gas turbine combustors by Jones and Priddin (1978) has shown overpredic-
tion of CO and H
2
levels in fuel-rich regions. This is caused by the fact
that the local turbulent and diffusion time scales in practical combustor
applications are much smaller than the time required to achieve equilib-
rium. Hence, predictions based on the built-in assumption that the minor
species reactions reach equilibrium tend to overestimate minor species
levels. As a general rule, the chemical equilibrium model should only be
used in situations where the residence time is sufficiently long. However,
it is useful since minor species can be predicted and implementation
is straightforward. An alternative to the equilibrium model is the partial
equilibrium model, which assumes partial equilibrium for some species and
non-equilibrium for others. The reader is referred to Eickhoff and Grethe
(1979) for further details.
Another simple and very efficient model used in combustion calculations is
the eddy break-up model due to Spalding (1971). In the eddy break-up
model, the rate of consumption of fuel is specified as a function of local flow
properties. The mixing-controlled rate of reaction is expressed in terms of
the turbulence time scale k/
ε
, where k is the turbulent kinetic energy and
ε
is the rate of dissipation of k. The reaction rate is equal to the turbulent
dissipation rate, which, for fuel, oxygen and products, may be expressed as
follows:
Ä
fu
=−C
R
4J
fu
(12.162)
Ä
ox
=−C
R
4 (12.163)
Ä
pr
=−C ′
R
4 (12.164)
Note that the above expressions describe the Favre-averaged reaction rates.
The eddy break-up model solves one transport equation for the mass frac-
tion of fuel J
fu
. The individual dissipation rates (12.162)–(12.164) of fuel,
oxygen and products are considered, and the model takes the actual reaction
rate of fuel to be equal to the slowest of these dissipation rates:
Ä
fu
=−4 min C
R
J
fu
,C
R
, C ′
R
(12.165)
J
K
L
J
pr
1 + s
J
ox
s
G
H
I
ε
k
ε
k
J
pr
(1 + s)
ε
k
J
ox
s
ε
k
Eddy break-up
model of
combustion
12.22
ANIN_C12.qxd 29/12/2006 04:44PM Page 385
where C
R
and C ′
R
are model constants. Using a different constant for the
products allows the existence of the burnt gases to be taken into account.
Typical values used in the literature are C
R
= 1.0 and C′
R
= 0.5. Since the
above expressions contain mass fractions of fuel, oxidant and products,
initialisation of these mass fractions is required to start the calculation. In
addition to the equation for J
fu
a transport equation for mixture fraction 7 is
also solved to deduce the product and oxygen mass fractions using relation-
ships such as (12.103), (12.104) and (12.106).
Figures 12.8a and b show the results of Magnussen and Hjertager (1976),
who obtained good predictions of the temperature field in furnace configura-
tions with the eddy break-up model. Figure 12.9 shows a further application
of the eddy break-up approach combined with the pdf method to account for
scalar fluctuations by Gosman et al. (1978) in furnace simulations, and again
the prediction compares very well with the experimental data.
386 CHAPTER 12 CFD MODELLING OF COMBUSTION
Figure 12.8 Results of the eddy break-up model (Magnussen and Hjertager, 1976): (a) comparison of experimental
(Lockwood and Odidi, 1975) and eddy break-up model predictions of local mean temperatures of a city gas diffusion
flame (Re = 24000); (b) experimental mean temperatures on the axis of the city gas diffusion flame (Re = 24000)
compared with prediction by Lockwood and Naguib (1975) and the prediction of the eddy break-up model
Source: Magnussen and Hjertager (1976)
The eddy break-up model can also accommodate kinetically controlled
reaction terms. When the combustion processes are kinetically controlled,
the fuel dissipation rate may be expressed by the Arhennius kinetic rate
expression
Ä
fu,kinetic
=−A
1
4
a
J
b
fu
J
c
ox
exp(−E
a
/R
u
H) (12.166)
where A
1
is the pre-exponential constant for the Arrhenius reaction rate; a,
b and c are model constants; T is the temperature in K; E
a
is the activation
energy; and R
u
is the universal gas constant.
ANIN_C12.qxd 29/12/2006 04:44PM Page 386