April 26, 2002 8:8
Proceedings of IMECE’02
2002 ASME International Mechanical Engineering Congress and Exhibition
November 17-22, 2002, New Orleans, Louisiana, USA
IMECE02/2-6-2-11
IMPORTANCE OF TURBULENCE-RADIATION INTERACTIONS IN TURBULENT
REACTING FLOWS
Genong Li
Fluent Incorporated
10 Cavendish Court
Lebanon, NH 03766
Email: gnl@fluent.com
Michael F. Modest
Department of Mechanical Engineering
The Pennsylvania State University
University Park, PA 16802
mfm@mara.me.psu.edu
ABSTRACT
Traditional modeling of radiative transfer in reacting flows
has ignored turbulence-radiation interactions (TRI). Radiative
fluxes, flux divergences and radiative properties have been based
on mean temperature and concentration fields. However, both
experimental and theoretical work have suggested that mean ra-
diative quantities may differ significantly from those predictions
based on the mean parameters because of their strongly non-
linear dependence on the temperature and concentration fields.
The composition PDF method is able to consider many nonlin-
ear interactions rigorously, and the method is used here to study
turbulence-radiation interactions. This paper tries to answer two
basic questions: (1) whether turbulence-radiation interactions are
important in turbulent flames or not; (2) if they are important,
then what correlations need to be considered in the simulation
to capture them. After conducting many flame simulations, it
was observed that, on average, TRI effects account for about 1/3
of the total drop in flame peak temperature caused by radiative
heat losses. In addition, this study shows that consideration of
the temperature self correlation alone is not sufficient to capture
TRI, but that the complete absorption coefficient–Planck func-
tion correlation must be considered.
NOMENCLATURE
a weight factor in FSK model
C
φ
model constant
C
µ
constant in turbulence modeling
d
j
jet diameter
f probability density function
f
rad
radiant fraction
G incident radiation
∆H
comb
heat of combustion
I radiative intensity
I
b
Planck function
J
α
i
molecular diffusive flux of α-th composition variable
k spectral dependence part of absorption coefficient
k turbulent kinetic energy
L flame length
˙m
fuel
mass flow rate of the fuel
q
R
radiative heat flux
˙
Q
em
radiative emission
˙
Q
net
net radiative heat loss
~
s directional vector
S
radiation
source due to radiation
T temperature
u spatial dependence part of absorption coefficient
u
c
coflow air velocity
u
j
jet velocity
u
p
pilot flow velocity
w numerical quadrature weight
x
i
space variable in ith direction
Y species concentration vector
Greek
κ absorption coefficient
κ
p
Planck mean absorption coefficient
ρ mixture density
Ω solid angle
1 Copyright 2002 by ASME
η wave number
φ composition variables
ψ sample space of the composition variables
Γ
T
turbulent diffusivity
ω turbulent mixing frequency, = ε/k
INTRODUCTION
In turbulent reacting flows the turbulent fluctuations of the
flow field cause fluctuations of species concentrations and tem-
perature. Consequently, the radiation field, which is determined
by species concentration and temperature fields, will fluctuate as
well. In a numerical simulation fluctuations in the radiation field
interact with the fluctuations of the flow field, causing the so-
called turbulent-radiation interactions (TRI). It has been a great
challenge to consider these interactions numerically simulation
because they are strongly nonlinear in nature.
For several decades, radiation and turbulence were treated
as independent phenomena, and radiative heat fluxes were com-
puted neglecting fluctuations in the radiative intensity and in ra-
diative properties [1]. Some early simple numerical analyses and
limited experimental data have indicated the importance of these
correlations. Through a Taylor series expansion of the Planck
function, Cox [2] estimated that the contribution from temper-
ature fluctuations to radiative emission may dominate the con-
tribution from the mean temperature field when the temperature
fluctuation intensity exceeds approximately 40%. Gore et al. [3]
showed through experiments that actual radiative fluxes may be
two times or more larger than would be expected based on the
mean values alone. In the late eighties, some researchers [4–8]
performed numerical simulations taking turbulence-radiation in-
teractions (TRI) into account in some simplified fashion, and
their predictions were observed to match better with experimen-
tal data. In these early studies, either correlations for the turbu-
lent medium or the shape of the PDF had to be assumed. As a
result, turbulence-radiation interactions could not be rigorously
considered and many claims that were made about TRI need to
be further examined.
Probability density function (PDF) methods have the unique
feature that many nonlinear interactions can be treated ex-
actly [9], and have been widely used in the modeling of react-
ing flows in the absence of radiation, in which the chemical
reactions, no matter how complicated they are, can be consid-
ered exactly [10, 11]. Such methods have been introduced to
the study of turbulence-radiation interactions by Mazumder and
Modest [12] and by Li and Modest [13]. Mazumder and Mod-
est [12] employed the velocity-composition joint PDF method in
their simulation of a bluff body combustor and found inclusion of
the absorption coefficient–temperature correlation alone may in-
crease radiative heat flux by 40-45%. The inclusion of velocities
and time scale information within the PDF, although allowing
closure of more terms, adds further mathematical complexities
to the modeling of the PDF equation as well as stability prob-
lem in the numerical simulations. For the purpose of capturing
TRI, the composition PDF method is as rigorous as the velocity-
composition joint PDF method, but computationally more robust
and more efficient. Its use in the study of TRI was demonstrated
by Li and Modest [13]. By employing the same method, this pa-
per aims to check the importance of turbulence-radiation interac-
tions, and the relative importance of the different contributions to
TRI. Since the Planck function is the most nonlinear function in
the radiation calculation, it has been hypothesized that considera-
tion of the temperature self correlation alone can capture most of
the TRI [14]. If this were the case, one could treat TRI with the
traditional Reynolds average approach, constructing the first few
higher moments of temperature. Such issues will be discussed in
this paper.
MATHEMATICAL FORMULATION
Turbulence-radiation coupling
In the presence of radiative heat transfer, the energy equation
needs to include a radiative source term,
S
radiation
= −∇ · q
R
=
Z
∞
0
κ
η
Z
4π
I
η
dΩ − 4πI
bη
!
dη, (1)
where q
R
denotes the radiative heat flux; κ
η
is the spectral ab-
sorption coefficient of the radiating gas, which may be a func-
tion of temperature T and species concentrations of the radiating
medium Y; here I
η
is the spectral intensity, I
bη
is the spectral
blackbody intensity (or Planck function), the subscript η is used
to indicate spectral dependence and Ω denotes solid angle. The
radiation intensity is governed by the radiative transfer equation
(RTE): for an absorbing-emitting but nonscattering gas, the in-
stantaneous radiant energy balance on a pencil of radiation prop-
agating in direction
~
s and confined to a solid angle dΩ is given
by [15],
(
~
s· ∇)I
η
= κ
η
(I
bη
− I
η
) (2)
where the first term on the right-hand side represents augmenta-
tion due to emission and the second term is attenuation due to
absorption.
To include radiation effects in conventional turbulence cal-
culations, Eqs. (1) and (2) need to be time-averaged, resulting
in
hS i
radiation
=
Z
∞
0
"
Z
4π
hκ
η
I
η
idΩ − 4πhκ
η
I
bη
i
#
dη, (3)
(
~
s· ∇)hI
η
i = hκ
η
I
bη
i − hκ
η
I
η
i. (4)
2 Copyright 2002 by ASME
Due to the strongly nonlinear dependence of radiative prop-
erties on temperature and species concentrations, hκ
η
(T,Y)I
η
i
does not equal κ
η
(hTi,hYi)hI
η
i and hκ
η
(T,Y)I
bη
(T)i does not
equal κ
η
(hTi,hYi)I
bη
(hTi), making these two terms unclosed.
hκ
η
I
η
i represents a correlation between the spectral absorption
coefficient and the spectral incident intensity, and hκ
η
I
bη
i rep-
resents a correlation between the spectral absorption coefficient
and the spectral blackbody intensity. Complete information of
the statistics among the composition variables is needed for
their determination. For the convenience of later discussion,
these two correlations are loosely defined as ‘spectral absorption
coefficient–spectral incident intensity correlation’ and ‘spectral
absorption coefficient–spectral blackbody intensity correlation’.
The time averaging procedure can be applied to any solu-
tion technique for radiation calculation and different unclosed
terms may arise for different spectral models and solution meth-
ods. However, all of them can be categorized as belonging to
two groups: (a) correlations that can be calculated from scalars
φ directly or indirectly, and (b) correlations that cannot. The set
of scalars φ is defined as
φ = (Y, T) = (φ
1
,φ
2
,· · · , φ
s
) (5)
where s is the total number of scalar variables (number of species
plus one) and the last scalar, φ
s
, is reserved for temperature (or
enthalpy). Variables in the set φ are often called the composition
variables, since they determine the composition of the mixture.
The unclosed term hκ
η
I
bη
i belongs to group (a), since both κ
η
and I
bη
are functions of variables in set φ only. The unclosed term
hκ
η
I
η
i belongs to group (b), because I
η
is not a local quantity, i.e.,
cannot be expressed in terms of the local scalar variables.
One of the most common approximations made in the open
literature on turbulence-radiation interactions is the optically thin
eddy approximation as described by Kabashnikov and Myas-
nikova [16]. Kabashnikov suggested that if the mean free path
for radiation is much larger than the turbulence length scale, then
the local radiative intensity is weakly correlated with the local ab-
sorption coefficient, i.e., hκ
η
I
η
i = hκ
η
ihI
η
i, in which hκ
η
i is loosely
defined as the ‘absorption coefficient self correlation.’ The ra-
tionale behind these assumptions is that the instantaneous local
intensity at a point is formed over a path traversing several tur-
bulent eddies. Therefore, the local intensity is weakly correlated
to the local radiative properties. The validity of this assumption
depends on the eddy size distribution and the radiation properties
of the absorbing gases. In a numerical simulation of combustion
chambers, Hartick et al. [8] showed that, although the thin eddy
assumption may not be valid over some highly absorbing parts
of the spectrum, these spectral zones affect the total radiation ex-
change only slightly, thus allowing straightforward application
of the thin eddy assumption in their simulation. The thin eddy
assumption is also employed in the current study. As a result, all
correlations needed to capture TRI belong to group (a).
Radiation Submodel
The radiative transfer equation is a spectrally, spatially and
directionally dependent integro-differential equation, and is ex-
tremely difficult to solve for general, multi-dimensional geome-
tries. Several approaches are available to reduce this equation to
a simpler form. Among them, one of the most popular methods
is the P
1
-approximation, in which the incident radiation is gov-
erned by a Helmholtz equation, which is relatively easy to solve.
For the vast majority of important engineering problems (i.e.,
in the absence of extreme anisotropy in the intensity field), the
method provides high accuracy at very reasonable computational
cost. Another challenge in gas radiation calculations comes
from the strong spectral dependence of radiation properties. Al-
though line-by-line calculations provide best accuracy, such cal-
culations are too time-consuming for any practical combustion
system. Global methods such as the Weighted-Sum-of-Gray-
Gases Model (WSGG) are commonly used [17]. Recently, the
Full-Spectrum k-Distribution method (FSK) developed by Mod-
est and Zhang [18] has been shown to be superior to the WSGG
model, to which it reduces in its crudest implementation. The
method is exact within its limitations [gray walls, gray scattering
properties, spectral absorption coefficient obeying the so-called
scaling approximation, i.e., the spectral and spatial dependence
of the absorption coefficient are separable as κ
η
(η,φ) = k
η
(η)u(φ)
where φ are the composition variables]. The P
1
-approximation
in conjunction with the FSK model will be used in this study.
Radiative properties and, consequently, the radiative inten-
sity change dramatically across spectral space. In the FSK
method the radiative quantities’ spectral dependence has been
transformed to a g-dependence, where g is the cumulative distri-
bution function of the absorption coefficient calculated over the
whole spectrum and weighted by the Planck function. For exam-
ple, the source term in the energy equation due to radiative heat
transfer is calculated as
S
radiation
= −
Z
∞
0
κ
η
(4πI
bη
− G
η
)dη
= −
Z
1
0
k
g
u(4πa
g
I
b
− G
g
)dg (6)
where u is the spatial dependence of the absorption coefficient
as mentioned before and a
g
is a weight factor introduced during
the transformation. The advantage of this transformation lies in
the fact that k
g
(g) is a smooth, monotonically increasing func-
tion of g, thus requiring only a few numerical quadrature points.
Readers are referred to [18] for the details about this method.
In practical calculations, the integration is replaced by numerical
3 Copyright 2002 by ASME
quadrature. If Gaussian quadrature is used, Eq. (6) becomes
S
radiation
≈ −
M
X
j=1
w
j
k
j
u(4πa
j
I
b
− G
j
), (7)
where M is the total number of quadrature points and the w
j
are
the quadrature weights. The incident radiation G
j
must be deter-
mined by solving the P
1
-equation, i.e. [15],
∇ ·
1
3k
g
u
∇G
g
!
= k
g
u[G
g
− 4πa
g
I
b
] (8)
subject to the boundary condition,
−
2(2− )
3
ˆn· ∇G
g
= k
g
u(4πa
g
I
b
− G
g
) (9)
where is surface emittance and ˆn is a unit normal at a boundary
surface. The values of k
g
,u(φ) and a
g
are obtained from a pre-
calculated FSK data base.
Reynolds averaging of the radiative source term and the P
1
-
equation leads to
hS i
radiation
= −
M
X
j=1
w
j
k
j
[4πhua
j
I
b
i − huihG
j
i] (10)
∇ ·
"
1
3k
j
1
hui
∇hG
j
i
#
= k
j
huihG
j
i − 4πk
j
hua
j
I
b
i,
j = 1, · · · , M (11)
where the optically thin-eddy approximation has been em-
ployed. As a result of turbulence-radiation interactions two
terms, k
j
hua
j
I
b
i and k
j
hui, representing correlations between de-
pendent variables, need to be modeled.
Composition PDF Methods
The philosophy of the PDF approach is to treat species con-
centration and temperature as random variables and consider the
transport of their PDFs rather than their finite moments. Once
that PDF is known, the mean of any quantity can be evaluated
exactly from the PDF, as long as it is a function of the species
concentrations or/and temperature. For example,
hu(φ)a
j
(φ)I
b
(φ)i=
Z
u(ψ)a
j
(ψ)I
b
(ψ)f (ψ)dψ (12)
hu(φ)i=
Z
u(ψ)f (ψ)dψ (13)
In these equations, ψ represents the composition space variable,
ψ ≡ (ψ
1
,ψ
2
,· · · , ψ
s
), and f (ψ) is defined to be the probability den-
sity of the compound event φ = ψ (i.e.,φ
1
= ψ
1
,φ
2
= ψ
2
,· · · , φ
s
=
ψ
s
), so that,
f(ψ)dψ = Probability(ψ ≤ φ ≤ ψ + dψ) (14)
In a general turbulent reacting flow, the composition PDF is also
a function of space, x, and time, t. The transport equation for the
composition PDF, f(ψ, x,t), can be derived from the conservation
laws of scalars, which is
∂
∂t
[ρ f ]+
∂
∂x
i
[˜u
i
ρ f ]+
∂
∂ψ
α
[S
α,reaction
(ψ)ρ f ]
−
M
X
j=1
4πw
j
k
j
∂
∂ψ
s
h
ua
j
I
b
f
i
= −
∂
∂x
i
[hu
00
i
|ψiρ f ]
+
∂
∂ψ
α
"*
1
ρ
∂J
α
i
∂x
i
ψ
+
ρ f
#
−
M
X
j=1
w
j
k
j
∂
∂ψ
s
h
uhG
j
i f
i
(15)
where i and α are summation indices in physical space and com-
position space, respectively and hA|Bi is the conditional proba-
bility of the event A, given that the event B occurs.
On the left-hand side of Eq. (15), the first two terms repre-
sent the rate of change of the PDF when following the Favre-
averaged mean flow. The third term is the divergence of the flux
of probability in the composition space due to chemical reaction
and radiative emission. The form of this term clearly shows the
advantage the PDF method has over moment methods: no matter
how complicated and nonlinear these source terms are, they re-
quire no modeling. In contrast, the terms on the right-hand side
of Eq. (15) need to be modeled. The first two terms represent
transport in physical space due to turbulent convection and trans-
port in scalar space due to molecular mixing, respectively. They
are usually modeled by using the gradient-diffusion hypothesis
and a simple mixing model such as Dopazo’s model [9], respec-
tively, leading to
−hu
00
i
|ψiρ f ≈ Γ
T
∂(ρ f )
∂x
i
, (16)
*
1
ρ
J
α
i
∂x
i
ψ
+
≈
1
2
C
φ
ω(ψ
α
−
˜
φ
α
) (17)
where Γ
T
= c
µ
hρiσ
−1
φ
k
2
/ is the turbulent diffusivity, and k,, c
µ
and σ
φ
are, respectively, the turbulent kinetic energy, dissipa-
tion rate of turbulent kinetic energy, a modeling coefficient in the
standard k- turbulence model, and turbulent Schmidt or Prandtl
numbers; finally, ω = /k is a turbulence ‘frequency’ and C
φ
is a
4 Copyright 2002 by ASME
model constant. The third term on the right-hand side of Eq. (15)
is closed by invoking the optically thin eddy approximation.
As a result, the modeled transport equation for the composi-
tion mass density PDF function is closed and contains all neces-
sary information about all scalars. The composition PDF trans-
port equation is a partial differential equation in (4+ s) dimen-
sions. Traditional finite volume or finite element methods are
very inefficient to solve an equation of such high dimensional-
ity. Instead, the Monte Carlo method is generally used, in which
the PDF is represented by a large number of computational par-
ticles. Each particle evolves in time and space according to a set
of stochastic equations and carries with it all composition vari-
ables. The PDF is then obtained approximately as a histogram
of the particles’ properties in sufficiently small neighborhoods in
physical space, and the mean quantities are deduced statistically
by sampling the particles.
Chemical Reaction Submodel
Although PDF methods allow the use of detailed chemical
reaction mechanisms in principle, computational intractability
has limited their application. In practice, reduced mechanisms
are often used in the PDF calculations. A wide range of reduced
mechanisms of chemical reactions for hydrocarbon fuels is avail-
able in the literature [19] and the simplest –a single-step skeletal
mechanism– is used in this study. It takes the form:
CH
4
+ 2O
2
→ CO
2
+ 2H
2
O. (18)
Westbrook and Dryer [20] provided an Arrhenius relationship for
the reaction rate of methane as:
d[CH
4
]
dt
= −Aexp(−E
a
/R
u
T)[CH
4
]
a
[O
2
]
b
, (19)
where the quantities within square brackets represent molar con-
centrations; R
u
is the universal gas constant, and E
a
is the ac-
tivation energy of the methane. A, a and b are constants in the
general Arrhenius equation, which may be obtained from West-
brook and Dryer.
FLAME SIMULATIONS
Flame optical thickness has an important impact on radiative
transfer. Three jet flames with different optical thickness have
been considered, where flame optical thickness has be defined as
τ = κ
P
L (20)
where κ
P
is an average Planck mean absorption coefficient of the
participating medium, i.e., the combustion products of H
2
O and
CO
2
, and L is the flame length. For turbulent jet flames, flame
length is approximately a linear function of jet diameter [21] and,
in this study, is estimated to be L = 40d
j
. The base flame is San-
dia’s Flame D [22]. The basic experimental setup of this flame
is summarized here. The fuel jet (d
j
= 7.2mm) with high ve-
locity (u
j
= 49.6m/s) is accompanied by an annular pilot flow
(d
p
= 18.4mm, u
p
= 11.4m/s), which is then surrounded by a
slow coflow of air (u
c
= 0.9m/s). The fuel is a mixture of air and
methane with a ratio of 3:1 by volume. A bank of measured data
is available for this flame. Detailed information of code valida-
tion in the simulation of this flame is reported elsewhere [23] and
not repeated in this paper. Flame optical thickness for this flame
is 0.237 by Eq. (20). The other two considered (artificial) flames
were derived from Flame D by doubling and quadrupling the jet
diameter, and their flame optical thickness is 0.474 and 0.948, re-
spectively. For future reference the three flames will be denoted
as κL.1, κL.2 and κL.3, respectively.
To simulate these flames, a rectangular axisymmetric com-
putational domain of 70d
j
×18d
j
was used, and a non-uniform
grid system of 60×70 was found to be fine enough to give grid-
independent solutions in the finite volume code. The global time
step used in the PDF/particle code was 2.0ms, and 4.0ms and
8.0ms for flames κL.1, κL.2 and κL.3, respectively. For each sim-
ulation a total of approximate 1100 iterations was required to get
to a statistically stationary result and about 58,000 particles were
used in the simulation, taking about 22 cpu hours on a four pro-
cessor Silicon Graphics O200 machine. The conventional way to
define residual error in finite volume methods is meaningless in
the hybrid FV/PDF Monte Carlo simulation because the statisti-
cal error is generally larger than the truncation error. In the cur-
rent study, the overall numerical error for a variable φ after the
j-th iteration is defined as err = 1/N
P
N
i=1
[φ
j
i
− φ
j−1
i
]
2
/[φ
j−1
i
]
2
,
where N is the total number of nodal points. This error never
converges to zero, but rather to a value representative of the sta-
tistical fluctuation of the solution when steady state is reached.
This level mainly depends on the number of particles in the sim-
ulation. If temperature is used to monitor the numerical error, a
value on the order of 10
−4
has been reached in the calculations.
Importance of TRI
In order to study turbulence-radiation interactions, three dif-
ferent scenarios were considered for each flame. In the first sce-
nario, radiation is completely ignored in order to study the im-
portance of radiation in flame simulations in general. In the sec-
ond and third scenarios, radiation is considered but turbulence-
radiation interactions are ignored and considered, respectively.
The importance of turbulence-radiation interactions can be as-
sessed by comparing numerical results from these two scenarios.
By ignoring turbulence-radiation interactions, it is implied that
the two unclosed terms hui and huaI
b
i are evaluated based on the
cell means; when by considering it, these two terms are treated
5 Copyright 2002 by ASME
