# Consistent Conditional Moment Closure Modelling of a Lifted Turbulent Jet Flame Using the Presumedβ-PDF Approach

Authors: Ahmad El Sayed; Roydon A. Fraser
Journal: Journal of Combustion
 (2014)

---

## Abstract

A liftedH
2
/
N
2 turbulent jet flame issuing into a vitiated coflow is investigated using the conditional moment closure. The conditional velocity (CV) and the conditional scalar dissipation rate (CSDR) submodels are chosen such that they are fully consistent with the moments of the presumed β probability density function (PDF). The CV is modelled using the PDF-gradient diffusion model. Two CSDR submodels based on the double integration of the homogeneous and inhomogeneous mixture fraction PDF transport equations are implemented. The effect of CSDR modelling is investigated over a range of coflow temperatures (T
c) and the stabilisation mechanism is determined from the analysis of the transport budgets and the history of radical build-up ahead of the stabilisation height. For all T
c, the balance between chemistry, axial convection, and micromixing, and the absence of axial diffusion upstream of the stabilisation height indicate that the flame is stabilized by autoignition. This conclusion is confirmed from the rapid build-up of H
O
2 ahead of H, O, and OH. The inhomogeneous CSDR modelling yields higher dissipation levels at the most reactive mixture fraction, which results in longer ignition delays and larger liftoff heights. The effect of the spurious sources arising from homogeneous modelling is found to be small but nonnegligible, mostly notably within the flame zone.


---

## Body

## 1. Introduction

The investigation of stabilisation mechanisms in lifted flames is an active area in combustion research. A substantial number of experimental and numerical studies have been dedicated to the understanding and analysis of the process by which lifted flames are stabilised. Several theories have been proposed in the literature. Vanquickenborne and van Tiggelen [1] suggest the premixed flame propagation stabilisation theory. Their experimental findings indicate that the fuel and the oxidiser are premixed ahead of the base of lifted diffusion flames and that stabilisation takes place at stoichiometric locations where the local mean velocity and the turbulent burning velocity of a stoichiometric premixed flame are equal. Peters and Williams [2] suggest that stabilisation is achieved by the local quenching of diffusion flamelets due to excessive straining. Cabra et al. [3] conduct experimental measurements for a lifted turbulent H
2
/
N
2 jet flame issuing into a coaxial vitiated coflow consisting of the hot combustion products of a lean premixed H
2/air flame. The coflow conditions allow for the possibility of stabilisation by means of autoignition. It is postulated that autoignition takes place as inert fuel parcels are convected downstream and become well mixed with the hot coflow. However, premixed flame propagation remains another possible stabilisation mechanism.The Cabra flame is attractive for the validation of combustion models due to its simple configuration and well-defined boundary conditions. It is extensively investigated using PDF methods [3–6]. Masri et al. [4] use the composition PDF approach coupled with the k
-
ε turbulence model and conclude that flame is largely controlled by chemical kinetics. They identify key reactions and show that the predicted liftoff height and the composition at the base of the flame are highly sensitive to the rate parameters of these reactions. Gordon et al. [6] use the same computational methodology in order to determine the underlying stabilisation mechanism. They develop two numerical indicators for the distinction between stabilisation via autoignition and premixed flame propagation. These indicators are based on the analysis of the transport budgets of convection, diffusion, and reaction and the history of radical build-up ahead of the reaction zone. A convection-reaction balance without significant contribution from axial diffusion, along with an early and rapid build-up of H
O
2 ahead of the remaining radicals, indicate that stabilisation occurs via autoignition, whereas an axial diffusion-convection balance with negligible reaction accompanied by a simultaneous build-up of all radicals points to stabilisation by means of premixed flame propagation. With the aid of these indicators, Gordon et al. [6] conclude that autoignition is the controlling mechanism for a wide range of coflow temperatures. Cao et al. [5] employ the joint velocity-turbulence frequency-composition PDF method and perform a detailed parametric study using several mixing models and chemical kinetic mechanisms, over a wide range of jet and coflow conditions. They show that although the modelling of mixing affects the prediction of liftoff height, the flame is mainly controlled by chemical kinetics. They also demonstrate that the liftoff height increases with increasing jet and coflow velocities and decreases with increasing jet and coflow temperatures. In all of the PDF studies described above, autoignition was found to be the controlling stabilisation mechanism. In the context of CMC, Stankovic and Merci [7] and Navarro-Martinez and Kronenburg [8] investigate this flame using Large Eddy Simulation (LES) with a first-order closure for the conditional chemical source. The main difference between the two LES-CMC studies lies in the treatment of the CSDR. The first employs the Amplitude Mapping Closure (AMC) [9], which produces symmetric profiles in mixture fraction space. The second employs a different conditioning technique described in [10], which allows for skewed profiles with the possibility of local minima and maxima in mixture fraction space. In both sets of calculations the coflow temperature ranges between 1020 and 1080 K. The first authors conclude that stabilisation is controlled by autoignition for all coflow temperature. The second authors draw the same conclusions for high coflow temperatures; however, they report a transition to flame stabilization by premixed flame propagation as the coflow temperature is reduced. The different conclusions are most likely attributed to the treatment of the CSDR and, to a lesser extent, to the employed numerical methods, discretization schemes, and grid resolutions in physical and mixture fraction spaces. Patwardhan et al. [11] also use the first-order CMC to study this flame. They employ the standard k
-
ε model for the computation of the turbulence and mixing fields. To a certain degree, their conclusions in regard to the controlling stabilisation mechanism are in line with those of Navarro-Martinez and Kronenburg [8]. Premixed flame propagations are found to be the controlling stabilisation mechanism of the lifted flame for low and intermediate coflow temperatures (1025 and 1045 K), whereas the flame is stabilised by autoignition for high coflow temperatures (1080 K). To close the CMC equations, they employ the β-distribution for the PDF, the linear velocity model [12] for the CV, and the AMC for the CSDR. The assumption of homogeneous Gaussian turbulence is inherent to the adopted CV and CSDR closures. The linear velocity model is exact only if the velocity and the mixture fraction are jointly Gaussian [13, 14], which is in general not the case for inhomogeneous flows. Though being consistent with the first moment of the PDF, this CV model is inconsistent with the second moment [15]. As for the AMC, this closure is a particular solution of the homogeneous PDF transport equation, wherein the employed mapping function is Gaussian. Patwardhan et al. [11] use a β-distribution instead. Further, the fact that the AMC does not account for inhomogeneous effects leads to an inconsistent CMC implementation. When such closure is employed, the integration of the PDF-weighted CMC equations over the mixture fractions space yields the unconditional equations in addition to extra spurious source terms [14, 16]. This term may lead to inaccurate results and misleading conclusions, depending on the magnitudes of these sources. As such, the CV and CSDR closures employed in [11] do not provide a fully consistent CMC implementation.In this work, the lifted hydrogen jet flame of Cabra et al. [3] is revisited and investigated using the two-dimensional first-order CMC and a modified version of the k
-
ε model. The conditional submodels required for the closure of the CMC equations are selected such that they are fully consistent with the moments of the presumed PDF. Additionally, a transport equation is solved for the Favre-averaged scalar dissipation rate rather employing the traditional “equal time scales” algebraic approach. The simulations are set up such that the settings of Patwardhan et al. [11] are reproduced as closely as possible. This study is aimed at determining (1) the effect of micromixing on liftoff via the modelling of the CSDR, (2) the response of the flame to changes in the coflow temperature, (3) the controlling stabilisation mechanism, (4) the influence of the spurious sources, and (5) the impact of chemical kinetics.


## 2. Experimental Configuration

The flame under investigation is the lifted turbulent jet flame of Cabra et al. [3]. A schematic of the experimental setup is shown in Figure 1 and the flow conditions are summarized in Table 1. The burner consists of a central H
2
/
N
2 jet surrounded with a coaxial flow composed of the hot combustion products of a lean premixed H
2/Air flame stabilized on a perforated disk. The perforated disk is surrounded by an exit collar for the purpose of delaying the entrainment of ambient air into the coflow. The nozzle exit is placed above the surface of the disk in order to allow the fuel to exit into the coflow with a uniform composition. The experimental criterion for the determination of the liftoff height (H
exp
⁡) is taken as the first location where the mass fraction of OH reaches 600 ppm. The measured height is equal to 10 nozzle diameters.Table 1
Conditions of the jet and coflow in the lifted flame of Cabra et al. [3].
Jet
Coflow
d [mm]
4.57
D [mm]
210
T
j [K]
305
T
c [K]
1045
U
j [m/s]
107
U
c [m/s]
3.5
X
H
2
0.2537
X
H
2
0.0005
X
O
2
0.0021
X
O
2
0.1474
X
N
2
0.7427
X
N
2
0.7534
X
H
2
O
0.0015
X
H
2
O
0.0989
ξ
st
0.474
ϕ
0.25
Nomenclature:d: jet diameter; D: coflow diameter; T: temperature; U: velocity; X: mole fraction; ξ: mixture fraction; ϕ: equivalence ratio. j: jet; c: coflow; st: stoichiometric.Figure 1
The liftedH
2
/
N
2 flame. Reproduced from [3].


## 3. Mathematical Model

### 3.1. Conditional Moment Closure

The modelling of the turbulence-chemistry interactions is carried out using the first-order CMC [14]. In CMC, reactive scalars such as the species mass fractions and the temperature are conditionally averaged with respect to the mixture fraction, ξ, and their conditional transport equations are solved. The conditional average of the mass fraction of a species κ is defined as
(1)
Q
κ
(
η
,
x
,
t
)
=
〈
Y
κ
(
x
,
t
)
∣
ξ
(
x
,
t
)
=
η
〉
=
〈
Y
κ
∣
η
〉
,
where Y
κ is the mass fraction of the species, 〈
·
∣
·
〉 denotes the conditional average of the quantity to the left of the vertical bar subject to the quantity to its right, and η is a sample variable of ξ, with 0
≤
η
≤
1. Following the decomposition approach [14], Y
κ is written as the sum of its conditional average Q
κ and a fluctuation y
κ
′′:
(2)
Y
κ
(
x
,
t
)
=
Q
κ
(
η
,
x
,
t
)
+
y
κ
′′
(
η
,
x
,
t
)
such that 〈
y
κ
′′
∣
η
〉
=
0. Substitution of (2) into the transport equation of species κ, followed by the conditional averaging of the resulting expression, leads to [14]:
(3)
∂
Q
κ
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
κ
∂
x
i
+
〈
χ
∣
η
〉
2
∂
2
Q
κ
∂
η
2
+
〈
ω
˙
κ
∣
η
〉
〈
ρ
∣
η
〉
+
e
Q
κ
+
e
y
κ
′′
〈
ρ
∣
η
〉
,
where
(4)
e
Q
κ
=
∂
∂
x
i
(
ρ
D
∂
Q
κ
∂
x
i
)
+
〈
ρ
D
∂
ξ
∂
x
i
∂
∂
ξ
(
∂
Q
κ
∂
x
i
)
∣
η
〉
,
(5)
e
y
κ
′′
=
-
〈
ρ
∂
y
κ
′′
∂
t
+
ρ
u
i
′′
∂
y
κ
′′
∂
x
i
-
∂
∂
x
i
(
ρ
D
∂
y
κ
′′
∂
x
i
)
∣
η
〉
,
(6)
χ
=
2
D
∂
ξ
∂
x
i
∂
ξ
∂
x
i
is the scalar dissipation rate. The species and mixture fraction diffusivities in (3)–(6) are assumed to be equal to the molecular diffusivity, D, and the Lewis number is set to unity. In (3), 〈
u
i
∣
η
〉 is the CV, 〈
χ
∣
η
〉 is the CSDR, and 〈
ω
˙
κ
∣
η
〉 is the conditional chemical source. The expressions of e
Q (see (4)) and e
y
κ
′′ (see (5)) are unclosed. They are modelled using the primary closure hypothesis. Klimenko and Bilger [14] show that, given finite Schmidt numbers, all the terms of e
Q scale as the inverse of the Reynolds number (Re). Therefore, this term is neglected in high Re applications. By neglecting the conditional density and diffusivity fluctuations, they further show that and e
y
κ
′′ reduces to
(7)
e
y
κ
′′
=
-
〈
ρ
∣
η
〉
ρ
¯
P
~
(
η
)
∂
∂
x
i
[
ρ
¯
〈
u
i
′′
y
κ
′′
∣
η
〉
P
~
(
η
)
]
,
where P
~
(
η
) is the Favre PDF of mixture fraction and 〈
u
i
′′
y
κ
′′
∣
η
〉 is the conditional turbulent flux. Using these approximations, the conditional species transport equation takes the form:(8)
∂
Q
κ
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
κ
∂
x
i
-
1
ρ
¯
P
~
(
η
)
∂
∂
x
i
[
ρ
¯
〈
u
i
′′
y
κ
′′
∣
η
〉
P
~
(
η
)
]
+
〈
χ
∣
η
〉
2
∂
2
Q
κ
∂
η
2
+
〈
ω
˙
κ
∣
η
〉
〈
ρ
η
〉
.
The left-hand side term of (8) is the rate of change of the conditional species mass fraction. The first term on the right-hand side (r.h.s.) represents transport by means of convection, the second corresponds to transport by conditional turbulent fluxes (diffusion in physical space), the third accounts for micromixing (diffusion in mixture fraction space), and the fourth is the chemical source. The conditional temperature equation is derived in a similar fashion and is given by
(9)
∂
Q
T
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
T
∂
x
i
︸
T
C
,
x
and
T
C
,
y
-
1
ρ
-
P
~
(
η
)
∂
∂
x
i
[
ρ
-
〈
u
i
′
′
T
′′
∣
η
〉
P
~
(
η
)
]
︸
T
D
,
x
and
T
D
,
y
+
〈
χ
∣
η
〉
2
{
∂
2
Q
T
∂
η
2
+
1
〈
c
p
∣
η
〉
×
[
∂
〈
c
p
∣
η
〉
∂
η
+
∑
κ
=
1
N
(
〈
c
p
,
κ
∣
η
〉
∂
Q
κ
∂
η
)
]
×
∂
Q
T
∂
η
}
T
MM
+
1
〈
c
p
∣
η
〉
〈
1
ρ
∂
p
∂
t
∣
η
〉
-
〈
ω
˙
h
∣
η
〉
〈
ρ
∣
η
〉
〈
c
p
∣
η
〉
︸
T
CS
-
〈
ω
˙
r
∣
η
〉
〈
ρ
∣
η
〉
〈
c
p
∣
η
〉
︸
T
RS
,
where 〈
ω
˙
h
∣
η
〉
=
∑
κ
=
1
N
〈
h
κ
∣
η
〉
〈
ω
˙
κ
∣
η
〉 is the chemical source term, N is the total number of species in the mixture, and 〈
ω
˙
r
∣
η
〉 is the conditional radiative source term. The unsteady pressure term in (9) is neglected due to the fact that the investigated flame is open to the atmosphere, and hence temporal pressure changes are expected very small.The unconditional (Favre) averages of the reactive scalars are calculated by integrating their conditional values weighted byP
~
(
η
) over the mixture fraction space:
(10)
ϕ
~
=
∫
0
1
〈
ϕ
∣
η
〉
P
~
(
η
)
d
η
,
ϕ
=
Y
κ
,
T
.


### 3.2. Conditional Submodels

The quantitiesP
~
(
η
), 〈
u
i
′′
y
κ
′′
∣
η
〉, 〈
u
i
′′
T
′′
∣
η
〉, 〈
u
i
∣
η
〉, 〈
χ
∣
η
〉, 〈
ω
˙
κ
∣
η
〉, and 〈
ω
˙
r
∣
η
〉 in (8) and (9) are unclosed and require further modelling. The submodels employed in this study are discussed in this section.


#### 3.2.1. Probability Density Function

The PDF is presumed using theβ-distribution. The Favre β-PDF is given by
(11)
P
~
(
η
;
v
,
w
)
=
η
v
-
1
(
1
-
η
)
w
-
1
B
(
v
,
w
)
,
where the parameters v and w are related to the mean and variance of the mixture fraction by v
=
γ
ξ
~ and w
=
γ
(
1
-
ξ
~
) with γ
=
[
ξ
~
(
1
-
ξ
~
)
/
ξ
′′
2
~
-
1
]
≥
0, and B
(
v
,
w
)
=
∫
0
1
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η is the beta function.


#### 3.2.2.  Conditional Turbulent Fluxes

The conditional turbulent fluxes are modelled using the gradient diffusion assumption:(12)
〈
ϕ
′
′
u
i
′′
∣
η
〉
=
-
D
t
∂
〈
ϕ
∣
η
〉
∂
x
i
,
ϕ
=
Y
κ
,
T
,
where
(13)
D
t
=
C
μ
S
c
t
k
~
2
ɛ
~
is the turbulent diffusivity. The constant C
μ is equal to 0.09 and the turbulent Schmidt number S
c
t is set to 0.7 following Jones and Whitelaw [17]. Richardson and Mastorakos [18] suggest adding a correction to (12) in order to account for counter-gradient transport effects. When applied to a lifted flame, they conclude that the inclusion of counter-gradient effects leads to a slight increase in liftoff height and to a decrease in flame thickness. This extension is not employed in the current study. However, it is worth investigating in future studies.


#### 3.2.3.  Conditional Velocity

The PDF-gradient model proposed by Pope [19] is employed to model the CV fluctuations. The resulting expression for the CV takes the form:
(14)
〈
u
i
∣
η
〉
=
u
~
-
D
t
P
~
(
η
)
∂
P
~
(
η
)
∂
x
i
=
u
~
i
-
D
t
{
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
~
∂
ξ
~
∂
x
i
+
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
}
.
As P
~
(
η
) tends to zero, 〈
u
i
∣
η
〉 diverges to ±
∞, depending on the sing of the gradient of P
~
(
η
). This behavior is of minor importance, because low-probability events have a negligible effect on mixing [20]. One important feature of the PDF-gradient model is its consistency with both the first and second moments of the mixture fraction [15, 20]. This does not hold for the commonly used linear model [12], which is only consistent with the first moment [15].


#### 3.2.4.  Conditional Scalar Dissipation Rate

Girimaji’s Model. Girimaji [21] derives a model for 〈
χ
∣
η
〉 by doubly integrating the homogeneous mixture fraction PDF transport equation. His formulation is based on the observation that a presumed β-PDF is capable of accurately characterizing the evolution of the scalar PDF over all stages of two-scalar, constant-density mixing in statistically stationary, isotropic turbulence. The model is given by:
(15)
〈
χ
η
〉
=
-
2
χ
~
ξ
~
(
1
-
ξ
~
)
(
ξ
′′
2
~
)
2
I
(
η
)
P
~
(
η
)
,
where χ
~ is the Favre-averaged scalar dissipation rate (discussed in more detail in Section 3.3) and I
(
η
) is given by the integral:
(16)
I
(
η
)
=
∫
0
η
{
ξ
~
[
ln
⁡
η
′
-
∫
0
1
ln
⁡
η
′′
P
~
(
η
′′
)
d
η
′′
]
+
(
1
-
ξ
~
)
×
[
ln
⁡
(
1
-
η
′
)
-
∫
0
1
ln
⁡
⁡
(
1
-
η
′′
)
P
~
(
η
′′
)
d
η
′′
]
}
×
P
~
(
η
′
)
(
η
-
η
′
)
d
η
′
.
As mentioned above, Girimaji’s model is valid for homogeneous turbulence. Therefore, its usage in the CMC of inhomogeneous flows is inconsistent and results in spurious source terms.Mortensen’s Model. In a similar fashion, Mortensen [22] derives an expression for 〈
χ
∣
η
〉 by doubly integrating the mixture fraction PDF transport equation. However, in his derivation, the inhomogeneous terms are retained. The CV fluctuations appearing in the transport equation are modelled using the PDF gradient model of Pope [19] and the PDF is presumed using a functional form described by the vector of mixture fraction moments. When the β-PDF is employed, this model reads as follows [23]:(17)
〈
χ
∣
η
〉
=
2
P
~
(
η
)
{
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
+
D
t
[
∂
2
I
I
(
η
)
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
+
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
~
×
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
}
,
where
(18)
I
I
(
η
;
v
,
w
)
=
∫
0
η
I
(
η
′
;
v
,
w
)
d
η
′
,
(19)
I
(
η
;
v
,
w
)
=
∫
0
η
P
~
(
η
;
v
,
w
)
d
η
=
∫
0
η
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η
B
(
v
,
w
)
is the incomplete β-function. By applying some identities and integral properties of the incomplete β-function, Mortensen [23] shows that I
I
(
η
;
v
,
w
) simplifies to
(20)
I
I
(
η
;
v
,
w
)
=
(
η
-
ξ
~
)
I
(
η
;
v
,
w
)
+
ξ
′′
2
~
P
~
(
η
;
v
+
1
,
w
+
1
)
.
The term S
ξ
′′
2
~ in (17) represents the source term of the mixture fraction variance transport equation (discussed in more detail in Section 3.3). It is given by
(21)
S
ξ
′′
2
~
=
2
D
t
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
-
χ
~
.
Mortensen’s model ensures a fully consistent CMC implementation as it accounts for the inhomogeneous terms of the PDF transport equation and employs a consistent closure for the CV fluctuations. As such, this model does not yield any spurious sources. It is important to note that Girimaji’s model is the exact equivalent of the homogeneous portion of Mortensen’s, 〈
χ
∣
η
〉
=
2
[
χ
~
/
P
~
(
η
)
]
∂
I
I
(
η
)
/
∂
ξ
′′
2
~.


#### 3.2.5.  Conditional Chemical Source

The conditional chemical source term,〈
ω
˙
κ
∣
η
〉, is modelled using a first-order closure. In this closure, it is assumed that the conditional fluctuations about the conditional averages of the reactive scalars are small. Accordingly, 〈
ω
˙
κ
∣
η
〉 is modelled as a function of the conditional density, mass fractions, and temperature:
(22)
〈
ω
˙
κ
(
ρ
,
T
,
Y
)
∣
η
〉
≈
〈
ω
˙
κ
(
〈
ρ
∣
η
〉
,
〈
T
∣
η
〉
,
〈
Y
∣
η
〉
)
∣
η
〉
=
ω
˙
κ
(
〈
ρ
∣
η
〉
,
Q
T
,
Q
)
,
where Y
=
{
Y
κ
∣
κ
=
1,2
,
…
,
N
} and Q
=
〈
Y
∣
η
〉. This closure has been successfully applied in the CMC of lifted flames [7, 8, 11, 24].


#### 3.2.6.  Conditional Radiative Source

The conditional radiative source term is modelled using the optically thin assumption withH
2
O being the predominantly participating species. To a first order approximation,
(23)
〈
ω
˙
r
∣
η
〉
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
〈
T
∣
η
〉
4
-
T
b
4
)
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
Q
T
4
-
T
b
4
)
,
where σ
=
5.669
×
1
0
8 W/m2K4 is the Stefan-Boltzmann constant, p
H
2
O and a
p
,
H
2
O are the partial pressure, and the Planck mean absorption coefficient of H
2
O, respectively, and T
b is the background temperature (300 K). A Curve fit for a
p
,
H
2
O is obtained from [25].


### 3.3.  Turbulent Flow Field Calculations

The Favre averages of the velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are obtained from an inert flow calculation using the commercial Computational Fluid Dynamics (CFD) software FLUENT [26]. This approach was successfully applied in the CMC calculations of lifted flames performed by Kim and Mastorakos [24]. They justify the validity of this simplification by the fact that the flow remains frozen until the stabilisation height. Nevertheless, they note that density changes caused by the flame are non-negligible and may affect the flow field.A two-dimensional axisymmetric computational domain is employed to perform the inert flow calculations. The domain extends 55 nozzle diameters above the nozzle exit in the axial direction and 15 diameters in the radial direction with the origin of the domain placed at the centre of the nozzle. The thickness of the wall of the nozzle is neglected [3, 27]. A quadratic unstructured mesh is generated using ANSYS ICEM CFD [28]. It consists of 1,000
×
265 nodes (axial × radial). The k
-
ε model is employed to perform the calculations. Default model constants are employed, except for C
ɛ
1 which is modified from the standard value of 1.44 to 1.6 in order to improve the prediction of the spreading rate of the jet [29]. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm is used for pressure-velocity coupling. The PREssure STaggering Option (PRESTO) pressure interpolation scheme is employed to compute the pressure at cell faces. Spatial discretization is performed using the second-order upwind scheme. The inlet velocity boundary condition is specified using the 1/7 power law with the centreline velocity set to 5/4 times the experimentally estimated mean value of 107 m/s as in [27]. The coflow velocity is assumed to be uniform and is set to 3.5 m/s. The turbulence intensity is taken to be 5% at both the inlet and the coflow. The boundary conditions of the temperature and species mass fractions are specified following Table 1. Transport equations for the mean and variance of the mixture fraction are added to the solver. These are given by(24)
ρ
¯
∂
ξ
~
∂
t
+
ρ
¯
u
~
i
∂
ξ
~
∂
x
i
=
-
∂
(
ρ
¯
u
j
′′
ξ
′′
~
)
∂
x
i
,
(25)
ρ
¯
∂
ξ
′′
2
~
∂
t
+
ρ
¯
u
~
i
∂
ξ
′′
2
~
∂
x
i
=
-
∂
(
ρ
¯
u
j
′′
ξ
′′
2
~
)
∂
x
i
-
2
ρ
¯
u
i
′′
ξ
′′
~
∂
ξ
~
∂
x
i
-
ρ
¯
χ
~
.
The gradient diffusion assumption is employed to close the turbulent fluxes as u
i
′′
ξ
′′
~
=
-
D
t
∂
ξ
~
/
∂
x
i and u
i
′′
ξ
′′
2
~
=
-
D
t
∂
ξ
′′
2
~
/
∂
x
i. The last term on the r.h.s. of (25), χ
~, is the Favre-averaged scalar dissipation rate. It is a standard practice to model this term as
(26)
χ
~
=
C
χ
ɛ
~
k
~
ξ
′′
2
~
,
where C
χ is a constant set to 2. In this algebraic model, it is assumed that the flow time (k
~
/
ɛ
~) is proportional to the time scale of scalar turbulence (ξ
′′
2
~
/
χ
~), with C
χ being the proportionality constant (the time scale ratio). This modelling approach was adopted in the early stages of this study. However, it was abandoned later because the calculation yielded very short liftoff heights for moderate and low T
c and nearly attached flames for sufficiently high T
c. This behavior was not observed in the lifted flames simulations of Patwardhan et al. [11] and Kim and Mastorakos [24]. Alternatively, the transport equation derived for χ
~ by Jones and Musonge [30] is solved. This equation reads
(27)
ρ
¯
∂
χ
~
∂
t
+
ρ
¯
u
~
i
∂
χ
~
∂
x
i
=
∂
∂
x
j
(
C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
ɛ
~
∂
χ
~
∂
x
l
)
-
C
1
ρ
¯
χ
~
2
2
ξ
′′
2
~
-
C
2
ρ
¯
ɛ
~
2
k
~
χ
~
-
2
C
3
ρ
¯
ɛ
~
k
~
u
i
′′
ξ
′′
~
∂
ξ
~
∂
x
i
-
C
4
ρ
¯
χ
~
k
~
u
i
′′
u
j
′′
~
∂
u
~
i
∂
x
j
,
where C
s
=
0.22, C
1
=
2, C
2
=
2
(
C
ɛ
2
-
1
)
=
1.84, C
3
=
1.70, and C
4
=
C
ɛ
1
=
1.6 are model constants. The turbulent flux u
i
′′
ξ
′′
~ is modelled using the gradient diffusion assumption. For consistency with the Boussinesq hypothesis, the production of turbulence kinetic energy appearing in the last term on the r.h.s of (27) is modelled as
(28)
-
ρ
¯
u
i
′′
u
j
′′
~
∂
u
~
i
∂
x
j
=
μ
t
S
2
,
where μ
t is the turbulent viscosity and S is the modulus of the mean strain rate tensor. The tensor C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
/
ɛ
~ appearing in the first term on the r.h.s. of (27) represents the anisotropic diffusivity. Initial trials showed numerical instabilities and suffered from severe convergence issues when accounting for anisotropic effects. Therefore, for simplicity, the diffusivity was assumed to be isotropic and standard modelling was employed (μ
t
/
S
c
t used instead of C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
/
ɛ
~). The inclusion of the transport equation of χ
~ (see (27)) in the context of CMC was first employed in the counterflow flames calculations of Kim and Mastorakos [31]. The authors solve an additional transport equation for u
i
′′
ξ
′′
~ rather than resorting to gradient diffusion modelling.


### 3.4.  The CMC Implementation

A smaller physical domain is employed in the CMC calculations for computational efficiency. The chosen domain is 30 nozzle diameters in length and 5 diameters in width. The mesh consists of70
×
40 (axial × radial) nonuniformly distributed nodes, with the mesh density being highest near the experimentally measured stabilization height (10
d) [3]. The mixture fraction grid consists of 80 points. These resolutions ensure grid independence. The CMC equations are discretized using the finite difference method. The first-order derivative appearing in the convective term is discretized using the second-order upwind difference scheme with the kappa flux limiter of Koren [32], following Patwardhan et al. [11]. The second-order derivatives in physical and mixture fraction spaces are discretized using the second-order central difference scheme. The steady-state solution is obtained from the transient CMC transport equations, (8) and (9), by relaxation of the time step, for which a value of 10−5 s is used. A three-step fractional method (Strang operator splitting [33]) is implemented in order to treat the stiff chemical source term separately. The first fractional step accounts for convection and diffusion in physical space over the first half of the time step. The second step handles the chemical source and micromixing in mixture fractions space over the whole time step. The third and last step accounts for convection and diffusion in physical space over the second half of the time step. Each step uses the solution of the previous one as initial conditions. Splitting errors are negligible due to the usage of a small time step. The Ordinary Differential Equation (ODE) solver VODPK [34–36] is used to solve the system of equations. The ODEs in the first and third steps are nonstiff and hence solved with the implicit Adams method. Conversely, the ODEs in the second step are stiff and therefore treated with the Backward Differentiation Formulas (BDF). The boundary conditions for Q
κ and Q
T at the inlet and the coflow are specified using the inert mixing solution. The composition of the fuel and coflow streams is obtained from the experiments [3]. Zero-gradient boundary conditions are specified at the remaining boundaries. The solution is initialized by igniting some nodes above the experimentally measured stabilization height with a steady laminar flamelet computed with a strain rate of 1000 s−1. The Inert mixing solution is set elsewhere. The Favre-averaged velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are transferred from FLUENT to the CMC solver by applying bilinear interpolation. This interpolation scheme is accurate enough given the high resolution of the CFD grid. The integrals appearing in Girimaji’s model are computed using the quadrature integration package QUADPACK [37]. In Mortensen’s model, the incomplete beta-function I
(
η
;
v
,
w
) (see (19)) is calculated using the method of continued fractions [38] and the partial derivatives of I
I
(
η
;
v
,
w
) (see (20)) with respect to the mean and variance of the mixture fraction are computed using Ridders’ method of polynomial extrapolations [38, 39]. Due to the division by the PDF in the PDF-gradient and Mortensen’s models, numerical instabilities may arise at low probabilities. To resolve this issue, 〈
u
i
∣
η
〉 is set equal to its unconditional value, u
~
i, and 〈
χ
∣
η
〉 is set to zero when η falls outside ξ
~
±
7
(
ξ
′′
2
~
)
1
/
2. Additionally, in the event where Moretnsen’s model yields unphysical negative values, 〈
χ
∣
η
〉 is set to zero.The chemical kinetics mechanisms developed by Mueller et al. [40], Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43] are considered in this study. The corresponding numbers of species and reactions are summarised in Table 2. All four mechanisms have a comparable level of detail, with those of Mueller and Burke being the least and most detailed, respectively. The Mueller mechanism is the most extensively used throughout this work. It is provided in Table 3.Table 2
Number of species and reactions in the considered chemical kinetics mechanisms.
Mechanism
Species
Reactions
Reference
Mueller
9
21
[40]
Li
13
25
[41]
Burke
13
27
[42]
Ó Conaire
10
21
[43]Table 3
Chemical kinetics mechanism of Mueller et al. [40].
R1
H + O2
⇌  O + OH
R2
O + H2
⇌  H + OH
R3
H2 + OH  ⇌  H2O + H
R4
O + H2O  ⇌  OH + OH
R5
H2 + M  ⇌  H + H + M
R6
O + O + M⇌  O2 + M
R7
O + H + M⇌  OH + M
R8
H + OH + M⇌  H2O + M
R9
H + O2 (+M) ⇌  HO2 (+M)
R10
HO2 + H  ⇌  H2 + O2
R11
HO2 + H  ⇌  OH + OH
R12
HO2 + O  ⇌  O2+ OH
R13
HO2 + OH  ⇌  H2O + O2
R14
HO2 + HO2
⇌  H2O2 + O
2
†
R15
H2O2 (+M) ⇌  OH + OH (+M)
R16
H2O2 + H  ⇌  H2O + OH
R17
H2O2 + H  ⇌  HO2 + H2
R18
H2O2 + O  ⇌  OH + HO2
R19
H2O2 + OH  ⇌  HO2+ H2O†
†R14 and R19 are expressed as the sum of two rate expressions. Hence, the mechanism has a total of 21 reactions.


## 3.1. Conditional Moment Closure

The modelling of the turbulence-chemistry interactions is carried out using the first-order CMC [14]. In CMC, reactive scalars such as the species mass fractions and the temperature are conditionally averaged with respect to the mixture fraction, ξ, and their conditional transport equations are solved. The conditional average of the mass fraction of a species κ is defined as
(1)
Q
κ
(
η
,
x
,
t
)
=
〈
Y
κ
(
x
,
t
)
∣
ξ
(
x
,
t
)
=
η
〉
=
〈
Y
κ
∣
η
〉
,
where Y
κ is the mass fraction of the species, 〈
·
∣
·
〉 denotes the conditional average of the quantity to the left of the vertical bar subject to the quantity to its right, and η is a sample variable of ξ, with 0
≤
η
≤
1. Following the decomposition approach [14], Y
κ is written as the sum of its conditional average Q
κ and a fluctuation y
κ
′′:
(2)
Y
κ
(
x
,
t
)
=
Q
κ
(
η
,
x
,
t
)
+
y
κ
′′
(
η
,
x
,
t
)
such that 〈
y
κ
′′
∣
η
〉
=
0. Substitution of (2) into the transport equation of species κ, followed by the conditional averaging of the resulting expression, leads to [14]:
(3)
∂
Q
κ
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
κ
∂
x
i
+
〈
χ
∣
η
〉
2
∂
2
Q
κ
∂
η
2
+
〈
ω
˙
κ
∣
η
〉
〈
ρ
∣
η
〉
+
e
Q
κ
+
e
y
κ
′′
〈
ρ
∣
η
〉
,
where
(4)
e
Q
κ
=
∂
∂
x
i
(
ρ
D
∂
Q
κ
∂
x
i
)
+
〈
ρ
D
∂
ξ
∂
x
i
∂
∂
ξ
(
∂
Q
κ
∂
x
i
)
∣
η
〉
,
(5)
e
y
κ
′′
=
-
〈
ρ
∂
y
κ
′′
∂
t
+
ρ
u
i
′′
∂
y
κ
′′
∂
x
i
-
∂
∂
x
i
(
ρ
D
∂
y
κ
′′
∂
x
i
)
∣
η
〉
,
(6)
χ
=
2
D
∂
ξ
∂
x
i
∂
ξ
∂
x
i
is the scalar dissipation rate. The species and mixture fraction diffusivities in (3)–(6) are assumed to be equal to the molecular diffusivity, D, and the Lewis number is set to unity. In (3), 〈
u
i
∣
η
〉 is the CV, 〈
χ
∣
η
〉 is the CSDR, and 〈
ω
˙
κ
∣
η
〉 is the conditional chemical source. The expressions of e
Q (see (4)) and e
y
κ
′′ (see (5)) are unclosed. They are modelled using the primary closure hypothesis. Klimenko and Bilger [14] show that, given finite Schmidt numbers, all the terms of e
Q scale as the inverse of the Reynolds number (Re). Therefore, this term is neglected in high Re applications. By neglecting the conditional density and diffusivity fluctuations, they further show that and e
y
κ
′′ reduces to
(7)
e
y
κ
′′
=
-
〈
ρ
∣
η
〉
ρ
¯
P
~
(
η
)
∂
∂
x
i
[
ρ
¯
〈
u
i
′′
y
κ
′′
∣
η
〉
P
~
(
η
)
]
,
where P
~
(
η
) is the Favre PDF of mixture fraction and 〈
u
i
′′
y
κ
′′
∣
η
〉 is the conditional turbulent flux. Using these approximations, the conditional species transport equation takes the form:(8)
∂
Q
κ
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
κ
∂
x
i
-
1
ρ
¯
P
~
(
η
)
∂
∂
x
i
[
ρ
¯
〈
u
i
′′
y
κ
′′
∣
η
〉
P
~
(
η
)
]
+
〈
χ
∣
η
〉
2
∂
2
Q
κ
∂
η
2
+
〈
ω
˙
κ
∣
η
〉
〈
ρ
η
〉
.
The left-hand side term of (8) is the rate of change of the conditional species mass fraction. The first term on the right-hand side (r.h.s.) represents transport by means of convection, the second corresponds to transport by conditional turbulent fluxes (diffusion in physical space), the third accounts for micromixing (diffusion in mixture fraction space), and the fourth is the chemical source. The conditional temperature equation is derived in a similar fashion and is given by
(9)
∂
Q
T
∂
t
=
-
〈
u
i
∣
η
〉
∂
Q
T
∂
x
i
︸
T
C
,
x
and
T
C
,
y
-
1
ρ
-
P
~
(
η
)
∂
∂
x
i
[
ρ
-
〈
u
i
′
′
T
′′
∣
η
〉
P
~
(
η
)
]
︸
T
D
,
x
and
T
D
,
y
+
〈
χ
∣
η
〉
2
{
∂
2
Q
T
∂
η
2
+
1
〈
c
p
∣
η
〉
×
[
∂
〈
c
p
∣
η
〉
∂
η
+
∑
κ
=
1
N
(
〈
c
p
,
κ
∣
η
〉
∂
Q
κ
∂
η
)
]
×
∂
Q
T
∂
η
}
T
MM
+
1
〈
c
p
∣
η
〉
〈
1
ρ
∂
p
∂
t
∣
η
〉
-
〈
ω
˙
h
∣
η
〉
〈
ρ
∣
η
〉
〈
c
p
∣
η
〉
︸
T
CS
-
〈
ω
˙
r
∣
η
〉
〈
ρ
∣
η
〉
〈
c
p
∣
η
〉
︸
T
RS
,
where 〈
ω
˙
h
∣
η
〉
=
∑
κ
=
1
N
〈
h
κ
∣
η
〉
〈
ω
˙
κ
∣
η
〉 is the chemical source term, N is the total number of species in the mixture, and 〈
ω
˙
r
∣
η
〉 is the conditional radiative source term. The unsteady pressure term in (9) is neglected due to the fact that the investigated flame is open to the atmosphere, and hence temporal pressure changes are expected very small.The unconditional (Favre) averages of the reactive scalars are calculated by integrating their conditional values weighted byP
~
(
η
) over the mixture fraction space:
(10)
ϕ
~
=
∫
0
1
〈
ϕ
∣
η
〉
P
~
(
η
)
d
η
,
ϕ
=
Y
κ
,
T
.


## 3.2. Conditional Submodels

The quantitiesP
~
(
η
), 〈
u
i
′′
y
κ
′′
∣
η
〉, 〈
u
i
′′
T
′′
∣
η
〉, 〈
u
i
∣
η
〉, 〈
χ
∣
η
〉, 〈
ω
˙
κ
∣
η
〉, and 〈
ω
˙
r
∣
η
〉 in (8) and (9) are unclosed and require further modelling. The submodels employed in this study are discussed in this section.


### 3.2.1. Probability Density Function

The PDF is presumed using theβ-distribution. The Favre β-PDF is given by
(11)
P
~
(
η
;
v
,
w
)
=
η
v
-
1
(
1
-
η
)
w
-
1
B
(
v
,
w
)
,
where the parameters v and w are related to the mean and variance of the mixture fraction by v
=
γ
ξ
~ and w
=
γ
(
1
-
ξ
~
) with γ
=
[
ξ
~
(
1
-
ξ
~
)
/
ξ
′′
2
~
-
1
]
≥
0, and B
(
v
,
w
)
=
∫
0
1
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η is the beta function.


### 3.2.2.  Conditional Turbulent Fluxes

The conditional turbulent fluxes are modelled using the gradient diffusion assumption:(12)
〈
ϕ
′
′
u
i
′′
∣
η
〉
=
-
D
t
∂
〈
ϕ
∣
η
〉
∂
x
i
,
ϕ
=
Y
κ
,
T
,
where
(13)
D
t
=
C
μ
S
c
t
k
~
2
ɛ
~
is the turbulent diffusivity. The constant C
μ is equal to 0.09 and the turbulent Schmidt number S
c
t is set to 0.7 following Jones and Whitelaw [17]. Richardson and Mastorakos [18] suggest adding a correction to (12) in order to account for counter-gradient transport effects. When applied to a lifted flame, they conclude that the inclusion of counter-gradient effects leads to a slight increase in liftoff height and to a decrease in flame thickness. This extension is not employed in the current study. However, it is worth investigating in future studies.


### 3.2.3.  Conditional Velocity

The PDF-gradient model proposed by Pope [19] is employed to model the CV fluctuations. The resulting expression for the CV takes the form:
(14)
〈
u
i
∣
η
〉
=
u
~
-
D
t
P
~
(
η
)
∂
P
~
(
η
)
∂
x
i
=
u
~
i
-
D
t
{
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
~
∂
ξ
~
∂
x
i
+
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
}
.
As P
~
(
η
) tends to zero, 〈
u
i
∣
η
〉 diverges to ±
∞, depending on the sing of the gradient of P
~
(
η
). This behavior is of minor importance, because low-probability events have a negligible effect on mixing [20]. One important feature of the PDF-gradient model is its consistency with both the first and second moments of the mixture fraction [15, 20]. This does not hold for the commonly used linear model [12], which is only consistent with the first moment [15].


### 3.2.4.  Conditional Scalar Dissipation Rate

Girimaji’s Model. Girimaji [21] derives a model for 〈
χ
∣
η
〉 by doubly integrating the homogeneous mixture fraction PDF transport equation. His formulation is based on the observation that a presumed β-PDF is capable of accurately characterizing the evolution of the scalar PDF over all stages of two-scalar, constant-density mixing in statistically stationary, isotropic turbulence. The model is given by:
(15)
〈
χ
η
〉
=
-
2
χ
~
ξ
~
(
1
-
ξ
~
)
(
ξ
′′
2
~
)
2
I
(
η
)
P
~
(
η
)
,
where χ
~ is the Favre-averaged scalar dissipation rate (discussed in more detail in Section 3.3) and I
(
η
) is given by the integral:
(16)
I
(
η
)
=
∫
0
η
{
ξ
~
[
ln
⁡
η
′
-
∫
0
1
ln
⁡
η
′′
P
~
(
η
′′
)
d
η
′′
]
+
(
1
-
ξ
~
)
×
[
ln
⁡
(
1
-
η
′
)
-
∫
0
1
ln
⁡
⁡
(
1
-
η
′′
)
P
~
(
η
′′
)
d
η
′′
]
}
×
P
~
(
η
′
)
(
η
-
η
′
)
d
η
′
.
As mentioned above, Girimaji’s model is valid for homogeneous turbulence. Therefore, its usage in the CMC of inhomogeneous flows is inconsistent and results in spurious source terms.Mortensen’s Model. In a similar fashion, Mortensen [22] derives an expression for 〈
χ
∣
η
〉 by doubly integrating the mixture fraction PDF transport equation. However, in his derivation, the inhomogeneous terms are retained. The CV fluctuations appearing in the transport equation are modelled using the PDF gradient model of Pope [19] and the PDF is presumed using a functional form described by the vector of mixture fraction moments. When the β-PDF is employed, this model reads as follows [23]:(17)
〈
χ
∣
η
〉
=
2
P
~
(
η
)
{
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
+
D
t
[
∂
2
I
I
(
η
)
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
+
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
~
×
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
}
,
where
(18)
I
I
(
η
;
v
,
w
)
=
∫
0
η
I
(
η
′
;
v
,
w
)
d
η
′
,
(19)
I
(
η
;
v
,
w
)
=
∫
0
η
P
~
(
η
;
v
,
w
)
d
η
=
∫
0
η
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η
B
(
v
,
w
)
is the incomplete β-function. By applying some identities and integral properties of the incomplete β-function, Mortensen [23] shows that I
I
(
η
;
v
,
w
) simplifies to
(20)
I
I
(
η
;
v
,
w
)
=
(
η
-
ξ
~
)
I
(
η
;
v
,
w
)
+
ξ
′′
2
~
P
~
(
η
;
v
+
1
,
w
+
1
)
.
The term S
ξ
′′
2
~ in (17) represents the source term of the mixture fraction variance transport equation (discussed in more detail in Section 3.3). It is given by
(21)
S
ξ
′′
2
~
=
2
D
t
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
-
χ
~
.
Mortensen’s model ensures a fully consistent CMC implementation as it accounts for the inhomogeneous terms of the PDF transport equation and employs a consistent closure for the CV fluctuations. As such, this model does not yield any spurious sources. It is important to note that Girimaji’s model is the exact equivalent of the homogeneous portion of Mortensen’s, 〈
χ
∣
η
〉
=
2
[
χ
~
/
P
~
(
η
)
]
∂
I
I
(
η
)
/
∂
ξ
′′
2
~.


### 3.2.5.  Conditional Chemical Source

The conditional chemical source term,〈
ω
˙
κ
∣
η
〉, is modelled using a first-order closure. In this closure, it is assumed that the conditional fluctuations about the conditional averages of the reactive scalars are small. Accordingly, 〈
ω
˙
κ
∣
η
〉 is modelled as a function of the conditional density, mass fractions, and temperature:
(22)
〈
ω
˙
κ
(
ρ
,
T
,
Y
)
∣
η
〉
≈
〈
ω
˙
κ
(
〈
ρ
∣
η
〉
,
〈
T
∣
η
〉
,
〈
Y
∣
η
〉
)
∣
η
〉
=
ω
˙
κ
(
〈
ρ
∣
η
〉
,
Q
T
,
Q
)
,
where Y
=
{
Y
κ
∣
κ
=
1,2
,
…
,
N
} and Q
=
〈
Y
∣
η
〉. This closure has been successfully applied in the CMC of lifted flames [7, 8, 11, 24].


### 3.2.6.  Conditional Radiative Source

The conditional radiative source term is modelled using the optically thin assumption withH
2
O being the predominantly participating species. To a first order approximation,
(23)
〈
ω
˙
r
∣
η
〉
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
〈
T
∣
η
〉
4
-
T
b
4
)
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
Q
T
4
-
T
b
4
)
,
where σ
=
5.669
×
1
0
8 W/m2K4 is the Stefan-Boltzmann constant, p
H
2
O and a
p
,
H
2
O are the partial pressure, and the Planck mean absorption coefficient of H
2
O, respectively, and T
b is the background temperature (300 K). A Curve fit for a
p
,
H
2
O is obtained from [25].


## 3.2.1. Probability Density Function

The PDF is presumed using theβ-distribution. The Favre β-PDF is given by
(11)
P
~
(
η
;
v
,
w
)
=
η
v
-
1
(
1
-
η
)
w
-
1
B
(
v
,
w
)
,
where the parameters v and w are related to the mean and variance of the mixture fraction by v
=
γ
ξ
~ and w
=
γ
(
1
-
ξ
~
) with γ
=
[
ξ
~
(
1
-
ξ
~
)
/
ξ
′′
2
~
-
1
]
≥
0, and B
(
v
,
w
)
=
∫
0
1
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η is the beta function.


## 3.2.2.  Conditional Turbulent Fluxes

The conditional turbulent fluxes are modelled using the gradient diffusion assumption:(12)
〈
ϕ
′
′
u
i
′′
∣
η
〉
=
-
D
t
∂
〈
ϕ
∣
η
〉
∂
x
i
,
ϕ
=
Y
κ
,
T
,
where
(13)
D
t
=
C
μ
S
c
t
k
~
2
ɛ
~
is the turbulent diffusivity. The constant C
μ is equal to 0.09 and the turbulent Schmidt number S
c
t is set to 0.7 following Jones and Whitelaw [17]. Richardson and Mastorakos [18] suggest adding a correction to (12) in order to account for counter-gradient transport effects. When applied to a lifted flame, they conclude that the inclusion of counter-gradient effects leads to a slight increase in liftoff height and to a decrease in flame thickness. This extension is not employed in the current study. However, it is worth investigating in future studies.


## 3.2.3.  Conditional Velocity

The PDF-gradient model proposed by Pope [19] is employed to model the CV fluctuations. The resulting expression for the CV takes the form:
(14)
〈
u
i
∣
η
〉
=
u
~
-
D
t
P
~
(
η
)
∂
P
~
(
η
)
∂
x
i
=
u
~
i
-
D
t
{
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
~
∂
ξ
~
∂
x
i
+
∂
ln
⁡
[
P
~
(
η
)
]
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
}
.
As P
~
(
η
) tends to zero, 〈
u
i
∣
η
〉 diverges to ±
∞, depending on the sing of the gradient of P
~
(
η
). This behavior is of minor importance, because low-probability events have a negligible effect on mixing [20]. One important feature of the PDF-gradient model is its consistency with both the first and second moments of the mixture fraction [15, 20]. This does not hold for the commonly used linear model [12], which is only consistent with the first moment [15].


## 3.2.4.  Conditional Scalar Dissipation Rate

Girimaji’s Model. Girimaji [21] derives a model for 〈
χ
∣
η
〉 by doubly integrating the homogeneous mixture fraction PDF transport equation. His formulation is based on the observation that a presumed β-PDF is capable of accurately characterizing the evolution of the scalar PDF over all stages of two-scalar, constant-density mixing in statistically stationary, isotropic turbulence. The model is given by:
(15)
〈
χ
η
〉
=
-
2
χ
~
ξ
~
(
1
-
ξ
~
)
(
ξ
′′
2
~
)
2
I
(
η
)
P
~
(
η
)
,
where χ
~ is the Favre-averaged scalar dissipation rate (discussed in more detail in Section 3.3) and I
(
η
) is given by the integral:
(16)
I
(
η
)
=
∫
0
η
{
ξ
~
[
ln
⁡
η
′
-
∫
0
1
ln
⁡
η
′′
P
~
(
η
′′
)
d
η
′′
]
+
(
1
-
ξ
~
)
×
[
ln
⁡
(
1
-
η
′
)
-
∫
0
1
ln
⁡
⁡
(
1
-
η
′′
)
P
~
(
η
′′
)
d
η
′′
]
}
×
P
~
(
η
′
)
(
η
-
η
′
)
d
η
′
.
As mentioned above, Girimaji’s model is valid for homogeneous turbulence. Therefore, its usage in the CMC of inhomogeneous flows is inconsistent and results in spurious source terms.Mortensen’s Model. In a similar fashion, Mortensen [22] derives an expression for 〈
χ
∣
η
〉 by doubly integrating the mixture fraction PDF transport equation. However, in his derivation, the inhomogeneous terms are retained. The CV fluctuations appearing in the transport equation are modelled using the PDF gradient model of Pope [19] and the PDF is presumed using a functional form described by the vector of mixture fraction moments. When the β-PDF is employed, this model reads as follows [23]:(17)
〈
χ
∣
η
〉
=
2
P
~
(
η
)
{
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
+
D
t
[
∂
2
I
I
(
η
)
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
ξ
′′
2
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
+
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
~
×
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
+
2
∂
2
I
I
(
η
)
∂
ξ
~
∂
ξ
′′
2
~
∂
ξ
~
∂
x
i
∂
ξ
′′
2
~
∂
x
i
]
-
∂
I
I
(
η
)
∂
ξ
′′
2
~
S
ξ
′′
2
~
}
,
where
(18)
I
I
(
η
;
v
,
w
)
=
∫
0
η
I
(
η
′
;
v
,
w
)
d
η
′
,
(19)
I
(
η
;
v
,
w
)
=
∫
0
η
P
~
(
η
;
v
,
w
)
d
η
=
∫
0
η
‍
η
v
-
1
(
1
-
η
)
w
-
1
d
η
B
(
v
,
w
)
is the incomplete β-function. By applying some identities and integral properties of the incomplete β-function, Mortensen [23] shows that I
I
(
η
;
v
,
w
) simplifies to
(20)
I
I
(
η
;
v
,
w
)
=
(
η
-
ξ
~
)
I
(
η
;
v
,
w
)
+
ξ
′′
2
~
P
~
(
η
;
v
+
1
,
w
+
1
)
.
The term S
ξ
′′
2
~ in (17) represents the source term of the mixture fraction variance transport equation (discussed in more detail in Section 3.3). It is given by
(21)
S
ξ
′′
2
~
=
2
D
t
∂
ξ
~
∂
x
i
∂
ξ
~
∂
x
i
-
χ
~
.
Mortensen’s model ensures a fully consistent CMC implementation as it accounts for the inhomogeneous terms of the PDF transport equation and employs a consistent closure for the CV fluctuations. As such, this model does not yield any spurious sources. It is important to note that Girimaji’s model is the exact equivalent of the homogeneous portion of Mortensen’s, 〈
χ
∣
η
〉
=
2
[
χ
~
/
P
~
(
η
)
]
∂
I
I
(
η
)
/
∂
ξ
′′
2
~.


## 3.2.5.  Conditional Chemical Source

The conditional chemical source term,〈
ω
˙
κ
∣
η
〉, is modelled using a first-order closure. In this closure, it is assumed that the conditional fluctuations about the conditional averages of the reactive scalars are small. Accordingly, 〈
ω
˙
κ
∣
η
〉 is modelled as a function of the conditional density, mass fractions, and temperature:
(22)
〈
ω
˙
κ
(
ρ
,
T
,
Y
)
∣
η
〉
≈
〈
ω
˙
κ
(
〈
ρ
∣
η
〉
,
〈
T
∣
η
〉
,
〈
Y
∣
η
〉
)
∣
η
〉
=
ω
˙
κ
(
〈
ρ
∣
η
〉
,
Q
T
,
Q
)
,
where Y
=
{
Y
κ
∣
κ
=
1,2
,
…
,
N
} and Q
=
〈
Y
∣
η
〉. This closure has been successfully applied in the CMC of lifted flames [7, 8, 11, 24].


## 3.2.6.  Conditional Radiative Source

The conditional radiative source term is modelled using the optically thin assumption withH
2
O being the predominantly participating species. To a first order approximation,
(23)
〈
ω
˙
r
∣
η
〉
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
〈
T
∣
η
〉
4
-
T
b
4
)
=
4
σ
p
H
2
O
a
p
,
H
2
O
(
Q
T
4
-
T
b
4
)
,
where σ
=
5.669
×
1
0
8 W/m2K4 is the Stefan-Boltzmann constant, p
H
2
O and a
p
,
H
2
O are the partial pressure, and the Planck mean absorption coefficient of H
2
O, respectively, and T
b is the background temperature (300 K). A Curve fit for a
p
,
H
2
O is obtained from [25].


## 3.3.  Turbulent Flow Field Calculations

The Favre averages of the velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are obtained from an inert flow calculation using the commercial Computational Fluid Dynamics (CFD) software FLUENT [26]. This approach was successfully applied in the CMC calculations of lifted flames performed by Kim and Mastorakos [24]. They justify the validity of this simplification by the fact that the flow remains frozen until the stabilisation height. Nevertheless, they note that density changes caused by the flame are non-negligible and may affect the flow field.A two-dimensional axisymmetric computational domain is employed to perform the inert flow calculations. The domain extends 55 nozzle diameters above the nozzle exit in the axial direction and 15 diameters in the radial direction with the origin of the domain placed at the centre of the nozzle. The thickness of the wall of the nozzle is neglected [3, 27]. A quadratic unstructured mesh is generated using ANSYS ICEM CFD [28]. It consists of 1,000
×
265 nodes (axial × radial). The k
-
ε model is employed to perform the calculations. Default model constants are employed, except for C
ɛ
1 which is modified from the standard value of 1.44 to 1.6 in order to improve the prediction of the spreading rate of the jet [29]. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm is used for pressure-velocity coupling. The PREssure STaggering Option (PRESTO) pressure interpolation scheme is employed to compute the pressure at cell faces. Spatial discretization is performed using the second-order upwind scheme. The inlet velocity boundary condition is specified using the 1/7 power law with the centreline velocity set to 5/4 times the experimentally estimated mean value of 107 m/s as in [27]. The coflow velocity is assumed to be uniform and is set to 3.5 m/s. The turbulence intensity is taken to be 5% at both the inlet and the coflow. The boundary conditions of the temperature and species mass fractions are specified following Table 1. Transport equations for the mean and variance of the mixture fraction are added to the solver. These are given by(24)
ρ
¯
∂
ξ
~
∂
t
+
ρ
¯
u
~
i
∂
ξ
~
∂
x
i
=
-
∂
(
ρ
¯
u
j
′′
ξ
′′
~
)
∂
x
i
,
(25)
ρ
¯
∂
ξ
′′
2
~
∂
t
+
ρ
¯
u
~
i
∂
ξ
′′
2
~
∂
x
i
=
-
∂
(
ρ
¯
u
j
′′
ξ
′′
2
~
)
∂
x
i
-
2
ρ
¯
u
i
′′
ξ
′′
~
∂
ξ
~
∂
x
i
-
ρ
¯
χ
~
.
The gradient diffusion assumption is employed to close the turbulent fluxes as u
i
′′
ξ
′′
~
=
-
D
t
∂
ξ
~
/
∂
x
i and u
i
′′
ξ
′′
2
~
=
-
D
t
∂
ξ
′′
2
~
/
∂
x
i. The last term on the r.h.s. of (25), χ
~, is the Favre-averaged scalar dissipation rate. It is a standard practice to model this term as
(26)
χ
~
=
C
χ
ɛ
~
k
~
ξ
′′
2
~
,
where C
χ is a constant set to 2. In this algebraic model, it is assumed that the flow time (k
~
/
ɛ
~) is proportional to the time scale of scalar turbulence (ξ
′′
2
~
/
χ
~), with C
χ being the proportionality constant (the time scale ratio). This modelling approach was adopted in the early stages of this study. However, it was abandoned later because the calculation yielded very short liftoff heights for moderate and low T
c and nearly attached flames for sufficiently high T
c. This behavior was not observed in the lifted flames simulations of Patwardhan et al. [11] and Kim and Mastorakos [24]. Alternatively, the transport equation derived for χ
~ by Jones and Musonge [30] is solved. This equation reads
(27)
ρ
¯
∂
χ
~
∂
t
+
ρ
¯
u
~
i
∂
χ
~
∂
x
i
=
∂
∂
x
j
(
C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
ɛ
~
∂
χ
~
∂
x
l
)
-
C
1
ρ
¯
χ
~
2
2
ξ
′′
2
~
-
C
2
ρ
¯
ɛ
~
2
k
~
χ
~
-
2
C
3
ρ
¯
ɛ
~
k
~
u
i
′′
ξ
′′
~
∂
ξ
~
∂
x
i
-
C
4
ρ
¯
χ
~
k
~
u
i
′′
u
j
′′
~
∂
u
~
i
∂
x
j
,
where C
s
=
0.22, C
1
=
2, C
2
=
2
(
C
ɛ
2
-
1
)
=
1.84, C
3
=
1.70, and C
4
=
C
ɛ
1
=
1.6 are model constants. The turbulent flux u
i
′′
ξ
′′
~ is modelled using the gradient diffusion assumption. For consistency with the Boussinesq hypothesis, the production of turbulence kinetic energy appearing in the last term on the r.h.s of (27) is modelled as
(28)
-
ρ
¯
u
i
′′
u
j
′′
~
∂
u
~
i
∂
x
j
=
μ
t
S
2
,
where μ
t is the turbulent viscosity and S is the modulus of the mean strain rate tensor. The tensor C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
/
ɛ
~ appearing in the first term on the r.h.s. of (27) represents the anisotropic diffusivity. Initial trials showed numerical instabilities and suffered from severe convergence issues when accounting for anisotropic effects. Therefore, for simplicity, the diffusivity was assumed to be isotropic and standard modelling was employed (μ
t
/
S
c
t used instead of C
s
ρ
¯
u
j
′′
u
l
′′
~
k
~
/
ɛ
~). The inclusion of the transport equation of χ
~ (see (27)) in the context of CMC was first employed in the counterflow flames calculations of Kim and Mastorakos [31]. The authors solve an additional transport equation for u
i
′′
ξ
′′
~ rather than resorting to gradient diffusion modelling.


## 3.4.  The CMC Implementation

A smaller physical domain is employed in the CMC calculations for computational efficiency. The chosen domain is 30 nozzle diameters in length and 5 diameters in width. The mesh consists of70
×
40 (axial × radial) nonuniformly distributed nodes, with the mesh density being highest near the experimentally measured stabilization height (10
d) [3]. The mixture fraction grid consists of 80 points. These resolutions ensure grid independence. The CMC equations are discretized using the finite difference method. The first-order derivative appearing in the convective term is discretized using the second-order upwind difference scheme with the kappa flux limiter of Koren [32], following Patwardhan et al. [11]. The second-order derivatives in physical and mixture fraction spaces are discretized using the second-order central difference scheme. The steady-state solution is obtained from the transient CMC transport equations, (8) and (9), by relaxation of the time step, for which a value of 10−5 s is used. A three-step fractional method (Strang operator splitting [33]) is implemented in order to treat the stiff chemical source term separately. The first fractional step accounts for convection and diffusion in physical space over the first half of the time step. The second step handles the chemical source and micromixing in mixture fractions space over the whole time step. The third and last step accounts for convection and diffusion in physical space over the second half of the time step. Each step uses the solution of the previous one as initial conditions. Splitting errors are negligible due to the usage of a small time step. The Ordinary Differential Equation (ODE) solver VODPK [34–36] is used to solve the system of equations. The ODEs in the first and third steps are nonstiff and hence solved with the implicit Adams method. Conversely, the ODEs in the second step are stiff and therefore treated with the Backward Differentiation Formulas (BDF). The boundary conditions for Q
κ and Q
T at the inlet and the coflow are specified using the inert mixing solution. The composition of the fuel and coflow streams is obtained from the experiments [3]. Zero-gradient boundary conditions are specified at the remaining boundaries. The solution is initialized by igniting some nodes above the experimentally measured stabilization height with a steady laminar flamelet computed with a strain rate of 1000 s−1. The Inert mixing solution is set elsewhere. The Favre-averaged velocity, mixture fraction mean and variance, turbulence kinetic energy, and turbulence eddy dissipation are transferred from FLUENT to the CMC solver by applying bilinear interpolation. This interpolation scheme is accurate enough given the high resolution of the CFD grid. The integrals appearing in Girimaji’s model are computed using the quadrature integration package QUADPACK [37]. In Mortensen’s model, the incomplete beta-function I
(
η
;
v
,
w
) (see (19)) is calculated using the method of continued fractions [38] and the partial derivatives of I
I
(
η
;
v
,
w
) (see (20)) with respect to the mean and variance of the mixture fraction are computed using Ridders’ method of polynomial extrapolations [38, 39]. Due to the division by the PDF in the PDF-gradient and Mortensen’s models, numerical instabilities may arise at low probabilities. To resolve this issue, 〈
u
i
∣
η
〉 is set equal to its unconditional value, u
~
i, and 〈
χ
∣
η
〉 is set to zero when η falls outside ξ
~
±
7
(
ξ
′′
2
~
)
1
/
2. Additionally, in the event where Moretnsen’s model yields unphysical negative values, 〈
χ
∣
η
〉 is set to zero.The chemical kinetics mechanisms developed by Mueller et al. [40], Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43] are considered in this study. The corresponding numbers of species and reactions are summarised in Table 2. All four mechanisms have a comparable level of detail, with those of Mueller and Burke being the least and most detailed, respectively. The Mueller mechanism is the most extensively used throughout this work. It is provided in Table 3.Table 2
Number of species and reactions in the considered chemical kinetics mechanisms.
Mechanism
Species
Reactions
Reference
Mueller
9
21
[40]
Li
13
25
[41]
Burke
13
27
[42]
Ó Conaire
10
21
[43]Table 3
Chemical kinetics mechanism of Mueller et al. [40].
R1
H + O2
⇌  O + OH
R2
O + H2
⇌  H + OH
R3
H2 + OH  ⇌  H2O + H
R4
O + H2O  ⇌  OH + OH
R5
H2 + M  ⇌  H + H + M
R6
O + O + M⇌  O2 + M
R7
O + H + M⇌  OH + M
R8
H + OH + M⇌  H2O + M
R9
H + O2 (+M) ⇌  HO2 (+M)
R10
HO2 + H  ⇌  H2 + O2
R11
HO2 + H  ⇌  OH + OH
R12
HO2 + O  ⇌  O2+ OH
R13
HO2 + OH  ⇌  H2O + O2
R14
HO2 + HO2
⇌  H2O2 + O
2
†
R15
H2O2 (+M) ⇌  OH + OH (+M)
R16
H2O2 + H  ⇌  H2O + OH
R17
H2O2 + H  ⇌  HO2 + H2
R18
H2O2 + O  ⇌  OH + HO2
R19
H2O2 + OH  ⇌  HO2+ H2O†
†R14 and R19 are expressed as the sum of two rate expressions. Hence, the mechanism has a total of 21 reactions.


## 4. Results and Discussion

### 4.1. Flow Field Results

The radial profiles of the Favre-averaged axial velocity are compared to the experimental measurements of Kent [44] in Figure 2(a). The predictions are in very good agreement with the experimental data up to x
/
d
=
11. The measurements are slightly underpredicted at x
/
d
=
14 for y
/
d
>
0.5. No velocity measurements are available at x
/
d
=
26. The profile at this location is provided for completeness. The radial profiles of the Favre-averaged mixture fraction mean are displayed in Figure 2(b) along with the experimental data of Cabra et al. [3]. The predictions are in very good agreement with the experiments at x
/
d = 1, 8, 9, 10, and 11. Slight underpredictions are obtained at x
/
d
=
14 and 26 for 0
≤
y
/
d
≤
1. Figure 2(c) shows a comparison of the radial Favre-averaged mixture fraction variance profiles with the experimental measurements. A very good agreement is obtained at all axial locations except at x
/
d
=
26. The discrepancies at this location are attributed to the underprediction of the mixture fraction mean (bottom pane in Figure 2(b)). The results obtained with the algebraic modelling of χ
~ (see (26)) tend to be underpredictive near the centreline and overpredictive around the peaks of the experimental data (not shown here). Therefore, it is concluded that the alternative approach of incorporating an additional transport equation for χ
~ (see (27)) provides more accurate mixture fraction variance predictions. Overall, the results presented in Figure 2 show that the inert flow calculations employed in this study provide reliable flow and mixing fields for the CMC calculations.Radial profiles of (a)u
~, (b) ξ
~, and (c) ξ
′′
2
~. Lines, numerical results; symbols, experimental data (u
~ obtained from [44] and ξ
~ and ξ
′′
2
~ obtained from [3]). No velocity measurements are available at x
/
d
=
26.
(a)
u
~
(b)
ξ
~
(c)
ξ
′′
2
~The centreline evolution of the time scale ratio,C
χ, is shown in Figure 3. Here, C
χ is computed using (26) with χ
~ obtained from the solution of (27). In the vicinity of the nozzle, C
χ is well above the commonly used constant value of 2. Away from the nozzle, C
χ decays gradually and tends asymptotically to 2.76. This suggests that the usage of (26) with C
χ
=
2 would generally underpredict the centreline scalar dissipation levels, especially near the nozzle area. The trend obtained in Figure 3 is in line with the experimental observations of Markides [45] who performs scalar dissipation measurements in autoignitive nitrogen-diluted fuel jets injected into coflowing heated air, a configuration similar to the one under investigation. His results show that the centreline values of C
χ decay from high levels near the injector to a constant value around 2 at downstream locations. Accordingly, in a subsequent CMC study in the same context, Markides et al. [46] tune C
χ in (26) in the vicinity of the injector in order to emulate the decay of the time scale ratio and set it to 2 downstream. In the absence of scalar dissipation measurements, as in the present case, such treatment isad hoc and tedious: C
χ must be tuned by trial and error until a reasonable agreement between the predicted and measured ξ
′′
2
~ profiles is achieved, and most likely, a nonconstant functional form for C
χ may be necessary. This process is easily avoided by means of solving (27).Figure 3
Evolution of the time scale ratio along the centreline.


### 4.2.  Effect of the CSDR Modelling for Various Coflow Temperatures

The effect of the modelling of the CSDR is investigated using the homogeneous model of Girimaji [21] and the inhomogeneous model of Mortensen [22]. The reader is reminded that the latter model degenerates to the former when the inhomogeneous terms in (17) are discarded. Calculations are performed for three coflow temperatures. In addition to the experimentally reported value T
c
=
1045 K [3], two other temperatures are investigated: T
c
=
1030 and 1060 K. Given the experimental uncertainty of 3% in the temperature measurements [3], these two values lie within the experimental error. The chemical kinetics mechanism of Mueller et al. [40] is employed throughout this section.


#### 4.2.1. Results in Physical Space

The radial profiles of the Favre-averaged reactive scalars are shown in Figures4(a)–4(e) at the axial locations x
/
d
=
8, 9, 10, 11, 14, and 26. For T
c
=
1045 K (thick black lines), the results obtained using the two CSDR models are almost identical, except for some minor differences at x
/
d
=
10 and 11. When T
c is decreased to 1030 K (thin black lines), the predictions remain in good agreement with the experiments up to x
/
d
=
10 and the two CSDR models yield very similar results. Farther downstream, the temperature and products mass fractions are underpredicted and the reactants mass fractions are overpredicted, most notably at x
/
d
=
11 and 14. The influence of the CSDR modelling becomes more apparent at these axial locations. Overall, the results remain in reasonable agreement with the experiments. When T
c is increased to 1060 K (grey lines), the CSDR models yield distinct results for x
/
d
≤
14. Using both models, the temperature and products mass fractions are grossly overpredicted and the reactants mass fractions are underpredicted. All three coflow temperatures result in very similar profiles at x
/
d
=
26, irrespective of the CSDR model. As a general observation in Figure 4, the usage of Girimaji’s model always yields lower temperatures and product mass fractions and higher reactant mass fractions when compared to Mortensen’s, which leads to a shorter liftoff height. This is attributed to the fact that Girimaji’s model results in earlier ignition delays and hence shorter ignition kernels for all coflow temperatures. This phenomenon will be discussed shortly in more detail. Before doing so, the liftoff heights obtained from the cases considered in Figure 4 are first presented. The contours of Y
~
OH
(
×
1
0
-
3
) are shown in Figure 5(a) for T
c
=
1030 K, Figure 5(b) for T
c
=
1045 K, and Figure 5(c) for T
c
=
1060 K. In each subfigure, the left and right panes correspond to the results obtained with Mortensen’s and Girimaji’s models, respectively. As shown in Figure 5(a), stabilization takes place at slightly distinct locations due to the very small Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1030 K). Mortensen’s model (left) results in H
/
d
=
11.91 (19.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
11.48 (14.8% in excess of H
exp
⁡
/
d). The same scenario is observed in Figure 5(b) where T
c
=
1045 K. However, in this case, the location of stabilisation is better predicted. Mortensen’s model (left) yields H
/
d
=
10.61 (6.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) results in H
/
d
=
10.50 (5% in excess of H
exp
⁡
/
d). The advantage of the inclusion of the effects of inhomogeneity Mortensen’s model is judged to be insignificant for these two coflow temperatures. Conversely, the differences in the stabilisation locations are more substantial for T
c
=
1060 K, as displayed in Figure 5(c). This is due to the much larger Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1060 K). Mortensen’s model (left) yields H
/
d
=
7.39 (26.1% below H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
6.61 (33.9% below H
exp
⁡
/
d). As for the radial location of stabilisation (the radial distance from the stabilisation point to the centreline normalised d, denoted here as W
/
d), the values obtained using the two CSDR models for T
c
=
1030 K and 1045 K are the same (1.58 and 1.47, resp.). As for T
c
=
1060 K, the two models yield slightly different values (Mortensen’s, 1.26; Girimaji’s, 1.16). It is obvious from the above results that the liftoff height becomes smaller and that flame base becomes narrower as T
c is decreased. This is due to the fact that at lower T
c, the mixture can autoignite in areas closer to the nozzle which are characterised by high scalar dissipation rate levels. This behaviour will be addressed in more detail in Section 4.2.3.Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
) (e), and Y
~
OH
(
×
1
0
3
) using Girimaji’s and Mortensen’s models for the modelling of the CSDR. In (a)–(e): solid lines, Girimaji’s model; dashed lines, Mortensen’s model; symbols, experimental data [3]. Colours: black (thin), T
c
=
1030K; black (thick), T
c
=
1045K; grey, T
c
=
1060K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
): (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K obtained using Mortensen’s model (left pane) and Girimaji’s model (right pane). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively.
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


#### 4.2.2.  Stabilisation Mechanism

The numerical indicators developed by Gordon et al. [6] for the distinction between stabilisation by autoignition and premixed flame propagation are employed here to determine the underlying controlling mechanism. For this purpose, the transport budget of the temperature and the history of radical build-up ahead of the stabilisation height are analysed.Budgets in Mixture Fraction Space. Figure 6 shows the transport budget of the steady-state conditional temperature equation, that is the individual contributions of the r.h.s. terms of (9). Only the results obtained using Girimaji’s model are presented here. Mortensen’s model produces similar trends. For each T
c, the budgets are displayed at several axial locations around the corresponding liftoff height. The radial locations are held fixed and set to the W
/
d value obtained with each T
c. As such, the chosen coordinates cover the preflame, liftoff, and postflame regions. Before analysing the balance in these budgets, the axial evolution of the chemical source term, T
CS, is first examined. As shown in Figures 6(a), 6(b), and 6(c), common trends are observed in the evolution of this term for all three coflow temperatures. Below the stabilisation height (first two locations in the preflame regions, Pre-F1 and Pre-F2), as x
/
d increases, T
CS shifts from very lean to less lean mixtures and its amplitude increases dramatically. At liftoff (middle panes), T
CS peaks at η
≈
0.24 with a significantly higher amplitude. Within approximately two to three nozzle diameters downstream the stabilisation height (first locations in the postflame regions, Post-F1) the peak of T
CS increases further and occurs around the stoichiometric mixture fraction. Further downstream (second locations in the postflame regions, Post-F2), the magnitude of T
c decreases substantially from the Post-F1 locations. The reaction zone becomes wider and presents two peaks. The first occurs around the stoichiometric mixture fraction and the second in rich mixtures. The second peak is most likely attributed to the propagation of a rich reaction zone towards stoichiometric mixtures [24]. All of the above trends are consistent with the direct numerical simulation Yoo et al. [47]. As in the evolution of T
CS, common trends are observed in the axial variation of the remaining terms contributing to the r.h.s. of (9). In the preflame region (locations Pre-F1 and Pre-F2), there is a clear balance in lean mixtures between T
CS, the axial convection term, T
C
,
x, and micromixing, T
MM. The axial and radial diffusion terms, T
D
,
x and T
D
,
y, and the radial convection term, T
C
,
y, are nonnegligible but have little contribution to the overall budget. This balance suggests that stabilisation by premixed flame propagation wherein T
D
,
x balances T
C
,
x in the preheat zone with negligible T
CS is not the case. Instead, stabilisation occurs via autoignition as can be seen from the T
CS-T
C
,
x-T
MM balance. At liftoff, the balance shifts to η
≈
0.24. The role of T
C
,
x in balancing T
CS diminishes but this term remains a major contributor to the overall balance, and T
MM prevails. As for the remaining terms, T
D
,
x becomes much weaker while T
D
,
y and T
C
,
y remain important. In the postflame regions (locations Post-F1 and Post-F2), the role of T
c
,
x diminishes further, T
MM becomes more dominant, and T
D
,
x is virtually null. The contributions of T
D
,
y and T
C
,
y are as significant as that of T
C
,
x at Post-F1 and supersede it at Post-F2. As T
CS increases in mixture fraction space, the contributions of T
D
,
y, T
C
,
y, and T
C
,
x remain smaller than that of T
MM, which acts as the major heat sink. In contrast, as T
CS decays, T
MM increases and acts as a source (notice its positive contribution), and both of these terms are counterbalanced by T
D
,
y, T
C
,
y, and T
C
,
x. Thus, beyond the stabilisation height, the flame budgets indicate the structure of a nonpremixed flame, which is largely characterised by a T
CS-T
MM balance. To be noted that the radiative source, T
RS, is negligible at all locations for all coflow temperatures due to the fact that hydrogen in the fuel stream is highly diluted with nitrogen. Therefore, this term can be safely neglected without loss of accuracy.Transport budget of the steady-stateQ
T equation (r.h.s. terms of (9)) at several axial locations around the stabilisation height: (a) T
c
=
1030K (y
/
d
=
1.58), (b) T
c
=
1045K (y
/
d
=
1.47), and (c) T
c
=
1060K (y
/
d
=
1.16). T
RS, radiative source; T
CS, chemical source; T
MM, micromixing; T
C
,
x, axial convection; T
C
,
y, radial convection; T
D
,
x, axial diffusion; T
D
,
y, radial diffusion. The vertical dashed line corresponds to the location of the stoichiometric mixture fraction. All terms are scaled down by a factor of 103 and the units are K/s. The CSDR is modelled using Girimaji’s model.
(a)
T
c
=
1030 K, y
/
d
=
1.58
(b)
T
c
=
1045 K, y
/
d
=
1.47
(c)
T
c
=
1060 K, y
/
d
=
1.16Budgets in Physical Space. The fact that autoignition is the controlling stabilisation mechanism can be further confirmed by analysing the transport budget of Q
T in physical space. To do so, the PDF-weighted integration of the r.h.s. terms of (9) is performed over the mixture fractions space. This yields I
T
RS, I
T
CS, I
T
MM, IT
C
,
x, I
T
C
,
y, IT
D
,
x, and IT
D
,
y (“IT” stands for “Integrated Transport”). The axial profiles of these terms are plotted in Figure 7(a) for Girimaji’s model and in Figure 7(b) for Mortensen’s model at the respective W
/
d locations obtained with each T
c. The parallels between the observation made in Figure 6 and this figure are evident. I
T
CS is mostly balanced by IT
C
,
x and I
T
MM in the preflame regions, with the former being more dominant. The remaining terms are small but nonnegligible. The absence of IT
D
,
x eliminates the possibility of stabilisation by premixed flame propagation and the current balance indicates that autoignition is the controlling mechanism. Right ahead of liftoff, IT
C
,
y and IT
D
,
y emerge and I
T
MM becomes more important. These three terms start to dominate at the expense of IT
C
,
x, which remains the main heat sink. Downstream of the stabilisation height, IT
C
,
x, IT
C
,
y, and IT
D
,
y diminish gradually and I
T
MM becomes more important. Further downstream I
T
CS is mostly counterbalanced by I
T
MM, with smaller contributions IT
C
,
x, IT
C
,
y, and IT
D
,
y, which indicates the structure of a nonpremixed flame.Axial profiles of the Integrated Transport (IT) budget of the steady-stateQ
T equation (integrated r.h.s. terms of (9)): (a) Girimaji’s model and (b) Mortensen’s model. The bottom, middle, and top panes correspond to T
c
=
1030K, T
c
=
1045K, and T
c
=
1060K, respectively. The IT subscripts follow the same notation used in Figure 6. All terms are scaled down by a factor of 106 and the units are K/s.
(a)
Girimaji's model
(b)
Mortensen's modelRadical History ahead of the Stabilisation Height. The transport budgets of the temperature in mixture fraction and physical spaces are both indicative of stabilisation via autoignition. In order to ascertain this conclusion, the history of radical build-up ahead of the reaction zone is now investigated. Figure 8 displays the axial profiles of the normalised Favre-averaged temperature ((
T
~
-
T
~
min
⁡
)
/
(
T
~
max
⁡
-
T
~
min
⁡
)) and mass fractions (Y
~
/
Y
~
max
⁡) of H, O, OH, and H
O
2. The subscripts “min” and “max” denote the minimum and maximum values of the reactive scalars at y
/
d
=
W
/
d. The results are reported for Girimaji’s and Mortensen’s models in Figures 8(a) and 8(b), respectively. Common trends are observed ahead of the stabilisation height for the different combinations of T
c values and CSDR models.(1)
The production of the intermediateH
O
2 initiates upstream of the axial locations where the radicals H, O, and OH emerge.
(2)
H, O, and OH start building up before the consumption of H
O
2 begins. As shown, the mass fractions of these three species start to increase before H
O
2 reaches its peak.
(3)
H
O
2 builds up rapidly prior to the runaway of H, O, and OH, as can be seen from the axial locations of the maximum slopes.Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of stabilisation corresponding to each T
c: (a) Girimaji’s model and (b) Mortensen’s model. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff height, respectively. The circles indicate the locations of the maximum slopes.
(a)
Girimaji's model
(b)
Mortensen's modelTherefore, radicals do not build up simultaneously as in premixed flame propagation, but ratherH
O
2 acts as a precursor to the production of H, O, and OH as in autoignition scenarios [6]. This confirms that the flame is stabilised by autoignition.


#### 4.2.3.  Impact of CSDR Modelling on Autoignition

To better understand the autoignition process, standalone zero-dimensional CMC (0DCMC) calculations are performed for eachT
c. In these calculations, the convective and (physical space) diffusive terms in (8) and (9) are turned off and the CSDR is modelled using the AMC [9]. This CSDR model is parametrised by its maximum value χ
∘ and it is given by
(29)
〈
χ
∣
η
〉
=
χ
∘
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
-
1
)
]
2
}
=
χ
∘
G
(
η
)
,
where er
f
-
1 is the inverse error function. The resulting homogeneous CMC equations are identical to those of the unsteady laminar flamelet model (with unity Lewis number) [48]. For a given χ
∘, the solution evolves from the inert mixing solution at t
=
0 up to the moment when autoignition occurs, t
=
τ. Autoignition is declared at the moment when the maximum OH mass fraction reaches 2
×
1
0
-
4 at any point in mixture fraction space, following Stanković and Merci [49]. Figure 9 displays the ignition delay (τ) as a function of χ
∘. These results are obtained using the mechanism of Mueller et al. [40] for T
c
=
1030 K, 1045 K, and 1060 K. It is evident from this figure that, for a given coflow temperature, τ reaches an asymptotic limit as χ
∘ is gradually increased. This means that the occurrence of autoignition becomes impossible beyond a critical value of χ
∘, denoted here by χ
∘
,
c. The approximate values of χ
∘
,
c are 227 s−1, 482 s−1, and 991 s−1 for T
c
=
1030 K, 1045 K, and 1060 K, respectively. Therefore, χ
∘
,
c increases with increasing T
c. In order to link these findings to the flow field of the flame under investigation, χ
∘ is first obtained by integrating (29) weighted by P
~
(
η
) over the mixture fraction space:
(30)
χ
∘
=
∫
0
1
‍
〈
χ
∣
η
〉
P
~
(
η
)
d
η
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
=
χ
~
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
which requires the usage of the readily available quantities χ
~, ξ
~, and ξ
′′
2
~. Figure 10 shows the Contours of χ
∘ computed from (30) for T
c
=
1045 K. The contours for T
c
=
1030 K and 1060 K are virtually the same and are not included here for brevity. The contour of Figure 10 will be used in all of the subsequent analysis irrespective of the value of T
c. The thick contours represent the levels 227 s−1, 482 s−1, and 991 s−1 which correspond to the χ
∘
,
c values determined in Figure 9. Also shown are the ignition kernel locations obtained from the different combinations of CSDR models and coflow temperatures. As expected, for all three coflow temperatures, ignition occurs at locations where χ
∘ is lower that χ
∘
,
c (the circles, squares, and triangles lie outside the 227 s−1, 482 s−1, and 991 s−1 contour levels, resp.). Since χ
~ decays away from the nozzle, autoignition does not take place until χ
∘ (which depends on χ
~ via (30)) falls below χ
∘
,
c. This explains why increasing T
c results in earlier ignition and shorter liftoff from a mixing point of view.Figure 9
Ignition delays as a function of the peak scalar dissipation rate obtained using the mechanism of Mueller et al. [40] for T
c
=
1030, 1045, and 1060 K. The vertical dashed lines represent the asymptotic limit beyond which ignition cannot occur (χ
∘
=
χ
∘
,
c).Figure 10
Contours ofχ
∘ (see (30)). The thick contours correspond to χ
∘
=
χ
∘
,
c for T
c
=
1030 K (227 s−1), 1045 K (482 s−1), and 1060 K (991 s−1). The symbols represent the ignition locations.In another set of 0DCMC calculations, the “most reactive” mixture fraction,η
mr, is computed for each T
c by setting χ
∘ in (29) to 10−20 s−1 (homogeneous reactor). The results are summarised in Table 4. Consistent with the findings of Stanković and Merci [49], η
mr decreases with increasing T
c. The objective behind determining η
mr is to compare the local CSDR values obtained using Mortensen’s and Girimaji’s models at η
=
η
mr, 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G, respectively. As shown in Figure 11, for all given coflow temperatures and radial locations, the axial profiles of 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G reveal that Mortensen’s models generally result in higher scalar dissipation rate levels and therefore yield longer ignition delays and higher liftoff heights. Also displayed in Figure 11 are the critical CSDR values computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr; that is, 〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
mr
-
1
)
]
2
} (horizontal grey lines). These values are also provided in Table 4 for all coflow temperatures. As expected, 〈
χ
∣
η
mr
〉
c
AMC increases with increasing T
c. Figure 11(a) (T
c
=
1030 K) shows that ignition cannot occur as long as 〈
χ
∣
η
mr
〉
M or 〈
χ
∣
η
mr
〉
G is below 〈
χ
∣
η
mr
〉
c
AMC. Indeed, when Mortensen’s model is used, ignition happens at the first axial and radial locations where 〈
χ
∣
η
mr
〉
M falls below 〈
χ
∣
η
mr
〉
c
AMC. This takes place at x
/
d
=
11.09 and y
/
d
=
1.69 (first pane from the top in Figure 11(a)). A similar scenario is observed when Girimaji’s model is used. Ignition occurs at the first axial and radial locations where 〈
χ
∣
η
mr
〉
G falls below 〈
χ
∣
η
mr
〉
c
AMC, which happens at x
/
d
=
10.81 and y
/
d
=
1.58 (second pane from the top in Figure 11(a)). Similar observation can be made in Figure 11(b) (T
c
=
1045 K) and Figure 11(c) (T
c
=
1060 K). However, in these two cases, although 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G fall quickly below 〈
χ
∣
η
mr
〉
c
AMC at y
/
d locations smaller than the ones where ignition is indicated, ignition occurs at the first axial location with lowest 〈
χ
∣
η
mr
〉. As shown in Figure 11(b) (T
c
=
1045 K), ignition takes place at x
/
d
=
9.97 and y
/
d
=
1.58 using Mortensen’s model and at x
/
d
=
9.69 and y
/
d
=
1.69 using Girimaji’s model. Similarly, in Figure 11(c) (T
c
=
1060 K), ignition happens at x
/
d
=
6.90 and y
/
d
=
1.26 using Mortensen’s model and at x
/
d
=
6.04 and y
/
d
=
1.26 using Girimaji’s model.Table 4
The most reactive mixture fraction,η
mr, and the CSDR values obtained from the AMC parametrised by χ
∘
,
c at η
mr.
T
c [K]
η
mr
〈
χ
∣
η
mr
〉
c
AMC
†
1030
0.0620
21.27
1045
0.0543
36.65
1060
0.0465
59.02
†
〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
erf
⁡
-
1
(
2
η
mr
-
1
)
]
2
}  [s−1].Axial profiles of〈
χ
∣
η
mr
〉 at several radial locations. (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K. Solid black lines, Girimaji’s model; dashed black lines, Mortensen’s model. The vertical black lines indicate the radial locations where ignition occurs. The grey horizontal lines represent the value of the CSDR computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr (〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1 (
2
η
mr
-
1
)
]
2
}).
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


#### 4.2.4.  Effect of Spurious Sources

Ideally, the integration of the PDF-weighted CMC equations overη-space should yield the unconditional set of equations without any additional source terms [14]. This outcome is guaranteed when fully consistent CSDR models (e.g., the inhomogeneous model of Mortensen) are employed. To be noted that in order for an inhomogeneous CSDR model to be fully consistent, the CV model employed in the derivation should be consistent with the moments of the presumed PDF. This condition is satisfied in Moretnsen’s model due to the usage of the PDF-gradient model. An example of an inhomogeneous yet inconsistent CSDR closure is the one proposed in [50]. The discrepancy arises from the modelling of the CV using the linear velocity model, which is not consistent with the second moment of the presumed β-PDF.When inconsistent CSDR models are employed, the integration of the CMC equations results in spurious (or false) source terms. In this study, this occurs when the CSDR model of Girimaji is used. The spurious source associated with a speciesκ is calculated through [14, 16]:
(31)
S
~
κ
=
1
2
∫
0
1
[
〈
χ
∣
η
〉
G
︸
Inonsistent
CSDR
model
-
〈
χ
∣
η
〉
M
︸
Consistent
CSDR
model
]
∂
2
Q
κ
,
G
∂
η
2
P
~
(
η
)
d
η
,
where Q
κ
,
G is the conditional mass fraction of species κ obtained in the inconsistent realisation (using Girimaji’s model).Figure12 shows the radial profiles of the mean chemical sources (I
T
CS) and spurious sources of O
2, H
2, H
2
O, and OH for T
c
=
1045 K. In comparison to the mean chemical sources, the magnitudes of the spurious sources are small but nonnegligible, in particular within the flame zone (x
/
d
=
14 and 26). The effect of the spurious sources ahead of the predicted stabilisation heights (x
/
d
=
9 and 10) is not sufficiently large to change the nature of the stabilisation mechanism, as the flame is stabilised by autoignition using both Mortensen’s and Girimaji’s models. However, the fact that Girimaji’s model yields a relatively earlier ignition and a smaller liftoff height indicates that the spurious sources are influential. Although this is not the case in the calculations of the current flame, inconsistent CMC implementations may lead erroneous results and conclusions.Radial profiles of the chemical and spurious sources forT
c
=
1045 K: (a) O
2, (b) H
2, (c) H
2
O, and (d) OH. Thin lines, chemical sources; thick lines, spurious sources. The units are s−1.
(a)
Y
~
O
2
(b)
Y
~
H
2
(c)
Y
~
H
2
O
(d)
Y
~
OH


### 4.3.  Effect of Chemical Kinetics

The mechanism of Mueller et al. [40] has been used in all of the results reported so far. The sensitivity of the flame to chemical kinetics is now assessed by testing the performance of the mechanisms of Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43]. The coflow temperature is set to the experimentally reported value of 1045 K.


#### 4.3.1. Comparison of the Chemical Kinetic Mechanisms

0DCMC A Priori Analysis. 0DCMC calculations are first performed using the three mechanisms in order to determine the corresponding χ
∘
,
c values. The results are shown in Figure 13, along with those obtained with the Mueller mechanism in Section 4.2.3 for T
c
=
1045K. By inspecting the calculated χ
∘
,
c values and referring to Figure 10, it can be postulated that the Li mechanism would yield the shortest ignition kernel and smallest liftoff height followed by the Mueller mechanism, then the Ó Conaire mechanism, and finally the Burke mechanism. This analysis provides a preliminary qualitative description of the sensitivity of liftoff height to the considered chemical kinetics.Figure 13
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the mechanisms of Mueller et al. [53], Ó Conaire et al. [43], Li et al. [41], and Burke et al. [42] for T
c
=
1045K. The inset shows the critical values of χ
° above which ignition cannot occur.Results in Physical Space. To ascertain these findings, the flame is simulated using the three newly considered mechanisms. The CSDR is closed using Girimaji’s model in all of the following calculations. The results are displayed in Figure 14 at the axial locations x
/
d = 9, 10, 11, 14, and 26, along with those obtained using Mueller mechanism in Section 4.2.1. Unlike the Burke mechanism, the other three kinetic schemes are capable of reproducing the experimental trends. The results of the Mueller and Ó Conaire mechanisms show best agreement with the experimental data. The Mueller mechanism yields better predictions at x
/
d
=
9 and 10, while Ó Conaire’s results in better agreement with the measurements at x
/
d
=
11. The differences between the two are small at x
/
d
=
14 and 26. When the Li mechanism is employed, the temperature rises at early axial locations owing to the ability of the mixture to ignite at high scalar dissipation rate levels (χ
∘
,
c
Li
=
562 s−1). This effect propagates downstream and results in the deviation of the numerical predictions from the experiments. The performance of the Burke mechanisms is unsatisfactory, as the mixture fails to achieve early ignition (χ
∘
,
c
Burke
=
169s−1) and remains inert up to approximately x
/
d
=
14. The contours of Y
~
OH (×
1
0
3) are displayed in Figure 15. As predicted from the 0DCMC calculations, the Li and Burke mechanisms yield the smallest (H
/
d
=
9.73) and largest (H
/
d
=
16.11) liftoff heights, as shown in Figures 15(a) and 15(b), respectively. Also being consistent with the 0DCMC findings, the liftoff heights obtained using the Mueller and Ó Conaire mechanisms fall between those of the Li and Burke mechanisms, as demonstrated in Figures 15(c) and 15(d). However, the Mueller mechanism results in a slightly larger liftoff height (H
/
d
=
10.50) in comparison to the Ó Conaire mechanism (H
/
d
=
10.26), which contradicts with the 0DCMC predictions. Apart from the Burke mechanism, the liftoff heights computed from the remaining mechanisms are similar and are in close agreement with the experimental value H
exp
⁡
/
d
=
10. As for the radial location of liftoff, the Mueller, Li, and Ó Conaire mechanisms result is W
/
d
=
1.47, whereas Burke’s yields a wider flame base with W
/
d
=
1.92. As in the Mueller mechanism calculations, the flame is also stabilised by autoignition when the Ó Conaire, Li, and Burke mechanisms are employed. This can be clearly seen in Figure 16 where H, O, OH, and H
O
2 exhibit the same trends observed in Figure 8(a).Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
) using different chemical kinetic mechanisms. In (a)–(e): thick solid lines, Mueller et al. [40] mechanism; thick dashed lines, Li et al. [41] mechanism; thin solid lines, Burke et al. [42] mechanism; thin dashed lines, Ó Conaire et al. [43] mechanism; Symbols, experimental data [3]. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
) obtained using the mechanisms of (a) Li et al. [41], (b) Burke et al. [42], (c) Mueller et al. [40], and (d) Ó Conaire et al. [43]. The thick contours correspond to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Li
(b)
Burke
(c)
Mueller
(d)
Ó ConaireFigure 16
Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of liftoff using the mechanisms of (from top to bottom) Burke et al. [42], Ó Conaire et al. [43], and Li et al. [41]. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff, respectively. The circles indicate the locations of the maximum slopes. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.Results in Mixture Fraction Space. For completeness, the conditional profiles are shown in Figure 17. The calculation of the conditional data from the experimental scatter is described in the Appendix. The numerical results are reported at x
/
d = 9, 10, 11, 14, and 26, and y
/
d
=
1.47 (W
/
d obtained with all mechanisms except Burke’s). Overall, the results of the Mueller and Ó Conaire mechanisms are in reasonable agreement with the experiments, especially within mixture fraction intervals where the PDF is large (see Figure 17(f)). This explains the good agreement these two mechanisms yield with the measurements in physical space (Figure 14). The Li mechanism predicts the highest OH levels up to x
/
d
=
11 (Figure 17(e)) which leads to the occurrence of the earliest ignition and consequently results in the smallest liftoff height among the considered kinetic schemes. In turn, this early ignition produces higher than expected temperature levels, an effect which propagates downstream resulting in the deviation of the conditional predictions from the experiments. In the three mechanisms discussed up to this point, the peak of Q
T in Figure 17(a) shifts gradually from lean to stoichiometric mixtures as x
/
d increases, which is consistent with the axial evolution of the conditional heat release term T
CS in Figure 6. The profiles of the Burke mechanism remain very close to the inert mixing solution up to x
/
d
=
14. This behaviour results in the most delayed ignition and leads to the largest liftoff height. The conditional profiles at y
/
d
=
1.92 where liftoff occurs undergo very similar trends (not shown here).Conditional profiles atx
/
d
=
9, 10, 11, 14 and 26 and y
/
d
=
1.47 using different chemical kinetic mechanisms: (a) Q
T, (b) Q
O
2, (c) Q
H
2
(
×
1
0
2
), (d) Q
H
2
O
(
×
1
0
2
), (e) Q
OH
(
×
1
0
3
), and (f) P
~
(
η
). In (a)–(e): thick solid lines, Mueller et al. [40]; thick dashed lines, Li et al. [41]; thin solid lines, Burke et al. [42]; thin dashed lines, Ó Conaire et al. [43]; symbols, experimental data (see the Appendix). The vertical dashed line indicates the location of the stoichiometric mixture fraction. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Q
T
[
K
]
(b)
Q
O
2
(c)
Q
H
2
×
1
0
2
(d)
Q
H
2
O
×
1
0
2
(e)
Q
OH
×
1
0
3
(f)
P
~
(
η
)


#### 4.3.2.  Diagnosis of the Burke Mechanism

In order to diagnose the unsatisfactory performance of the Burke mechanism, the rate parameters of two of the key reactions are investigated. The recombination reaction(
R
9
)
H
+
O
2
(
+
M
)
⇌
H
O
2
(
+
M
)
and the chain-branching reaction
(
R
1
)
H
+
O
2
⇌
O
+
OH
are two of the most important reactions in H
2 oxidation chemistry. (
R
9
) and (
R
1
) (labelled here as in [40–43]) compete for H atoms and therefore control the overall branching ratio and determine the second explosion limit in homogeneous H
2
/
O
2 systems [42]. The treatment of (
R
9
) and (
R
1
) in the Li and Burke mechanisms differs significantly from the Ó Conaire and Mueller mechanisms. The rate parameters of (
R
9
) are provided in Table 5. The high-pressure limit rate parameters are identical in the Ó Conaire, Mueller, and Li mechanisms [51]. Burke’s mechanism uses an updated set of parameters [52]. In the low-pressure limit, the Ó Conaire and Mueller mechanisms use the same parameters [53] except for the third-body efficiency factors of H
2 and H
2
O. The Li mechanism employs a completely different set of parameters [54]. The Burke parameters differ from Li’s by the broadening factor (decreased by a factor of 1.6 and equal to the ones employed in the Ó Conaire and Mueller mechanisms) and the third-body efficiency factor for H
2
O (increased by a factor of 1.3) [42]. As for (
R
1
), the corresponding rate parameters are listed in Table 6. The parameters employed in the Ó Conaire and Mueller mechanisms are identical [55]. Two distinct sets of updated parameters are employed in the mechanisms of Li [56] and Burke [57]. In light of the above description, the variability of the predictions in Figure 14 and the close agreement between the results of the Mueller and Ó Conaire mechanisms may be explained in part by the treatment of (
R
9
) and (
R
1
). This can be demonstrated, for instance, by adjusting the rate parameters of these two reactions in the Burke mechanism, in an attempt to improve its predictions. Since this mechanism is an updated version of Li’s, two separate experiments are carried out as follows. (1)
The parameters of(
R
9
) are substituted by those of the Li mechanism, with the parameters of  (
R
1
) left unchanged.
(2)
The parameters of both(
R
9
) and (
R
1
) are substituted by those of the Li mechanism.Table 5
Rate parameters of reaction(
R
9
), H + O2 (+M) ⇌ HO2 (+M).
Ó Conaire
Mueller
Li
Burke
High-pressure limit
A
1.475
×
10
12
1.475
×
10
12
1.475
×
10
12
4.65084
×
10
12
n
0.6
0.6
0.6
0.44
E
a
0
0
0
0
Low-pressure limit
A
3.482
×
10
16
3.482
×
10
16
6.366
×
10
20
6.366
×
10
20
n
−0.411
−0.411
−1.72
−1.72
E
a
-
1.115
×
10
3
-
1.115
×
10
3
5.248
×
10
2
5.248
×
10
2
F
c
0.5
0.5
0.8
0.5
ε
H
2
1.3
2.5
2
2
ε
H
2
O
14
12
11
14
ε
O
2
1
1
0.78
0.78
Units are cm3-mole-sec-kcal-K; k
=
A
T
n
exp
⁡
(
-
E
a
/
R
T
).
Nomenclature:A: preexponential coefficient; n: temperature exponent; E
a: activation energy; F
c: broadening factor; ε: third-body efficiency factor.Table 6
Rate parameters of reaction(
R
1
), H + O2
⇌ O + OH.
Ó Conaire
Mueller
Li
Burke
A
1.915
×
10
14
1.915
×
10
14
3.547
×
10
15
1.040
×
10
14
n
0
0
−0.406
0
E
a
1.6439
×
10
4
1.6439
×
10
4
1.6599
×
10
4
1.5286
×
10
4
Same units and nomenclature as in Table5.The preliminary 0DCMC calculations shown in Figure18 reveal that χ
∘
,
c increases significantly when these modifications are applied. Therefore, the mixture is now anticipated to ignite earlier and the liftoff height is expected to decrease. Two additional CMC runs are performed using these modifications while keeping the same settings for the CSDR and the coflow temperature. The results are shown in Figure 19 at x
/
d = 10, 11, and 14. It is evident that the modification of (
R
9
) alone brings a substantial improvement over the rate parameters of the original mechanism. The temperature (Figure 19(a)) and the species mass fractions (Figures 19(b)–19(e)) are much better predicted, especially at x
/
d
=
11 and 14. This further proves the importance of (
R
9
) as a precursor to autoignition reactions. Better agreement with the experimental data is obtained when both (
R
9
) and (
R
1
) are modified. As shown in Figure 19(f), the liftoff height decreases by 26.69% when (
R
9
) is modified and by 29.24% when both (
R
9
) and (
R
1
) are modified. In both cases the computed liftoff heights are substantially closer to the experimentally measured value. Figure 19(f) also shows a remarkable decrease in W
/
d as liftoff occurs closer to the centreline. With the modified settings, both the axial and radial positions of liftoff become in line with those obtained using the Ó Conaire, Mueller, and Li mechanisms. There should be further room for improvement by applying additional modifications to the rate parameters of other key reactions. However, extreme care must be taken while doing so because chemical kinetic mechanisms are usually optimised as a whole, rather than on an individual reaction basis. A sensitivity analysis is beyond the scope of this work and was not performed here. The modifications applied above aread hoc.Figure 18
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the original and modified mechanisms of Burke et al. [42]. The inset shows the critical values of  χ
° above which ignition cannot occur.Radial profiles obtained using the original and modified versions of the Burke mechanism: (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
). In (a)–(e): dashed lines, original mechanism; thin solid lines, modified (
R
9
); thick solid lines, modified (
R
9
) and (
R
1
); symbols, experimental data [3]. Subfigure (f) shows the contours of Y
~
OH
(
×
1
0
3
). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3
(f)
Contours ofY
~
OH
(
×
1
0
3
)
n


## 4.1. Flow Field Results

The radial profiles of the Favre-averaged axial velocity are compared to the experimental measurements of Kent [44] in Figure 2(a). The predictions are in very good agreement with the experimental data up to x
/
d
=
11. The measurements are slightly underpredicted at x
/
d
=
14 for y
/
d
>
0.5. No velocity measurements are available at x
/
d
=
26. The profile at this location is provided for completeness. The radial profiles of the Favre-averaged mixture fraction mean are displayed in Figure 2(b) along with the experimental data of Cabra et al. [3]. The predictions are in very good agreement with the experiments at x
/
d = 1, 8, 9, 10, and 11. Slight underpredictions are obtained at x
/
d
=
14 and 26 for 0
≤
y
/
d
≤
1. Figure 2(c) shows a comparison of the radial Favre-averaged mixture fraction variance profiles with the experimental measurements. A very good agreement is obtained at all axial locations except at x
/
d
=
26. The discrepancies at this location are attributed to the underprediction of the mixture fraction mean (bottom pane in Figure 2(b)). The results obtained with the algebraic modelling of χ
~ (see (26)) tend to be underpredictive near the centreline and overpredictive around the peaks of the experimental data (not shown here). Therefore, it is concluded that the alternative approach of incorporating an additional transport equation for χ
~ (see (27)) provides more accurate mixture fraction variance predictions. Overall, the results presented in Figure 2 show that the inert flow calculations employed in this study provide reliable flow and mixing fields for the CMC calculations.Radial profiles of (a)u
~, (b) ξ
~, and (c) ξ
′′
2
~. Lines, numerical results; symbols, experimental data (u
~ obtained from [44] and ξ
~ and ξ
′′
2
~ obtained from [3]). No velocity measurements are available at x
/
d
=
26.
(a)
u
~
(b)
ξ
~
(c)
ξ
′′
2
~The centreline evolution of the time scale ratio,C
χ, is shown in Figure 3. Here, C
χ is computed using (26) with χ
~ obtained from the solution of (27). In the vicinity of the nozzle, C
χ is well above the commonly used constant value of 2. Away from the nozzle, C
χ decays gradually and tends asymptotically to 2.76. This suggests that the usage of (26) with C
χ
=
2 would generally underpredict the centreline scalar dissipation levels, especially near the nozzle area. The trend obtained in Figure 3 is in line with the experimental observations of Markides [45] who performs scalar dissipation measurements in autoignitive nitrogen-diluted fuel jets injected into coflowing heated air, a configuration similar to the one under investigation. His results show that the centreline values of C
χ decay from high levels near the injector to a constant value around 2 at downstream locations. Accordingly, in a subsequent CMC study in the same context, Markides et al. [46] tune C
χ in (26) in the vicinity of the injector in order to emulate the decay of the time scale ratio and set it to 2 downstream. In the absence of scalar dissipation measurements, as in the present case, such treatment isad hoc and tedious: C
χ must be tuned by trial and error until a reasonable agreement between the predicted and measured ξ
′′
2
~ profiles is achieved, and most likely, a nonconstant functional form for C
χ may be necessary. This process is easily avoided by means of solving (27).Figure 3
Evolution of the time scale ratio along the centreline.


## 4.2.  Effect of the CSDR Modelling for Various Coflow Temperatures

The effect of the modelling of the CSDR is investigated using the homogeneous model of Girimaji [21] and the inhomogeneous model of Mortensen [22]. The reader is reminded that the latter model degenerates to the former when the inhomogeneous terms in (17) are discarded. Calculations are performed for three coflow temperatures. In addition to the experimentally reported value T
c
=
1045 K [3], two other temperatures are investigated: T
c
=
1030 and 1060 K. Given the experimental uncertainty of 3% in the temperature measurements [3], these two values lie within the experimental error. The chemical kinetics mechanism of Mueller et al. [40] is employed throughout this section.


### 4.2.1. Results in Physical Space

The radial profiles of the Favre-averaged reactive scalars are shown in Figures4(a)–4(e) at the axial locations x
/
d
=
8, 9, 10, 11, 14, and 26. For T
c
=
1045 K (thick black lines), the results obtained using the two CSDR models are almost identical, except for some minor differences at x
/
d
=
10 and 11. When T
c is decreased to 1030 K (thin black lines), the predictions remain in good agreement with the experiments up to x
/
d
=
10 and the two CSDR models yield very similar results. Farther downstream, the temperature and products mass fractions are underpredicted and the reactants mass fractions are overpredicted, most notably at x
/
d
=
11 and 14. The influence of the CSDR modelling becomes more apparent at these axial locations. Overall, the results remain in reasonable agreement with the experiments. When T
c is increased to 1060 K (grey lines), the CSDR models yield distinct results for x
/
d
≤
14. Using both models, the temperature and products mass fractions are grossly overpredicted and the reactants mass fractions are underpredicted. All three coflow temperatures result in very similar profiles at x
/
d
=
26, irrespective of the CSDR model. As a general observation in Figure 4, the usage of Girimaji’s model always yields lower temperatures and product mass fractions and higher reactant mass fractions when compared to Mortensen’s, which leads to a shorter liftoff height. This is attributed to the fact that Girimaji’s model results in earlier ignition delays and hence shorter ignition kernels for all coflow temperatures. This phenomenon will be discussed shortly in more detail. Before doing so, the liftoff heights obtained from the cases considered in Figure 4 are first presented. The contours of Y
~
OH
(
×
1
0
-
3
) are shown in Figure 5(a) for T
c
=
1030 K, Figure 5(b) for T
c
=
1045 K, and Figure 5(c) for T
c
=
1060 K. In each subfigure, the left and right panes correspond to the results obtained with Mortensen’s and Girimaji’s models, respectively. As shown in Figure 5(a), stabilization takes place at slightly distinct locations due to the very small Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1030 K). Mortensen’s model (left) results in H
/
d
=
11.91 (19.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
11.48 (14.8% in excess of H
exp
⁡
/
d). The same scenario is observed in Figure 5(b) where T
c
=
1045 K. However, in this case, the location of stabilisation is better predicted. Mortensen’s model (left) yields H
/
d
=
10.61 (6.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) results in H
/
d
=
10.50 (5% in excess of H
exp
⁡
/
d). The advantage of the inclusion of the effects of inhomogeneity Mortensen’s model is judged to be insignificant for these two coflow temperatures. Conversely, the differences in the stabilisation locations are more substantial for T
c
=
1060 K, as displayed in Figure 5(c). This is due to the much larger Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1060 K). Mortensen’s model (left) yields H
/
d
=
7.39 (26.1% below H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
6.61 (33.9% below H
exp
⁡
/
d). As for the radial location of stabilisation (the radial distance from the stabilisation point to the centreline normalised d, denoted here as W
/
d), the values obtained using the two CSDR models for T
c
=
1030 K and 1045 K are the same (1.58 and 1.47, resp.). As for T
c
=
1060 K, the two models yield slightly different values (Mortensen’s, 1.26; Girimaji’s, 1.16). It is obvious from the above results that the liftoff height becomes smaller and that flame base becomes narrower as T
c is decreased. This is due to the fact that at lower T
c, the mixture can autoignite in areas closer to the nozzle which are characterised by high scalar dissipation rate levels. This behaviour will be addressed in more detail in Section 4.2.3.Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
) (e), and Y
~
OH
(
×
1
0
3
) using Girimaji’s and Mortensen’s models for the modelling of the CSDR. In (a)–(e): solid lines, Girimaji’s model; dashed lines, Mortensen’s model; symbols, experimental data [3]. Colours: black (thin), T
c
=
1030K; black (thick), T
c
=
1045K; grey, T
c
=
1060K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
): (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K obtained using Mortensen’s model (left pane) and Girimaji’s model (right pane). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively.
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


### 4.2.2.  Stabilisation Mechanism

The numerical indicators developed by Gordon et al. [6] for the distinction between stabilisation by autoignition and premixed flame propagation are employed here to determine the underlying controlling mechanism. For this purpose, the transport budget of the temperature and the history of radical build-up ahead of the stabilisation height are analysed.Budgets in Mixture Fraction Space. Figure 6 shows the transport budget of the steady-state conditional temperature equation, that is the individual contributions of the r.h.s. terms of (9). Only the results obtained using Girimaji’s model are presented here. Mortensen’s model produces similar trends. For each T
c, the budgets are displayed at several axial locations around the corresponding liftoff height. The radial locations are held fixed and set to the W
/
d value obtained with each T
c. As such, the chosen coordinates cover the preflame, liftoff, and postflame regions. Before analysing the balance in these budgets, the axial evolution of the chemical source term, T
CS, is first examined. As shown in Figures 6(a), 6(b), and 6(c), common trends are observed in the evolution of this term for all three coflow temperatures. Below the stabilisation height (first two locations in the preflame regions, Pre-F1 and Pre-F2), as x
/
d increases, T
CS shifts from very lean to less lean mixtures and its amplitude increases dramatically. At liftoff (middle panes), T
CS peaks at η
≈
0.24 with a significantly higher amplitude. Within approximately two to three nozzle diameters downstream the stabilisation height (first locations in the postflame regions, Post-F1) the peak of T
CS increases further and occurs around the stoichiometric mixture fraction. Further downstream (second locations in the postflame regions, Post-F2), the magnitude of T
c decreases substantially from the Post-F1 locations. The reaction zone becomes wider and presents two peaks. The first occurs around the stoichiometric mixture fraction and the second in rich mixtures. The second peak is most likely attributed to the propagation of a rich reaction zone towards stoichiometric mixtures [24]. All of the above trends are consistent with the direct numerical simulation Yoo et al. [47]. As in the evolution of T
CS, common trends are observed in the axial variation of the remaining terms contributing to the r.h.s. of (9). In the preflame region (locations Pre-F1 and Pre-F2), there is a clear balance in lean mixtures between T
CS, the axial convection term, T
C
,
x, and micromixing, T
MM. The axial and radial diffusion terms, T
D
,
x and T
D
,
y, and the radial convection term, T
C
,
y, are nonnegligible but have little contribution to the overall budget. This balance suggests that stabilisation by premixed flame propagation wherein T
D
,
x balances T
C
,
x in the preheat zone with negligible T
CS is not the case. Instead, stabilisation occurs via autoignition as can be seen from the T
CS-T
C
,
x-T
MM balance. At liftoff, the balance shifts to η
≈
0.24. The role of T
C
,
x in balancing T
CS diminishes but this term remains a major contributor to the overall balance, and T
MM prevails. As for the remaining terms, T
D
,
x becomes much weaker while T
D
,
y and T
C
,
y remain important. In the postflame regions (locations Post-F1 and Post-F2), the role of T
c
,
x diminishes further, T
MM becomes more dominant, and T
D
,
x is virtually null. The contributions of T
D
,
y and T
C
,
y are as significant as that of T
C
,
x at Post-F1 and supersede it at Post-F2. As T
CS increases in mixture fraction space, the contributions of T
D
,
y, T
C
,
y, and T
C
,
x remain smaller than that of T
MM, which acts as the major heat sink. In contrast, as T
CS decays, T
MM increases and acts as a source (notice its positive contribution), and both of these terms are counterbalanced by T
D
,
y, T
C
,
y, and T
C
,
x. Thus, beyond the stabilisation height, the flame budgets indicate the structure of a nonpremixed flame, which is largely characterised by a T
CS-T
MM balance. To be noted that the radiative source, T
RS, is negligible at all locations for all coflow temperatures due to the fact that hydrogen in the fuel stream is highly diluted with nitrogen. Therefore, this term can be safely neglected without loss of accuracy.Transport budget of the steady-stateQ
T equation (r.h.s. terms of (9)) at several axial locations around the stabilisation height: (a) T
c
=
1030K (y
/
d
=
1.58), (b) T
c
=
1045K (y
/
d
=
1.47), and (c) T
c
=
1060K (y
/
d
=
1.16). T
RS, radiative source; T
CS, chemical source; T
MM, micromixing; T
C
,
x, axial convection; T
C
,
y, radial convection; T
D
,
x, axial diffusion; T
D
,
y, radial diffusion. The vertical dashed line corresponds to the location of the stoichiometric mixture fraction. All terms are scaled down by a factor of 103 and the units are K/s. The CSDR is modelled using Girimaji’s model.
(a)
T
c
=
1030 K, y
/
d
=
1.58
(b)
T
c
=
1045 K, y
/
d
=
1.47
(c)
T
c
=
1060 K, y
/
d
=
1.16Budgets in Physical Space. The fact that autoignition is the controlling stabilisation mechanism can be further confirmed by analysing the transport budget of Q
T in physical space. To do so, the PDF-weighted integration of the r.h.s. terms of (9) is performed over the mixture fractions space. This yields I
T
RS, I
T
CS, I
T
MM, IT
C
,
x, I
T
C
,
y, IT
D
,
x, and IT
D
,
y (“IT” stands for “Integrated Transport”). The axial profiles of these terms are plotted in Figure 7(a) for Girimaji’s model and in Figure 7(b) for Mortensen’s model at the respective W
/
d locations obtained with each T
c. The parallels between the observation made in Figure 6 and this figure are evident. I
T
CS is mostly balanced by IT
C
,
x and I
T
MM in the preflame regions, with the former being more dominant. The remaining terms are small but nonnegligible. The absence of IT
D
,
x eliminates the possibility of stabilisation by premixed flame propagation and the current balance indicates that autoignition is the controlling mechanism. Right ahead of liftoff, IT
C
,
y and IT
D
,
y emerge and I
T
MM becomes more important. These three terms start to dominate at the expense of IT
C
,
x, which remains the main heat sink. Downstream of the stabilisation height, IT
C
,
x, IT
C
,
y, and IT
D
,
y diminish gradually and I
T
MM becomes more important. Further downstream I
T
CS is mostly counterbalanced by I
T
MM, with smaller contributions IT
C
,
x, IT
C
,
y, and IT
D
,
y, which indicates the structure of a nonpremixed flame.Axial profiles of the Integrated Transport (IT) budget of the steady-stateQ
T equation (integrated r.h.s. terms of (9)): (a) Girimaji’s model and (b) Mortensen’s model. The bottom, middle, and top panes correspond to T
c
=
1030K, T
c
=
1045K, and T
c
=
1060K, respectively. The IT subscripts follow the same notation used in Figure 6. All terms are scaled down by a factor of 106 and the units are K/s.
(a)
Girimaji's model
(b)
Mortensen's modelRadical History ahead of the Stabilisation Height. The transport budgets of the temperature in mixture fraction and physical spaces are both indicative of stabilisation via autoignition. In order to ascertain this conclusion, the history of radical build-up ahead of the reaction zone is now investigated. Figure 8 displays the axial profiles of the normalised Favre-averaged temperature ((
T
~
-
T
~
min
⁡
)
/
(
T
~
max
⁡
-
T
~
min
⁡
)) and mass fractions (Y
~
/
Y
~
max
⁡) of H, O, OH, and H
O
2. The subscripts “min” and “max” denote the minimum and maximum values of the reactive scalars at y
/
d
=
W
/
d. The results are reported for Girimaji’s and Mortensen’s models in Figures 8(a) and 8(b), respectively. Common trends are observed ahead of the stabilisation height for the different combinations of T
c values and CSDR models.(1)
The production of the intermediateH
O
2 initiates upstream of the axial locations where the radicals H, O, and OH emerge.
(2)
H, O, and OH start building up before the consumption of H
O
2 begins. As shown, the mass fractions of these three species start to increase before H
O
2 reaches its peak.
(3)
H
O
2 builds up rapidly prior to the runaway of H, O, and OH, as can be seen from the axial locations of the maximum slopes.Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of stabilisation corresponding to each T
c: (a) Girimaji’s model and (b) Mortensen’s model. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff height, respectively. The circles indicate the locations of the maximum slopes.
(a)
Girimaji's model
(b)
Mortensen's modelTherefore, radicals do not build up simultaneously as in premixed flame propagation, but ratherH
O
2 acts as a precursor to the production of H, O, and OH as in autoignition scenarios [6]. This confirms that the flame is stabilised by autoignition.


### 4.2.3.  Impact of CSDR Modelling on Autoignition

To better understand the autoignition process, standalone zero-dimensional CMC (0DCMC) calculations are performed for eachT
c. In these calculations, the convective and (physical space) diffusive terms in (8) and (9) are turned off and the CSDR is modelled using the AMC [9]. This CSDR model is parametrised by its maximum value χ
∘ and it is given by
(29)
〈
χ
∣
η
〉
=
χ
∘
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
-
1
)
]
2
}
=
χ
∘
G
(
η
)
,
where er
f
-
1 is the inverse error function. The resulting homogeneous CMC equations are identical to those of the unsteady laminar flamelet model (with unity Lewis number) [48]. For a given χ
∘, the solution evolves from the inert mixing solution at t
=
0 up to the moment when autoignition occurs, t
=
τ. Autoignition is declared at the moment when the maximum OH mass fraction reaches 2
×
1
0
-
4 at any point in mixture fraction space, following Stanković and Merci [49]. Figure 9 displays the ignition delay (τ) as a function of χ
∘. These results are obtained using the mechanism of Mueller et al. [40] for T
c
=
1030 K, 1045 K, and 1060 K. It is evident from this figure that, for a given coflow temperature, τ reaches an asymptotic limit as χ
∘ is gradually increased. This means that the occurrence of autoignition becomes impossible beyond a critical value of χ
∘, denoted here by χ
∘
,
c. The approximate values of χ
∘
,
c are 227 s−1, 482 s−1, and 991 s−1 for T
c
=
1030 K, 1045 K, and 1060 K, respectively. Therefore, χ
∘
,
c increases with increasing T
c. In order to link these findings to the flow field of the flame under investigation, χ
∘ is first obtained by integrating (29) weighted by P
~
(
η
) over the mixture fraction space:
(30)
χ
∘
=
∫
0
1
‍
〈
χ
∣
η
〉
P
~
(
η
)
d
η
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
=
χ
~
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
which requires the usage of the readily available quantities χ
~, ξ
~, and ξ
′′
2
~. Figure 10 shows the Contours of χ
∘ computed from (30) for T
c
=
1045 K. The contours for T
c
=
1030 K and 1060 K are virtually the same and are not included here for brevity. The contour of Figure 10 will be used in all of the subsequent analysis irrespective of the value of T
c. The thick contours represent the levels 227 s−1, 482 s−1, and 991 s−1 which correspond to the χ
∘
,
c values determined in Figure 9. Also shown are the ignition kernel locations obtained from the different combinations of CSDR models and coflow temperatures. As expected, for all three coflow temperatures, ignition occurs at locations where χ
∘ is lower that χ
∘
,
c (the circles, squares, and triangles lie outside the 227 s−1, 482 s−1, and 991 s−1 contour levels, resp.). Since χ
~ decays away from the nozzle, autoignition does not take place until χ
∘ (which depends on χ
~ via (30)) falls below χ
∘
,
c. This explains why increasing T
c results in earlier ignition and shorter liftoff from a mixing point of view.Figure 9
Ignition delays as a function of the peak scalar dissipation rate obtained using the mechanism of Mueller et al. [40] for T
c
=
1030, 1045, and 1060 K. The vertical dashed lines represent the asymptotic limit beyond which ignition cannot occur (χ
∘
=
χ
∘
,
c).Figure 10
Contours ofχ
∘ (see (30)). The thick contours correspond to χ
∘
=
χ
∘
,
c for T
c
=
1030 K (227 s−1), 1045 K (482 s−1), and 1060 K (991 s−1). The symbols represent the ignition locations.In another set of 0DCMC calculations, the “most reactive” mixture fraction,η
mr, is computed for each T
c by setting χ
∘ in (29) to 10−20 s−1 (homogeneous reactor). The results are summarised in Table 4. Consistent with the findings of Stanković and Merci [49], η
mr decreases with increasing T
c. The objective behind determining η
mr is to compare the local CSDR values obtained using Mortensen’s and Girimaji’s models at η
=
η
mr, 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G, respectively. As shown in Figure 11, for all given coflow temperatures and radial locations, the axial profiles of 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G reveal that Mortensen’s models generally result in higher scalar dissipation rate levels and therefore yield longer ignition delays and higher liftoff heights. Also displayed in Figure 11 are the critical CSDR values computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr; that is, 〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
mr
-
1
)
]
2
} (horizontal grey lines). These values are also provided in Table 4 for all coflow temperatures. As expected, 〈
χ
∣
η
mr
〉
c
AMC increases with increasing T
c. Figure 11(a) (T
c
=
1030 K) shows that ignition cannot occur as long as 〈
χ
∣
η
mr
〉
M or 〈
χ
∣
η
mr
〉
G is below 〈
χ
∣
η
mr
〉
c
AMC. Indeed, when Mortensen’s model is used, ignition happens at the first axial and radial locations where 〈
χ
∣
η
mr
〉
M falls below 〈
χ
∣
η
mr
〉
c
AMC. This takes place at x
/
d
=
11.09 and y
/
d
=
1.69 (first pane from the top in Figure 11(a)). A similar scenario is observed when Girimaji’s model is used. Ignition occurs at the first axial and radial locations where 〈
χ
∣
η
mr
〉
G falls below 〈
χ
∣
η
mr
〉
c
AMC, which happens at x
/
d
=
10.81 and y
/
d
=
1.58 (second pane from the top in Figure 11(a)). Similar observation can be made in Figure 11(b) (T
c
=
1045 K) and Figure 11(c) (T
c
=
1060 K). However, in these two cases, although 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G fall quickly below 〈
χ
∣
η
mr
〉
c
AMC at y
/
d locations smaller than the ones where ignition is indicated, ignition occurs at the first axial location with lowest 〈
χ
∣
η
mr
〉. As shown in Figure 11(b) (T
c
=
1045 K), ignition takes place at x
/
d
=
9.97 and y
/
d
=
1.58 using Mortensen’s model and at x
/
d
=
9.69 and y
/
d
=
1.69 using Girimaji’s model. Similarly, in Figure 11(c) (T
c
=
1060 K), ignition happens at x
/
d
=
6.90 and y
/
d
=
1.26 using Mortensen’s model and at x
/
d
=
6.04 and y
/
d
=
1.26 using Girimaji’s model.Table 4
The most reactive mixture fraction,η
mr, and the CSDR values obtained from the AMC parametrised by χ
∘
,
c at η
mr.
T
c [K]
η
mr
〈
χ
∣
η
mr
〉
c
AMC
†
1030
0.0620
21.27
1045
0.0543
36.65
1060
0.0465
59.02
†
〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
erf
⁡
-
1
(
2
η
mr
-
1
)
]
2
}  [s−1].Axial profiles of〈
χ
∣
η
mr
〉 at several radial locations. (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K. Solid black lines, Girimaji’s model; dashed black lines, Mortensen’s model. The vertical black lines indicate the radial locations where ignition occurs. The grey horizontal lines represent the value of the CSDR computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr (〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1 (
2
η
mr
-
1
)
]
2
}).
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


### 4.2.4.  Effect of Spurious Sources

Ideally, the integration of the PDF-weighted CMC equations overη-space should yield the unconditional set of equations without any additional source terms [14]. This outcome is guaranteed when fully consistent CSDR models (e.g., the inhomogeneous model of Mortensen) are employed. To be noted that in order for an inhomogeneous CSDR model to be fully consistent, the CV model employed in the derivation should be consistent with the moments of the presumed PDF. This condition is satisfied in Moretnsen’s model due to the usage of the PDF-gradient model. An example of an inhomogeneous yet inconsistent CSDR closure is the one proposed in [50]. The discrepancy arises from the modelling of the CV using the linear velocity model, which is not consistent with the second moment of the presumed β-PDF.When inconsistent CSDR models are employed, the integration of the CMC equations results in spurious (or false) source terms. In this study, this occurs when the CSDR model of Girimaji is used. The spurious source associated with a speciesκ is calculated through [14, 16]:
(31)
S
~
κ
=
1
2
∫
0
1
[
〈
χ
∣
η
〉
G
︸
Inonsistent
CSDR
model
-
〈
χ
∣
η
〉
M
︸
Consistent
CSDR
model
]
∂
2
Q
κ
,
G
∂
η
2
P
~
(
η
)
d
η
,
where Q
κ
,
G is the conditional mass fraction of species κ obtained in the inconsistent realisation (using Girimaji’s model).Figure12 shows the radial profiles of the mean chemical sources (I
T
CS) and spurious sources of O
2, H
2, H
2
O, and OH for T
c
=
1045 K. In comparison to the mean chemical sources, the magnitudes of the spurious sources are small but nonnegligible, in particular within the flame zone (x
/
d
=
14 and 26). The effect of the spurious sources ahead of the predicted stabilisation heights (x
/
d
=
9 and 10) is not sufficiently large to change the nature of the stabilisation mechanism, as the flame is stabilised by autoignition using both Mortensen’s and Girimaji’s models. However, the fact that Girimaji’s model yields a relatively earlier ignition and a smaller liftoff height indicates that the spurious sources are influential. Although this is not the case in the calculations of the current flame, inconsistent CMC implementations may lead erroneous results and conclusions.Radial profiles of the chemical and spurious sources forT
c
=
1045 K: (a) O
2, (b) H
2, (c) H
2
O, and (d) OH. Thin lines, chemical sources; thick lines, spurious sources. The units are s−1.
(a)
Y
~
O
2
(b)
Y
~
H
2
(c)
Y
~
H
2
O
(d)
Y
~
OH


## 4.2.1. Results in Physical Space

The radial profiles of the Favre-averaged reactive scalars are shown in Figures4(a)–4(e) at the axial locations x
/
d
=
8, 9, 10, 11, 14, and 26. For T
c
=
1045 K (thick black lines), the results obtained using the two CSDR models are almost identical, except for some minor differences at x
/
d
=
10 and 11. When T
c is decreased to 1030 K (thin black lines), the predictions remain in good agreement with the experiments up to x
/
d
=
10 and the two CSDR models yield very similar results. Farther downstream, the temperature and products mass fractions are underpredicted and the reactants mass fractions are overpredicted, most notably at x
/
d
=
11 and 14. The influence of the CSDR modelling becomes more apparent at these axial locations. Overall, the results remain in reasonable agreement with the experiments. When T
c is increased to 1060 K (grey lines), the CSDR models yield distinct results for x
/
d
≤
14. Using both models, the temperature and products mass fractions are grossly overpredicted and the reactants mass fractions are underpredicted. All three coflow temperatures result in very similar profiles at x
/
d
=
26, irrespective of the CSDR model. As a general observation in Figure 4, the usage of Girimaji’s model always yields lower temperatures and product mass fractions and higher reactant mass fractions when compared to Mortensen’s, which leads to a shorter liftoff height. This is attributed to the fact that Girimaji’s model results in earlier ignition delays and hence shorter ignition kernels for all coflow temperatures. This phenomenon will be discussed shortly in more detail. Before doing so, the liftoff heights obtained from the cases considered in Figure 4 are first presented. The contours of Y
~
OH
(
×
1
0
-
3
) are shown in Figure 5(a) for T
c
=
1030 K, Figure 5(b) for T
c
=
1045 K, and Figure 5(c) for T
c
=
1060 K. In each subfigure, the left and right panes correspond to the results obtained with Mortensen’s and Girimaji’s models, respectively. As shown in Figure 5(a), stabilization takes place at slightly distinct locations due to the very small Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1030 K). Mortensen’s model (left) results in H
/
d
=
11.91 (19.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
11.48 (14.8% in excess of H
exp
⁡
/
d). The same scenario is observed in Figure 5(b) where T
c
=
1045 K. However, in this case, the location of stabilisation is better predicted. Mortensen’s model (left) yields H
/
d
=
10.61 (6.1% in excess of H
exp
⁡
/
d) while Girimaji’s model (right) results in H
/
d
=
10.50 (5% in excess of H
exp
⁡
/
d). The advantage of the inclusion of the effects of inhomogeneity Mortensen’s model is judged to be insignificant for these two coflow temperatures. Conversely, the differences in the stabilisation locations are more substantial for T
c
=
1060 K, as displayed in Figure 5(c). This is due to the much larger Y
~
OH differentials (see the radial profiles in Figure 4(e) for T
c
=
1060 K). Mortensen’s model (left) yields H
/
d
=
7.39 (26.1% below H
exp
⁡
/
d) while Girimaji’s model (right) gives H
/
d
=
6.61 (33.9% below H
exp
⁡
/
d). As for the radial location of stabilisation (the radial distance from the stabilisation point to the centreline normalised d, denoted here as W
/
d), the values obtained using the two CSDR models for T
c
=
1030 K and 1045 K are the same (1.58 and 1.47, resp.). As for T
c
=
1060 K, the two models yield slightly different values (Mortensen’s, 1.26; Girimaji’s, 1.16). It is obvious from the above results that the liftoff height becomes smaller and that flame base becomes narrower as T
c is decreased. This is due to the fact that at lower T
c, the mixture can autoignite in areas closer to the nozzle which are characterised by high scalar dissipation rate levels. This behaviour will be addressed in more detail in Section 4.2.3.Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
) (e), and Y
~
OH
(
×
1
0
3
) using Girimaji’s and Mortensen’s models for the modelling of the CSDR. In (a)–(e): solid lines, Girimaji’s model; dashed lines, Mortensen’s model; symbols, experimental data [3]. Colours: black (thin), T
c
=
1030K; black (thick), T
c
=
1045K; grey, T
c
=
1060K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
): (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K obtained using Mortensen’s model (left pane) and Girimaji’s model (right pane). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively.
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


## 4.2.2.  Stabilisation Mechanism

The numerical indicators developed by Gordon et al. [6] for the distinction between stabilisation by autoignition and premixed flame propagation are employed here to determine the underlying controlling mechanism. For this purpose, the transport budget of the temperature and the history of radical build-up ahead of the stabilisation height are analysed.Budgets in Mixture Fraction Space. Figure 6 shows the transport budget of the steady-state conditional temperature equation, that is the individual contributions of the r.h.s. terms of (9). Only the results obtained using Girimaji’s model are presented here. Mortensen’s model produces similar trends. For each T
c, the budgets are displayed at several axial locations around the corresponding liftoff height. The radial locations are held fixed and set to the W
/
d value obtained with each T
c. As such, the chosen coordinates cover the preflame, liftoff, and postflame regions. Before analysing the balance in these budgets, the axial evolution of the chemical source term, T
CS, is first examined. As shown in Figures 6(a), 6(b), and 6(c), common trends are observed in the evolution of this term for all three coflow temperatures. Below the stabilisation height (first two locations in the preflame regions, Pre-F1 and Pre-F2), as x
/
d increases, T
CS shifts from very lean to less lean mixtures and its amplitude increases dramatically. At liftoff (middle panes), T
CS peaks at η
≈
0.24 with a significantly higher amplitude. Within approximately two to three nozzle diameters downstream the stabilisation height (first locations in the postflame regions, Post-F1) the peak of T
CS increases further and occurs around the stoichiometric mixture fraction. Further downstream (second locations in the postflame regions, Post-F2), the magnitude of T
c decreases substantially from the Post-F1 locations. The reaction zone becomes wider and presents two peaks. The first occurs around the stoichiometric mixture fraction and the second in rich mixtures. The second peak is most likely attributed to the propagation of a rich reaction zone towards stoichiometric mixtures [24]. All of the above trends are consistent with the direct numerical simulation Yoo et al. [47]. As in the evolution of T
CS, common trends are observed in the axial variation of the remaining terms contributing to the r.h.s. of (9). In the preflame region (locations Pre-F1 and Pre-F2), there is a clear balance in lean mixtures between T
CS, the axial convection term, T
C
,
x, and micromixing, T
MM. The axial and radial diffusion terms, T
D
,
x and T
D
,
y, and the radial convection term, T
C
,
y, are nonnegligible but have little contribution to the overall budget. This balance suggests that stabilisation by premixed flame propagation wherein T
D
,
x balances T
C
,
x in the preheat zone with negligible T
CS is not the case. Instead, stabilisation occurs via autoignition as can be seen from the T
CS-T
C
,
x-T
MM balance. At liftoff, the balance shifts to η
≈
0.24. The role of T
C
,
x in balancing T
CS diminishes but this term remains a major contributor to the overall balance, and T
MM prevails. As for the remaining terms, T
D
,
x becomes much weaker while T
D
,
y and T
C
,
y remain important. In the postflame regions (locations Post-F1 and Post-F2), the role of T
c
,
x diminishes further, T
MM becomes more dominant, and T
D
,
x is virtually null. The contributions of T
D
,
y and T
C
,
y are as significant as that of T
C
,
x at Post-F1 and supersede it at Post-F2. As T
CS increases in mixture fraction space, the contributions of T
D
,
y, T
C
,
y, and T
C
,
x remain smaller than that of T
MM, which acts as the major heat sink. In contrast, as T
CS decays, T
MM increases and acts as a source (notice its positive contribution), and both of these terms are counterbalanced by T
D
,
y, T
C
,
y, and T
C
,
x. Thus, beyond the stabilisation height, the flame budgets indicate the structure of a nonpremixed flame, which is largely characterised by a T
CS-T
MM balance. To be noted that the radiative source, T
RS, is negligible at all locations for all coflow temperatures due to the fact that hydrogen in the fuel stream is highly diluted with nitrogen. Therefore, this term can be safely neglected without loss of accuracy.Transport budget of the steady-stateQ
T equation (r.h.s. terms of (9)) at several axial locations around the stabilisation height: (a) T
c
=
1030K (y
/
d
=
1.58), (b) T
c
=
1045K (y
/
d
=
1.47), and (c) T
c
=
1060K (y
/
d
=
1.16). T
RS, radiative source; T
CS, chemical source; T
MM, micromixing; T
C
,
x, axial convection; T
C
,
y, radial convection; T
D
,
x, axial diffusion; T
D
,
y, radial diffusion. The vertical dashed line corresponds to the location of the stoichiometric mixture fraction. All terms are scaled down by a factor of 103 and the units are K/s. The CSDR is modelled using Girimaji’s model.
(a)
T
c
=
1030 K, y
/
d
=
1.58
(b)
T
c
=
1045 K, y
/
d
=
1.47
(c)
T
c
=
1060 K, y
/
d
=
1.16Budgets in Physical Space. The fact that autoignition is the controlling stabilisation mechanism can be further confirmed by analysing the transport budget of Q
T in physical space. To do so, the PDF-weighted integration of the r.h.s. terms of (9) is performed over the mixture fractions space. This yields I
T
RS, I
T
CS, I
T
MM, IT
C
,
x, I
T
C
,
y, IT
D
,
x, and IT
D
,
y (“IT” stands for “Integrated Transport”). The axial profiles of these terms are plotted in Figure 7(a) for Girimaji’s model and in Figure 7(b) for Mortensen’s model at the respective W
/
d locations obtained with each T
c. The parallels between the observation made in Figure 6 and this figure are evident. I
T
CS is mostly balanced by IT
C
,
x and I
T
MM in the preflame regions, with the former being more dominant. The remaining terms are small but nonnegligible. The absence of IT
D
,
x eliminates the possibility of stabilisation by premixed flame propagation and the current balance indicates that autoignition is the controlling mechanism. Right ahead of liftoff, IT
C
,
y and IT
D
,
y emerge and I
T
MM becomes more important. These three terms start to dominate at the expense of IT
C
,
x, which remains the main heat sink. Downstream of the stabilisation height, IT
C
,
x, IT
C
,
y, and IT
D
,
y diminish gradually and I
T
MM becomes more important. Further downstream I
T
CS is mostly counterbalanced by I
T
MM, with smaller contributions IT
C
,
x, IT
C
,
y, and IT
D
,
y, which indicates the structure of a nonpremixed flame.Axial profiles of the Integrated Transport (IT) budget of the steady-stateQ
T equation (integrated r.h.s. terms of (9)): (a) Girimaji’s model and (b) Mortensen’s model. The bottom, middle, and top panes correspond to T
c
=
1030K, T
c
=
1045K, and T
c
=
1060K, respectively. The IT subscripts follow the same notation used in Figure 6. All terms are scaled down by a factor of 106 and the units are K/s.
(a)
Girimaji's model
(b)
Mortensen's modelRadical History ahead of the Stabilisation Height. The transport budgets of the temperature in mixture fraction and physical spaces are both indicative of stabilisation via autoignition. In order to ascertain this conclusion, the history of radical build-up ahead of the reaction zone is now investigated. Figure 8 displays the axial profiles of the normalised Favre-averaged temperature ((
T
~
-
T
~
min
⁡
)
/
(
T
~
max
⁡
-
T
~
min
⁡
)) and mass fractions (Y
~
/
Y
~
max
⁡) of H, O, OH, and H
O
2. The subscripts “min” and “max” denote the minimum and maximum values of the reactive scalars at y
/
d
=
W
/
d. The results are reported for Girimaji’s and Mortensen’s models in Figures 8(a) and 8(b), respectively. Common trends are observed ahead of the stabilisation height for the different combinations of T
c values and CSDR models.(1)
The production of the intermediateH
O
2 initiates upstream of the axial locations where the radicals H, O, and OH emerge.
(2)
H, O, and OH start building up before the consumption of H
O
2 begins. As shown, the mass fractions of these three species start to increase before H
O
2 reaches its peak.
(3)
H
O
2 builds up rapidly prior to the runaway of H, O, and OH, as can be seen from the axial locations of the maximum slopes.Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of stabilisation corresponding to each T
c: (a) Girimaji’s model and (b) Mortensen’s model. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff height, respectively. The circles indicate the locations of the maximum slopes.
(a)
Girimaji's model
(b)
Mortensen's modelTherefore, radicals do not build up simultaneously as in premixed flame propagation, but ratherH
O
2 acts as a precursor to the production of H, O, and OH as in autoignition scenarios [6]. This confirms that the flame is stabilised by autoignition.


## 4.2.3.  Impact of CSDR Modelling on Autoignition

To better understand the autoignition process, standalone zero-dimensional CMC (0DCMC) calculations are performed for eachT
c. In these calculations, the convective and (physical space) diffusive terms in (8) and (9) are turned off and the CSDR is modelled using the AMC [9]. This CSDR model is parametrised by its maximum value χ
∘ and it is given by
(29)
〈
χ
∣
η
〉
=
χ
∘
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
-
1
)
]
2
}
=
χ
∘
G
(
η
)
,
where er
f
-
1 is the inverse error function. The resulting homogeneous CMC equations are identical to those of the unsteady laminar flamelet model (with unity Lewis number) [48]. For a given χ
∘, the solution evolves from the inert mixing solution at t
=
0 up to the moment when autoignition occurs, t
=
τ. Autoignition is declared at the moment when the maximum OH mass fraction reaches 2
×
1
0
-
4 at any point in mixture fraction space, following Stanković and Merci [49]. Figure 9 displays the ignition delay (τ) as a function of χ
∘. These results are obtained using the mechanism of Mueller et al. [40] for T
c
=
1030 K, 1045 K, and 1060 K. It is evident from this figure that, for a given coflow temperature, τ reaches an asymptotic limit as χ
∘ is gradually increased. This means that the occurrence of autoignition becomes impossible beyond a critical value of χ
∘, denoted here by χ
∘
,
c. The approximate values of χ
∘
,
c are 227 s−1, 482 s−1, and 991 s−1 for T
c
=
1030 K, 1045 K, and 1060 K, respectively. Therefore, χ
∘
,
c increases with increasing T
c. In order to link these findings to the flow field of the flame under investigation, χ
∘ is first obtained by integrating (29) weighted by P
~
(
η
) over the mixture fraction space:
(30)
χ
∘
=
∫
0
1
‍
〈
χ
∣
η
〉
P
~
(
η
)
d
η
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
=
χ
~
∫
0
1
‍
G
(
η
)
P
~
(
η
)
d
η
which requires the usage of the readily available quantities χ
~, ξ
~, and ξ
′′
2
~. Figure 10 shows the Contours of χ
∘ computed from (30) for T
c
=
1045 K. The contours for T
c
=
1030 K and 1060 K are virtually the same and are not included here for brevity. The contour of Figure 10 will be used in all of the subsequent analysis irrespective of the value of T
c. The thick contours represent the levels 227 s−1, 482 s−1, and 991 s−1 which correspond to the χ
∘
,
c values determined in Figure 9. Also shown are the ignition kernel locations obtained from the different combinations of CSDR models and coflow temperatures. As expected, for all three coflow temperatures, ignition occurs at locations where χ
∘ is lower that χ
∘
,
c (the circles, squares, and triangles lie outside the 227 s−1, 482 s−1, and 991 s−1 contour levels, resp.). Since χ
~ decays away from the nozzle, autoignition does not take place until χ
∘ (which depends on χ
~ via (30)) falls below χ
∘
,
c. This explains why increasing T
c results in earlier ignition and shorter liftoff from a mixing point of view.Figure 9
Ignition delays as a function of the peak scalar dissipation rate obtained using the mechanism of Mueller et al. [40] for T
c
=
1030, 1045, and 1060 K. The vertical dashed lines represent the asymptotic limit beyond which ignition cannot occur (χ
∘
=
χ
∘
,
c).Figure 10
Contours ofχ
∘ (see (30)). The thick contours correspond to χ
∘
=
χ
∘
,
c for T
c
=
1030 K (227 s−1), 1045 K (482 s−1), and 1060 K (991 s−1). The symbols represent the ignition locations.In another set of 0DCMC calculations, the “most reactive” mixture fraction,η
mr, is computed for each T
c by setting χ
∘ in (29) to 10−20 s−1 (homogeneous reactor). The results are summarised in Table 4. Consistent with the findings of Stanković and Merci [49], η
mr decreases with increasing T
c. The objective behind determining η
mr is to compare the local CSDR values obtained using Mortensen’s and Girimaji’s models at η
=
η
mr, 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G, respectively. As shown in Figure 11, for all given coflow temperatures and radial locations, the axial profiles of 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G reveal that Mortensen’s models generally result in higher scalar dissipation rate levels and therefore yield longer ignition delays and higher liftoff heights. Also displayed in Figure 11 are the critical CSDR values computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr; that is, 〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1
(
2
η
mr
-
1
)
]
2
} (horizontal grey lines). These values are also provided in Table 4 for all coflow temperatures. As expected, 〈
χ
∣
η
mr
〉
c
AMC increases with increasing T
c. Figure 11(a) (T
c
=
1030 K) shows that ignition cannot occur as long as 〈
χ
∣
η
mr
〉
M or 〈
χ
∣
η
mr
〉
G is below 〈
χ
∣
η
mr
〉
c
AMC. Indeed, when Mortensen’s model is used, ignition happens at the first axial and radial locations where 〈
χ
∣
η
mr
〉
M falls below 〈
χ
∣
η
mr
〉
c
AMC. This takes place at x
/
d
=
11.09 and y
/
d
=
1.69 (first pane from the top in Figure 11(a)). A similar scenario is observed when Girimaji’s model is used. Ignition occurs at the first axial and radial locations where 〈
χ
∣
η
mr
〉
G falls below 〈
χ
∣
η
mr
〉
c
AMC, which happens at x
/
d
=
10.81 and y
/
d
=
1.58 (second pane from the top in Figure 11(a)). Similar observation can be made in Figure 11(b) (T
c
=
1045 K) and Figure 11(c) (T
c
=
1060 K). However, in these two cases, although 〈
χ
∣
η
mr
〉
M and 〈
χ
∣
η
mr
〉
G fall quickly below 〈
χ
∣
η
mr
〉
c
AMC at y
/
d locations smaller than the ones where ignition is indicated, ignition occurs at the first axial location with lowest 〈
χ
∣
η
mr
〉. As shown in Figure 11(b) (T
c
=
1045 K), ignition takes place at x
/
d
=
9.97 and y
/
d
=
1.58 using Mortensen’s model and at x
/
d
=
9.69 and y
/
d
=
1.69 using Girimaji’s model. Similarly, in Figure 11(c) (T
c
=
1060 K), ignition happens at x
/
d
=
6.90 and y
/
d
=
1.26 using Mortensen’s model and at x
/
d
=
6.04 and y
/
d
=
1.26 using Girimaji’s model.Table 4
The most reactive mixture fraction,η
mr, and the CSDR values obtained from the AMC parametrised by χ
∘
,
c at η
mr.
T
c [K]
η
mr
〈
χ
∣
η
mr
〉
c
AMC
†
1030
0.0620
21.27
1045
0.0543
36.65
1060
0.0465
59.02
†
〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
erf
⁡
-
1
(
2
η
mr
-
1
)
]
2
}  [s−1].Axial profiles of〈
χ
∣
η
mr
〉 at several radial locations. (a) T
c
=
1030K, (b) T
c
=
1045K, and (c) T
c
=
1060K. Solid black lines, Girimaji’s model; dashed black lines, Mortensen’s model. The vertical black lines indicate the radial locations where ignition occurs. The grey horizontal lines represent the value of the CSDR computed using the AMC model parametrised by χ
∘
,
c at η
=
η
mr (〈
χ
∣
η
mr
〉
c
AMC
=
χ
∘
,
c
exp
⁡
{
-
2
[
er
f
-
1 (
2
η
mr
-
1
)
]
2
}).
(a)
T
c
=
1030 K
(b)
T
c
=
1045 K
(c)
T
c
=
1060 K


## 4.2.4.  Effect of Spurious Sources

Ideally, the integration of the PDF-weighted CMC equations overη-space should yield the unconditional set of equations without any additional source terms [14]. This outcome is guaranteed when fully consistent CSDR models (e.g., the inhomogeneous model of Mortensen) are employed. To be noted that in order for an inhomogeneous CSDR model to be fully consistent, the CV model employed in the derivation should be consistent with the moments of the presumed PDF. This condition is satisfied in Moretnsen’s model due to the usage of the PDF-gradient model. An example of an inhomogeneous yet inconsistent CSDR closure is the one proposed in [50]. The discrepancy arises from the modelling of the CV using the linear velocity model, which is not consistent with the second moment of the presumed β-PDF.When inconsistent CSDR models are employed, the integration of the CMC equations results in spurious (or false) source terms. In this study, this occurs when the CSDR model of Girimaji is used. The spurious source associated with a speciesκ is calculated through [14, 16]:
(31)
S
~
κ
=
1
2
∫
0
1
[
〈
χ
∣
η
〉
G
︸
Inonsistent
CSDR
model
-
〈
χ
∣
η
〉
M
︸
Consistent
CSDR
model
]
∂
2
Q
κ
,
G
∂
η
2
P
~
(
η
)
d
η
,
where Q
κ
,
G is the conditional mass fraction of species κ obtained in the inconsistent realisation (using Girimaji’s model).Figure12 shows the radial profiles of the mean chemical sources (I
T
CS) and spurious sources of O
2, H
2, H
2
O, and OH for T
c
=
1045 K. In comparison to the mean chemical sources, the magnitudes of the spurious sources are small but nonnegligible, in particular within the flame zone (x
/
d
=
14 and 26). The effect of the spurious sources ahead of the predicted stabilisation heights (x
/
d
=
9 and 10) is not sufficiently large to change the nature of the stabilisation mechanism, as the flame is stabilised by autoignition using both Mortensen’s and Girimaji’s models. However, the fact that Girimaji’s model yields a relatively earlier ignition and a smaller liftoff height indicates that the spurious sources are influential. Although this is not the case in the calculations of the current flame, inconsistent CMC implementations may lead erroneous results and conclusions.Radial profiles of the chemical and spurious sources forT
c
=
1045 K: (a) O
2, (b) H
2, (c) H
2
O, and (d) OH. Thin lines, chemical sources; thick lines, spurious sources. The units are s−1.
(a)
Y
~
O
2
(b)
Y
~
H
2
(c)
Y
~
H
2
O
(d)
Y
~
OH


## 4.3.  Effect of Chemical Kinetics

The mechanism of Mueller et al. [40] has been used in all of the results reported so far. The sensitivity of the flame to chemical kinetics is now assessed by testing the performance of the mechanisms of Li et al. [41], Burke et al. [42], and Ó Conaire et al. [43]. The coflow temperature is set to the experimentally reported value of 1045 K.


### 4.3.1. Comparison of the Chemical Kinetic Mechanisms

0DCMC A Priori Analysis. 0DCMC calculations are first performed using the three mechanisms in order to determine the corresponding χ
∘
,
c values. The results are shown in Figure 13, along with those obtained with the Mueller mechanism in Section 4.2.3 for T
c
=
1045K. By inspecting the calculated χ
∘
,
c values and referring to Figure 10, it can be postulated that the Li mechanism would yield the shortest ignition kernel and smallest liftoff height followed by the Mueller mechanism, then the Ó Conaire mechanism, and finally the Burke mechanism. This analysis provides a preliminary qualitative description of the sensitivity of liftoff height to the considered chemical kinetics.Figure 13
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the mechanisms of Mueller et al. [53], Ó Conaire et al. [43], Li et al. [41], and Burke et al. [42] for T
c
=
1045K. The inset shows the critical values of χ
° above which ignition cannot occur.Results in Physical Space. To ascertain these findings, the flame is simulated using the three newly considered mechanisms. The CSDR is closed using Girimaji’s model in all of the following calculations. The results are displayed in Figure 14 at the axial locations x
/
d = 9, 10, 11, 14, and 26, along with those obtained using Mueller mechanism in Section 4.2.1. Unlike the Burke mechanism, the other three kinetic schemes are capable of reproducing the experimental trends. The results of the Mueller and Ó Conaire mechanisms show best agreement with the experimental data. The Mueller mechanism yields better predictions at x
/
d
=
9 and 10, while Ó Conaire’s results in better agreement with the measurements at x
/
d
=
11. The differences between the two are small at x
/
d
=
14 and 26. When the Li mechanism is employed, the temperature rises at early axial locations owing to the ability of the mixture to ignite at high scalar dissipation rate levels (χ
∘
,
c
Li
=
562 s−1). This effect propagates downstream and results in the deviation of the numerical predictions from the experiments. The performance of the Burke mechanisms is unsatisfactory, as the mixture fails to achieve early ignition (χ
∘
,
c
Burke
=
169s−1) and remains inert up to approximately x
/
d
=
14. The contours of Y
~
OH (×
1
0
3) are displayed in Figure 15. As predicted from the 0DCMC calculations, the Li and Burke mechanisms yield the smallest (H
/
d
=
9.73) and largest (H
/
d
=
16.11) liftoff heights, as shown in Figures 15(a) and 15(b), respectively. Also being consistent with the 0DCMC findings, the liftoff heights obtained using the Mueller and Ó Conaire mechanisms fall between those of the Li and Burke mechanisms, as demonstrated in Figures 15(c) and 15(d). However, the Mueller mechanism results in a slightly larger liftoff height (H
/
d
=
10.50) in comparison to the Ó Conaire mechanism (H
/
d
=
10.26), which contradicts with the 0DCMC predictions. Apart from the Burke mechanism, the liftoff heights computed from the remaining mechanisms are similar and are in close agreement with the experimental value H
exp
⁡
/
d
=
10. As for the radial location of liftoff, the Mueller, Li, and Ó Conaire mechanisms result is W
/
d
=
1.47, whereas Burke’s yields a wider flame base with W
/
d
=
1.92. As in the Mueller mechanism calculations, the flame is also stabilised by autoignition when the Ó Conaire, Li, and Burke mechanisms are employed. This can be clearly seen in Figure 16 where H, O, OH, and H
O
2 exhibit the same trends observed in Figure 8(a).Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
) using different chemical kinetic mechanisms. In (a)–(e): thick solid lines, Mueller et al. [40] mechanism; thick dashed lines, Li et al. [41] mechanism; thin solid lines, Burke et al. [42] mechanism; thin dashed lines, Ó Conaire et al. [43] mechanism; Symbols, experimental data [3]. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
) obtained using the mechanisms of (a) Li et al. [41], (b) Burke et al. [42], (c) Mueller et al. [40], and (d) Ó Conaire et al. [43]. The thick contours correspond to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Li
(b)
Burke
(c)
Mueller
(d)
Ó ConaireFigure 16
Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of liftoff using the mechanisms of (from top to bottom) Burke et al. [42], Ó Conaire et al. [43], and Li et al. [41]. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff, respectively. The circles indicate the locations of the maximum slopes. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.Results in Mixture Fraction Space. For completeness, the conditional profiles are shown in Figure 17. The calculation of the conditional data from the experimental scatter is described in the Appendix. The numerical results are reported at x
/
d = 9, 10, 11, 14, and 26, and y
/
d
=
1.47 (W
/
d obtained with all mechanisms except Burke’s). Overall, the results of the Mueller and Ó Conaire mechanisms are in reasonable agreement with the experiments, especially within mixture fraction intervals where the PDF is large (see Figure 17(f)). This explains the good agreement these two mechanisms yield with the measurements in physical space (Figure 14). The Li mechanism predicts the highest OH levels up to x
/
d
=
11 (Figure 17(e)) which leads to the occurrence of the earliest ignition and consequently results in the smallest liftoff height among the considered kinetic schemes. In turn, this early ignition produces higher than expected temperature levels, an effect which propagates downstream resulting in the deviation of the conditional predictions from the experiments. In the three mechanisms discussed up to this point, the peak of Q
T in Figure 17(a) shifts gradually from lean to stoichiometric mixtures as x
/
d increases, which is consistent with the axial evolution of the conditional heat release term T
CS in Figure 6. The profiles of the Burke mechanism remain very close to the inert mixing solution up to x
/
d
=
14. This behaviour results in the most delayed ignition and leads to the largest liftoff height. The conditional profiles at y
/
d
=
1.92 where liftoff occurs undergo very similar trends (not shown here).Conditional profiles atx
/
d
=
9, 10, 11, 14 and 26 and y
/
d
=
1.47 using different chemical kinetic mechanisms: (a) Q
T, (b) Q
O
2, (c) Q
H
2
(
×
1
0
2
), (d) Q
H
2
O
(
×
1
0
2
), (e) Q
OH
(
×
1
0
3
), and (f) P
~
(
η
). In (a)–(e): thick solid lines, Mueller et al. [40]; thick dashed lines, Li et al. [41]; thin solid lines, Burke et al. [42]; thin dashed lines, Ó Conaire et al. [43]; symbols, experimental data (see the Appendix). The vertical dashed line indicates the location of the stoichiometric mixture fraction. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Q
T
[
K
]
(b)
Q
O
2
(c)
Q
H
2
×
1
0
2
(d)
Q
H
2
O
×
1
0
2
(e)
Q
OH
×
1
0
3
(f)
P
~
(
η
)


### 4.3.2.  Diagnosis of the Burke Mechanism

In order to diagnose the unsatisfactory performance of the Burke mechanism, the rate parameters of two of the key reactions are investigated. The recombination reaction(
R
9
)
H
+
O
2
(
+
M
)
⇌
H
O
2
(
+
M
)
and the chain-branching reaction
(
R
1
)
H
+
O
2
⇌
O
+
OH
are two of the most important reactions in H
2 oxidation chemistry. (
R
9
) and (
R
1
) (labelled here as in [40–43]) compete for H atoms and therefore control the overall branching ratio and determine the second explosion limit in homogeneous H
2
/
O
2 systems [42]. The treatment of (
R
9
) and (
R
1
) in the Li and Burke mechanisms differs significantly from the Ó Conaire and Mueller mechanisms. The rate parameters of (
R
9
) are provided in Table 5. The high-pressure limit rate parameters are identical in the Ó Conaire, Mueller, and Li mechanisms [51]. Burke’s mechanism uses an updated set of parameters [52]. In the low-pressure limit, the Ó Conaire and Mueller mechanisms use the same parameters [53] except for the third-body efficiency factors of H
2 and H
2
O. The Li mechanism employs a completely different set of parameters [54]. The Burke parameters differ from Li’s by the broadening factor (decreased by a factor of 1.6 and equal to the ones employed in the Ó Conaire and Mueller mechanisms) and the third-body efficiency factor for H
2
O (increased by a factor of 1.3) [42]. As for (
R
1
), the corresponding rate parameters are listed in Table 6. The parameters employed in the Ó Conaire and Mueller mechanisms are identical [55]. Two distinct sets of updated parameters are employed in the mechanisms of Li [56] and Burke [57]. In light of the above description, the variability of the predictions in Figure 14 and the close agreement between the results of the Mueller and Ó Conaire mechanisms may be explained in part by the treatment of (
R
9
) and (
R
1
). This can be demonstrated, for instance, by adjusting the rate parameters of these two reactions in the Burke mechanism, in an attempt to improve its predictions. Since this mechanism is an updated version of Li’s, two separate experiments are carried out as follows. (1)
The parameters of(
R
9
) are substituted by those of the Li mechanism, with the parameters of  (
R
1
) left unchanged.
(2)
The parameters of both(
R
9
) and (
R
1
) are substituted by those of the Li mechanism.Table 5
Rate parameters of reaction(
R
9
), H + O2 (+M) ⇌ HO2 (+M).
Ó Conaire
Mueller
Li
Burke
High-pressure limit
A
1.475
×
10
12
1.475
×
10
12
1.475
×
10
12
4.65084
×
10
12
n
0.6
0.6
0.6
0.44
E
a
0
0
0
0
Low-pressure limit
A
3.482
×
10
16
3.482
×
10
16
6.366
×
10
20
6.366
×
10
20
n
−0.411
−0.411
−1.72
−1.72
E
a
-
1.115
×
10
3
-
1.115
×
10
3
5.248
×
10
2
5.248
×
10
2
F
c
0.5
0.5
0.8
0.5
ε
H
2
1.3
2.5
2
2
ε
H
2
O
14
12
11
14
ε
O
2
1
1
0.78
0.78
Units are cm3-mole-sec-kcal-K; k
=
A
T
n
exp
⁡
(
-
E
a
/
R
T
).
Nomenclature:A: preexponential coefficient; n: temperature exponent; E
a: activation energy; F
c: broadening factor; ε: third-body efficiency factor.Table 6
Rate parameters of reaction(
R
1
), H + O2
⇌ O + OH.
Ó Conaire
Mueller
Li
Burke
A
1.915
×
10
14
1.915
×
10
14
3.547
×
10
15
1.040
×
10
14
n
0
0
−0.406
0
E
a
1.6439
×
10
4
1.6439
×
10
4
1.6599
×
10
4
1.5286
×
10
4
Same units and nomenclature as in Table5.The preliminary 0DCMC calculations shown in Figure18 reveal that χ
∘
,
c increases significantly when these modifications are applied. Therefore, the mixture is now anticipated to ignite earlier and the liftoff height is expected to decrease. Two additional CMC runs are performed using these modifications while keeping the same settings for the CSDR and the coflow temperature. The results are shown in Figure 19 at x
/
d = 10, 11, and 14. It is evident that the modification of (
R
9
) alone brings a substantial improvement over the rate parameters of the original mechanism. The temperature (Figure 19(a)) and the species mass fractions (Figures 19(b)–19(e)) are much better predicted, especially at x
/
d
=
11 and 14. This further proves the importance of (
R
9
) as a precursor to autoignition reactions. Better agreement with the experimental data is obtained when both (
R
9
) and (
R
1
) are modified. As shown in Figure 19(f), the liftoff height decreases by 26.69% when (
R
9
) is modified and by 29.24% when both (
R
9
) and (
R
1
) are modified. In both cases the computed liftoff heights are substantially closer to the experimentally measured value. Figure 19(f) also shows a remarkable decrease in W
/
d as liftoff occurs closer to the centreline. With the modified settings, both the axial and radial positions of liftoff become in line with those obtained using the Ó Conaire, Mueller, and Li mechanisms. There should be further room for improvement by applying additional modifications to the rate parameters of other key reactions. However, extreme care must be taken while doing so because chemical kinetic mechanisms are usually optimised as a whole, rather than on an individual reaction basis. A sensitivity analysis is beyond the scope of this work and was not performed here. The modifications applied above aread hoc.Figure 18
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the original and modified mechanisms of Burke et al. [42]. The inset shows the critical values of  χ
° above which ignition cannot occur.Radial profiles obtained using the original and modified versions of the Burke mechanism: (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
). In (a)–(e): dashed lines, original mechanism; thin solid lines, modified (
R
9
); thick solid lines, modified (
R
9
) and (
R
1
); symbols, experimental data [3]. Subfigure (f) shows the contours of Y
~
OH
(
×
1
0
3
). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3
(f)
Contours ofY
~
OH
(
×
1
0
3
)
n


## 4.3.1. Comparison of the Chemical Kinetic Mechanisms

0DCMC A Priori Analysis. 0DCMC calculations are first performed using the three mechanisms in order to determine the corresponding χ
∘
,
c values. The results are shown in Figure 13, along with those obtained with the Mueller mechanism in Section 4.2.3 for T
c
=
1045K. By inspecting the calculated χ
∘
,
c values and referring to Figure 10, it can be postulated that the Li mechanism would yield the shortest ignition kernel and smallest liftoff height followed by the Mueller mechanism, then the Ó Conaire mechanism, and finally the Burke mechanism. This analysis provides a preliminary qualitative description of the sensitivity of liftoff height to the considered chemical kinetics.Figure 13
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the mechanisms of Mueller et al. [53], Ó Conaire et al. [43], Li et al. [41], and Burke et al. [42] for T
c
=
1045K. The inset shows the critical values of χ
° above which ignition cannot occur.Results in Physical Space. To ascertain these findings, the flame is simulated using the three newly considered mechanisms. The CSDR is closed using Girimaji’s model in all of the following calculations. The results are displayed in Figure 14 at the axial locations x
/
d = 9, 10, 11, 14, and 26, along with those obtained using Mueller mechanism in Section 4.2.1. Unlike the Burke mechanism, the other three kinetic schemes are capable of reproducing the experimental trends. The results of the Mueller and Ó Conaire mechanisms show best agreement with the experimental data. The Mueller mechanism yields better predictions at x
/
d
=
9 and 10, while Ó Conaire’s results in better agreement with the measurements at x
/
d
=
11. The differences between the two are small at x
/
d
=
14 and 26. When the Li mechanism is employed, the temperature rises at early axial locations owing to the ability of the mixture to ignite at high scalar dissipation rate levels (χ
∘
,
c
Li
=
562 s−1). This effect propagates downstream and results in the deviation of the numerical predictions from the experiments. The performance of the Burke mechanisms is unsatisfactory, as the mixture fails to achieve early ignition (χ
∘
,
c
Burke
=
169s−1) and remains inert up to approximately x
/
d
=
14. The contours of Y
~
OH (×
1
0
3) are displayed in Figure 15. As predicted from the 0DCMC calculations, the Li and Burke mechanisms yield the smallest (H
/
d
=
9.73) and largest (H
/
d
=
16.11) liftoff heights, as shown in Figures 15(a) and 15(b), respectively. Also being consistent with the 0DCMC findings, the liftoff heights obtained using the Mueller and Ó Conaire mechanisms fall between those of the Li and Burke mechanisms, as demonstrated in Figures 15(c) and 15(d). However, the Mueller mechanism results in a slightly larger liftoff height (H
/
d
=
10.50) in comparison to the Ó Conaire mechanism (H
/
d
=
10.26), which contradicts with the 0DCMC predictions. Apart from the Burke mechanism, the liftoff heights computed from the remaining mechanisms are similar and are in close agreement with the experimental value H
exp
⁡
/
d
=
10. As for the radial location of liftoff, the Mueller, Li, and Ó Conaire mechanisms result is W
/
d
=
1.47, whereas Burke’s yields a wider flame base with W
/
d
=
1.92. As in the Mueller mechanism calculations, the flame is also stabilised by autoignition when the Ó Conaire, Li, and Burke mechanisms are employed. This can be clearly seen in Figure 16 where H, O, OH, and H
O
2 exhibit the same trends observed in Figure 8(a).Radial profiles of (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
) using different chemical kinetic mechanisms. In (a)–(e): thick solid lines, Mueller et al. [40] mechanism; thick dashed lines, Li et al. [41] mechanism; thin solid lines, Burke et al. [42] mechanism; thin dashed lines, Ó Conaire et al. [43] mechanism; Symbols, experimental data [3]. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
[
K
]
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3Contours ofY
~
OH
(
×
1
0
3
) obtained using the mechanisms of (a) Li et al. [41], (b) Burke et al. [42], (c) Mueller et al. [40], and (d) Ó Conaire et al. [43]. The thick contours correspond to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Li
(b)
Burke
(c)
Mueller
(d)
Ó ConaireFigure 16
Axial profiles of the normalised Favre-averaged temperature and mass fractions ofH, O, OH, and H
O
2 obtained at the radial locations of liftoff using the mechanisms of (from top to bottom) Burke et al. [42], Ó Conaire et al. [43], and Li et al. [41]. The vertical dash-dotted and dashed lines correspond to the axial locations of maximum H
O
2 and liftoff, respectively. The circles indicate the locations of the maximum slopes. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.Results in Mixture Fraction Space. For completeness, the conditional profiles are shown in Figure 17. The calculation of the conditional data from the experimental scatter is described in the Appendix. The numerical results are reported at x
/
d = 9, 10, 11, 14, and 26, and y
/
d
=
1.47 (W
/
d obtained with all mechanisms except Burke’s). Overall, the results of the Mueller and Ó Conaire mechanisms are in reasonable agreement with the experiments, especially within mixture fraction intervals where the PDF is large (see Figure 17(f)). This explains the good agreement these two mechanisms yield with the measurements in physical space (Figure 14). The Li mechanism predicts the highest OH levels up to x
/
d
=
11 (Figure 17(e)) which leads to the occurrence of the earliest ignition and consequently results in the smallest liftoff height among the considered kinetic schemes. In turn, this early ignition produces higher than expected temperature levels, an effect which propagates downstream resulting in the deviation of the conditional predictions from the experiments. In the three mechanisms discussed up to this point, the peak of Q
T in Figure 17(a) shifts gradually from lean to stoichiometric mixtures as x
/
d increases, which is consistent with the axial evolution of the conditional heat release term T
CS in Figure 6. The profiles of the Burke mechanism remain very close to the inert mixing solution up to x
/
d
=
14. This behaviour results in the most delayed ignition and leads to the largest liftoff height. The conditional profiles at y
/
d
=
1.92 where liftoff occurs undergo very similar trends (not shown here).Conditional profiles atx
/
d
=
9, 10, 11, 14 and 26 and y
/
d
=
1.47 using different chemical kinetic mechanisms: (a) Q
T, (b) Q
O
2, (c) Q
H
2
(
×
1
0
2
), (d) Q
H
2
O
(
×
1
0
2
), (e) Q
OH
(
×
1
0
3
), and (f) P
~
(
η
). In (a)–(e): thick solid lines, Mueller et al. [40]; thick dashed lines, Li et al. [41]; thin solid lines, Burke et al. [42]; thin dashed lines, Ó Conaire et al. [43]; symbols, experimental data (see the Appendix). The vertical dashed line indicates the location of the stoichiometric mixture fraction. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
Q
T
[
K
]
(b)
Q
O
2
(c)
Q
H
2
×
1
0
2
(d)
Q
H
2
O
×
1
0
2
(e)
Q
OH
×
1
0
3
(f)
P
~
(
η
)


## 4.3.2.  Diagnosis of the Burke Mechanism

In order to diagnose the unsatisfactory performance of the Burke mechanism, the rate parameters of two of the key reactions are investigated. The recombination reaction(
R
9
)
H
+
O
2
(
+
M
)
⇌
H
O
2
(
+
M
)
and the chain-branching reaction
(
R
1
)
H
+
O
2
⇌
O
+
OH
are two of the most important reactions in H
2 oxidation chemistry. (
R
9
) and (
R
1
) (labelled here as in [40–43]) compete for H atoms and therefore control the overall branching ratio and determine the second explosion limit in homogeneous H
2
/
O
2 systems [42]. The treatment of (
R
9
) and (
R
1
) in the Li and Burke mechanisms differs significantly from the Ó Conaire and Mueller mechanisms. The rate parameters of (
R
9
) are provided in Table 5. The high-pressure limit rate parameters are identical in the Ó Conaire, Mueller, and Li mechanisms [51]. Burke’s mechanism uses an updated set of parameters [52]. In the low-pressure limit, the Ó Conaire and Mueller mechanisms use the same parameters [53] except for the third-body efficiency factors of H
2 and H
2
O. The Li mechanism employs a completely different set of parameters [54]. The Burke parameters differ from Li’s by the broadening factor (decreased by a factor of 1.6 and equal to the ones employed in the Ó Conaire and Mueller mechanisms) and the third-body efficiency factor for H
2
O (increased by a factor of 1.3) [42]. As for (
R
1
), the corresponding rate parameters are listed in Table 6. The parameters employed in the Ó Conaire and Mueller mechanisms are identical [55]. Two distinct sets of updated parameters are employed in the mechanisms of Li [56] and Burke [57]. In light of the above description, the variability of the predictions in Figure 14 and the close agreement between the results of the Mueller and Ó Conaire mechanisms may be explained in part by the treatment of (
R
9
) and (
R
1
). This can be demonstrated, for instance, by adjusting the rate parameters of these two reactions in the Burke mechanism, in an attempt to improve its predictions. Since this mechanism is an updated version of Li’s, two separate experiments are carried out as follows. (1)
The parameters of(
R
9
) are substituted by those of the Li mechanism, with the parameters of  (
R
1
) left unchanged.
(2)
The parameters of both(
R
9
) and (
R
1
) are substituted by those of the Li mechanism.Table 5
Rate parameters of reaction(
R
9
), H + O2 (+M) ⇌ HO2 (+M).
Ó Conaire
Mueller
Li
Burke
High-pressure limit
A
1.475
×
10
12
1.475
×
10
12
1.475
×
10
12
4.65084
×
10
12
n
0.6
0.6
0.6
0.44
E
a
0
0
0
0
Low-pressure limit
A
3.482
×
10
16
3.482
×
10
16
6.366
×
10
20
6.366
×
10
20
n
−0.411
−0.411
−1.72
−1.72
E
a
-
1.115
×
10
3
-
1.115
×
10
3
5.248
×
10
2
5.248
×
10
2
F
c
0.5
0.5
0.8
0.5
ε
H
2
1.3
2.5
2
2
ε
H
2
O
14
12
11
14
ε
O
2
1
1
0.78
0.78
Units are cm3-mole-sec-kcal-K; k
=
A
T
n
exp
⁡
(
-
E
a
/
R
T
).
Nomenclature:A: preexponential coefficient; n: temperature exponent; E
a: activation energy; F
c: broadening factor; ε: third-body efficiency factor.Table 6
Rate parameters of reaction(
R
1
), H + O2
⇌ O + OH.
Ó Conaire
Mueller
Li
Burke
A
1.915
×
10
14
1.915
×
10
14
3.547
×
10
15
1.040
×
10
14
n
0
0
−0.406
0
E
a
1.6439
×
10
4
1.6439
×
10
4
1.6599
×
10
4
1.5286
×
10
4
Same units and nomenclature as in Table5.The preliminary 0DCMC calculations shown in Figure18 reveal that χ
∘
,
c increases significantly when these modifications are applied. Therefore, the mixture is now anticipated to ignite earlier and the liftoff height is expected to decrease. Two additional CMC runs are performed using these modifications while keeping the same settings for the CSDR and the coflow temperature. The results are shown in Figure 19 at x
/
d = 10, 11, and 14. It is evident that the modification of (
R
9
) alone brings a substantial improvement over the rate parameters of the original mechanism. The temperature (Figure 19(a)) and the species mass fractions (Figures 19(b)–19(e)) are much better predicted, especially at x
/
d
=
11 and 14. This further proves the importance of (
R
9
) as a precursor to autoignition reactions. Better agreement with the experimental data is obtained when both (
R
9
) and (
R
1
) are modified. As shown in Figure 19(f), the liftoff height decreases by 26.69% when (
R
9
) is modified and by 29.24% when both (
R
9
) and (
R
1
) are modified. In both cases the computed liftoff heights are substantially closer to the experimentally measured value. Figure 19(f) also shows a remarkable decrease in W
/
d as liftoff occurs closer to the centreline. With the modified settings, both the axial and radial positions of liftoff become in line with those obtained using the Ó Conaire, Mueller, and Li mechanisms. There should be further room for improvement by applying additional modifications to the rate parameters of other key reactions. However, extreme care must be taken while doing so because chemical kinetic mechanisms are usually optimised as a whole, rather than on an individual reaction basis. A sensitivity analysis is beyond the scope of this work and was not performed here. The modifications applied above aread hoc.Figure 18
Ignition delay as a function of the peak scalar dissipation rateχ
° obtained using the original and modified mechanisms of Burke et al. [42]. The inset shows the critical values of  χ
° above which ignition cannot occur.Radial profiles obtained using the original and modified versions of the Burke mechanism: (a)T
~, (b) Y
~
O
2, (c) Y
~
H
2
(
×
1
0
2
), (d) Y
~
H
2
O
(
×
1
0
2
), and (e) Y
~
OH
(
×
1
0
3
). In (a)–(e): dashed lines, original mechanism; thin solid lines, modified (
R
9
); thick solid lines, modified (
R
9
) and (
R
1
); symbols, experimental data [3]. Subfigure (f) shows the contours of Y
~
OH
(
×
1
0
3
). The thick contour corresponds to 600 ppm. The solid and dashed horizontal lines correspond to the numerical and experimental liftoff heights, respectively. The CSDR is modelled using Girimaji’s model and T
c
=
1045K.
(a)
T
~
(b)
Y
~
O
2
(c)
Y
~
H
2
×
1
0
2
(d)
Y
~
H
2
O
×
1
0
2
(e)
Y
~
OH
×
1
0
3
(f)
Contours ofY
~
OH
(
×
1
0
3
)
n


## 5. Conclusions

A liftedH
2
/
N
2 jet flame issuing into a vitiated coflow was investigated using the first-order CMC. The flow and mixing fields were obtained using a modified version of the k
-
ɛ model. The calculations included an additional transport equation for the Favre-averaged scalar dissipation rate. Two formulations were implemented for the CSDR, the homogeneous model of Girimaji, and the inhomogeneous model of Mortensen. The CV was modelled using the PDF-gradient model. Calculations were performed for three coflow temperatures, 1030, 1045, and 1060 K, and four chemical kinetic mechanisms were assessed. In light of the previous results, the following conclusions are drawn.(1)
The modification of the constantC
ɛ
1 in the k
-
ɛ model from the standard value of 1.42 to 1.6 improves the predictions of the spreading rate of the jet.
(2)
The solution of a transport equation for the scalar dissipation rate in lieu of using the traditional algebraic modelling approach provides a more reliable mixing field.
(3)
The flame is very sensitive to small changes in the coflow temperature. A variation of roughly±
1.44
% in the experimentally reported value of 1045 K results in substantial changes in the predictions. T
c
=
1045K yields best agreement with the experiments. The calculated liftoff height at the 1030 K level remains close to the experimentally measured value. A drastic decrease occurs when the coflow temperature is increased to 1060 K.
(4)
For all coflow temperatures and using both CSDR models, the transport budgets in mixture fraction and physical spaces and the history of radical build-up ahead of the reaction zone indicate that the flame is stabilised by autoignition rather than premixed flame propagation. These findings fully agree with previous PDF calculations [4–6], with the CMC results of Stankovic and Merci [7], and in part with those of Navarro-Martinez and Kronenburg [8] and Patwardhan et al. [11].
(5)
Standalone zero-dimensional CMC calculation indicates that the mixture is capable of igniting at higher scalar dissipation levels as the coflow temperature is increased. This provides an explanation for the occurrence of ignition at locations closer to the nozzle exit and the decrease in liftoff height with increasing coflow temperature.
(6)
In comparison to Mortensen’s model, Girimaji’s model always results in lower CSDR at the “most reactive” mixture fraction and therefore results in earlier ignition and smaller liftoff heights for all coflow temperatures.
(7)
The spurious source terms resulting from the modelling of the CSDR are in general small but nonnegligible, mostly notably within the flame zone.
(8)
The flame shows high sensitivity to chemical kinetics. The results obtained with the Mueller, Li, and Ó Conaire mechanisms reproduce the experimental trends with varying degrees of accuracy. The predictions are in general in good agreement with the experimental measurements, particularly those of the Mueller and Ó Conaire mechanisms. The three kinetic schemes predict very similar liftoff heights. Conversely, the performance of the Burke mechanism is unsatisfactory due to the inability of the mixture to achieve early ignition. The temperature and species mass fractions are not well predicted and the liftoff height is grossly overestimated. The modification of the rate parameters of some of the key reactions improves the predictions substantially.Further improvements are possible by modelling the conditional chemical sources with a second-order closure and by inclusion of counter-gradient and differential diffusion effects.



---

*Source: 507459-2014-05-19.xml*