DNS of inhomogeneous reactants
premixed combustion
Kian Min Lim
Fitzwilliam College
University of Cambridge
A thesis submitted for the degree of
Doctor of Philosophy
September 2014
Declaration
This dissertation is the result of my own work and includes noth-
ing which is the outcome of work done in collaboration except where
specifically indicated in the text. This dissertation has not been sub-
mitted in any form for another qualification to this or any other uni-
versity. Information derived from the published and unpublished work
of others has been referenced in the text.
This dissertation contains approximately 58,000 words1 and 94 fig-
ures, meeting the requirements set by the Degree Committee for the
Department of Engineering.
Kian Min Lim
30th September 2014
1Using LibreOffice
To Chai Ling and my parents.
Acknowledgements
First and foremost, I would like to thank my supervisor, Professor
Stewart Cant for his constant support, advice and guidance during
my studies. His patience, continuous encouragement and confidence
in me kept me motivated throughout my studies.
I am grateful to my advisor, Professor Mastorakos for his useful dis-
cussions and Professor Nilanjan Chakraborty who has provided useful
insights and discussions on analysing the DNS data. Dr. Tom Dun-
stan has in numerous occasions helped me in DNS coding as well as
sharing his insight on managing HECToR.
The generous financial support from the Malaysian government through
the Biasiswa Yang Di-Pertuan Agong scholarship has made the pur-
sue of my research possible. Computational resources from the UK
supercomputer HECToR has been provided by the EPSRC and it is
gratefully acknowledged.
I am also grateful to my tutor Dr. John Cleaver and the Fitzwilliam
graduate secretary, Sue Free for helping me in numerous administra-
tive matters.
Many thanks to Ryan Griffiths and Girish Nivarti for their help in
visualisation using MATLAB and Python. Thanks is also due to my
friends and colleagues from the CFD laboratory, namely Dr. Peter
Cocks, Dr. Adam Comer, and Dr. Mohd. Fairuz, for the useful
discussions and the good time we had together.
Last but not least, I like to thank my family for their continuous
support and encouragement that has helped kept me motivated in
pursuing my passion.
DNS of inhomogeneous reactants premixed
combustion
Kian Min Lim
The search for clean and efficient combustors is motivated by the increasingly
stringent emissions regulations. New gas turbine engines are designed to operate
under lean conditions with inhomogeneous reactants to ensure cleanliness and
stability of the combustion. This ushers in a new mode of combustion, called the
inhomogeneous reactants premixed combustion.
The present study investigates the effects of inhomogeneous reactants on pre-
mixed combustion, specifically on the interactions of an initially planar flame with
a field of inhomogeneous reactants. Unsteady and unstrained laminar methane-
air flames are studied in one- and two-dimensional simulations to investigate the
effects of normally and tangentially (to the flame surface) stratified reactants. A
three-dimensional DNS of turbulent inhomogeneous reactants premixed combus-
tion is performed to extend the investigation into turbulent flames. The methane-
air combustion is represented by a complex chemical reaction mechanism with 18
species and 68 steps.
The flame surface density (FSD) and displacement speed Sd have been used
as the framework to analyse the inhomogeneous reactants premixed flame. The
flames are characterised by an isosurface of reaction progress variable. The un-
steady flames are compared to the steady laminar unstrained reference case. An
equivalence ratio dip is observed in all simulations and it can serve as a marker for
the premixed flame. The dip is attributed to the preferential diffusion of carbon-
and hydrogen- containing species.
Hysteresis of Sd is observed in the unsteady and unstrained laminar flames
that propagate into normally stratified reactants. Stoichiometric flames prop-
agating into lean mixture have a larger Sd than lean flames propagating into
stoichiometric mixtures. The cross-dissipation term contribution to Sd is small
(≈ 10%) but its contribution to the hysteresis of Sd is not (≈ 50%).
Differential propagation of the flame surface is observed in the laminar flame
that propagates into tangentially stratified reactants. Stretch on the flame sur-
face is induced by the differential propagation, which in turn increases the flame
surface area.
In the DNS, inhomogeneous reactants were observed to have interacted with
the turbulent premixed flame in the preheat zone. Results from DNS show similar
FSD and Sd statistics as previous thin reaction zones regime DNS of turbulent
premixed flames. Effects of inhomogeneous reactants as observed in the laminar
flame simulations are not pervalent in the turbulent flames.
Keywords: Inhomogeneous reactants premixed combustion, Premixed combus-
tion, Direct Numerical Simulation (DNS), Progress variable, Flame surface den-
sity (FSD), Displacement speed.
Contents
Contents vi
List of Figures xi
List of Tables xix
1 Introduction 1
1.1 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Summary of thesis structure . . . . . . . . . . . . . . . . . . . . . 4
2 Background 6
2.1 Premixed combustion . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 DNS in turbulent premixed combustion . . . . . . . . . . . . . . . 9
2.3 Composition variables for the inhomogeneous reactants premixed
combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.1 Mixture fraction . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.2 Progress variable . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.3 Equivalence ratio . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Numerical studies of inhomogeneous reactants premixed combus-
tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4.1 Multidimensional turbulent inhomogeneous reactants pre-
mixed combustion numerical simulations . . . . . . . . . . 12
2.4.2 Reported results for DNS of turbulent inhomogeneous re-
actants premixed combustion . . . . . . . . . . . . . . . . 14
vi
CONTENTS
2.4.3 One-dimensional laminar inhomogeneous reactants premixed
combustion numerical simulations . . . . . . . . . . . . . . 17
2.4.4 Reported results for one-dimensional simulations of laminar
inhomogeneous reactants premixed combustion . . . . . . 17
2.5 Experiments on inhomogeneous reactants premixed combustion . 19
2.6 Research approach . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 Mathematical Background and Numerical Implementations 24
3.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1.1 Thermochemical properties . . . . . . . . . . . . . . . . . 26
3.1.2 Transport properties . . . . . . . . . . . . . . . . . . . . . 27
3.1.3 Diffusion correction velocity . . . . . . . . . . . . . . . . . 27
3.2 Spatial discretisation and temporal advancement schemes . . . . . 28
3.3 Determination of the chemical time-step and grid resolution re-
quirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Turbulence Initialisation . . . . . . . . . . . . . . . . . . . . . . . 31
3.4.1 Initial turbulent velocity flow field in a rectangular domain 32
3.5 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.5.1 Locally One-Dimensional Inviscid (LODI) relations . . . . 34
3.5.2 Inflow Boundary Condition . . . . . . . . . . . . . . . . . 35
3.5.3 Outflow Boundary Conditions . . . . . . . . . . . . . . . . 36
3.6 The progress variable . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.6.1 The element mixture fraction and the local equivalence ratio 37
3.7 The progress variable transport equations . . . . . . . . . . . . . . 38
3.7.1 Mixture fraction transport equation . . . . . . . . . . . . . 40
3.8 Displacement speed components . . . . . . . . . . . . . . . . . . . 40
3.9 Surface density function and fine-grained flame surface density . . 43
4 One-dimensional laminar flame simulations 45
4.1 One-dimensional steady laminar flame calculations at different equiv-
alence ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2 Equivalence ratio profile across laminar flame . . . . . . . . . . . 48
vii
CONTENTS
4.3 The properties of progress variable . . . . . . . . . . . . . . . . . 53
4.4 Unsteady One-dimensional Laminar Flame Simulations . . . . . . 55
4.4.1 Unsteady laminar flame with sinusoidal inlet equivalence
ratio time history . . . . . . . . . . . . . . . . . . . . . . . 56
4.4.2 Unsteady laminar flame simulation with complimentary er-
ror function inlet equivalence ratio time history . . . . . . 58
4.5 Temporal analysis of case S00 displacement speed . . . . . . . . . 60
4.5.1 Displacement speed Sd analysis for case S00 . . . . . . . . 62
4.5.2 Displacement speed deviation ∆Sd analysis for case S00 . . 64
4.6 Temporal analysis of R cases displacement speed . . . . . . . . . . 67
4.7 Displacement speed deviation ∆Sd analysis for cases R00, R01,
R02, and R03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.7.1 The contributors to the deviation of normal displacement
speed component Sn . . . . . . . . . . . . . . . . . . . . . 73
4.7.2 The contributors to the deviation of reaction displacement
speed component Sr . . . . . . . . . . . . . . . . . . . . . 74
4.7.3 The contributors to the deviation of cross dissipation dis-
placement speed component Scd . . . . . . . . . . . . . . . 76
4.7.4 Dominant terms that influence Sn, Sr, and Scd deviation . 78
4.7.5 Total heat release rate Q˙hr . . . . . . . . . . . . . . . . . . 85
4.7.6 Heat flux on the c = 0.90 iso-surface q˙ . . . . . . . . . . . 86
4.7.7 Species that influence the deviation of carbon and hydrogen
elements (zC , zH) and their gradients. . . . . . . . . . . . . 88
4.7.8 Relevant species mass fractions deviations and diffusive mass
flux deviations . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.8 The displacement speed deviation ∆Sd analysis for cases R01,
R01b, R01c, and R01d . . . . . . . . . . . . . . . . . . . . . . . . 93
4.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5 Two-dimensional laminar flame simulations 99
5.1 Two-dimensional inhomogeneous laminar flame simulation config-
urations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.1.1 Test matrix of two-dimensional simulations . . . . . . . . . 102
viii
CONTENTS
5.1.2 Initial solution for the simulation . . . . . . . . . . . . . . 103
5.2 The methane species mass fraction transport budget . . . . . . . . 105
5.2.1 Equivalence ratio and progress variable profile for Case A00 107
5.3 Case B00 flame evolution . . . . . . . . . . . . . . . . . . . . . . . 109
5.3.1 Time evolution of the equivalence ratio (Φ0.90) profile . . . 111
5.3.2 Evolution of the surface-weighted averaged of equivalence
ratio (Φ0.90) profile . . . . . . . . . . . . . . . . . . . . . . 112
5.3.3 Time evolution of the displacement speed profile . . . . . . 113
5.3.4 Time evolution of the displacement speed components . . . 114
5.4 The two-dimensional flame evolution stages . . . . . . . . . . . . 118
5.4.1 Normal stratification stage . . . . . . . . . . . . . . . . . . 119
5.4.2 Flame stretching stage . . . . . . . . . . . . . . . . . . . . 123
5.4.3 Cusp formation stage . . . . . . . . . . . . . . . . . . . . . 127
5.5 SDF scatter against progress variable across the flame brush for
case B00 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.6 The budget of the progress variable transport equations . . . . . . 129
5.7 Scatter plot of displacement speed and its components coloured by
local equivalence ratio . . . . . . . . . . . . . . . . . . . . . . . . 134
5.8 The surface density function transport equation . . . . . . . . . . 136
5.9 Scatter of the surface density function transport equation across
the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
5.10 Mixture fraction profile along the c=0.90 isosurface . . . . . . . . 141
5.11 Total flame surface area and the flame surface area weighted aver-
aged local equivalence ratio . . . . . . . . . . . . . . . . . . . . . 148
5.12 Flame surface area weighted averaged displacement speed and its
components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5.13 Global flame response . . . . . . . . . . . . . . . . . . . . . . . . . 153
5.14 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
6 Direct numerical simulations of turbulent premixed flame with
inhomogeneous reactants 159
6.1 Initial solution for the DNS . . . . . . . . . . . . . . . . . . . . . 159
6.2 Contour plots of equivalence ratio across the flame . . . . . . . . . 164
ix
CONTENTS
6.2.1 Vorticity colour contour plot for t/tL = 0.87 . . . . . . . . 166
6.3 Reactants segregation across the flame brush . . . . . . . . . . . . 169
6.4 Scatter plot of the SDF across the flame brush . . . . . . . . . . . 171
6.5 Mean behaviour of the source terms of the SDF transport equation 172
6.6 Scatter plot of displacement speed and its components against
flame brush . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
6.7 Strain rate and curvature statistics . . . . . . . . . . . . . . . . . 175
6.7.1 Joint PDF of curvature and strain rate . . . . . . . . . . . 179
6.8 Strain rate and curvature effect on SDF . . . . . . . . . . . . . . . 181
6.9 Mean behaviour of displacement speed across the flame brush . . 183
6.10 Scatter plot of displacement speed and its components against
flame brush . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
6.11 PDFs of displacement speed and its components . . . . . . . . . . 188
6.12 Mean variation of displacement speed against curvature on differ-
ent iso-surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
6.13 Mean variation of displacement speed against tangential strain rate
on different iso-surfaces. . . . . . . . . . . . . . . . . . . . . . . . 193
6.14 The local equivalence ratio profile across the turbulent flame brush 195
6.14.1 Mixture fraction scatter plot across the flame brush . . . . 197
6.14.2 Species that contributes to the element mixture fraction
profile across the flame brush . . . . . . . . . . . . . . . . 201
6.15 Selected species scatter plot . . . . . . . . . . . . . . . . . . . . . 204
6.16 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
7 Conclusions 210
7.1 Primary conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . 210
7.2 Future Investigations . . . . . . . . . . . . . . . . . . . . . . . . . 213
A Derivation of the balance equation for the progress variable 216
References 218
x
List of Figures
1.1 Venn diagram depicting the combustion modes. The inhomoge-
neous reactants premixed combustion is a subset of premixed com-
bustion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Modified Borghi diagram . . . . . . . . . . . . . . . . . . . . . . . 7
3.1 The laminar flame species chemical time tc,α and length lc,α scales
for flame with lean to unity global equivalence ratio. . . . . . . . . 30
3.2 Turbulent velocity flow field initialisation in a rectangular domain.The
contour plot is coloured by vorticity magnitude. . . . . . . . . . . 33
4.1 Comparison between calculated laminar flame speed and experi-
mentally measured values by Vagelopoulos et. al. [112]. The devi-
ations are plotted in green. . . . . . . . . . . . . . . . . . . . . . . 48
4.2 Profile of local equivalence ratio normalised by the global equiva-
lence ratio. Spatial variable is normalised by laminar flame thickness. 49
4.3 Normalised carbon element mixture fraction. . . . . . . . . . . . . 50
4.4 Normalised hydrogen element mixture fraction. . . . . . . . . . . 50
4.5 Comparison of the local equivalence ratio measured at c = 0.90
with the global equivalence ratio. . . . . . . . . . . . . . . . . . . 52
4.6 Three instantaneous local equivalence ratio profiles of the unsteady
flame calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.7 Position of different flame properties plotted in the progress vari-
able space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.8 Local instantaneous equivalence ratio case S00. . . . . . . . . . . 55
xi
LIST OF FIGURES
4.9 Local instantaneous equivalence ratio for cases R00, R01a, R02,
and R03. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.10 Sd and Φ0.90 against elapsed time te for case S00. . . . . . . . . . 60
4.11 Displacement speed and its components for case S00 . . . . . . . . 63
4.12 Deviations of the displacement speed and its components from the
steady state counterpart for case S00 . . . . . . . . . . . . . . . . 65
4.13 The time histories of Sd and Φ0.90 for cases R00, R01, R02, and R03 67
4.14 The time histories of Sd and Φ0.90 for cases R01, R01b, R01c, and
R01d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.15 Displacement speed and its components for cases R00, R01, R02,
and R03. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.16 Deviations of the displacement speed and its components from the
steady state counterpart for cases R00, R01, R02, and R03 . . . . 72
4.17 The contributions towards the deviation of normal component of
the displacement speed Sn. . . . . . . . . . . . . . . . . . . . . . . 75
4.18 The contribution of relevant terms towards the deviation of the
reaction rate component of the displacement speed Sr. . . . . . . 76
4.19 The contribution of relevant terms towards the deviation of the
cross dissipation component of the displacement speed Scd. . . . . 77
4.20 Density ρ, surface density function |∇c|, and surface density func-
tion spatial gradient ∂|∇c|/∂x profiles and their deviations for
cases R00, R01, R02, and R03. . . . . . . . . . . . . . . . . . . . . 79
4.21 The profiles and deviations of carbon and hydrogen element mix-
ture fractions (zC , zH). . . . . . . . . . . . . . . . . . . . . . . . . 82
4.22 The spatial gradient and the deviation of the gradient of carbon
and hydrogen element mixture fractions (∂zC/∂x, ∂zH/∂x) plotted
against local equivalence ratio Φ0.90 . . . . . . . . . . . . . . . . . 83
4.23 Methane reaction rate ω˙CH4 . . . . . . . . . . . . . . . . . . . . . 85
4.24 Heat release rate Q˙hr and the deviation from its reference values
∆Q˙hr/Q˙refhr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.25 The profiles and deviations of total heat flux q˙, heat conduction
term q˙T , and species enthalpy term q˙s for cases R00, R01, R02,
and R03. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
xii
LIST OF FIGURES
4.26 The deviation of species contribution towards carbon and hydrogen
elements for case R01. . . . . . . . . . . . . . . . . . . . . . . . . 89
4.27 The deviation of species contribution towards carbon and hydrogen
elements spatial gradient for case R01. . . . . . . . . . . . . . . . 91
4.28 The displacement speed Sd against local equivalence ratio Φ0.90 for
cases R01, R01b, R01c, and R01d. . . . . . . . . . . . . . . . . . . 94
4.29 Displacement speed and its components for cases R00, R01, R02,
and R03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.1 Equivalence ratio Φ specified at the inlet boundary. . . . . . . . . 101
5.2 The initial equivalence ratio profile for case B00. Φ = 0.50 in blue
to Φ = 1.00 in red. . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3 Methane species mass fraction . . . . . . . . . . . . . . . . . . . . 106
5.4 Field plots of Φ and c for case A00 at t/tL = 2.2. . . . . . . . . . 108
5.5 Progress variable location (c = 0.90) for case B00 . . . . . . . . . 110
5.6 Normalised total flame surface area (AΣ/AL) at c = 0.90 for case
B00. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.7 Local equivalence ratio on the c = 0.90 iso-surface for case B00 at
different elapsed times. . . . . . . . . . . . . . . . . . . . . . . . . 112
5.8 The surface averaged local equivalence ratio Φ0.90 against t/tL for
case B00 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5.9 Displacement speed on the c = 0.90 iso-surface for case B00 at
different t/tL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.10 The snap-shots of displacement speed and its components for case
B00 at different elapsed times. . . . . . . . . . . . . . . . . . . . . 115
5.11 Normalised difference of displacement speed ∆Sd/Srefd on the c =
0.90 iso-surface for case B00 at different t/tL. . . . . . . . . . . . 117
5.12 The displacement speed Sd against local equivalence ratio Φ for
case B00 at t/tL = 0.50 to 2.50 with an increment of 0.50. The
black dashed lin is the steady laminar flame reference displacement
speed Srefd as a function of the local equivalence ratio on c = 0.90
isosurface Φ0.90. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
xiii
LIST OF FIGURES
5.13 The displacement speed Sd against local equivalence ratio Φ for
case B00 at t/tL = 0.50 and 1.00. The black solid line is the steady
laminar flame reference displacement speed Srefd as a function of
the local equivalence ratio on c = 0.90 iso-surface Φ0.90. . . . . . . 120
5.14 The normalised Sd deviation (∆r,n(Sd)) against the normalised
∂Φ/∂n deviation (∆r,n(∂nΦ)) for Case B00 at t/tL = 0.5 and 1.0. . 122
5.15 The normalised flame curvature κ×δL against transverse direction
y/Ly for case B00 at five different instances. . . . . . . . . . . . . 123
5.16 The surface averaged rate of flame surface stretch ∂(lnA)/∂t against
t/tL for case B00. . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.17 The surface averaged strain rate aT against t/tL for case B00. . . 126
5.18 The surface average of the curvature component of the surface
stretch 2κSd against t/tL for case B00. . . . . . . . . . . . . . . . 127
5.19 Scatter plot of SDF σ = |∇c| across the flame brush coloured
by local equivalence ratio Φ. The mean of SDF conditioned on
progress variable is plotted in red asterisks. . . . . . . . . . . . . . 128
5.20 Budget for the r.h.s of progress variable transport equation for case
B00 in the flame stretching stage t/tL = 1.00 (1 of 2). Note that
the SDF is σ = |∇c| (refer to equation 3.46). . . . . . . . . . . . . 131
5.21 Budget for the r.h.s of progress variable transport equation for case
B00 in the flame stretching stage t/tL = 1.00 (2 of 2). . . . . . . . 132
5.22 Scatter plot of r.h.s progress variable transport equation terms
against progress variable. . . . . . . . . . . . . . . . . . . . . . . . 135
5.23 Surface density function transport equations budget . . . . . . . . 137
5.24 Scatter plot of r.h.s SDF transport equation terms against progress
variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
5.25 Carbon zC and hydrogen zH element mixture fractions profile for
case B00 at t/tL = 1.50 (flame stretching stage) plotted against
transverse direction y/Ly. . . . . . . . . . . . . . . . . . . . . . . 141
5.26 The species contribution towards carbon element mixture fraction
for case B00 at t/tL = 1.50 on the c=0.90 isosurface. . . . . . . . 143
5.27 The species contribution towards hydrogen element mixture frac-
tion for case B00 at t/tL = 1.50 on the c=0.90 isosurface. . . . . 143
xiv
LIST OF FIGURES
5.28 The normalised difference of carbon and hydrogen element mixture
fractions for case B00 at t/tL = 1.50 on c=0.90 isosurface. . . . . 144
5.29 The species contributions towards the deviations of carbon element
mixture fraction on c=0.90 isosurface for case B00 at t/tL = 1.50. 145
5.30 The species contributions towards the deviations of hydrogen ele-
ment mixture fraction on c=0.90 isosurface for case B00 at t/tL =
1.50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.31 The deviation of H2, OH, H, H2O2, HCO, and CH2OH from lami-
nar flame values with corresponding local equivalence ratio on the
c=0.90 isosurface for case B00 at t/tL = 1.50 across the transverse
direction y/Ly. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.32 The deviation of temperature from laminar flame values with cor-
responding local equivalence ratio on the c=0.90 isosurface for case
B00 at t/tL = 1.50 across the transverse direction y/Ly. . . . . . . 147
5.33 Normalised total flame surface are AΣ/AL on c=0.90 isosurface for
all six cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
5.34 The surface area weighted averaged local equivalence ratio Φ0.90
against t/tL for all six cases. . . . . . . . . . . . . . . . . . . . . . 149
5.35 Area-averaged displacement speed components for A and B cases,
plotted against elapsed time. . . . . . . . . . . . . . . . . . . . . . 151
5.36 Global flame response plotted against t/tL for all A and B cases. . 155
6.1 Initial local equivalence ratio of the DNS on x−y plane (z/Lz = 0.50).160
6.2 DNS configuration in the regime diagram . . . . . . . . . . . . . . 162
6.3 Initial temperature contour of the DNS on x − y plane (z/Lz =
0.50), combined with the initial velocity vector field. . . . . . . . . 163
6.4 The colour legend varies from blue Φ = 0.50 to red Φ = 1.00. Five
different progress variable isosurfaces are marked on the plots. . . 165
xv
LIST OF FIGURES
6.5 Vorticity contour plot on x-y plane (y/Ly = 0.50, y/δL ≈ 4.88)
and x-z plane (z/Lz = 0.50, z/δL ≈ 4.88). Local equivalence ratio
contours and progress variable contours are superimposed in white
and black respectively. Blue colour indicates vorticity magnitude
of √ωkωk × δL = 5 × 10−3 while the red colour indicates vorticity
magnitude of √ωkωk × δL = 2.88. . . . . . . . . . . . . . . . . . . 168
6.6 (a) The intensity of segregation of local equivalence ratio I(ΦL)
against progress variable c (b) The difference between maximum
and minimum local equivalence ratio max(∆ΦL) against progress
variable c for five different times. . . . . . . . . . . . . . . . . . . 170
6.7 Scatter plot of SDF σ = |∇c| across the flame brush coloured
by local equivalence ratio Φ. The mean of SDF conditioned on
progress variable is plotted in red asterisks. . . . . . . . . . . . . . 172
6.8 Mean of the source terms for SDF transport equation (equation 3.47):
strain rate term aTσ, curvature term Sdκσ, and the propagation
term −∇.(Sd~nσ) conditioned on progress variable. All terms have
been normalised by SL/δL. . . . . . . . . . . . . . . . . . . . . . . 173
6.9 Scatter plot of r.h.s SDF transport equation terms against progress
variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
6.10 The surfaces of c = 0.90 isosurface coloured by normalised strain
rate aT × δL/SL in (a) and normalised curvature κ× δL in (b). . . 176
6.11 The pdfs of normalised strain rate aT × δL/SL and normalised
curvature κ × δL at t/tL = 0.87 on c = 0.10, 0.30, 0.50, 0.70, 0.90
across the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . 177
6.12 Mean value of normalised curvature κ×δL conditioned on: (a) nor-
malised tangential strain rate aT×δL/SL, (b) normalised dilatation
rate ∇.~u× δL/SL. . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
6.13 The mean of normalised SDF |∇c| × δL conditioned on (a) nor-
malised tangential strain rate aT× δLSL and (b) normalised curvature
κ× δL for c = 0.5, 0.7, and 0.9 isosurface across the flame brush. . 182
xvi
LIST OF FIGURES
6.14 Mean behaviour of displacement speed components conditioned on
progress variable across the flame brush. All terms have been nor-
malised by inlet density ρ0, laminar flame speed SL, and laminar
flame thickness δL . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
6.15 Scatter plot of r.h.s progress variable transport equation terms
against progress variable. . . . . . . . . . . . . . . . . . . . . . . . 186
6.16 The pdfs of normalised displacement speed and its components for
five c isosurfaces from c = 0.10 to c = 0.90 with increment of 0.20. 189
6.17 Mean variation of the normalised displacement speed and its com-
ponent S/SL conditioned on the normalised curvature κ×δL on five
different iso-surfaces (c = 0.50 0.70 0.90) at t/tL = 0.87, t/tt = 3.19. 192
6.18 Mean variation of the normalised displacement speed and its com-
ponent S/SL conditioned on the normalised tangential strain rate
aT × δL/SL on five different iso-surfaces (c = 0.50, 0.70, 0.90) at
t/tL = 0.87, t/tt = 3.19. . . . . . . . . . . . . . . . . . . . . . . . 194
6.19 Scatter plot of Φ across the flame brush. The mean is plotted in
red asterisks. The black solid line and the dashed line represent the
Φ and c relation of laminar flame with ΦG = 0.55 and ΦG = 0.50
respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
6.20 PDF of Φ on five different c isosurfaces. . . . . . . . . . . . . . . . 197
6.21 Scatter plot of z = zC + zH against the progress variable c at
t/tL = 0.87, t/tt = 3.19. The mean of the scatter is plotted in red
asterisks. The black solid line and the dashed line are values from
the laminar flame with Φlam,1G = 0.55 and Φlam,1G = 0.50 respectively. 198
6.22 Scatter plot of zC , and zH against the progress variable c at t/tL =
0.87, t/tt = 3.19. The mean of the scatter is plotted in red aster-
isks. The black solid line and the dashed line are values from the
laminar flame with Φlam,1G = 0.55 and Φ
lam,1
G = 0.50 respectively. . 199
6.23 Mean variation of zC,α/z conditioned on progress variable c across
the flame brush and laminar flame variation of zC,α/z across the
flame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
xvii
LIST OF FIGURES
6.24 Mean variation of zH,α/z conditioned on progress variable c across
the flame brush and laminar flame variation of zH,α/z across the
flame. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
6.25 Scatter plot of YCH4 , YH2O, YCO, and YCO2 against progress variable
c across the flame brush at t/tL = 0.87, t/tt = 3.19. The red line
represents the mean conditioned on progress variable c. . . . . . . 205
6.26 Scatter plot of YH, YH2 , and YOH against progress variable c across
the flame brush at t/tL = 0.87, t/tt = 3.19. The red line represents
the mean conditioned on progress variable c. . . . . . . . . . . . . 206
xviii
List of Tables
3.1 The chemical time tc,α and length lc,α scales for flame with unity
global equivalence ratio ΦG = 1.00. . . . . . . . . . . . . . . . . . 31
4.1 Summary of the steady laminar flame calculations . . . . . . . . . 47
4.2 The inlet boundary parameters prescribed for the unsteady lam-
inar flame simulation with sinusoidal incoming reactant mixture
strength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.3 S00 section parameters. The starting value of Φ0.90 for each section
is denoted by (). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.4 The simulation parameters prescribed for unsteady laminar flame
simulations with complementary error function inlet equivalence
ratio time history. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.5 Time averaged properties for case S00. . . . . . . . . . . . . . . . 61
4.6 The equivalence ratio and displacement speed transition time for
all the R cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.7 Percentage enhancement of key species in mass fractions and dif-
fusive mass flux across the c = 0.90 iso-surface for case R01. . . . 92
5.1 The inlet boundary conditions of the two-dimensional simulations. 102
6.1 Summary of initial set-up parameters for the DNS. . . . . . . . . 161
6.2 Statistics of aT × δL/SL on five different isosurfaces. . . . . . . . . 178
xix
LIST OF TABLES
6.3 Correlation coefficient between tangential strain rate and curvature
aT −κ , dilatation against curvature ∇.~u−κ, surface density func-
tion and tangential strain rate |∇c| − aT for different c isosurfaces
across the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . 181
6.4 Correlation coefficient between SDF and curvature |∇c| − κ, SDF
and tangential strain rate |∇c|−aT for different c isosurfaces across
the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
6.5 The mean µ, and standard deviation σ of normalised displace-
ment speed Sd/SL and its components on five different c isosurfaces
across the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . 188
6.6 The correlation coefficient of displacement speed and its compo-
nents with curvature at c = 0.5, 0.7, 0.9 across the flame brush. . 193
6.7 The correlation coefficient of the magnitude of the mixture fraction
gradient |∇z|−aT with tangential strain rate aT at c = 0.5, 0.7, 0.9
across the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . 193
6.8 The correlation coefficient of displacement speed and its compo-
nents with tangential strain rate aT at c = 0.5, 0.7 0.9 across the
flame brush. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
6.9 The correlation coefficient of the magnitude of the mixture fraction
gradient |∇z|−aT with tangential strain rate aT at c = 0.5, 0.7, 0.9
across the flame brush. . . . . . . . . . . . . . . . . . . . . . . . . 195
6.10 The mean µ, standard deviation σ, and coefficient of variation σ/µ
of the local equivalence ratio Φ on five different c isosurfaces. . . . 197
xx
Chapter 1
Introduction
At present, 90% of global energy supply is provided by combustion [118]. It is
an important method of power generation and this trend will continue for the
foreseeable future. It is estimated by Airbus1 that in year 2032, the number of
passenger planes in the world will double, reaching a total of 33,000 passenger
aircraft fleet in the world. With the growing awareness of environmental preser-
vation, pollution control regulations have become increasingly stringent.
Combustion can be coarsely divided into two categories; premixed and non-
premixed combustion. The first has its fuel and oxidiser premixed prior to com-
bustion and the latter introduces fuel and oxidiser from different streams. Both
combustion modes have their advantages and disadvantages, most notably, the
premixed flame has cleaner burnt products but it has a smaller range of oper-
ating conditions, making this mode of flame difficult to control. Also, having
premixed reactants, a premixed flame poses a danger of flashback, which may
lead to catastrophic consequences. On the contrary, having a wider operating
range, a non-premixed flame is easier to control and it will not pose any dan-
ger of flashback as fuel and air were not mixed prior to combustion. However,
this combustion mode often experiences incomplete combustion locally, forming
pollutants like soot and NOx.
In order to reduce harmful emissions, industrial combustors are designed to op-
erate under lean condition [34, 41] and with inhomogeneous mixtures, e.g. richer
1http://www.airbus.com/company/market/forecast/
1
pilot flames [4] were introduced to provide stability to the fuel-lean combustor.
This leads to a new mode of combustion, where the flame interacts with the inho-
mogeneous reactants: the is called inhomogeneous reactants premixed combus-
tion. Occurrence of inhomogeneous reactants premixed combustion is prevalent
in real world combustion systems [1, 39]. Many practical combustion systems
operate in inhomogeneous reactants mode to gain the advantages of a spatially
varying mixture fraction field, including lean premixed prevaporised (LPP) gas
turbine combustors [77], afterburners, and direct-injection spark-ignition internal
combustion engines [62].
Figure 1.1: Venn diagram depicting the combustion modes. The inhomogeneous
reactants premixed combustion is a subset of premixed combustion.
Inhomogeneous reactants premixed combustion is considered a subset of pre-
mixed combustion (refer to figure 1.1), where a premixed flame front travels
through mixture of varying stoichiometry which may be either all lean or all rich
and the stoichiometry variation is independent of the direction of flame propaga-
tion. In the literature, this is also referred to as stratified premixed [9], premixed
with variable stoichiometry, stratified mixture, imperfect premixing and other
variations.
Direct numerical simulations (DNS), a computational fluid dynamics (CFD)
tool that resolves all flow features explicitly is often used in combustion research.
2
Cant [15] presented the feasibility and challenges of turbulent premixed combus-
tion. Turbulent reacting flow is inherently a multi-scale problem and will impose
restrictive resolution requirements onthe computational mesh for DNS therefore,
DNS depends on the availability of very large computing power. However, DNS is
highly accurate and provides an enormous amount of detail for research, allowing
several turbulent flame interaction mechanisms to be studied in detail.
1.1 Research Objectives
Two questions regarding inhomogeneous reactants premixed combustion are: how
will an initially planar premixed flame respond when it propagates into inhomo-
geneous reactants? Since the mixture strength varies spatially, will the flame
exhibit different propagation speed?
The aim of the present work is to utilise the DNS to study the effect of
inhomogeneous reactants on laminar and turbulent premixed combustion under
lean condition. The objectives of this study are:
1. to assess and adapt a framework to analyse the effect of inhomogeneous
reactants on the premixed flame.
2. to investigate the effects of inhomogeneous reactants on laminar premixed
flames, specifically on two orientations:
(a) normally stratified reactants where the mixture spatial gradient is nor-
mal to the flame surface, or parallel to the direction of flame propaga-
tion.
(b) tangentially stratified reactants where the mixture spatial gradient is
parallel to the flame surface, or perpendicular to the direction of flame
propagation.
3. to investigate if differential propagation speed occurs.
4. to indicate the implications of differential propagation speed on premixed
flames.
3
5. to determine if the observations in laminar flames will be manifested also
in turbulent premixed flames.
6. to elucidate the underlying causes of these effects.
1.2 Summary of thesis structure
The outline of this thesis is as follows. In chapter 2, existing work on inhomo-
geneous reactants premixed combustion is reviewed along with other practices
available for the numerical and theoretical analysis of inhomogeneous reactants
premixed combustion.
Prior to discussing the results, chapter 3 provides the necessary mathematical
background for the present work. The basic assumptions made in the simulations
are presented as well. Along with the mathematical background, the numerical
methods adopted for the present study together with the boundary conditions
and the DNS turbulent initial conditions are discussed.
Results and discussion is presented in the subsequent three chapters, organised
according to the number of spatial dimensions considered in the simulations. The
initial flame solution of each simulations are incorporated at the beginning of the
chapter while a summary is provided at the end of each chapter.
In Chapter 4 results from the one-dimensional laminar flame simulation is
presented together with discussions concerning the normal gradient effects on the
laminar premixed flame.
Chapter 5 presents and discusses the results from the two-dimensional laminar
flame simulations. The effect of tangential mixture gradient is explored here and
the effects of differential flame propagation are presented.
Once the effects of inhomogeneous reactants on laminar premixed flames have
been established, chapter 6 presents the results of three-dimensional DNS of a
turbulent inhomogeneous reactants premixed flame. The DNS is an extension
from the laminar two-dimensional simulations that shows differential propagation
along the flame. This extends the investigation to include the effects of turbulent
straining on an inhomogeneous reactants premixed flame.
Lastly, chapter 7 draws conclusion from these observations and presents a
4
number of suggestions for future investigations into the effects of nonhomogeneous
reactants on premixed flames.
5
Chapter 2
Background
A brief introduction to turbulent premixed combustion is presented along with
relevant DNS findings. Next, a review of the definition of compositional vari-
ables used in inhomogeneous reactants premixed combustion is presented. This
is followed by key studies and findings in both the experimental and numeri-
cal simulation of inhomogeneous reactants premixed combustion are presented.
Lastly, the research approach is outlined and explained.
2.1 Premixed combustion
A premixed flame is self-propagating and it moves in the direction normal to
itself to consume the available reactant mixture [18]. It has propagation speed
and thickness determined by the thermochemical properties of the mixture. This
propagation speed and thickness are also known as the laminar flame speed SL
and laminar flame thickness δL respectively. The laminar flame speed quantifies
the rate of which the flame can consume the reactants, hence the rate of heat
release.
In turbulent premixed combustion, turbulent eddies with a range of length
scales can interact with the premixed flame. Central to turbulent premixed com-
bustion modelling is the laminar flamelet theory, where the inner structure of the
premixed flame is unaffected by the underlying turbulence flow field, providing a
flamelet-like structure. This conveniently separates the combustion phenomenon
6
from the turbulent flow field as the reaction is assumed fast, and takes place in a
thin sheet, confining the flame into thin layers.
Turbulent eddies will interact with the thin flamelets, enlarging the available
flame surface are. The larger flame surface area provides a larger rate of re-
actant consumption, and therefore a larger associated (turbulent) flame speed.
Damko¨hler [37] hypothesised that the turbulent to laminar flame speed is pro-
portional to the area ratio.
The ratio of turbulent r.m.s fluctuating velocity u to laminar flame speed,
and the turbulent integral length scale lt to laminar flame thickness is used to
characterise the turbulent premixed flame in a diagram called the Borghi diagram
[11]. This diagram has been refined over the years and a modified Borghi diagram
is proposed by Peters [82] (shown in figure 2.1).
10−1 100 101 102 103 104
10−1
100
101
102
103
u/
S L
lt/δL
laminar flames wrinkled flamelets
corrugated flamelets
thin reaction zones
broken reaction zones
ηk = δL
ηk = δr
 
 
Ret = 1
u/SL = 1
Ka = 1
Kaδ = 1
Figure 2.1: Modified Borghi diagram
The dimensionless numbers in figure 2.1 are the Karlovitz number Ka, Damko¨hler
number Da, and Reynolds number Re. The Karlovitz number Ka is the ratio be-
7
tween chemical time scale τc and the Kolmogorov time scale τη [82].
Ka = τcτη
(2.1)
The Damko¨hler number Da is another time scale ratio with the flow time scale
being the turbulent eddy turnover time τt
Da = τtτc
(2.2)
The chemical time scale τc and turbulent eddy flow time τt are defined as
τc =
D
S2L
= δLSL
(2.3a)
τt =
lt
u (2.3b)
where D is the diffusivity. The Reynolds number Re is related to Ka and Da
through the following expression [18, 32]:
Re ∼ Da2Ka2 (2.4)
The thin reaction zones regime is important in the present work. This regime
is located above the Klimov-Williams criterion (red line, where Ka=1)[119]. The
flame thickness is comparable to the turbulent length scale while the flame reac-
tion zone thickness δr is still smaller than the turbulent length scale. Therefore,
the turbulent eddies are expected to disturb the preheat zone of the flame while
the reaction zone remains thin. However the thin reaction zone is well-preserved
and is the marker that delineates the burnt and unburned gases. Thus, flamelet
theories are still applicable in the thin reaction zone regime. Poinsot et al. [89]
and Helie and Trouve´ [61] have noted that the flamelet assumption provided a
good framework to analyse turbulent inhomogeneous reactants premixed com-
bustion.
The other zones shown in figure 2.1 are discussed in Peters [82] and Cant and
Mastorakos [18].
8
2.2 DNS in turbulent premixed combustion
DNS has been used to investigate the physics of turbulent premixed flame. DNS is
computationally expensive as it explicitly resolves all flow features [15, 87]. DNS
of combustion is often done in either two-dimensions with reduced chemistry
[30, 35, 36, 42, 43, 49, 57, 64, 83] or three-dimensions with simple chemistry
[21, 24, 53, 60, 66, 72, 99, 111, 122]. With the availability of high performance
computing, several three-dimensional DNS with reduced chemistry have been
performed [7, 100, 108, 121]. Note that this is not an exhaustive list, for example
DNS that uses the flamelet-generated manifold method [79] is not considered
here.
Predicting available flame surface area in a volume is important in flamelet-
based model as it is proportional to the local mean burning rate [19, 116]. Flame
surface density Σ is a measure of the flame surface area per unit volume. In
turbulent premixed combustion modelling, Σ can be modelled using an algebraic
expression, or obtained using a differential equation. The latter has advatages
over the former for cases where flame surface wrinkling is high and flame propaga-
tion is highly unsteady [54] as it captures the ‘memory’ effect of the evolution of
FSD [40]. FSD has been studied theoretically [13, 19, 90]. The effects of unequal
mass and heat diffusivity, characterised by Lewis number Le, on FSD have been
studied by DNS [29, 30, 53, 60, 111].
In analysing the DNS database for turbulent premixed flame, the location of
the flame is tracked using a scalar: either using a passive scalar G as used in the
level-set approach [86] or progress variable isosurface c∗ [28] in the FSD approach.
In the thin reaction zones regime, eddies can enter the flame, the flamelet
becomes thinner and more wrinkled and making the flame more susceptible to
curvature effects, making curvature effects on the Sd important. Gran et al. [49],
Peters et al. [83] and Echekki et al. [42] decompose the displacement speed
into reaction, normal and tangential diffusion components. Using scaling anlysis,
Peters [82] showed that dependence of displacement speed becomes the leading
order contributor in the thin reaction zone. Hawkes and Chen [57] concur that
inclusion of the tangential diffusion term improves the modelling of Sd in thin
9
reaction zone. Recent DNS [24, 31, 66, 72] shows that the reaction and normal
diffusion term will also have dependency on curvature in thin reaction zones
regime. In regions with large positive curvature (flame curved convex towards
reactants), Gran et al. [49] observed that negative Sd is possible and this is
supported by many recent DNS, e.g. [22] where non-zero probability of negative
Sd is observed. Other DNS [23, 25–28, 30] have looked into effects of local strain
rate and curvature effects on the displacement speed Sd.
2.3 Composition variables for the inhomogeneous
reactants premixed combustion
Inhomogeneous reactants premixed combustion may be described in terms of
two composition variables [12]: one variable to describe the mixture compo-
sition and another to describe the progress of the premixed reaction [61]. A
modified progress variable is necessary to accommodate the variation in mixture
strength [61, 78, 89]. However, the modification of progress variable is often linked
with the definition of mixture fraction. An additional composition variable in the
form of local equivalence ratio is also defined in simulations [59, 67, 85, 123] of
inhomogeneous reactants premixed combustion. The equivalence ratio is often
defined as a function of mixture fraction.
2.3.1 Mixture fraction
The mixture fraction is often expressed using both fuel and oxidiser mass frac-
tion [8]. For DNS using simple chemistry, the mixture fraction is defined as [61,
75, 89]:
Z = YF − YO/rs + 1/(rs + rsb)1 + 1/(rs + rsb)
(2.5)
where YF and YO are the fuel and oxidiser mass fraction respectively, rs is the
stoichiometric oxidiser-fuel mass ratio, and b is the oxidiser nitrogen-oxygen mass
ratio.
In complex chemistry simulations, equation 2.5 is modified to include the
radicals. Element mixture fraction is often used to define the mixture fraction of
10
a multiple species system [59, 67, 85, 93, 96]. Haworth et al. [59] and Jimenez et
al. [67] define their mixture fraction as:
z = zC + zH (2.6)
where zC and zH are the carbon and hydrogen element mixture fractions respec-
tively. The element mixture fraction is obtained from:
zβ =
Ns∑
α=1
nβαYα
Wβ
Wα
(2.7)
where β refers to the elements, Ns is the total number of species considered, Wβ
is the molecular weight of element β and nβα is the number of atoms of element
β in species α. A similar definition of element mixture fraction is also given by
Pires Da Cruz et al. [85].
2.3.2 Progress variable
The progress variable is often used as the second variable for inhomogeneous
reactants premixed combustion. As mentioned previously, a modified progress
variable, which is a function of mixture fraction is defined to account for the
variation in mixture strength. Below is the fuel based progress variable adopted
by [52, 59, 61, 67, 75, 85, 89, 96]
c(z) = YF,∞(z)− YF (z)YF,∞(z)− YF,−∞(z)
= 1− YFzYF,∞
(2.8)
Note that an oxidiser based progress variable has also been used in the one-
dimensional simulation by Richardson [93]. This is because a fuel based progress
variable for methane-air flames has a very high progress variable value at the
location of maximum heat release [58]. Grout [50] introduced the aging progress
variable to extend the view of the flame structure.
Haworth et al. [59] and Jimenez et al. [67] used a modified version of equa-
tion 2.8:
c(z) = 1− YF/z (2.9)
11
The above formulation is closely coupled with the definition of z. The simplifica-
tion is possible because z represents the combined element mixture fractions of
the fuel (refer to equation 2.6).
In the absence of element mixture fraction measurements, a progress variable
definition based on termperature is often used in the experiments [105, 107].
2.3.3 Equivalence ratio
Mu¨ller et al. [78] relates the definition of equivalence ratio Φ with mixture fraction
z:
φ = rs
YF
YO
= z(1− zst)zst(1− z)
(2.10)
where zst is the stoichiometric mixture fraction. This definition is adopted by
Pires Da Cruz et al. [85], Haworth et al. [59] and Jimenez et al. [67] in their
definition of local equivalence ratio.
Note that Sweeney et al. [105] and Barlow et al. [3] used alternative definition
of equivalence ratio (and therefore mixture fraction) that consists of only the
experimentally measurable species mole fractions of CH4, O2, CO2, H2O, CO,
and H2.
2.4 Numerical studies of inhomogeneous reac-
tants premixed combustion
Inhomogeneous premixed combustion has been studied numerically in one-dimensional
laminar flame simulations and in two- and three-dimensional DNS of turbulent
flames. These simulations will be presented in two parts: the multi-dimensional
numerical simulation, and the one-dimensional numerical simulations.
2.4.1 Multidimensional turbulent inhomogeneous reactants
premixed combustion numerical simulations
Three-dimensional DNS of a flame propagating into stratified reactants has been
investigated by Poinsot et al. [89], Helie and Trouve´ [61], Grout et al. [51, 52],
12
and Malkeson and Chakraborty [74, 75]. Limited by the computational resources,
the calculations by Poinsot et al. and Helie and Trouve´ were restricted to sim-
ple one-step irreversible reaction mechanism representing propane-air chemistry.
Similarly, Malkeson and Chakraborty used a simple one-step irreversible reaction
mechanism representing methane-air chemistry while Grout et al. used a 2-step
4-species reduced mechanism for methane-air combustion.
Aside from DNS using simple chemistry, DNS of partially premixed flames
using realistic chemistry and transport were also performed. Haworth et al. [59]
and Jime´nez et al. [67] used this approach to perform a two-dimensional DNS
with detailed chemistry and transport to investigate the effect of stratified re-
actants towards both primary and secondary flames. The reduction to a two-
dimensional domain is necessary to compensate for the additional cost incurred
in complex chemistry calculations. Both simulations used propane-air reactants
with chemical kinetics based on the 29-species reaction mechanism published by
Warnatz [117]. Subsequent modifications were made by Haworth et al. to reduce
the stiffness of the mechanism. Jimene´z et al. extended the modified chemistry
to incorporate the Zeldovich mechanism for thermal NO production.
All simulations mentioned above employed a freely propagating turbulent pre-
mixed flame that interacts with an oncoming inhomogeneous reactant stream in-
side a decaying turbulent flow field, with the exception for Grout et al. [52] who
employ a turbulent velocity flow field supplied to the inlet boundary in the man-
ner pioneered by Rutland and Cant [99]. The simulations can be broadly divided
into two categories, depending on the methodology in which the inhomogeneous
reactants are introduced to the turbulent premixed flame.
In the first category, the equivalence ratio Φ variation in the reactants is
initialised using a pseudo-spectral method proposed by Eswaran and Pope et al.
[45]. This yields a random distribution of Φ following a bi-modal distribution
with a prescribed mean Φ, r.m.s. fluctuation Φ′, and scalar length scale lΦ that
is comparable to that of the turbulent eddies lt. Poinsot et al. [89], Helie and
Trouve´ [61], Grout et al. [52], and Malkeson and Chakraborty [74] have adopted
this method to initialise the inhomogeneous reactants.
The second category make use of a larger scalar length scale lΦ > lt of the in-
homogeneous reactants. The equivalence ratio variation is initialised using a sinu-
13
soidal function extending across the span-wise direction of the simulation domain.
This provides an equivalence ratio variation parallel to the initial planar flame.
Haworth et al. [59], Jime´nez et al. [67], and Malkeson and Chakraborty [75]
follow this approach.
2.4.2 Reported results for DNS of turbulent inhomoge-
neous reactants premixed combustion
Global flame properties like the total reaction rate R, flame wrinkling factor 〈W 〉,
and local mass burning rate rˆ were reported by Poinsot et al. [89] and Helie and
Trouve´ [61]. Poinsot et al. [89] used the same numerical configuration as adopted
by Trouve´ and Poinsot [111], placing the flame in the thin reaction zones. The
initial scalar bi-modal distribution is set to have equivalence ratio values of 0.5
and 1.1 at the two peaks with averaged equivalence ratio of 0.8. They concluded
that the flame stretch increases in the presence of reactants stratification, thus
slight enhancement of the heat release is observed in the form of production of
flame surface area. However, the flame surface enhancement by stratification is
negligible when compared to the flame surface enhancement due to fluid motion.
Helie and Trouve´ [61] used the same numerical approach as Poinsot et al..
The key difference in their approach is that the peaks of the initial scalar bi-
modal distribution are set to equivalence ratios equal to 0.4 and 1.6 with averaged
equivalence ratio set to be unity, giving the emphasis on a configuration featuring
both rich and lean conditions. Also, the scalar r.m.s fluctuation is set higher by
Helie and Trouve´ [61].
The stratified case has higher total reaction rate compared to the homogeneous
case initially. However, after about one eddy-turn over time, the total reaction
rate for the homogeneous case is reported to be higher than that of stratified
case. A similar change in wrinkling factor is also observed. The stratified case
has significantly higher wrinkling factor while having almost equal mean mass
burning rate before one eddy-turnover time. Since total reaction rate is the
product of wrinkling factor and mean mass burning rate, it can be concluded that
stratification will incur more flame wrinkling (higher wrinkling factor), at least
14
at the initial state as observed by Helie and Trouve´ [61]. This non-sustainable
gain in terms of total reaction rate and wrinkling factor in the stratified case
can be explained either as the effect of initial transients or due to the effect
of decaying velocity and scalar fluctuations in the reactants. Given that the
numerical configuration used was inherently unsteady, the contribution of the
initial gain cannot be unambiguously determined.
Helie and Trouve´ [61] reported that the reduced mass burning rate per unit
flame area or the normalised flame speed is seen to correlate strongly with the
local mixture composition. The correlation follows closely the variation of mass
burning rate and it is found to be nonlinear. However, if only the lean to stoichio-
metric region were considered, the z-distribution is close to symmetric. Therefore
at non-rich regions of partially premixed combustion, departures of r from unity
will cancel in the mean when averaged over the whole flame.
Haworth et al. [59] and Jime´nez et al. [67] also used similar global variables
to quantify the partially premixed flame. In the studies of Haworth et al. [59],
the reported results showed that the partially premixed flame has no significant
departure from the homogeneous counterpart in terms of its flame length and
heat release for the first two flame time τf . It is argued by Haworth et al. [59]
that the global heat release for the inhomogeneous reactants flame will be lower
if longer time integration was performed. However, the turbulence intensity is
not sustained throughout the simulation, giving a decay in turbulence with time.
Hence an ambiguity arises as the turbulence effect and the partially premixed
reactants effect cannot be clearly separated.
Jime´nez et al. [67] extended from Haworth et al. to study the mixture dis-
tribution effect in inhomogeneous reactants premixed combustion. It was shown
that for a given mixture composition pdf that includes a high probability of
finding lean reactants, the strongly stratified mixture distribution which is char-
acterised by lower unmixedness length scale lφ will produce a higher heat release
rate and flame length. In fact, with high lφ the heat release will be smaller than
the homogeneous counterpart. This is explained by the residence time of the lean
pockets found in the partially premixed flame. A high value of lφ will enhance
the probability of finding lean pockets, therefore local mixing of the lean and rich
15
pockets becomes important which leads to a lower local flame speed. However if
lφ is small, the reactants were stratified and the flame at the locally lean mixtures
will be back supported by heat from the rich pockets thus increasing the locally
lean mixture flame speed, resulting in an increase of global flame area.
The simulations performed by Grout et al. [52] fall into the wrinkled flamelet
regime. Under these condition, the velocity fluctuation is comparable to the
flame speed fluctuation (u′/SL ≈ S ′L/SL) where S ′L is the laminar flame speed
fluctuations arising due to compositional variations in the reactant field and SL
is the averaged value. The surface area is reported to have an additional 20%
increase when compared to the homogeneous counterpart. Evidence of differential
propagation was observed as mixture fraction and progress variable gradients
align strongly at the leading edge of the flame but the alignment grows weaker
towards the trailing edge of the flame. Grout et al. suggested that with larger
turbulence intensity, the effects of stratification will be less pronounced.
Malkeson and Chakraborty [75] performed several simulations that belong to
the thin reaction zone regime. Displacement speed statistics were emphasised.
The mean of the cross-dissipation term Scd is found to be negligible compared
to the mean of the diffusion Sn and reaction Sn terms. Flames with unity mean
global equivalence ratio Φ are found to have lower displacement speed when
compared to homogeneous counterpart. Enhancement of the displacement speed
Sd against homogeneous counterpart is observed for flames with lean Φ at high
turbulence intensity.
The statistics of displacement speed Sd are observed to be consistent with pre-
vious DNS of turbulent premixed combustion [22, 27]. Moreover, the curvature
κ and strain rate aT dependence of displacement speed Sd in the case of turbu-
lent stratified flames are quanlitatively similar to those observed in homogeneous
premixed flames and mixture inhomogeneity is reported to have no significant
impact on these dependencies.
16
2.4.3 One-dimensional laminar inhomogeneous reactants
premixed combustion numerical simulations
One dimensional simulations of laminar inhomogeneous reactants premixed com-
bustion were performed to investigate the effects of inhomogeneous reactants that
have composition gradient normal to the direction of flame propagation.
Pires da Cruz et al. [85] have performed a numerical calculation on one-
dimensional laminar methane-air freely-propagating flames with inhomogeneous
reactants using the GRI 2.11 mechanism.
Marzouk et al. [76] and Richardson et al. [93] performed numerical calcula-
tions on one-dimensional laminar methane-air flames in a counterflow configura-
tion. Marzouk et al. used a C1 kinetic model for methane-air combustion that
consists of 46 reactions and 16 species [103]. Richardson et al. adopted the GRI
3.0 natural gas mechanism.
2.4.4 Reported results for one-dimensional simulations of
laminar inhomogeneous reactants premixed combus-
tion
Pires da Cruz [85] performed unsteady one-dimensional numerical simulations
of both rich and lean methane flames, each propagating into leaner or richer
mixtures, conveniently described as back or front supported flames. The variation
in equivalence ratio is introduced in a step-change manner. However, due to flame
attenuation, the lean flame cases covered the range of equivalence ratio between
unity and 0.43 within an elapsed time of 60 ms. This gives a temporal equivalence
gradient of about 9.5 s−1. It is found that the premixed flame propagating from
stoichiometric to lean conditions (back support) is controlled by high burned gas
temperature, providing a faster flame relative to the homogeneous equivalent.
This back supported premixed flame is found to have higher concentrations of
molecular and atomic hydrogen in the burned gas. In general, the acceleration
effect due to the higher burned gas temperature provides 10-20% boost of local
flame speed.
Marzouk et al. [76] subjected the methane-air counterflow flame to a gradi-
17
ent of equivalence ratio and strain. They observed that the rate of heat release
is higher for lean flames than for their homogeneous equivalent, consistent with
Pires Da Cruz et al. [85]. They attributed this to support from the burned gas
via diffusion of heat and radicals. When comparing to a homogeneous flame with
equal heat release rate, the stoichiometric to lean flame has a broader flame thick-
ness. In other words, the back supported flame is observed to become thinner.
Richardson et al. [93] systematically investigated the behaviour of steady,
stabilized strained flames with fresh reactants on one side and burned products
on the other. The opposing streams have different equivalence ratios, creating an
equivalence ratio gradient across the flame.
Agreeing with Pires Da Cruz et al. and Marzouk et al., Richardson et al.
showed that under lean conditions, the back supported flames exhibit higher
propagation speed. The flame thickness is also observed to have decreased, con-
sistent with Marzouk et al. [76]. They showed that back supported flames display
enhanced concentrations and upstream fluxes of H, H2 and OH through the reac-
tion zone. The enhanced concentrations of highly reactive radical species in the
reaction zone causes back supported flames to propagate faster than equivalent
premixed flames. This shows that the effect of equivalence ratio gradients on
steady flames is primarily due to radical pool flux rather than heat flux.
Definition of stratified flame propagation based upon the displacement speed
of a mixture fraction dependent progress variable is advocated. Using displace-
ments speed components, Richardson et al. found that the normal diffusion
component of the displacement speed responds rapidly to strain variation (from
the counterflow), causing the changes of the reaction term and subsequently al-
tering the location of the flame. The cross dissipation term has a slow response
to strain and its response is determined by the flow and flame time scales.
18
2.5 Experiments on inhomogeneous reactants pre-
mixed combustion
Kang and Kyritsis [69] measured the laminar flame speed for a back supported
flame that propagates into lean mixtures. The results indicate that the back
supported flame propagates at a higher speed compared to an equivalent pre-
mixed flame. They indicated that the elevated heat flux from the products is the
dominant mechanism of increasing the flame speed.
Similarly, Galizzi and Escudie´ [46] performed laminar V-flame experiments
and found that laminar flame speeds are larger than the expected values based
on the local equivalence ratio.
The increased flame speed measured in both laminar inhomogeneous reactants
premixed combustion experiments above are consistent with the numerical studies
[76, 85, 93].
The low turbulence (u/SL ∼ O(1)) and lean inhomogeneous reactants premixed
combustion has been studied using different experimental set ups. In chronolog-
ical order, they are presented below.
Pasquier et al. [81] studied the propagation of a spark ignited flame through
a low turbulence inhomogeneous propane mixture in a constant volume chamber
(combustion bomb). Instantaneous velocity and equivalence ratio are measured.
The averaged flame propagation is enhanced due to inhomogeneous reactants.
Local flame propagation is shown to be influenced by the local equivalence ratio,
where flame burning towards lean pockets is found to have enhanced propagation
speed.
Robin et al. [94] studied a low turbulence rod-stabilised methane-air flame in
a V-configuration. The inhomogeneous mixture has a Gaussian distribution of
mean mixture fraction (equivalence ratio). The flame thickness is measured and
found to be relatively thinner than the equivalent homogeneous flame. This is
attributed to additional stretch arising from equivalence ratio fluctuation induced
differential propagation.
Anselmo-Filho et al. [2] used a low turbulence, rod-stabilised slot burner to
investigate the effect of inhomogeneous reactants on the geometric properties of
19
the flame. The flame is supplied with lean mean equivalence ratio of 0.77 and
the low turbulence intensity placed the flame in the corrugated flamelet regime.
The FSD and curvature of the flame was measured and analysed. The reaction
rate and FSD are found to be higher compared to the homogeneous counterpart.
Analysis of curvature indicates that the distribution of curvature is broadened by
inhomogeneous reactants. Differential propagation is suggested as the mechanism
for these effects.
Extending from their previous work [46], Galizzi and Escudie´ [47] studied
the low turbulence inhomogeneous reactants methane-air V-flame. Inhomoge-
neous reactants are reported to alter the velocity field, increasing the flame brush
thickness, and to accelerate the burned gas towards the flame front opposite the
stratification zone.
Sweeney et al. [106] extended the investigation of Anselmo-Filho et al. [2]
to include an axisymmetric co-annular burner for comparison. This enables a
comparison between swirling and non-swirling lean inhomogeneous reactants pre-
mixed combustion. FSD is reduced in the non-swirl burner, but a similar reduc-
tion is not observed in the swirl burner.
Vena et al. [114] studied equivalence ratio gradient effects on the flame front
topology of low-turbulence, iso-octane/air V-flames. Contrary to Anselmo-Filho
et al., the variation in curvature distribution is reported to be negligible and the
FSD increase is modest in the inhomogeneous reactants premixed combustion
due to a decrease in flame length and an increase in flame brush thickness.
Extending the analysis from Anselmo-Filho et al. [2] and Sweeney et al. [106],
Sweeney et al. [105] measured temperatures and major species concentrations of
the slot burner. In agreement with Anselmo-Filho, the curvature distribution for
flames with inhomogeneous reactants is broader. The turbulent flame brush is
found to be thicker, in agreement with DNS [59, 67]. FSD is increased compared
to the homogeneous case and the flame structure is well represented by laminar
flamelets, in agreement with DNS [52, 61, 89]. A deviation of mixture fraction is
observed across the reaction zone. The deviation is attributed to the differential
diffusion of carbon- and hydrogen-containing species.
20
2.6 Research approach
Present research aims to investigate the displacement speed Sd of unsteady pre-
mixed flame propagating into inhomogeneous reactants. The progress variable
transport equation proposed by Bray et al. [12] is adapted for the complex chem-
istry context and also extended to include the non-equal diffusivity term arises
from the difference in diffusivity between mixture fractions and progress variable
(refer to equation 3.34).
One-dimensional unstrained laminar premixed flame is simulated in both
steady and unsteady configurations. The first will provide baselines or refer-
ence cases which are used to compare with the latter. From the reference cases,
the mixture profile across the flame is examined and shown to have a dip across
the flame (refer to section 4.2)
The unsteady premixed flame response towards inhomogeneous reactants is
studied using the unstrained laminar flame configuration. This configuration is
similar to Pires Da Cruz et al. [85], however the analyses are performed using
the flamelet theory framework. The unsteady flame performance is compared
with the corresponding steady or reference flame. This analysis extend from
Richardson et al. [93] to investigate the Sd behaviour of the premixed flames that
are freely propagating into inhomogeneous reactants. This eliminates any effects
from the enhanced mixture gradient across the flame due to strain rate associated
with counterflow configuration. The behaviour of Sd is explored by studying the
contributions of its components and the time response to the inhomogeneous
reactants. This provides the understanding of the occurrence of Sd hysteresis in
the unsteady flame which will be discussed in detail in chapter 4.
Extending from the one-dimensional flame analysis, the two-dimensional un-
steady laminar flame simulations will explicitly show the effect of differential flame
propagation. Moreover, the evolution of two-dimensional flame propagating into
reactants with mixture fraction parallel to the direction of flame propagation is
captured and studied, showing that the two-dimensional flame is not an almaga-
tion of differnt equivalence ratio one-dimensional laminar flames. Note that this
is done with the absence of turbulent straining.
Finally, DNS is performed on inhomogeneous reactants turbulent premixed
21
combustion with configurations similar to that of two-dimensional case. This
is done to examine the turbulent straining effect competitions with the mixture
fraction gradient effects. This extends from the study performed by Malkeson
and Charkraborty [75] to include complex chemistry. Similarly this DNS extends
from Haworth et al. [59] and Jimenez et al. [67] from two-dimensional DNS to
a fully compressible three-dimensional DNS. By using similar configuration from
the two-dimensional laminar flame simulations, the evolution of the turbulent
flame can be compared to the laminar counterpart, providing direct comparison
between flame in temporal and spatial evolution.
2.7 Summary
Background material on the experimental works and numerical simulations re-
garding lean inhomogeneous reactants premixed combustion and premixed com-
bustion are provided. The research approach is also outlined and explained.
Two compositional variables are required for characterising inhomogeneous
reactants premixed combustion [12]; mixture fraction z and progress variable c
are typically used in both experiments and numerical studies.
Both experiments and multi-species numerical simulations used element mass
fractions to define the mixture fractions. Progress variable is then linked with
mixture fraction to account for the variation of equivalence ratio.
Deviation of mixture fraction across the reaction zone is observed in experi-
ments [3, 105] and this deviation is attributed to preferential diffusion.
There is a consensus that a lean inhomogeneous reactants premixed flame that
is propagating from stoichiometric to lean mixtures in a laminar or low turbulence
state has a larger propagation speed than the equivalent homogeneous premixed
counterpart.
The species H, H2 and OH are observed to be the key species that influence
the propagation speed of inhomogeneous reactants premixed flames shown in
numerical simulations [76, 85, 93].
There is a disagreement between experimental measurements of flame thick-
ness in lean inhomogeneous reactants premixed flames in low turbulence. How-
ever, later experiments [105] showed agreement with DNS [59, 67]. The FSD
22
for inhomogeneous reactants premixed flames is larger in value compared to the
homogeneous counterpart.
23
Chapter 3
Mathematical Background and
Numerical Implementations
The mathematical formulations used in this dissertation and the numerical imple-
mentations are introduced here. The governing equations for the simulations are
outlined, followed by a brief introduction on the numerical schemes and chemical
scales. Initialisation of the turbulent flow field for three-dimensional DNS and
the specification of boundary conditions are outlined. Lastly, equations relevant
for post-processing analysis are presented.
3.1 Governing Equations
The DNS code SENGA2 developed by Cant [17] is used to perform all the simula-
tions presented in this disseration. It solves the fully compressible Navier-Stokes
equations on a regular finite difference grid. A brief description of the governing
equations is presented next.
The governing equations that describe the gaseous reacting flow consist of
continuity, momentum, energy, and species transport equations. There are 5+Ns
governing equations used in SENGA2 [17]; five compressible flow equations and
an additional Ns chemical species mass fraction balance equations. The five
compressible flow equations are: mass continuity (3.1), three for the Navier-Stokes
momentum (3.2), and the internal energy (3.3). The extra Ns equations (3.4) are
24
solved to account for the transport and reaction of the chemical species involved in
the combustion process. Note that the these equations are presented in Cartesian
coordinates.
∂
∂tρ +
∂
∂xk
ρuk = 0 (3.1)
∂
∂tρui +
∂
∂xk
ρukui = −
∂
∂xi
P + ∂∂xk
τki (3.2)
∂
∂tρE +
∂
∂xk
ρukE = −
∂
∂xk
Puk −
∂
∂xk
qk +
∂
∂xk
τkmum (3.3)
∂
∂tρYα +
∂
∂xk
ρukYα = ω˙α −
∂
∂xk
ρVα,kYα α = 1, 2, . . . Ns (3.4)
where ρ is the mixture density, ui is the i-th component of the velocity vector, P
is the pressure, E is the stagnation internal energy and Yα and ω˙α are the mass
fraction and reaction rate of species α of N total species in the reacting mixture.
The viscous stress tensor in equation 3.2 and 3.3 is given by:
τki = µ
(∂uk
∂xi
+ ∂ui∂xk
)
− 23µ
∂um
∂xm
δki, (3.5)
where µ is the dynamic viscosity of the mixture. The stagnation internal energy
in equation 3.3 is defined as:
E =
N∑
α=1
Yαhα −
P
ρ +
1
2ukuk (3.6)
where the enthalpy of species α is defined as:
hα =
∫ T
T0
CpαdT + h0α, (3.7)
in which Cpα is the mass-based specific heat capacity of species α and h0α is the
species enthalpy at the reference temperature T0. The heat flux vector is given
25
as:
qk = −λ
∂T
∂xk
+
Ns∑
α=1
ρVα,kYαhα (3.8)
where λ is the mixture thermal conductivity and Vα,k is the k-th component of
the diffusion velocity vector of species α relative to the mixture.
3.1.1 Thermochemical properties
An ideal gas is assumed, hence the thermal equation of state for the mixture is
given by
P = ρR0T
Ns∑
α=1
Yα
Wα
(3.9)
where Wα is the molar mass of species α. The chemical production rate ω˙α for
species α in equation 3.4 is expressed as:
ω˙α = Wα
M∑
m=1
(
ν ′′α,m − ν ′α,m
)
km(T )
Ns∏
β=1
cν
′
β,m
β . (3.10)
with ν ′α,m and ν ′′α,m represent the reactant and product stoichiometric coefficients.
The specific reaction rate coefficient km(T ) from equation 3.10 is given by the
Arrhenius expression
km(T ) = AmT nm exp
(
− EmR0T
)
(3.11)
where Am, nm and Em are the frequency factor, temperature exponent and ac-
tivation energy. In the presence of a reversible chemical reaction, the forward
reaction rate coefficient kf,m and the reverse reaction rate coefficient kr,m are
related by:
ln kr,m = ln kf,m +
Ns∑
α=1
(
ν ′′α,m − ν ′α,m
)
( g¯α
R0T + ln
P0
R0 − lnT
)
(3.12)
where P0 is the reference pressure used for nondimensionalising the partial pres-
sure equilibrium constant and g¯α is the molar Gibbs function for species α.
26
3.1.2 Transport properties
The molecular transport coefficients and mixture viscosity are represented as:
µ
Pr =
λ
Cp
= Aλ
( T
T0
)r
(3.13)
where µ is mixture dynamic viscosity, Pr is the mixture Prandtl number which
is taken constant at Pr=0.7, λ is the thermal conductivity, Cp is the constant
pressure specific heat capacity of the mixture and Aλ, r and T0 are constants
with values of Aλ = 2.58 × 10−5kg m−1s−1, r = 0.7 and T0 = 298K. The species
diffusive mass flux assumes agreement with the Fick’s law:
ρVα,kYα = −ρDα
∂Yα
∂xk
(3.14)
The diffusion coefficient Dα for each species is calculated using its prescribed
Lewis number Leα as given by:
Dα =
λ
ρCP
1
Leα
(3.15)
The Soret and the Dufour effects are the secondary diffusion effects that account
for the diffusion of mass by the temperature gradient and the diffusion of heat
by the concentration gradient. These secondary diffusion effects are assumed to
be negligible in the current simulations.
3.1.3 Diffusion correction velocity
Since all the Ns species transport equations are considered, a correction velocity
V (c)k is used to ensure that the continuity equation 3.1 is satisfied.
ρVα,kYα = −ρDα
∂Yα
∂xk
+ ρV (c)k Yα
= −ρDα
∂Yα
∂xk
+
( Ns∑
α=1
ρDα
∂Yα
∂xk
)
Yα (3.16)
27
A full account of the SENGA2 governing equations can be found in Cant [17].
3.2 Spatial discretisation and temporal advance-
ment schemes
The high fidelity numerical solver requires high accuracy spatial discretisation
schemes. The SENGA2 uses the tenth order central difference scheme for the
interior nodes. The central difference scheme formulations for the first derivatives
are given as:
f ′i =
m/2
∑
j=1
aj
2jh (fi+j − fi−j) (3.17)
where m=10 for the tenth order scheme. The constants aj are obtained by ex-
panding the Taylor series and h is the mesh spacing.
The tenth order interior scheme is reduced in order to 8th, 6th, 4th, 4th skew
and 4th one-sided as the boundary is approached.
The two registers five-stage explicit Runge-Kutta time integration method
derived by Kennedy, Carpenter and Lewis [71] is adopted for time advancement.
This method is fourth order accurate and requires five sub-steps for each time
step. Only two storage locations are required per equation per grid point using
the storage reduction methodology introduced by van der Houwen [113] and Wray
[120], one for the time derivative and one for the dependent variable.
The adaptive time stepping strategy, where the time step magnitude varies
from one time step to another, is used on the one- and two-dimensional flame. The
magnitude of this time step is determined by a proportional-integral-derivative
(PID) controller. A constant time step is used for the three-dimensional DNS
(refer to chapter 6).
28
3.3 Determination of the chemical time-step and
grid resolution requirements
The C1 methane-air mechanism using 64 steps and 18 species proposed by War-
natz [118] and modified by Gran et al. [49] is used for the present study. This
mechanism has been successfully used in previous DNS [35, 42, 44, 49, 73].
In a complex chemistry reacting flow, the species in the chemical reaction
introduce additional time and length scales. These scales should be captured for
accurate numerical simulations. To assess these characteristic chemical scales,
the chemical time scale (tc,α) and length scale (lc,α) for each species is defined
as [59]:
tc,α =
Yα,max
|Y˙ |α,max
(3.18)
lc,α =
Yα,max − Yα,min
|∇Yα|max
(3.19)
where Yα,max and Yα,min are the maximum and minimum mass fraction of species
α, the value of |∇Yα|max is the maximum spatial gradient of the species α mass
fraction and |Y˙ |α,max = ω˙α/ρ is the maximum molar chemical production rate for
species α.
The characteristic chemical time and length scales are obtained based on
steady one-dimensional laminar premixed flame calculations. The time and length
scales of the chemical species for laminar flame with global equivalence ratio from
ΦG = 0.50 to ΦG = 1.00 are presented in three-dimensional bar chart in figure 3.1a
and figure 3.1b respectively. The vertical axes denotes the global equivalence
ratio, while the x and y axes denote the time or length scales, and the species.
From figure 3.1, the time and length scales of the chemical species are observed
to increase with decreasing global equivalence ratio. The smallest time and length
scales in the equivalence ratio range presented is found in laminar flame with
ΦG = 1.00. The species chemical time and length scales for the flame with
ΦG = 1.00 are tabulated in table 3.1.
The hydrogen atom H is found to have the smallest time scale with the value of
5.3076×10−5 s while carbon dioxide CO2 has the longest time scale. The methyl
29
0
0.02
0.04
1.00
0.95
0.90
0.85
0.80
0.75
0.70
0.65
0.60
0.55
0.50
CH
3 O
CH
3
CH
3 OH
CH
2 OH
HCOCH2 OCOH2 O
2
HO
2
H
Yα
OHOH
2H2 OCO
2
O
2CH
4tc,α
Φ G
(a) The laminar flame species chemical time scale tc,α.
0
2
4x 10
−3
1.00
0.95
0.90
0.85
0.80
0.75
0.70
0.65
0.60
0.55
0.50
CH
3 O
CH
3
CH
3 OH
CH
2 OH
HCOCH2 OCOH2 O
2
HO
2
H
Yα
OHOH
2H2 OCO
2
O
2CH
4lc,α
Φ G
(b) The laminar flame species chemical length scale lc,α.
Figure 3.1: The laminar flame species chemical time tc,α and length lc,α scales for
flame with lean to unity global equivalence ratio.
30
species tc,α [s] lc,α [m]
CH4 3.7532 ×10−4 3.0181 ×10−4
O2 5.6150 ×10−4 3.7014 ×10−4
CO2 7.8510 ×10−4 4.5202 ×10−4
H2O 4.8870 ×10−4 3.9695 ×10−4
H2 1.2306 ×10−4 4.2889 ×10−4
O 1.0822 ×10−4 1.6096 ×10−4
OH 1.4792 ×10−4 1.8466 ×10−4
H 5.3076 ×10−5 1.9735 ×10−4
HO2 1.9941 ×10−4 1.4704 ×10−4
H2O2 1.2422 ×10−4 1.1147 ×10−4
CO 2.5253 ×10−4 2.6423 ×10−4
CH2O 1.2509 ×10−4 1.2902 ×10−4
HCO 7.5642 ×10−5 1.0609 ×10−4
CH2OH 6.1104 ×10−5 1.1767 ×10−4
CH3OH 1.2341 ×10−4 1.4089 ×10−4
CH3 6.9015 ×10−5 1.0579 ×10−4
CH3O 2.2065 ×10−4 1.5211 ×10−4
Table 3.1: The chemical time tc,α and length lc,α scales for flame with unity global
equivalence ratio ΦG = 1.00.
radical CH3 contributes the smallest length scale with a value of 1.0579×10−4 m.
Carbon dioxide CO2 again has the longest length scale. The long time and length
scales from CO2 are consistent with the thick H2-CO reaction zone [84]. Note
that the formyl radical HCO also registers a short length scale with value of
1.0609× 10−4 m.
The acoustic time scale with unity Courant-Friedrichs-Lewy (CFL) number is
calculated as ∆t = ∆x/amax = 2.74 × 10−8s (amax = 913.90 m/s). This value is
smaller than the methyl radical CH3 time scale, therefore the stability criterion
(CFL ≤ 1) is the bottleneck for a longer simulation time step.
3.4 Turbulence Initialisation
In the absence of a closed form mathematical solution, the turbulent initial con-
dition is numerically constructed. The present work uses the DNS code SENGA2
31
written by Cant [17] to initialise the turbulent flow field. It generates a peri-
odic velocity field of three-dimensional homogeneous isotropic turbulence having
a prescribed energy spectrum and satisfying the continuity constraint of an in-
compressible flow. The velocity field is first generated in the Fourier space and an
inverse Fourier transform is performed to obtain the real space velocity compo-
nents. The inverse Fourier transform is performed using the Fast Fourier Trans-
form (FFT) algorithm by Temperton [109]. Full details are given by Cant [14].
This method is based on the pseudo-spectral method pioneered by Orszag and
Patterson [80] and refined by Rogallo [95].
Batchelor and Townsend’s spectrum [102] which represents the final period of
decay of grid turbulence [5] is adopted in the present work:
E(κ) = 323
√
2
πk
κ4
κ50
exp
[
−
( κ
κ0
)2
]
(3.20)
where k is the prescribed kinetic energy and κ and κ0 = 1/(
√
2πlt) is the
wavenumber and the largest wavenumber respectively. This energy spectrum
conforms to the final period of decay of turbulence, where the inertial effects are
diminished relative to viscous process [91]. A detailed description of the imple-
mentation can be found in [16, 33, 54].
3.4.1 Initial turbulent velocity flow field in a rectangular
domain
A rectangular domain is used for the current DNS in order to capture the antic-
ipated flame wrinkling.
Using the turbulence generator from SENGA2, a cubic turbulent flow field is
produced and the vorticity contour plot of an x-y plane is presented in figure 3.2a.
The plot is coloured by the vorticity magnitude and this turbulent flow field has
turbulent velocity of u=1.25m/s and integral length scale of lt =6.64×10−4m.
This turbulent flow field is split into two sections along the stream-wise direction,
namely: section A and section B as labelled in the figure.
Utilising the periodic boundary condition of the initially cubic turbulent flow
32
01
2 3
4 5
x ×10
−3
 [m]
0
1
2
3
4
5
y
 
×
1
0
−
3
 
[
m
]
A B
√
ω
k
ω
k
[s
−1
]
0
750
1500
2250
3000
3750
4500
5250
6000
(a) Vorticity magnitude √ωkωk contour of the cubic tur-
bulent velocity flow field generator on x-y plane. The
contour plot is coloured by vorticity magnitude.
0 2 4 6 8 10
x ×10
−3
 [m]
0
1
2
3
4
5
y
 
×
1
0
−
3
 
[
m
]
B A B A
(b) Vorticity magnitude √ωkωk contour of the initial turbulent ve-
locity flow field on x-y plane.
Figure 3.2: Turbulent velocity flow field initialisation in a rectangular domain.The
contour plot is coloured by vorticity magnitude.
33
field, the initial solution is elongated by copying section B to precede section A,
and an additional section A is copied after section B, as shown in figure 3.2b. This
produces the required turbulent flow field for the DNS. This elongated turbulent
flow field in figure 3.2b is allowed to decay isotropically as recommended by
Rogallo [95] before it is used as the initial turbulent velocity flow field solution.
3.5 Boundary conditions
The Navier-Stokes Characteristic Boundary Condition (NSCBC) formulation is
used to specify the boundary conditions of the current DNS study. It is an ex-
tension from the Euler Characteristic Boundary Conditions (ECBC) technique
by taking into account the additional viscous relations of Navier-Stokes systems.
Both NSCBC and ECBC formulations are derived from the analysis of the charac-
teristic waves crossing at the boundary. This is an appealing method of specifying
the DNS boundary conditions as it provides accurate control of wave reflections
from the boundaries and the possibility of specifying non-reflecting boundary
conditions while assuring stability [48].
NSCBC boundary formulation here follows the recommendations of Suther-
land and Kennedy [104]. The formulation is an improved version of that proposed
by Poinsot and Lele [88] by taking into account the reactive source term treat-
ment at the boundary and also expanded to be used for multicomponent gases.
The implementation in SENGA2 is described by Cant [17].
3.5.1 Locally One-Dimensional Inviscid (LODI) relations
Treatment of ECBC boundary conditions for a one-dimensional problem is pre-
sented by Thompson [110] and uses the Local One-Dimensional Inviscid (LODI)
relations. The objective is to recast the evolution equations in characteristic form
[104]. The LODI equations are written in conserved variables in equation 3.21,
primitive variables in equation 3.22, and characteristic variables in equation 3.23
are shown here. For a detailed derivation of the LODI equations, the reader is
referred to [17, 104].
∂Ui
∂t +PijS
(n)
jk Lk = si (3.21)
34
∂Ui
∂t + S
(n)
ij Lj = si (3.22)
(
S(n)
)−1
ji
∂Ui
∂t + Lj = si (3.23)
where Ui and Ui are the conservation and primitive variables with the former
chosen as U = {ρui, ρ, ρe0, ρYα}T and U = {ui, ρ, P, Yα}T . P = ∂U∂U , S
(n) is the
matric of right eigenvectors, and si and si are the conservation and primitive
source terms. Note that the source terms are set to zero; it is not intended that
the flame should cross the boundary. The Lj are the wave amplitudes for the
characteristic variables and λj are the associated characteristic velocities. Waves
propagating across the boundary can be directly controlled by manipulating the
wave amplitudes Lj .
3.5.2 Inflow Boundary Condition
Hard inflow conditions are prescribed for all the simulations. The primitive
variables UBC = {ui, T, Yα} are specified at the inflow. For two- and three-
dimensional simulations, the time derivatives of these primitive variables are set to
zero (∂UBC/∂t = {0}) while for one-dimensional simulations, the time derivatives
are prescribed according to the desired time dependent function. Note that UBC
does not include density to prevent over specification of the inlet boundary condi-
tion. The last condition to be specified is the viscous condition. Following Poinsot
and Lele [88] and Sutherland and Kennedy [104], boundary normal-viscous stress,
which is one of the 4 + N boundary equations required for well-posedness is set
to:
∂τ11
∂x1
= 0 (3.24)
Variable inlet species mass fraction
The reactant inhomogeneity is introduced by varying the mass fraction of methane
CH4 and air (O2 and N2) in the reactants (at the boundary plane, x = 0). This is
done in a manner such that the total quantity of fuel and oxidiser was maintained
at a prescribed equivalence ratio and the ratio of O2 and N2 is kept uniform
35
(uniform air composition). Having a uniform air composition, the inhomogeneous
reactants can be expressed in terms of equivalence ratio, Φ.
With the desired Φ profile in mind, the mass fractions for methane, oxygen,
and nitrogen in a methane fuel-air mixture are given as a function of equivalence
ratio:
YCH4,in =
nCH4WCH4
WCH4 + 2φ
(
WO2 + 0.790.21WN2
) (3.25a)
YO2,in =
nO2WO2
φWCH4 + 2
(
WO2 + 0.790.21WN2
) (3.25b)
YN2,in =1− YCH4,in − YO2,in (3.25c)
The above equations are used to specify Yα as required by UBC .
3.5.3 Outflow Boundary Conditions
The subsonic partially non-reflecting outflow boundary condition is used here.
There is only one incoming wave at the outflow boundary. The incoming eigen-
value (λ1 = u− c) has its assorted wave amplitude L1 set as:
L1 = [σc(1−Ma2)(p− p∞)]/2Lx (3.26)
where σ = 0.287 is used for current study as suggested by Rudy and Strikw-
erda [97] and tested by Chakraborty and Cant [23], c is the speed of sound,
Ma = 1 × 10−3 and p∞ = 1 × 105Pa is the reference pressure which the pressure
at the boundary is expected to track, and Lx is the domain length normal to
the boundary. The constant σ is not set to zero as expected of a non-reflection
criteria in order to eliminate pressure drift [88]. The performance of this partially
non-reflecting boundary condition has been tested by Chakraborty [33].
3.6 The progress variable
In the analysis of premixed flames, whether laminar or turbulent, it is convenient
to track the evolution of various quantities using the reaction progress variable.
36
The progress variable, which will be a function of both space and time, mono-
tonically tracks the premixed flame evolution from zero in the reactants to unity
in the products. To retain the physical meaning of the progress variable for
quantitative analysis, it should be defined using physical quantities like species
mass fraction or temperature. There are several definitions of progress variable in
the literature, which are mainly based on normalised temperature or normalised
species mass fraction.
Since the focus of the present work is on the variation of reactant mixture
strength, the progress variable based on species mass fraction is preferred over the
temperature based counterpart. A further requirement for the progress variable
definition is to be independent towards the variation of reactant mixture strength.
Species in the unburnt gas are preferred over those in burnt gases. The fuel mass
fraction is used because it is the deficient species in the current simulations.
Following Haworth [59], the progress variable adopted in the current analysis is
defined as:
c = 1− YCH4z (3.27)
where z = zC + zH, a conserved scalar, is the sum of mixture fractions of carbon
and hydrogen elements. In the unburnt region, z = YCH4 and leads to zero
progress variable. In regions where all methane has been consumed YCH4 = 0, and
the progress variable becomes unity. This definition of progress variable ensures
that the progress variable throughout the domain will be bounded between zero
and unity.
3.6.1 The element mixture fraction and the local equiva-
lence ratio
The element mixture fraction is obtained from [59]:
zβ =
Ns∑
α=1
nβαYα
Wβ
Wα
(3.28)
where β refers to the elements (C, H, and O), Ns is the total number of species
considered (18 for the current chemical scheme), Wβ is the molecular weight of
37
element β and nβα is the number of atoms of element β in species α. With the
mixture fractions defined, the local equivalence ratio is given as [59, 78, 89]:
Φ = z1− z
1− zst
zst
(3.29)
where subscript st denotes a value at unity equivalence ratio. The value for zst is
obtained by summing zC,st = 4.1275×10−2 and zH,st = 1.3758×10−2, which gives
zst = zC,st + zH,st = 5.5033 × 10−2. Both local equivalence ratio Φ and element
mixture fractions zβ are functions of space and time.
This particular set of definitions is used because all of them can be derived
from the primary variables and the concise relationship of c, YCH4 and z provides
the means to derive the progress variable transport equation.
3.7 The progress variable transport equations
The thermochemistry of inhomogeneous reactants premixed combustion is de-
scribed with two composition variables. These two variables are mixture fraction
z and the reaction progress variable c. When the progress variable is defined
using the species mass fraction Yα, the species mass fraction can be defined as:
Yα(xi, t) = Yα [c(xi, t), z(xi, t)] (3.30)
Using the above relations, the progress variable transport equation becomes:
ρDcDt =
∂
∂xk
(
ρDc
∂c
∂xk
)
+ ω˙Y∂Yα/∂c
+ ρ∂
2Yα/∂c2
∂Yα/∂c
χc + ρ
∂2Yα/∂z2
∂Yα/∂c
χz
+ 2ρ∂
2Yα/∂c∂z
∂Yα/∂c
χz,c +
∂Yα/∂z
∂Yα/∂c
[ ∂
∂xk
ρ(Dc −Dz)
∂z
∂xk
] (3.31)
where the diffusion coefficient for the progress variable Dc is the same as the
diffusion coefficient of the species DYα (DC = DYα), while Dz is the diffusion
38
coefficient for the mixture fraction. The three χ terms are given as:
χc = Dc
∂c
∂xk
∂c
∂xk
χz = Dc
∂z
∂xk
∂z
∂xk
χz,c = Dc
∂z
∂xk
∂c
∂xk
(3.32)
The diffusion coefficient of these three scalar dissipation terms is the same as the
diffusion coefficient of the species mass fraction. Their role as either a sink or
source depends on the second order partial derivative of species mass fraction.
The progress variable transport equation provided in equation 3.31 is similar
to that derived by Bray et. al. [12] and used by Ruan et. al. [96]. This equation
has also been extended to accomodate the non-equal diffusivity between reaction
progress variable and the mixture fraction. The last term on the right hand
side of equation 3.31 arises because of the non-equal diffusivity. Note that the c
transport equation expressed in equation 3.31 is applicable in both premixed and
non-premixed combustion. The full derivation can be found in the Appendix.
The progress variable definition has been chosen in equation 3.27, therefore
the following partial derivatives of the progress variable apply:
∂YF
∂c = −z
∂YF
∂z = 1− c
∂2YF
∂c∂z = −1
∂2YF
∂c2 =
∂2YF
∂z2 = 0
(3.33)
Using equation 3.31 and 3.33, the progress variable transport equation be-
comes:
ρ∂c∂t + ρuk
∂c
∂xk
= ∂∂xk
(
ρDc
∂c
∂xk
)
+ ω˙c
+2ρDc
∂ ln z
∂xk
∂c
∂xk
+ c− 1z
[ ∂
∂xk
(
ρ(Dc −Dz)
∂z
∂xk
)] (3.34)
The extended progress variable transport equation shown in equation 3.34
39
has two extra terms on the right hand side that are absent in the homogeneous
premixed case; These two terms are the cross dissipation term and the extra
diffusion term. Note that the scalar dissipation terms in equation 3.31 have
vanished. Similar extended progress variable transport equations have been used
without the extra diffusion term by Richardson et. al. [93] and Poinsot et. al. [89].
3.7.1 Mixture fraction transport equation
For completeness, the mixture fraction transport equation is given as:
∂ρz
∂t +
∂ρukz
∂xk
= ∂∂xk
(
ρDz
∂z
∂xk
)
= ∂∂xk
(
ρDzC
∂zC
∂xk
)
+ ∂∂xk
(
ρDzH
∂zH
∂xk
)
(3.35)
where Dz is the diffusion coefficient of mixture fraction and DzC and DzH are the
diffusion coefficient for the carbon and hydrogen element mixture fractions. The
diffusion term for mixture fraction can be separated into two terms consisting of
carbon element and hydrogen element diffusion terms, however the diffusion coef-
ficient for mixture fraction is not the linear combination of carbon and hydrogen
element diffusion coefficients, Dz 6= DzC + DzH .
3.8 Displacement speed components
The displacement speed Sd is commonly used in measuring the speed of a reactive
diffusive scalar. The displacement speed can be used when a properly defined
scalar field is available.
Using level set analysis, the kinematic form of the deficient species (methane)
mass fraction equation at a fixed isosurface can be expressed as:
∂c
∂t + Vf,k
∂c
∂xk
= 0 (3.36)
where Vf,k is the absolute velocity of the iso-contour, which is the sum of both
the convective component uk and its propagation speed relative to the gas Sd nk:
Vf,k = uk + Sd nk (3.37)
40
Here, ni is the unit normal vector of the iso-contour of c directed towards the
unburnt gas, given as:
nk =
−1
|∇c|
∂c
∂xk
(3.38)
By expressing the molecular diffusion term in terms of its tangential and normal
components [43, 82]:
∂
∂xk
(
ρDc
∂c
∂xk
)
= −nk
∂
∂xk
(ρDc|∇c|)− ρD|∇c|
∂nk
∂xk
(3.39)
and using the definition from equation 3.35, an expression for the flame displace-
ment speed can thne be obtained by substituting equation 3.34 into equation 3.36
and 3.37:
Sd =
Sn
︷ ︸︸ ︷
−nk
ρ|∇c|
∂ (ρDc|∇c|)
∂xk
St
︷ ︸︸ ︷
−Dc
∂nk
∂xk
+
Sr
︷ ︸︸ ︷
ω˙c
ρ|∇c|
Scd
︷ ︸︸ ︷
2Dc
z|∇c|
∂z
∂xk
∂c
∂xk
+
Sdiff,z
︷ ︸︸ ︷
c− 1
ρ|∇c|z
∂
∂xk
(
ρ(Dc −Dz)
∂z
∂xk
)
(3.40)
It is shown that the displacement speed Sd consists of five terms. The first
three terms, which are also found in the homogeneous premixed case, are the
normal, tangential and reaction components Sn, St, and Sr respectively. The
last two terms arise from the coupling with mixture fraction, they are the cross
dissipation term Scd and the extra diffusion term Sdiff,z.
The normal component Sn can be recast to:
Sn = −nk
[∂Dc
∂xk
+ Dc
∂ ln ρ
∂xk
+ Dc
∂ ln |∇c|
∂xk
]
(3.41)
This form of Sn is useful when analysing the one dimensional simulations. The
tangential component is related to the local curvature κ as:
St = −2Dcκ (3.42)
41
where the local mean curvature κ is defined as the arithmetic average of the two
principal curvatures for the iso-surface [90]:
κ = −tr(H(c))2|∇c| =
1
2
∂nk
∂xk
(3.43)
The cross dissipation term can be expressed as:
Scd =
2Dc
z|∇c|
∂z
∂xk
∂c
∂xk
= −2Dcnk
∂ ln z
∂xk
(3.44)
In the one-dimensional case, the cross dissipation term reduces to 2Dc∂ ln z/∂xi,
which is governed by the gradient of the mixture fraction.
The extra diffusion term is related to the mixture fraction molecular diffusion
term by:
Sdiff,z =
1− c
z
∂
∂xk
(
ρ (Dz −Dc)
∂z
∂xk
)
= 1− cz
[
− ∂∂xk
(
ρDc
∂z
∂xk
)
+ ρDzDt
] (3.45)
This shows that the extra diffusion term is coupled with the total time derivative
of the mixture fraction. The extra diffusion term contribution to Sd becomes
small when it is evaluated on an iso-surface that is situated close to the products.
Should the mixture fraction is defined such that Dc − Dz is large or its spatial
gradient becomes large or both, this extra diffusion term becomes non-trivial.
3.9 Surface density function and fine-grained flame
surface density
The surface density function is introduced as the amount of surface area per unit
of volume [115]. It is defined as [73]:
σ = |∇c| (3.46)
42
The balance equation for the surface density function can be obtained by mul-
tiplying the kinematic equation 3.36 with ∂c/∂xk, then dividing by the surface
density function |∇c|. The balance equation for σ is given as:
∂σ
∂t +
∂ukσ
∂xk
=(δij − ninj)
∂ui
∂xj
σ + Sd
∂nk
∂xk
σ − ∂Sdnkσ∂xk
(3.47)
where nk = −∇c/|∇c| is the normal vector and δij is the Kronecker delta.
The fine-grained flame surface density Σ′ is defined as the surface area per
unit volume [13, 20, 90]:
Σ′ = |∇c|δ(c = c∗) (3.48)
where c∗ is the prescribed iso-surface value. The balance equation for Σ′ has
been derived by Pope [90] using the differential geometry approach. An alterna-
tive derivation is presented by Candel and Poinsot [13] using the material flow
transport theorem [6] of an elementary flame surface to volume ratio (δA/δV ).
The balance equation the Σ′ is given by:
∂Σ′
∂t +
∂ukΣ′
∂xk
=
(
δij − n∗in∗j
) ∂ui
∂xj
Σ′ + Sd
∂n∗k
∂xk
Σ′ − ∂Sdn
∗
kΣ′
∂xk
(3.49)
where the ∗ denotes values evaluated on the c = c∗ iso-surface, i.e. n∗i = − ∇c|∇c| |(c =
c∗) = nk|(c = c∗). Although the FSD and SDF transport equation are analogous,
the fine-grained FSD balance equation is valid on the c = c∗ iso-surface while the
SDF balance equation is a field equation for regions where |∇c| 6= 0.
The three terms on the right hand side of equations 3.47 and 3.49 are the:
1. Strain rate term: production of σ or Σ′ by straining.
2. Curvature term: production or destruction of σ or Σ′ by curvature.
3. Propagation term: production or destruction of σ or Σ′ by propagation.
43
Chapter 4
One-dimensional laminar flame
simulations
One-dimensional laminar flame simulations are discussed in this chapter. Two
types of laminar premixed flame simulations were performed, for the steady lam-
inar premixed flame and the unsteady laminar premixed flame.
A series of steady premixed flame simulations with different equivalence ratio
ranging from lean to stoichiometric mixtures was performed. These simulations
will form a library of reference simulations for use in comparing unsteady flame
simulations with equal or higher dimensions.
Unsteady laminar premixed flame simulations were performed to investigate
the effect of unsteady reactant mixture strength towards an initially steady lami-
nar premixed flame. The time dependence of mixture strength, measured in terms
of equivalence ratio, produces a gradient of equivalence ratio that is parallel to
the direction of laminar flame propagation. This type of mixture stratification is
called the normally stratified mixture and this type of unsteady flame is named
as the normally-stratified unstrained premixed flame. Results of normal stratifi-
cation will be used in the next chapter on discussion of two-dimensional unsteady
laminar premixed flame simulations.
44
4.1 One-dimensional steady laminar flame cal-
culations at different equivalence ratio
The present investigation will see the flame propagating into a varying equivalence
ratio field, therefore it is necessary to assess the chemical mechanism performance
on producing accurate flame behaviour for a desired range of equivalence ratio.
One-dimensional calculations of steady laminar flames were performed for global
equivalence ratio ΦG ranging from ΦG = 0.50 to 1.00 with increments of 0.05.
The inlet temperature and pressure for all simulations were prescribed at T=300K
and P=1.0×105Pa.
The laminar flame speed, which is the speed at which the laminar flame is
propagating with respect to the unburnt gas in a one-dimensional geometry, is
defined as:
SL =
−1
ρuYCH4,u
∫ ∞
−∞
ω˙CH4dx (4.1)
where ρu and YCH4,u denote reactant density and methane mass-fraction respec-
tively. The laminar flame thickness based on temperature is defined as:
δL =
Tb − Tu
max |∇T | (4.2)
where Tu and Tb denotes unburnt and burnt gas temperature and max |∇T | is
the maximum temperature gradient across the flame.
The chemical time scale is calculated using:
τc =
δL
SL
(4.3)
The laminar flame speed SL, laminar flame thickness δL, and chemical time scale
τc for each steady laminar flame calculations are tabulated in table 4.1. The
calculated laminar flame speed is also compared with the experimentally obtained
value by Vagelopoulos et. al. [112] in figure 4.1. The flame speed obtained
from the numerical calculation clearly captures the effect of equivalence ratio
on the laminar flame speed and agrees fairly well with the experimental values.
The numerical flame speed is systematically overestimated compared with the
45
Table 4.1: Summary of the steady laminar flame calculations
ΦG SL ×10−1 [m/s] δL ×10−4 [m] τc ×10−3 [s]
1.00 3.8350 4.3009 1.1215
0.95 3.6637 4.4385 1.2115
0.90 3.4069 4.6518 1.3654
0.95 3.0908 4.9256 1.5936
0.80 2.7353 5.3162 1.9436
0.75 2.3351 5.8609 2.5099
0.70 1.9284 6.6584 3.4528
0.65 1.5183 7.8880 5.1953
0.60 1.1234 9.8989 8.8116
0.55 0.7614 13.5636 17.8140
0.50 0.4649 20.8587 44.8371
experiments with mean deviation of 0.0239 m/s.
46
0.5 0.6 0.7 0.8 0.9 1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
S L
[m
/s
]
ΦG
 
 
numerical
experiment
deviation
Figure 4.1: Comparison between calculated laminar flame speed and experimen-
tally measured values by Vagelopoulos et. al. [112]. The deviations are plotted
in green.
4.2 Equivalence ratio profile across laminar flame
The local equivalence ratio ΦL does not remain equal to the global equivalence
ratio ΦG across the flame and it shows a dip within the flame. This dip is consis-
tently manifested in all the laminar flame calculations with different equivalence
ratio that are above the lower flammability limit. The normalised local equiv-
alence ratio profile across the flame is plotted in figure 4.2. The normalised
equivalence ratio Φn is the local equivalence ratio divided by its global value,
Φn = ΦL/ΦG. The normalised equivalence ratio profile for each flame has been
shifted to have its local equivalence ratio minimum situated at zero and the spatial
coordinate is normalised by the relevant flame thickness xΦ,n = [x−x|min(Φ)]/δL.
The value of min(Φn) is the lowest for the ΦG = 0.50 flame, which is also
the most lean flame considered here. The value of min(Φn) increases for the
higher ΦG flame. The dip in the carbon element mixture fraction zC across the
flame is the main cause of the local equivalence ratio dip. Figure 4.3 shows the
47
−10 −5 0 5 10 15 20
0.9
0.92
0.94
0.96
0.98
1
1.02
Φ n
xΦ,n
 
 
ΦG = 1.0
ΦG = 0.9
ΦG = 0.8
ΦG = 0.7
ΦG = 0.6
ΦG = 0.5
Figure 4.2: Profile of local equivalence ratio normalised by the global equivalence
ratio. Spatial variable is normalised by laminar flame thickness.
normalised carbon element mixture fraction across the flame. The normalised
carbon element mixture fraction is the carbon element mixture fraction divided
by its reactant value zC,n = zC/zC,u. The spatial axis has been shifted to have
zero situated at the minimum Φn and normalised by the laminar flame thickness,
similar to figure 4.2. It can be observed that the spatial position of min zC,n is
the same as that of min(Φn).
An increase of Φn towards the unburnt gas (−3 < xΦ,n ≤ 0 and c > 0.01)
is observed for flames with ΦG ≥ 0.7. Also, the value of Φn is seen to return to
unity (Φn = 1) towards the burnt region after the minimum at a different spatial
gradient depending on the value of ΦG. These two observations can be attributed
to the distribution of hydrogen element mixture fraction zH across the flame.
The normalised hydrogen element mixture fraction zH,n = zH/zH,u is plotted in
figure 4.4. The value of zH,n is above unity (xΦ,n ≈ −1) in the same region where
the increase of Φn is observed (refer to figure 4.2). The highest spike of zH,n
happens in the stoichiometric flame and decreases for leaner flames. The value
of zH,n is seen dipping below unity for xΦ,n > 0 and slowly returning to zH,n = 1
48
−10 −5 0 5 10 15 20
0.85
0.9
0.95
1
1.05
z C
,n
xΦ,n
 
 
ΦG = 1.0
ΦG = 0.9
ΦG = 0.8
ΦG = 0.7
ΦG = 0.6
ΦG = 0.5
Figure 4.3: Normalised carbon element mixture fraction.
−10 −5 0 5 10 15 20
0.95
1
1.05
1.1
z H
,n
xΦ,n
 
 
ΦG = 1.0
ΦG = 0.9
ΦG = 0.8
ΦG = 0.7
ΦG = 0.6
ΦG = 0.5
Figure 4.4: Normalised hydrogen element mixture fraction.
49
towards the burned gas.
Figure 4.3 and 4.4 indicate that the preferential diffusion of species produces
the equivalence ratio dip in figure 4.2. This equivalence ratio dip has also been
reported in experiments by Barlow et. al. [3].
Unsteady flames, which is discussed in section 4.4 pose the ambiguity of flame
speed definition (refer to figure 4.6) specifically for the denominator of equa-
tion 4.1. This can be alleviated by using the displacement speed evaluated at
a specific isosurface to relate to the instantaneous behaviour of the flame. An
isosurface corresponding to the maximum reaction rate in the laminar flame has
been taken to represent the flame [35, 36, 58]. The current definition of progress
variable uses the methane mass fraction. For the range of equivalence ratio con-
sidered (refer to section 4.3), the isosurface closest to the methane maximum
reaction rate is located at c = 0.90.
The existence of the local equivalence ratio dip makes it necessary to estab-
lish the relationship between local equivalence ratio measured at progress variable
equal to 0.90, Φ0.90 and the global equivalence ratio ΦG. This is plotted in fig-
ure 4.5 and the graph shows a linear relationship between Φ0.90 and ΦG even with
the existence of the local equivalence ratio dip across the flame.
The linear relationship between Φ0.90 and ΦG provides for straightforward eval-
uation of ΦG for a given Φ0.90. This method of mapping a local equivalence ratio
onto the corresponding global equivalence ratio alleviates the need for numerical
interpolation for evaluating an instantaneous flame equivalence ratio. The ad-
vantage of this linear relationship is apparent when the need to compare between
steady and unsteady flame arises.
For demonstration, three local equivalence ratio snapshots from a unsteady
laminar flame calculation case S00 (refer to section 4.5) are plotted in figure 4.6.
The spatial axis xn,0.90 has been shifted so that the c = 0.90 iso-surface is situated
at zero and the profile is normalised with the laminar flame thickness, according
to: xn,0.90 = (x − x|c = 0.90)/δL. An equivalence ratio dip is observed in all the
three snapshots and being an unsteady flame, the reactants and products have
a non-zero local equivalence ratio gradient ∇ΦL 6= 0. This has made evaluation
of the boundaries of the flame ambiguous and interpolation between these two
boundaries to obtain the expected or equivalent global equivalence ratio will give
50
0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
ΦG
Φ 0
.9
0
Figure 4.5: Comparison of the local equivalence ratio measured at c = 0.90 with
the global equivalence ratio.
rise to numerical interpolation error.
With an established relationship between local equivalence ratio at c = 0.90,
Φ0.90 and its expected global equivalence ratio ΦG, comparison between steady
state and unsteady laminar flames can be made definitively. This thesis will
present results in terms of Φ0.90 unless mentioned otherwise.
51
−10 −5 0 5 10 15 20 25 30
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
xn,0.90
Φ L
 
 
t0
t1 = t0 + δt
t2 = t1 + δt
Figure 4.6: Three instantaneous local equivalence ratio profiles of the unsteady
flame calculation.
4.3 The properties of progress variable
The progress variable is a monotonic function that describes the progress from
reactants to products across the laminar premixed flame. The current progress
variable definition is based on fuel species mass fraction (refer to equation 3.27).
In order to examine the position of maximum methane reaction rate max|ω˙CH4|
in terms of progress variable c for global equivalence ratio ranging from ΦG=0.50
to 1.00, the locus of max|ω˙CH4 | is plotted on ΦG against c in figure 4.7. Also
plotted are the position of local equivalence ratio minimum min(ΦL), tempera-
ture gradient maximum max(∇T ), progress variable gradient maximum max(∇c),
heat release rate maximum magnitude, max |Qhr|, and reduced temperature at
Tn=0.60, 0.65, 0.70, and 0.75. The reduced temperature is defined as Tn =
(T − Tu)/(Tb − Tu), which is also the temperature-based progress variable, and
Qhr =
∑
α hαω˙α [W/m3] is the heat release rate.
It is observed from figure 4.7 that the local equivalence ratio minimum is
52
0.5 0.6 0.7 0.8 0.9 1
0.4
0.5
0.6
0.7
0.8
0.9
1
ΦG
c
 
 
min(ΦL)
max |∇T |
max |ω˙CH4 |
max |∇c|
Tn = 0.60
Tn = 0.65
Tn = 0.70
Tn = 0.75
max(Qhr)
Figure 4.7: Position of different flame properties plotted in the progress variable
space.
located at between c = 0.48 to c = 0.60 for flames with ΦG = 0.5 to ΦG = 1.0.
The locus moves forwards towards the reactants for leaner flames. This behaviour
is due to the effect of preferential diffusion, where a leaner flame will produce a
higher proportion of hydrogen radicals and these will be preferentially diffused
towards the reactants due to the lower Lewis number of the H radical.
Since the progress variable defined here is fuel (methane) based, the progress
variable profile is generally shifted slightly towards the reactants, as observed
from the Tn loci. The locus for progress variable gradient maximum, max |∇c| is
located around c = 0.65 with small variation across the different global equiva-
lence ratio ΦG.
The loci for max |Qhr| and max |ω˙CH4| lie in the region of c ≥ 0.90, and tend
towards c = 1.00 at lower ΦG. The locus for max |Qhr| is always larger than that
of max |ω˙CH4|. Consequently, the position of max |Qhr| is situated closer to the
product side of the flame than the max |ω˙CH4| counterpart. From figure 4.7, it is
observed that the positions of max |Qhr| and max |ω˙CH4| are close to c = 0.90 for
the range of ΦG=0.50 to 1.00.
53
4.4 Unsteady One-dimensional Laminar Flame
Simulations
Unsteady one-dimensional laminar flame simulations were performed to investi-
gate the influence of normally stratified reactant on laminar premixed flames.
An initially homogeneous premixed laminar flame is allowed to propagate into
reactants with time varying mixture strength, measured using the equivalence
ratio.
Two types of time varying equivalence ratio were prescribed, one that has
a time-history of a sinusoidal function, while another has a time-history of the
complementary error function. For the sinusoidal case, the time history for local
equivalence ratio at c = 0.90 is plotted against normalised time tn in figure 4.8.
Similarly, the time history of the local equivalence ratio at c = 0.90 for all the
complementary error function cases is plotted against tn in figure 4.9. Note that
the normalised time tn is the simulation elapsed time te normalised by relevant
laminar flame time scale τc.
0 10 20 30 40 50 60
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
tn
Φ 0
.9
0
Figure 4.8: Local instantaneous equivalence ratio case S00.
54
0 5 10 15 20 25 30 35 40
0.5
0.6
0.7
0.8
0.9
1
tn
Φ 0
.9
0
 
 
R00
R01
R02
R03
Figure 4.9: Local instantaneous equivalence ratio for cases R00, R01a, R02, and
R03.
4.4.1 Unsteady laminar flame with sinusoidal inlet equiv-
alence ratio time history
Case S00 is an unsteady laminar flame simulation that has its inlet equivalence
ratio given as:
Φin(tn) = ΦG + Φ′in sin (2πfinτctn) (4.4)
where ΦG and Φ′in are the time-averaged global equivalence ratio and the inlet
equivalence ratio fluctuation amplitude. The reactant mixture strength pertur-
bation frequency fin is prescribed at the inlet and is given as:
fin =
1
6π
1
τc
(4.5)
During the convection of the fluctuating mixture to the flame, the amplitude
of the equivalence ratio perturbation that the flames experience is attenuated
because of diffusion. This frequency ensures that the prescribed sinusoidal sig-
55
nal of equivalence ratio is retained when convecting towards the flame from the
inlet boundary. Note that the frequency is about one twentieth of the inverse
characteristic timescale τc.
The normalised time tn is defined as the elapsed time te normalised by the
time scale τc of a laminar flame with equivalence ratio similar to the time-averaged
global equivalence ratio ΦG = 0.80. Simulation parameters for case S00 are
tabulated in tables 4.2.
Table 4.2: The inlet boundary parameters prescribed for the unsteady laminar
flame simulation with sinusoidal incoming reactant mixture strength.
Case ΦG Φ′G,in fin [Hz] τc [s]
S00 0.80 0.20 27.27 1.944×10−3
In figure 4.8, although the equivalence ratio at the inlet was set to behave
sinusoidally in time, the local equivalence ratio at c = 0.90 is not sinusoidal. This
is because the flame is not forced to respond to the inlet changes of Φ, but it
is allowed to propagate freely into the time-varying mixture, hence the time-lag
between inlet and c = 0.90 is expected.
Case S00 is divided into four sections as colour coded in figure 4.8: blue,
green, red, and cyan sections. These sections are defined by the trajectory of the
instantaneous local equivalence ratio at iso-surface c = 0.90. The blue and red
sections have negative ∂Φ0.90/∂t while the green and cyan sections have positive
∂Φ0.90/∂t. The blue section will be referred to as section 00, the red as section
01, the green as section 02, and the cyan as section 03 (refer to table 4.3).
The sinusoidal nature of the incoming mixture strength time history has pro-
vided a response in the flame that encompasses flame history with both stoichio-
metric to lean and vice versa. Section 00 and section 02 are sections where the
flame is propagating into a leaner mixture while section 01 and section 03 are
sections where the flame is propagating into a richer mixture.
The simulation reference frame has a constant underlying flow speed while
the unsteady flame is allowed to propagate freely into the time-varying mixture
strength.
The maximum max (Φ0.90) and minimum min (Φ0.90) equivalence ratio mea-
56
Table 4.3: S00 section parameters. The starting value of Φ0.90 for each section is
denoted by ().
section colour ∆te/τc max (Φ0.90) min (Φ0.90) max
(∂Φ0.90
∂t
)
[s−1]
00 blue 10.908 (0.9422) 0.5842 -36.04
01 green 11.834 0.9293 (0.5842) 30.98
02 red 10.717 (0.9293) 0.5949 -33.19
03 cyan 11.484 0.9191 (0.5949) 28.98
sured on the c = 0.90 iso-surface for each section are tabulated in table 4.3. The
residence time ∆te for each section was normalised by τc. The deviation of the
adjacent peaks or troughs of Φ0.90 is less than 2% of the time-averaged Φ0.90. The
time spent in section 02 is about 2% shorter than in section 00. Similarly, the
time spent in section 03 is about 3% less than in section 01. The reduction of
∆te is expected as the amplitude of each section decays with time.
Analysis of case S00 in terms of displacement speed will be discussed in sec-
tion 4.5.
4.4.2 Unsteady laminar flame simulation with complimen-
tary error function inlet equivalence ratio time his-
tory
The R cases, denoting inverse ramp, have an inlet equivalence ratio prescribed
to adhere to a complementary error function in time. The inlet equivalence ratio
for the R cases are given as:
Φin(tn) = ∆Φ
(
1− erf
(
A0tn −
A0
2τc
))
+ ΦfinalG (4.6)
where ∆Φ is the difference between initial ΦinitG and final ΦfinalG global equivalence
ratio: ∆Φ = ΦinitG − ΦfinalG . The time factor A0 controls the steepness of the
inverse ramp. The normalised time tn is the simulation elapsed time normalised
by the laminar flame time corresponding to unity equivalence ratio τc = τc|(ΦG =
57
1.00) = 1.1215× 10−3s. The simulations parameters are tabulated in table 4.4.
Table 4.4: The simulation parameters prescribed for unsteady laminar flame sim-
ulations with complementary error function inlet equivalence ratio time history.
Case ΦinitG ΦfinalG ∆Φ A0 τc [s] max
(
∂Φ0.90
∂t
)
[s−1]
R00 1.00 0.60 0.40 4 1.1215 × 10−3 -115.81
R01 1.00 0.70 0.30 4 1.1215 × 10−3 -84.624
R01b 1.00 0.70 0.30 2 1.1215 × 10−3 -85.022
R01c 1.00 0.70 0.30 1 1.1215 × 10−3 -77.755
R01d 1.00 0.70 0.30 8 1.1215 × 10−3 -83.383
R02 1.00 0.80 0.20 4 1.1215 × 10−3 -56.441
R03 1.00 0.90 0.10 4 1.1215 × 10−3 -28.471
All the R cases are stoichiometric flames that were allowed to propagate freely
into leaner reactants. Two sets of R cases are presented: Cases R00, R01, R02,
and R03 have similar time factor A0 but a different change of equivalence ratio
magnitude ∆Φ. Cases R01, R01b, R01c, and R01d have different time factor A0
but the same of change of equivalence ratio magnitude ∆Φ.
The maximum temporal gradient of equivalence ratio Φ0.90 for all the R cases
is tabulated in table 4.4. For cases R00, R01, R02, and R03, it is observed that
an increase in ∆Φ also increases the magnitude of max (∂Φ0.90/∂t) even when A0,
which controls the inlet equivalence ratio temporal gradient ∂Φin/∂t of the time
varying signal, is held constant.
The effect of inlet equivalence ratio temporal gradient ∂Φin/∂t towards the
laminar flame is investigated in cases R01, R01b, R01c, and R01d, where they
have similar ∆Φ = 0.30 but different time factors A0. Case R01b and R01c
have value of A0 that are half and quarter of that of case R01 respectively. The
value of A0 for case R01d is double that of case R01. It is observed that the
maximum temporal gradient of equivalence ratio max (∂Φ0.90/∂t) for case R01,
R01b, and R01d are similar with less than 1% difference between their values.
However, there is a significant difference, about 10%, in max (∂Φ0.90/∂t) value for
case R01c compared to the other three cases.
58
4.5 Temporal analysis of case S00 displacement
speed
The instantaneous time history of Sd and Φ0.90 for case S00 is plotted in fig-
ure 4.10. The instantaneous Sd has the same temporal behaviour as the local
instantaneous Φ0.90 and they are both cyclic and in phase. The time averaged
displacement speed, Sd = 1∆te
∫
Sddt as well as the time averaged local equiv-
alence ratio, Φ0.90 = 1∆te
∫
Φ0.90dt for each section of case S00 is tabulated in
table 4.5. The laminar flame displacement speed Sd
(
Φ0.90
)
corresponding to the
time averaged equivalence ratio is used to normalise the time averaged unsteady
laminar flame displacement speed. The normalised difference between Sd and
Sd
(
Φ0.90
)
and the normalised residence time ∆te/τc is tabulated in table 4.5.
0 20 40 60
0.5
1
1.5
2
2.5
tn
S d
[m
/s
]o
rΦ
0.
90
 
 
Sd(00)
Sd(01)
Sd(02)
Sd(03)
Φ0.90(00)
Φ0.90(01)
Φ0.90(02)
Φ0.90(03)
Figure 4.10: Sd and Φ0.90 against elapsed time te for case S00.
It can be observed that the residence time for a flame that propagates into
leaner mixture i.e. sections 00 and 02, is about 10% smaller than that for a flame
59
Table 4.5: Time averaged properties for case S00.
section ∆te/τc Sd [m/s] Φ0.90 Sd
(
Φ0.90
)
[m/s] ∆nSd
00 10.908 1.2290 0.7167 1.2449 -0.01281
01 11.834 1.2014 0.7120 1.2211 -0.01618
02 10.717 1.2593 0.7221 1.2705 -0.01003
03 11.484 1.2337 0.7184 1.2531 -0.01555
All 44.943 1.2301 0.7172 1.2472 -0.01372
that propagates into richer mixture i.e. sections 01 and 03. The value of Sd for
section 00 is higher than for section 01. Similarly, the value of Sd for section 02
is higher than for section 03. This shows that the Sd value is higher for a flame
that is propagating into a leaner mixture compared to the opposite.
The observed time averaged pattern of Sd in each section was also observed in
the Φ0.90 value, where the Φ0.90 value for section 00 is larger than for section 01,
and that for section 02 is larger than for section 03. This particular observation
of Φ0.90 is due to the shorter flame residence time at higher Φ0.90 compared to a
longer residence time at lower Φ0.90.
The normalised difference of Sd and Sd
(
Φ0.90
)
is given as
∆nSd =
Sd
Sd
(
Φ0.90
) − 1 (4.7)
These normalised differences have negative values for all sections. Having a higher
Sd in section 00 and 02 compared to section 01 and 03, the normalised difference
between Sd and Sd
(
Φ0.90
)
is larger in section 00 and 02 compared to section 01
and 03.
Consistently, the value of ∆Sd over all four sections, or over these two cycles,
also has a small negative value, about -1.4% difference. This indicates that the
unsteady flame in case S00 has a lower effective displacement speed compared to
the displacement speed that corresponds to its time averaged equivalence ratio.
60
4.5.1 Displacement speed Sd analysis for case S00
The displacement speed is a linear combination of five components as shown in
equation 3.40. In the one-dimensional geometry, the tangential component St is
zero and the displacement speed can be reduced to:
Sd = Sn + Sr + Scd + Sdiff,z (4.8)
where Sn, Sr, Scd, and Sdiff,z are the normal, reaction, cross dissipation and extra
diffusive components of the displacement speed Sd.
The displacement speed for case S00 over the two periods of varying equiv-
alence ratio is plotted against local instantaneous equivalence ratio Φ0.90 in fig-
ure 4.11. Steady laminar flame displacement speeds are used as the reference and
they are also plotted as a dashed line in figure 4.11.
The sum of Sn and Sr is the main contributor to Sd in case S00. The cross
dissipation term Scd contributes about 10% to Sd while the sum of the extra
diffusive term contributes about 2% towards Sd. Unsteady laminar flame Sd, Sn,
and Sr trajectories have small deviations from the reference values. However,
the cross dissipation term Scd is observed to have significant deviation from the
reference value. The deviation of the extra diffusive term Sdiff,z from the reference
value remains small.
It is interesting to note that the small deviation curves observed in Sd (refer
to figure 4.11a) form a hysteresis loop. The hysteresis loop becomes clear when
the deviation of Sd is plotted as shown later.
The unsteady flame displacement speed of case S00 can be compared to the
laminar flame counterpart. Case S00 has time varying Φ0.90 as shown in fig-
ure 4.10. Therefore, the value of Sd for case S00 (solid line in figure 4.11) is
compared with the interpolated values of the steady laminar flame counterpart
for a similar instantaneous value of Φ0.90 (dashed line in figure 4.11). For example
at time t0, case S00 has local equivalence ratio of Φ0.90,0 and an instantaneous
displacement speed of Sd,0. The deviation of displacement speed ∆Sd for case
S00 at t0 is given as:
∆Sd = Sd,0 − Srefd
∣
∣
∣Φ0.90,0 (4.9)
61
0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
Φ0.90
S d
 
 
00
01
02
03
steady
(a) Sd
0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
Φ0.90
S n
+
S r
 
 
00
01
02
03
steady
(b) Sn + Sr
0.5 0.6 0.7 0.8 0.9 1
−6
−5.5
−5
−4.5
−4
−3.5
−3
−2.5
−2
−1.5
Φ0.90
S n
 
 
00
01
02
03
steady
(c) Sn
0.5 0.6 0.7 0.8 0.9 1
2
3
4
5
6
7
8
Φ0.90
S r
 
 
00
01
02
03
steady
(d) Sr
0.5 0.6 0.7 0.8 0.9 1
0.08
0.1
0.12
0.14
0.16
0.18
0.2
Φ0.90
S c
d
 
 
00
01
02
03
steady
(e) Scd
0.5 0.6 0.7 0.8 0.9 1
−6
−5
−4
−3
−2
−1
0
1
2
x 10−3
Φ0.90
S d
iff
,z
 
 
00
01
02
03
steady
(f) Sdiff,z
Figure 4.11: Displacement speed and its components for case S00
62
where Srefd |Φ0.90,0 is the displacement speed of a steady laminar flame corre-
sponding to Φ0.90,0 that was obtained through interpolation of the laminar flame
speed library (refer to table 4.1). A similar comparison method is used for other
quantities e.g. deviation of normal component ∆Sn = Sn − Srefn .
A positive deviation of ∆Sd indicates that the unsteady flame has higher
Sd at the given value of Φ0.90 compared to the reference steady laminar flame
counterpart. Similarly, a negative ∆Sd indicates that the unsteady flame has
lower Sd relative to the reference.
4.5.2 Displacement speed deviation ∆Sd analysis for case
S00
The deviation of displacement speed defined in equation 4.9 is normalised by the
corresponding Srefd giving:
∆nSd =
Sd − Srefd
Srefd
= ∆Sd
Srefd
(4.10)
The deviations of Sd and its components for case S00 are plotted against the local
equivalence ratio in figure 4.12. The displacement speed deviation for case S00 is
observed to be non-zero; moreover, the trajectory forms a continuous closed loop.
This is consistent with figure 4.11a where the hysteresis of displacement speed is
observed.
From figure 4.11a, the normalised deviations of displacement speed ∆nSd for
section 00 and 02 are observed to have positive values. Sections 00 and 02 are
sections where the unsteady flame propagates into a leaner mixture. On the other
hand, in sections 01 and 03, where the unsteady flame propagates into a richer
mixture, ∆nSd is observed to have negative values. This is in agreement with the
findings of Pires Da Cruz et. al. [85], Kang and Kyritsis [70], and Richardson et.
al. [93].
The major contributors to the Sd deviations are the normal ∆nSn, reaction
∆nSr, and cross dissipation ∆nScd terms as shown in figure 4.12c, 4.12d, and 4.12e
respectively. The contributions from ∆nSdiff,z are about two orders smaller than
the major contributors. Although Scd has only about 10% contribution towards
63
0.5 0.6 0.7 0.8 0.9 1
−0.03
−0.02
−0.01
0
0.01
0.02
0.03
Φ0.90
∆
S d
/S
re
f
d
 
 
00
01
02
03
(a) ∆Sd/Sd,ref
0.5 0.6 0.7 0.8 0.9 1
−0.04
−0.03
−0.02
−0.01
0
0.01
0.02
0.03
0.04
Φ0.90
∆
(S
n
+
S r
)/
Sr
ef d
 
 
00
01
02
03
(b) ∆(Sn + Sr)/Sd,ref
0.5 0.6 0.7 0.8 0.9 1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
Φ0.90
∆
S n
/S
re
f
d
 
 
00
01
02
03
(c) ∆Sn/Sd,ref
0.5 0.6 0.7 0.8 0.9 1
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
Φ0.90
∆
S r
/S
re
f
d
 
 
00
01
02
03
(d) ∆Sr/Sd,ref
0.5 0.6 0.7 0.8 0.9 1
−0.03
−0.02
−0.01
0
0.01
0.02
0.03
Φ0.90
∆
S c
d/
Sr
ef d
 
 
00
01
02
03
(e) ∆Scd/Sd,ref
0.5 0.6 0.7 0.8 0.9 1
−3
−2
−1
0
1
2
3
4
x 10−4
Φ0.90
∆
S d
iff
,z
/S
re
f
d
 
 
00
01
02
03
(f) ∆Sdiff,z/Sd,ref
Figure 4.12: Deviations of the displacement speed and its components from the
steady state counterpart for case S00
64
the Sd value, ∆nScd has similar magnitude to ∆nSd.
The two-loop pattern of ∆nSd in figure 4.12a, where one loop is found at 0.6 <
Φ0.90 < 0.7, and the another loop at 0.7 < Φ0.90 < 0.95, is directly influenced
by the ∆n(Sn + Sr) pattern in figure 4.12b. From figure 4.12b, the ∆n(Sn + Sr)
loop at 0.6 < Φ0.90 < 0.7 can be attributed to ∆nSr, where its values within
0.6 < Φ0.90 < 0.7 for section 00 and 02 are observed to decrease rapidly while
∆nSr is also observed not to have a monotonic relation with Φ0.90 in the range of
0.6 < Φ0.90 < 0.7 for section 01 and 03. The hysteresis loops observed in ∆nScd
and ∆nSdiff,z are observed to have larger magnitudes in the low equivalence ratio
regions.
It is observed that the displacement speed trajectories for flames propagating
into leaner mixture (sections 00 and 02) are not the same as for flames propagating
into richer mixtures (sections 01 and 03). This difference is consistent with the
difference in residence time observed between these two groups (c.f. table 4.5).
Sections 00 and 02 have relatively smaller residence time compared to sections 01
and 03. The maximum deviation magnitude for Sd is larger for sections 00 and
02 compared to sections 01 and 03. This indicates that a more abrupt change
in equivalence ratio will result in a greater departure from the reference steady
laminar values.
The normalised deviation of the cross dissipation term ∆nScd for sections 01
and 03 are observed to be the mirror image of that in section 00 and 02 about
∆nScd = 0. The trajectrories of ∆nScd is not significantly affected by the residence
time difference between the two groups of sections.
Since the cross dissipation term Scd in the one-dimensional context is inde-
pendent of |∇c|, and the normal component Sn and the reaction component Sr
are related to |∇c| (refer to equation 3.40), the expansion and contraction rate of
the flame thickness (δF ∼ 1/|∇c|) influences the non-symmetric trajectory loop
of ∆nSd shown in figure 4.11a.
65
4.6 Temporal analysis of R cases displacement
speed
The instantaneous displacement speed Sd and local equivalence ratio time history
for cases R00, R01, R02, and R03 are plotted in figure 4.13 and those for cases
R01, R01b, R01c, and R01d are plotted in figure 4.14. From figure 4.13, the
displacement speeds Sd for cases R00, R01, R02, and R03 are observed to have
similar behaviour to their local equivalence ratio Φ0.90, where Sd is observed to
adhere to the complementary error function trend. The displacement speed Sd for
cases R01, R01b, R01c, R01d also behaves like their Φ0.90, where a complementary
error function trend is observed in figure 4.14.
0 10 20 30 40
0.5
1
1.5
2
2.5
tn
S d
[m
/s
]o
rΦ
0.
90
 
 
Sd(R00)
Sd(R01)
Sd(R02)
Sd(R03)
Φ0.90(R00)
Φ0.90(R01)
Φ0.90(R02)
Φ0.90(R03)
Figure 4.13: The time histories of Sd and Φ0.90 for cases R00, R01, R02, and R03
The equivalence ratio transition time ∆tn(∆Φ0.90) is the time taken for the
laminar flame equivalence ratio to change from its initial to final value and is
measured by the time taken for the flame to change from 99% of initial equiv-
alence ratio (0.99×Φinit0.90) to 1% above the final equivalence ratio (1.01×Φfinal0.90 ).
Similarly, the displacement speed transition time ∆tn(∆Sd) is the time taken for
66
0 10 20 30 40
0.5
1
1.5
2
2.5
tn
S d
[m
/s
]o
rΦ
0.
90
 
 
Sd(R01)
Sd(R01b)
Sd(R01c)
Sd(R01d)
Φ0.90(R01)
Φ0.90(R01b)
Φ0.90(R01c)
Φ0.90(R01d)
Figure 4.14: The time histories of Sd and Φ0.90 for cases R01, R01b, R01c, and
R01d
the laminar flame displacement speed to change from its initial to final value, and
it is measured by the time taken for the flame Sd to change from 0.99×Sinitd to
1.01×Sfinald .
Along with the transition times, the maximum temporal gradient of the pre-
scribed inlet equivalence ratio Φin, the local equivalence ratio Φ0.90 and the corre-
sponding displacement speed Sd for each of the R cases are tabulated in table 4.6.
For cases R00, R01, R02, and R03, the maximum inlet equivalence ratio
temporal gradient max (∂Φin/∂t) is observed to be at least twice the magnitude of
the local equivalence ratio temporal gradient max (∂Φ0.90/∂t). This indicates that
the freely propagating flame is responding to the time-varying reactant mixture
strength at its limiting frequency. Although R00 has the steepest max (∂Φ0.90/∂t),
which implies the most rapid changes of Φ0.90, the Φ0.90 transition time is found
to be the longest. Case R03 has the smallest max (∂Φ0.90/∂t), and it has the
shortest transition time ∆tn(∆Φ0.90).
The displacement speed Sd response for cases R00, R01, R02, and R03 is
similar to the Φ0.90 response, where case R00 has the largest max (∂Sd/∂t) while
67
Table 4.6: The equivalence ratio and displacement speed transition time for all
the R cases.
Case R00 R01 R01b R01c R01d R02 R03
A0 4 4 2 1 8 4 4
∆Φ 0.4 0.3 0.3 0.3 0.3 0.2 0.1
max
(
∂Φin
∂t
)
[s−1] -312.82 -312.41 -158.48 -79.90 -610.48 -260.13 -156.01
max
(
∂Φ0.90
∂t
)
[s−1] -115.81 -84.62 -85.00 -77.76 -83.38 -56.44 -28.47
∆tn(∆Φ0.90) 12.30 7.14 7.03 6.69 7.08 5.11 3.55
max
(
∂Sd
∂t
)
[s−1] -558.96 -385.41 -390.07 -363.99 -379.19 -232.94 -98.00
∆tn(∆Sd) 13.39 8.72 8.25 7.53 8.49 5.97 4.04
case R03 has the smallest max (∂Sd/∂t) which reflects the max (∂Φ0.90/∂t) trend.
Case R00 has the largest ∆tn(∆Φ0.90) and also the longest Sd transition time
∆tn(∆Sd), while case R03 has the smallest ∆tn(∆Φ0.90) and the shortest duration
of ∆tn(∆Sd).
For cases R01, R01b, R01c, and R01d, ∆Φ was fixed but A0 was varied. It
can be observed that max (∂Φin/∂t) has a direct correlation with A0 as expected.
However, the response of max (∂Φ0.90/∂t) has been limited to a maximum value,
at about max (∂Φ0.90/∂t) ≈ −84.50 s−1 as observed in cases R01, R01b, and
R01d. Since the equivalence ratio temporal gradient of case R01c is lower than
this limiting maximum max (∂Φin/∂t) < −84.50 s−1, its local equivalence ratio
temporal gradient response is similar to the inlet equivalence ratio temporal gra-
dient. i.e. max (∂Φin/∂t) ≈ max (∂Φ0.90/∂t). A similar trend is observed for the
displacement speed temporal gradient, where cases R01, R01b, and R01d have
similar values but the case R01c value of max (∂Sd/∂t) is observed to be 15%
smaller. This indicates that for a given change of equivalence ratio ∆Φ, there is
a limiting rate of change of Φ0.90 or Sd for a flame that propagates freely into the
inhomogeneous reactants.
The transition time for Φ0.90 is again observed to have a negative relation
with max (∂Φ0.90/∂t). Cases R01, R01b, and R01d have larger max (∂Φ0.90/∂t)
than case R01c, and they also have a larger tn(∆Φ0.90) than case R01c. A similar
68
relation of tn(∆Sd) to max (∂Sd/∂t) is observed for the R01 cases.
Note that the max (∂Φ0.90/∂t) for the R cases is at least two times larger than
that for the S cases counterpart.
4.7 Displacement speed deviation ∆Sd analysis
for cases R00, R01, R02, and R03
The normal component of the displacement speed Sn can be written as:
Sn =
∂Dc
∂x + Dc
∂ ln ρ
∂x + Dc
∂ ln |∇c|
∂x (4.11)
and in the one-dimensional context, the first extra term of the displacement speed
can be written as:
Scd = 2Dc
∂ ln z
∂x (4.12)
where Dc = DCH4 is the effective diffusivity of progress variable.
The behaviour of displacement speed Sd and its components across the equiv-
alence ratio space Φ0.90 for cases R00, R01, R02, and R03 is presented in fig-
ure 4.15. The contributions of the deviation of each term towards the deviation
of displacement speed ∆Sd are shown in figure 4.16. Note that the deviation
of each component is normalised with the corresponding value of Srefd (refer to
equation 4.9). The normalised deviations are given by:
(∆nQ) =
Q−Qref
Srefd
(4.13)
where Q is the displacement speed component. A positive value of ∆nSd denotes
that the unsteady flame has a higher Sd relative to the reference value.
The normal and reaction components Sn + Sr are the major contributors to
the displacement speed as observed in figure 4.15b. Indeed the sum of Sn + Sr
contributes about 90% of the total Sd. The cross dissipation displacement speed
component Scd (refer to figure 4.15e) contributes about 9% to Sd and the last 1%
comes from the extra diffusion term Sdiff,z (refer to figure 4.15f).
From figure 4.15c, Sn is negative for all Φ0.90 while Sr (refer to figure 4.15d)
69
0.4 0.6 0.8 1
0
0.5
1
1.5
2
2.5
Φ0.90
S d
[m
/s
]
 
 
R00
R01
R02
R03
steady
(a) Sd
0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
Φ0.90
S n
+
S r
[m
/s
]
 
 
R00
R01
R02
R03
steady
(b) Sn + Sr
0.5 0.6 0.7 0.8 0.9 1
−6
−5
−4
−3
−2
−1
Φ0.90
S n
[m
/s
]
 
 
R00
R01
R02
R03
steady
(c) Sn
0.5 0.6 0.7 0.8 0.9 1
0
2
4
6
8
Φ0.90
S r
[m
/s
]
 
 
R00
R01
R02
R03
steady
(d) Sr
0.5 0.6 0.7 0.8 0.9 1
0.05
0.1
0.15
0.2
0.25
Φ0.90
S c
d
[m
/s
]
 
 
R00
R01
R02
R03
steady
(e) Scd
0.5 0.6 0.7 0.8 0.9 1
−6
−4
−2
0
2
x 10−3
Φ0.90
S d
iff
,z
[m
/s
]
 
 
R00
R01
R02
R02
steady
(f) Sdiff,z
Figure 4.15: Displacement speed and its components for cases R00, R01, R02,
and R03.
70
0.5 0.6 0.7 0.8 0.9 1
−0.05
0
0.05
Φ0.90
∆
n
S d
 
 
R00
R01
R02
R03
(a) ∆nSd
0.5 0.6 0.7 0.8 0.9 1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
Φ0.90
∆
n
(S
n
+
S r
)
 
 
R00
R01
R02
R03
(b) ∆n(Sn + Sr)
0.5 0.6 0.7 0.8 0.9 1
−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1
Φ0.90
∆
n
S n
 
 
R00
R01
R02
R03
(c) ∆nSn
0.5 0.6 0.7 0.8 0.9 1
−0.1
−0.05
0
0.05
0.1
0.15
0.2
0.25
Φ0.90
∆
n
S r
 
 
R00
R01
R02
R03
(d) ∆nSr
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
Φ0.90
∆
n
S c
d
 
 
R00
R01
R02
R03
(e) ∆nScd
0.5 0.6 0.7 0.8 0.9 1
−1
0
1
2
3
4
5
6
x 10−4
Φ0.90
∆
n
(S
di
ff
,z
)
 
 
R00
R01
R02
R02
(f) ∆n(Sdiff,z)
Figure 4.16: Deviations of the displacement speed and its components from the
steady state counterpart for cases R00, R01, R02, and R03
71
is positive for all Φ0.90.
The difference between the values of Sn and Sr becomes smaller for leaner
flames.
The cross dissipation component Scd, which is shown plotted against Φ0.90 in
figure 4.15e, has a positive value for all Φ0.90. Hence the Scd contribution to Sd
is always positive. This is in agreement with Richardson et. al. [93] where it is
reported that a back supported flame (flame that is propagating towards leaner
mixture) has a positive contribution from the cross dissipation term.
The deviation of displacement speed ∆nSd as plotted in figure 4.16 is shown
to have positive values in all cases. This indicates that the R case flames are
registering a higher displacement speed relative to their steady laminar flame
counterpart. This is consistent with the finding from case S00 in section 4.5.2.
Similarly, Pires Da Cruz et. al. [85], Kang and Kyritsis [70], and Richardson et.
al. [93] also reported similar findings where flames that propagates into a leaner
mixture is observed to have higher propagation speed.
In terms of the contribution towards the deviation of displacement speed ∆nSd
(figure 4.16), the contribution of ∆nScd is comparable to that of ∆n(Sn + Sr),
but the ∆nSdiff,z contribution has remained small (≈ 1% of Sd). Individually,
∆nSn and ∆nSr have values that are an order higher than ∆nSd. Since the ∆nScd
contribution to ∆nSd is always positive (as shown in figure 4.16e), the main reason
why ∆nSd drops below zero for cases R00 and R01 is due to the negative values of
∆n(Sn+Sr) which occur at Φ0.90 = 0.71 for case R00 and at Φ0.90 = 0.73 for case
R01. The negative ∆n(Sn +Sr) values stem from the large negative ∆nSn values
that were not compensated by positive values of ∆nSr. Moreover, the values of
∆nSr also drop below zero when approaching the end of the transition to lower
ΦG (Φ0.90 ≈ Φfinal0.90 ), compounding the negative ∆n(Sn + Sr) values.
4.7.1 The contributors to the deviation of normal dis-
placement speed component Sn
The normal displacement speed Sn is the linear combination of three contribu-
tions (c.f. equation 4.11): the gradient of the effective progress variable diffusion
72
coefficient,
Sn,00 =
∂Dc
∂x (4.14)
the logarithmic gradient of
Sn,01 = Dc
∂ ln ρ
∂x (4.15)
and the surface density
Sn,02 = Dc
∂ ln |∇c|
∂x (4.16)
The deviation of each term contributes to the net deviation of normal displace-
ment speed component ∆Sn. The deviation of each contributing term normalised
by the reference normal displacement speed component Srefn plotted in figure 4.17.
The deviation of the diffusion coefficient gradient term ∆Sn,00 (refer to fig-
ure 4.17a) is contributing negatively towards ∆Sn (refer to figure 4.17d). How-
ever, the deviation of the density gradient term ∆Sn,01 (figure 4.17b) and the de-
viation of the surface density function gradient term ∆Sn,02 (figure 4.17c) serve to
enhance ∆Sn. The term ∆Sn,02 is the dominant contributing term to ∆Sn. The
deviations of the diffusion coefficient term ∆Sn,00 and the density gradient term
∆Sn,01 are small and they have opposite signs. The sum of (∆Sn,00 + ∆Sn,01)
produces a small net negative contribution towards ∆Sn, supporting ∆Sn,02 as
the main contributor towards ∆Sn.
4.7.2 The contributors to the deviation of reaction dis-
placement speed component Sr
The reaction component of displacement speed Sr is influenced by four quantities:
density ρ, surface density function |∇c|, carbon and hydrogen element mixture
fraction zC + zH , and the reaction rate of methane ω˙CH4 . Since the comparison
is made at the local equivalence ratio on the c = 0.90 isosurface, the deviation
of the sum of carbon and hydrogen element mixture fraction is expected to be
zero. This is because there is only one value of z for a given value of Φ0.90.
The deviations of each of the other three quantities are plotted and compared in
figure 4.18.
Since Sr is the product of the inverse of density 1/ρ, the inverse of surface den-
73
0.5 0.6 0.7 0.8 0.9 1
−0.03
−0.025
−0.02
−0.015
−0.01
−0.005
0
0.005
Φ0.90
∆
S n
,0
0
/S
re
f
n
 
 
R00
R01
R02
R03
(a)
(
∆∂Dc∂x
)
/Srefn
0.5 0.6 0.7 0.8 0.9 1
−5
0
5
10
15
20
x 10−3
Φ0.90
∆
S n
,0
1
/S
re
f
n
 
 
R00
R01
R02
R03
(b)
(
∆Dc ∂ ln ρ∂x
)
/Srefn
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
0.1
Φ0.90
∆
S n
,0
2
/S
re
f
n
 
 
R00
R01
R02
R03
(c)
(
∆Dc ∂ ln |∇c|∂x
)
/Srefn
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
0.1
Φ0.90
∆
S n
/S
re
f
n
 
 
R00
R01
R02
R03
(d) (∆Sn) /Srefn
Figure 4.17: The contributions towards the deviation of normal component of the
displacement speed Sn.
sity function 1/|∇c|, and the progress variable reaction rate term ω˙c = −ω˙CH4/z,
the deviation of each term was normalised with the corresponding steady state
counterpart value. The contribution from the deviation of the inverse den-
sity term ∆(1/ρ) to ∆Sr is small compared to the deviation of the other two
terms. The magnitude of the deviation of ∆(1/|∇c|)/(1/|∇c|ref) is half that of
∆ω˙c/ω˙c,ref , and they are opposite in sign. The deviations of the two terms ω˙c
and |∇c| are the main contributors to ∆Sr.
74
0.5 0.6 0.7 0.8 0.9 1
−3
−2
−1
0
1
2
3
x 10−3
Φ0.90
∆
(1
/ρ
)/
(1
/ρ
)r
ef
 
 
R00
R01
R02
R03
(a) ∆(1/ρ)(1/ρ)ref
0.5 0.6 0.7 0.8 0.9 1
−0.06
−0.05
−0.04
−0.03
−0.02
−0.01
0
0.01
Φ0.90
∆
(1
/|
∇
c|)
/(
1/
|∇
c|)
re
f
 
 
R00
R01
R02
R03
(b) ∆(1/|∇c|)(1/|∇c|)ref
0.5 0.6 0.7 0.8 0.9 1
−0.05
0
0.05
0.1
0.15
Φ0.90
∆
(ω˙
C
H
4
/z
)/
(ω˙
C
H
4
/z
)re
f
 
 
R00
R01
R02
R03
(c) ∆(ω˙CH4/z)(ω˙CH4/z)ref
0.5 0.6 0.7 0.8 0.9 1
−0.04
−0.02
0
0.02
0.04
0.06
0.08
Φ0.90
∆
S r
/S
re
f
r
 
 
R00
R01
R02
R03
(d) ∆Sr
Srefr
Figure 4.18: The contribution of relevant terms towards the deviation of the
reaction rate component of the displacement speed Sr.
4.7.3 The contributors to the deviation of cross dissipa-
tion displacement speed component Scd
In the one-dimensional context, the cross dissipation displacement speed term
Scd is the product of diffusivity Dc and normalised gradient of mixture fraction
∂ ln z/∂x. The deviations of these two contributing terms, each normalised by its
corresponding steady state value, and the contributions towards the deviation of
Scd are plotted in figure 4.19.
The deviation of the diffusivity ∆Dc plotted in figure 4.19a is small compared
to the deviation of the normalised gradient of mixture fraction ∆∂ ln z/∂x, plotted
75
0.5 0.6 0.7 0.8 0.9 1
−4
−2
0
2
4
x 10−3
Φ0.90
∆
D
c/
D
re
f
c
 
 
R00
R01
R02
R03
(a) ∆Dc
Drefc
0.5 0.6 0.7 0.8 0.9 1
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Φ0.90
∆
(∂
x
ln
z)
/(
∂ x
ln
z)
re
f
 
 
R00
R01
R02
R03
(b) ∆
∂ ln z
∂x
∂ ln z
∂x
ref
0.5 0.6 0.7 0.8 0.9 1
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Φ0.90
∆
S c
d/
Sr
ef cd
 
 
R00
R01
R02
R03
(c) ∆Scd
Srefcd
Figure 4.19: The contribution of relevant terms towards the deviation of the cross
dissipation component of the displacement speed Scd.
in figure 4.19b. The small value of ∆Dc implies that the main contributing term
towards ∆Sn stems from ∆|∇c| (refer to section 4.7.1). Similar observations were
reported by Richardson et. al. [93].
The magnitude of the normalised value of ∆∂ ln z/∂x is comparable to that
of the normalised ∆Scd. This indicates that ∆Scd is directly linked to the de-
viation of the normalised gradient of element mixture fraction ∂ ln z/∂x (refer
to figures 4.19b and 4.19c). Note that Scd is unaffected by the change of sur-
face density function |∇c|. Also, Srefcd is not zero because there is an equivalence
ratio dip as observed in figure 4.2 that ensure non-zero gradient of z (refer to
equation 4.12).
76
4.7.4 Dominant terms that influence Sn, Sr, and Scd devi-
ation
It is shown in figure 4.17 that the the normalised gradient of surface density
function ∂ ln |∇c|/∂x is the main contributor towards the deviation of Sn. The
dynamics of |∇c| and ω˙c contribute towards the net deviation of Sr (figure 4.18).
Lastly, the deviation of the cross dissipation displacement speed component Scd
is directly affected by the deviations of normalised gradient of element mixture
fraction ∂ ln z/∂x (shown in figure 4.19). The behaviour of these four terms is
presented next.
Density ρ
The instantaneous density ρ is plotted against local equivalence ratio Φ0.90 in
figure 4.20a. The instantaneous density for an unsteady flame propagating into
leaner reactants (Cases R00 to R03) deviates slightly from the steady state coun-
terpart at a given instantaneous local equivalence ratio (Φ0.90). The deviation of
unsteady flame density on the c = 0.90 iso-surface from the steady state coun-
terpart is plotted against Φ0.90 in figure 4.20b. The density deviation is not
expected to be large as the instantaneous density is taken from the same iso-
surface at c = 0.90. However the existence of density deviations in the unsteady
flame points to the occurrence of a change in species mass fraction composition
at the iso-surface.
The density deviation is mainly negative except for case R00 and R01 at
lean Φ0.90. A negative deviation suggests that the mixture composition at the
c = 0.90 iso-surface consists of a larger proportion of low density molecules rather
than their heavier counterparts. Carbon element containing species are relatively
heavier in terms of molecular weight compared to non-carbon element containing
species like H radicals. Therefore, a negative density deviation ∆ρ is accompanied
by a negative deviation of zC . This is confirmed in the next section.
Surface density function |∇c| and its gradient ∂|∇c|∂x
The surface density function |∇c| is a measure of the surface area in a given
volume. In the one-dimensional context, the surface density function can be
77
0.5 0.6 0.7 0.8 0.9 1
0.19
0.2
0.21
0.22
0.23
0.24
0.25
Φ0.90
ρ
[k
g/
m
3 ]
 
 
R00
R01
R02
R03
steady
(a) ρ
0.5 0.6 0.7 0.8 0.9 1
−3
−2
−1
0
1
2
3
x 10−3
Φ0.90
∆
ρ/
ρr
ef
 
 
R00
R01
R02
R03
(b) (∆ρ)/ρref
0.5 0.6 0.7 0.8 0.9 1
800
1000
1200
1400
1600
1800
2000
2200
Φ0.90
|∇
c|
[1
/m
]
 
 
R00
R01
R02
R03
steady
(c) |∇c|
0.5 0.6 0.7 0.8 0.9 1
−0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
Φ0.90
∆
|∇
c|/
|∇
c|r
ef
 
 
R00
R01
R02
R03
(d) (∆|∇c|)/|∇c|ref
0.5 0.6 0.7 0.8 0.9 1
−3
−2.5
−2
−1.5
−1
−0.5
0
x 107
Φ0.90
∂|
∇
c|/
∂x
[1
/m
2 ]
 
 
R00
R01
R02
R03
steady
(e) ∂|∇c|∂x
0.5 0.6 0.7 0.8 0.9 1
−0.05
0
0.05
0.1
0.15
Φ0.90
∆
(∂
x
|∇
c|)
/(
∂ x
|∇
c|)
re
f
 
 
R00
R01
R02
R03
(f)
(
∆∂|∇c|∂x
)
/
(
∂|∇c|
∂x
)ref
Figure 4.20: Density ρ, surface density function |∇c|, and surface density function
spatial gradient ∂|∇c|/∂x profiles and their deviations for cases R00, R01, R02,
and R03.
78
perceived as an inverse of a length scale, and this length scale is generally taken
as the flame thickness. Although the maximum of the surface density gradient
does not reside on the c = 0.90 iso-surface, the changes of |∇c| at this iso-surface
reflect the change in the flame thickness. The instantaneous surface density
function |∇c|0.90 is plotted against Φ0.90 in figure 4.20c. The positive value of
∂|∇c|0.90/∂Φ0.90 reflects the fact that leaner flame has a larger flame thickness.
As the unsteady flame propagates into a leaner mixture, the deviation of
surface density function from its steady state counterpart is observed to increase
to a maximum at equivalence ratio of about Φ0.90 = Φfinal0.90 + 12∆Φ (refer to
figure 4.20d), then to return to zero as the flame approaches the final equivalence
ratio ΦfinalG . The increase in surface density function translates to a contraction
of flame thickness. All cases show a contraction of flame thickness during the
transition from stoichiometric to lean mixture.
The maximum of surface density function occurs between 0.64 < c < 0.69 for
the current range of equivalence ratio, therefore the spatial gradient of surface
density function is observed to be negative. A more negative value of ∂|∇c|/∂x
at c = 0.90 gives a higher value of max |∇c|, which translates to a narrower flame
thickness. This is consistent with the |∇c| findings presented earlier. A flame with
more negative instantaneous ∂|∇c|/∂x is a flame with a smaller flame thickness.
In the one-dimensional geometry, the deviation of the surface density spatial
gradient is identical in trend to the deviation of surface density function. However
the magnitude of the percentage deviation differs, with the spatial gradient of
surface density function having a higher deviation.
Carbon and hydrogen element mixture fractions, (zC, zH)
The carbon element mixture fraction (zC) and the hydrogen element mixture frac-
tion (zH) profiles are plotted against Φ0.90 in figure 4.21a and 4.21c respectively.
79
The deviations of these element mixture fractions are given as:
∆nzC =
zC − zrefC
zref =
∆zC
zref
∆nzH =
zH − zrefH
zref =
∆zH
zref
(4.17)
where ∆nzC and ∆nzH are the normalised deviations of carbon and hydrogen ele-
ment mixture fraction respectively. The superscript denotes reference values, e.g.
zref = zrefC + zrefH is the steady laminar flame mixture fraction that corresponds
to the unsteady flame instantaneous value of Φ0.90. Figures 4.21b and 4.21d show
the plots of ∆nzC and ∆nzH against Φ0.90 respectively.
Both zC and zH have positive values and positive gradient in the Φ0.90 space.
This is because Φ0.90 is a function of z = zC + zH where Φ0.90 ∝ z1−z . The
deviation of mixture fraction z for a specific local equivalence ratio on an iso-
surface will be zero because there is only one value for Φ0.90 for a given value of
z. However, this does not restrict the changes in individual carbon and hydrogen
elements. The carbon and hydrogen element mixture fraction for the unsteady
flame can individually deviate from the steady counterpart as long as the sum of
these deviations is zero for all cases.
From figures 4.21b and 4.21d, the deviations of hydrogen and carbon element
mixture fraction are small but not zero and their profiles are the mirror image of
one another. The noticeable change in zC and zH is consistent with the finding
of small density deviation presented previously (figure 4.20b), where the mixture
composition is found to differ in the unsteady cases from the steady counter-
part. The negative deviation in zC and positive deviation of zH point to lower
percentage of carbon containing species at Φ0.90 compared to the reference value.
Spatial gradients of carbon and hydrogen element mixture fractions
The spatial gradient of carbon element mixture fraction ∂zC/∂x and the spatial
gradient of hydrogen element mixture fraction ∂zH/∂x are plotted against Φ0.90
in figures 4.22a and 4.22c respectively. The ∂zC/∂x and ∂zH/∂x deviations nor-
malised by the steady state total mixture fraction ((∂z/∂x)ref ) are plotted in
80
0.5 0.6 0.7 0.8 0.9 1
0.02
0.025
0.03
0.035
0.04
Φ0.90
z C
 
 
R00
R01
R02
R03
steady
(a) zC
0.5 0.6 0.7 0.8 0.9 1
−3
−2
−1
0
1
2
3
4
x 10−3
Φ0.90
∆
z C
/z
re
f
 
 
R00
R01
R02
R03
(b) ∆nzC
0.5 0.6 0.7 0.8 0.9 1
0.008
0.009
0.01
0.011
0.012
0.013
0.014
Φ0.90
z H
 
 
R00
R01
R02
R03
steady
(c) zH
0.5 0.6 0.7 0.8 0.9 1
−4
−3
−2
−1
0
1
2
3
x 10−3
Φ0.90
∆
z H
/z
re
f
 
 
R00
R01
R02
R03
(d) ∆nzH
Figure 4.21: The profiles and deviations of carbon and hydrogen element mixture
fractions (zC , zH).
81
figure 4.22b and 4.22d.
0.5 0.6 0.7 0.8 0.9 1
4
6
8
10
12
14
16
18
Φ0.90
∂z
C
/∂
x
[1
/m
]
 
 
R00
R01
R02
R03
steady
(a) ∂zC∂x
0.5 0.6 0.7 0.8 0.9 1
−0.1
0
0.1
0.2
0.3
0.4
Φ0.90
∆
(∇
z C
)/
(∇
z)
re
f
 
 
R00
R01
R02
R03
(b) ∆n ∂zC∂x
0.5 0.6 0.7 0.8 0.9 1
−6
−5
−4
−3
−2
−1
Φ0.90
∂z
H
/∂
x
[1
/m
]
 
 
R00
R01
R02
R03
steady
(c) ∂zH∂x
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
0.1
0.12
Φ0.90
∆
(∇
z H
)/
(∇
z)
re
f
 
 
R00
R01
R02
R03
(d) ∆n ∂zH∂x
Figure 4.22: The spatial gradient and the deviation of the gradient of carbon
and hydrogen element mixture fractions (∂zC/∂x, ∂zH/∂x) plotted against local
equivalence ratio Φ0.90
The c = 0.90 iso-surface is located between the position where the minimum of
local equivalence ratio min(ΦL) is situated (c.f. figure 4.7) and the products. At
this iso-surface, the value of ∂zC/∂x is always positive and the gradient is larger
at higher ΦL, while the value of ∂zH/∂x is always negative and the gradient
is more negative at higher ΦL. This can be observed from figures 4.3 and 4.4
respectively.
82
The normalised element mixture fraction gradient is given by:
∆n
∂zβ
∂x =
∂zβ
∂x −
∂zβ
∂x
ref
∂z
∂x
ref =
∆∂zβ∂x
∂z
∂x
ref (4.18)
where β can be C or H.
Both ∆n∂zC/∂x and ∆n∂zH/∂x are positive and they both have positive
contributions towards the enhancement of ∂z/∂x. The carbon element mixture
fraction gradient ∂zC/∂x for the unsteady cases is observed to have been en-
hanced during the unsteady transition from the initial to final equivalence ratio.
However, the hydrogen element mixture fraction gradient ∂zH/∂x has positive de-
viation when it has negative values, resulting in a less steep spatial gradient of zH
during the unsteady transition. The quantity ∆∂z/∂x = (∆∂zC/∂x + ∆∂zC/∂x)
governs the deviation of the cross dissipation term of the displacement speed
component Scd. It is observed that the value of ∆n∂zC/∂x is about four times
larger than the value of ∆n∂zH/∂x in the enhancement of ∂z/∂x.
Methane reaction rate (ω˙c)
The effective progress variable reaction rate (ω˙c) is a function of methane reaction
rate and mixture fraction. The deviation of ω˙c is one of the major contributors
towards ∆Sr. It has been established that the mixture fraction deviation ∆z dur-
ing the unsteady transition is zero, hence the deviation of ω˙c stems solely from
the deviation of ω˙CH4. The methane reaction rate profile for the unsteady and
steady flame is plotted in figure 4.23a. The interpolation of the laminar steady
flame value of ω˙CH4 is superimposed for reference. The unsteady transition devi-
ation of methane reaction rate was normalised with the steady state counterpart
and plotted against Φ0.90 in figure 4.23b.
The deviation of the normalised methane reaction rate is given as:
∆nω˙CH4 =
ω˙CH4 − ω˙refCH4
ω˙refCH4
= ∆ω˙CH4
ω˙refCH4
(4.19)
The methane reaction rate magnitude |ω˙CH4| is larger for flames with higher
Φ0.90 as expected and it was enhanced during the unsteady transition from the
83
0.5 0.6 0.7 0.8 0.9 1
−200
−150
−100
−50
0
ω˙ C
H
4
kg
/(
m
3 s
)
Φ0.90
 
 
R00
R01
R02
R03
steady
(a) ω˙CH4
0.5 0.6 0.7 0.8 0.9 1
−0.05
0
0.05
0.1
0.15
Φ0.90
(∆
ω˙ C
H
4
)/
ω˙r
ef C
H
4
 
 
R00
R01
R02
R03
(b) (∆ω˙CH4)/ω˙refCH4
Figure 4.23: Methane reaction rate ω˙CH4
initial to final equivalence ratio. This results in the unsteady flame instantaneous
ω˙CH4 on c = 0.90 having a larger magnitude compared to the steady state coun-
terpart. Also, a negative ∆nω˙CH4 coincides with the change of sign in ∆nzC and
∆nzH . This implies that the change in sign of ∆nω˙CH4 is due to the change in
instantaneous mixture composition.
4.7.5 Total heat release rate Q˙hr
The enhancement of reaction rate also increases the heat release of the flame. The
total heat release rate for the unsteady and steady flame are plotted against Φ0.90
is figure 4.24a. The heat release is negative and the heat release rate magnitude
is larger at higher Φ0.90 as expected. The deviations of the heat release rate nor-
malised by its reference value ∆Q˙hr/Q˙refhr is plotted against Φ0.90 in figure 4.24b.
The heat release rate normalised deviation has the same profile as the nor-
malised reaction rate deviations. From figure 4.24b, heat release enhancement
of about 14% is observed for case R00. The enhanced heat release for flame
that propagates into leaner mixture (back supported flame) is in agreement with
Richardson et. al. [93]
84
0.5 0.6 0.7 0.8 0.9 1
−5
−4
−3
−2
−1
0
x 109
Φ0.90
Q˙
hr
[W
/m
3 ]
 
 
R00
R01
R02
R03
steady
(a) Q˙hr
0.5 0.6 0.7 0.8 0.9 1
−0.05
0
0.05
0.1
0.15
Φ0.90
(∆
Q˙
hr
)/
Q˙
re
f
hr
 
 
R00
R01
R02
R03
(b) (∆Q˙hr)/Q˙refhr
Figure 4.24: Heat release rate Q˙hr and the deviation from its reference values
∆Q˙hr/Q˙refhr .
4.7.6 Heat flux on the c = 0.90 iso-surface q˙
The heat flux across the c = 0.90 iso-surface is plotted against Φ0.90 in figure 4.25a.
The heat flux is negative across all values of Φ0.90, showing that the heat flux
vector is pointing towards the reactants. Deviations of heat flux are observed
during the unsteady transition. The deviations of heat flux for the unsteady cases
were normalised by the steady state counterpart. These normalised deviations
∆q˙/q˙ref are plotted in figure 4.25b.
The heat fluxes under unsteady conditions were generally higher in magnitude
compared to the steady state counterpart. The heat flux is given as:
q˙ = q˙T + q˙s = −λ
∂T
∂x +
Ns∑
α=1
ρVα,xYαhα (4.20)
The heat flux due to heat conduction (Fourier’s Law) q˙T = −λ∂T/∂x and the flux
of the species enthalpy q˙s =
∑Ns
α=1 ρVα,xYαhα are plotted in figure 4.25c and 4.25e
respectively.
The heat conduction term q˙T is negative for all Φ0.90. Similar to the heat flux
vector, it is pointing towards the reactants. The deviation of the normalised heat
85
0.5 0.6 0.7 0.8 0.9 1
−2
−1.5
−1
−0.5
0
x 106
Φ0.90
q˙
 
 
R00
R01
R02
R03
steady
(a) q˙
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
0.1
Φ0.90
(∆
q˙)
/q˙
re
f
 
 
R00
R01
R02
R03
(b) (∆q˙)/(q˙)ref
0.5 0.6 0.7 0.8 0.9 1
−2
−1.5
−1
−0.5
0
x 106
Φ0.90
q˙ T
 
 
R00
R01
R02
R03
steady
(c) q˙T = −λ∂T∂x
0.5 0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
0.08
0.1
0.12
Φ0.90
(∆
q˙ T
)/
q˙r
ef
 
 
R00
R01
R02
R03
(d) (∆q˙T )/(q˙)ref
0.5 0.6 0.7 0.8 0.9 1
1
1.5
2
2.5
3
3.5
x 105
Φ0.90
q˙ s
 
 
R00
R01
R02
R03
steady
(e) q˙s =
∑N
α=1 ρVαYαhα
0.5 0.6 0.7 0.8 0.9 1
−0.025
−0.02
−0.015
−0.01
−0.005
0
0.005
Φ0.90
(∆
q˙ s
)/
q˙r
ef
 
 
R00
R01
R02
R03
(f) (∆q˙s)/(q˙)ref
Figure 4.25: The profiles and deviations of total heat flux q˙, heat conduction
term q˙T , and species enthalpy term q˙s for cases R00, R01, R02, and R03.
86
conduction term is given as:
∆nq˙T =
q˙T − q˙refT
q˙ref (4.21)
and is plotted in figure 4.25d against Φ0.90.
The species enthalpy flux term q˙s is positive for all Φ0.90 and is transport-
ing heat in the opposite direction compared to the heat conduction term. The
deviation of the normalised species enthalpy flux is given as:
∆nq˙s =
q˙s − q˙refs
q˙ref (4.22)
and is plotted in figure 4.25f against Φ0.90.
During the unsteady transition, the instantaneous value of q˙T becomes more
negative, enhancing the overall heat flux vector. However q˙s becomes more posi-
tive during the unsteady process, counteracting the enhancement brought by q˙T .
The magnitude of ∆nq˙T is seven times higher than ∆nq˙s, and hence q˙ is main-
tained strongly negative. The quantity ∆q˙T is dominant in the deviation of total
heat flux. Note that the temperature deviation is in the order of 0.1% which is
relatively small.
4.7.7 Species that influence the deviation of carbon and
hydrogen elements (zC, zH) and their gradients.
The species contribution towards carbon and hydrogen element mixture fraction
provides insights into which species are responsible for the deviation of zC and zH .
The contribution of individual species Yα towards the element mixture fraction
zβ is given by:
zβ,α = nβ,α
Wβ
Wα
Yα (4.23)
where nβ,α is the number of atoms of element β in species α, and Wβ and Wα
are the molecular weight for element β and species α respectively. For example,
the contribution of methane mass fraction to hydrogen element mixture fraction
87
is given by zh,CH4 = 4× YCH4/16. The deviation of zβ,α is therefore given by:
∆zβ,α = zβ,α − zrefβ,α (4.24)
where zrefβ,α is obtained from interpolation from the one-dimensional flame library.
The deviation of species contribution is normalised with the local mixture fraction
as ∆zβ,α/z. The carbon element and hydrogen element species contribution time
history for case R01 is plotted against instantaneous Φ0.90 in figure 4.26a and 4.26b
respectively.
0.7 0.75 0.8 0.85 0.9 0.95
−8
−6
−4
−2
0
2
4
6
x 10−3
c
∆
z c
,α
/∆
z
 
 
CH4
CO2
CO
CH2O
HCO
CH2OH
CH3OH
CH3
CH3O
(a) ∆zc,α/zref
0.7 0.75 0.8 0.85 0.9 0.95
−5
0
5
10
x 10−4
c
∆
z h
,α
/∆
z
 
 
CH4
H2O
H2
OH
H
HO2
H2O2
CH2O
HCO
CH2OH
CH3OH
CH3
CH3O
(b) ∆zh,α/zref
Figure 4.26: The deviation of species contribution towards carbon and hydrogen
elements for case R01.
88
From figure 4.26a, it is observed that the two dominant species that con-
tribute toward ∆zC are carbon dioxide CO2 and carbon monoxide CO. The
carbon monoxide contributes positively towards ∆zC while the carbon dioxide
contributes negatively towards ∆zC .
From figure 4.26b, water molecule H2O, hydrogen molecule H2, hydroxyl rad-
ical OH, and hydrogen atom H are found to influence ∆zH . The first two species
each have an order of magnitude larger contribution compared to the latter two
species.
The six species that influence ∆zC and ∆zH are key species for the water-
gas shift reaction and the oxygen-hydrogen reactions [84]. These reactions are
summarised as:
CO + H2O ⇋ CO2 + H2 (4.25)
O2 + 3H2 ⇋ 2H2O + 2H (4.26)
The deviation of the spatial gradient of mixture fraction ∆n∂z/∂x dominantly
influences the deviation of the cross dissipation term ∆nScd (refer to figure 4.19).
The species contribution towards the deviation of carbon and hydrogen element
spatial gradient can be obtained using equation 4.23 by differentiating with re-
spect to x
∂zβ,α
∂x = nβ,α
Wβ
Wα
∂Yα
∂x (4.27)
Similarly, the deviations of ∂zβ,α∂x from the reference value is given as:
∆∂zβ,α∂x =
∂zβ,α
∂x −
∂zβ,α
∂x
ref
(4.28)
The species contribution towards the spatial gradient of element mixture frac-
tion time history is plotted against local instantaneous Φ0.90 for case R01 in fig-
ure 4.27a for carbon element and figure 4.27b for hydrogen element.
The species that are responsible for the deviations of ∂zc,α/∂x and ∂zh,α/∂x
are the same species that are responsible for the deviations of zc,α and zh,α. Car-
bon dioxide CO2 and carbon monoxide CO both contribute positively towards
the deviations of ∂zc,α/∂x (refer to figure 4.27a). Water molecule H2O, hydroxyl
89
0.7 0.75 0.8 0.85 0.9 0.95
−0.15
−0.1
−0.05
0
0.05
0.1
0.15
0.2
c
∆
(∇
x
z c
,α
)/
∇
x
z
 
 
CH4
CO2
CO
CH2O
HCO
CH2OH
CH3OH
CH3
CH3O
(a) ∆∂zc,α∂x / ∂z∂x
ref
0.7 0.75 0.8 0.85 0.9 0.95
−0.05
0
0.05
0.1
c
(∆
∇
x
z h
,α
)/
∇
x
z
 
 
CH4
H2O
H2
OH
H
HO2
H2O2
CH2O
HCO
CH2OH
CH3OH
CH3
CH3O
(b) ∆∂zh,α∂x / ∂z∂x
ref
Figure 4.27: The deviation of species contribution towards carbon and hydrogen
elements spatial gradient for case R01.
radical OH, and hydrogen atom H contribute positively towards the deviation
of ∂zh,α/∂x while hydrogen molecule H2 has a negative contribution (refer to
figure 4.27b). In addition, the methane species is found to have a negative con-
tribution towards the deviations of both carbon and hydrogen element spatial
gradient as observed in figure 4.27.
In section 4.19, the deviation of ∂zc/∂x is found to be the main contributor
90
towards ∆∂z/∂x which in turn dominantly influences the deviation of the cross
dissipation displacement speed term ∆Scd. Hence carbon dioxide CO2 and carbon
monoxide CO are identified as the major species that influence ∆Scd.
4.7.8 Relevant species mass fractions deviations and dif-
fusive mass flux deviations
The species spatial gradients, CH4, CO2, H2O, H2, OH, H, and CO all con-
tribute towards ∂z/∂x. Since ∆∂z/∂x is the main contributor towards ∆Scd,
these species will affect the value of ∆Sd. The maximum deviations of species
mass fraction and the deviations of diffusive mass fluxes for these key species in
case R01 are tabulated in table 4.7.
Table 4.7: Percentage enhancement of key species in mass fractions and diffusive
mass flux across the c = 0.90 iso-surface for case R01.
species α max
(
∆Yα
Y refα
)
× 100% max
(
∆m˙d,α
m˙refα
)
× 100%
CH4 ∼ O(-6) 2.789
CO2 -3.902 -22.34
H2O 0.375 -10.76
H2 6.240 7.137
OH 6.437 9.016
H 9.472 11.18
CO 2.254 -4.287
The deviation of the diffusive mass flux of these species ∆m˙d,α is normalised
by its corresponding total reference mass flux m˙α. The definition of m˙α is the
sum of convective mass flux m˙c,α and diffusive mass flux m˙d,α. It is given as:
m˙α = m˙c,α + m˙d,α
= ρSdYα + ρDα
∂Yα
∂x
(4.29)
The methane mass fraction at c = 0.90 is shown to have negligible changes
91
because the definition of progress variable c is based on YCH4. However, the diffu-
sive flux for methane is shown to have a small increase of 2% across the c = 0.90
isosurface during the transition of the unsteady flame from unity equivalence ra-
tio to ΦG = 0.70 in case R01. Carbon dioxide has a reduction in mass fraction
on the c = 0.90 isosurface. This is accompanied by an increase of carbon monox-
ide mass fraction. However, both carbon dioxide and carbon monoxide have a
reduction in flux across the c = 0.90 isosurface. The water molecule is observed
to have negligible change in mass fraction on the c = 0.90 iso-surface but its flux
is reduced by 10%. Both CO2 and H2O have a large reduction of flux across
the c = 0.90 isosurface with changes of about -20% and -10% respectively. The
hydrogen molecule, hydroxyl radical, and hydrogen atom are observed to have
increases in both mass fraction and diffusive flux across the c = 0.90 iso-surface.
The mass fraction of these key radical species for the unsteady flame is clearly
different from that of the steady state counterpart. These key species are also the
species involved in the water-gas shift reaction and the system of oxygen-hydrogen
reactions as reported by Peters and Williams [84]. The Arrhenius temperature
dependence of the equilibrium constants of these two reactions [84] modifies the
flux of highly reactive species such as OH through the reaction zone [93]. From
table 4.7, the H2, OH, and H radical species have their mass fraction enhanced
by 6%, 6% and 9% respectively. The diffusive fluxes for these species were also
enhanced by 7%, 9% and 11%. This finding is supports the observations of
Richardson et. al. [93].
4.8 The displacement speed deviation ∆Sd anal-
ysis for cases R01, R01b, R01c, and R01d
The behaviour of displacement speed Sd across the equivalence ratio space Φ0.90
for cases R00, R01, R02, and R03 is presented in figure 4.28. It is observed that
the trend for all R01 cases is similar, however, there is a small deviation observed
in case R01c. The differences in Sd for the R01 cases can be illustrated through
deviation of displacement speed ∆Sd against instantaneous local equivalence ratio
Φ0.90 as plotted in figure 4.29a.
92
0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
0.5
1
1.5
2
2.5
Φ0.90
S d
[m
/s
]
 
 
R01
R01b
R01c
R01d
steady
Figure 4.28: The displacement speed Sd against local equivalence ratio Φ0.90 for
cases R01, R01b, R01c, and R01d.
From figure 4.29a, it is observed that case R01c has a smaller maximum ∆nSd
compared to cases R01, R01b, and R01d. In fact the trajectories of ∆nSd for
cases R01, R01b and R01d are similar and the maximum value of ∆nSd is the
same. This highlights the fact that case R01c, which has max (∂Φ0.90/∂t) <
84.50s−1, has shown that ∆nSd is also smaller than the maximum attainable
value of ∆nSd ≈ 0.06 in cases R01, R01b, and R01d.
The deviation of the displacement speed components is plotted against the lo-
cal instantaneous equivalence ratio Φ0.90 in figure 4.29. It is shown that max (∆nScd)
(refer to figure 4.29e) for case R01c shows no difference compared to cases R01,
R01b, and R01d. The difference of max (∆nSd) between case R01c and the
rest can be attributed to ∆n(Sn + Sr) (refer to figure 4.29b). The maximum
of ∆n(Sn + Sr) that occurs at Φ0.90 ≈ 0.85 is shown to be about 15% smaller for
case R01c compared to the rest.
93
0.6 0.7 0.8 0.9 1
−0.02
0
0.02
0.04
0.06
Φ0.90
∆
S d
/S
re
f
d
(a) ∆nSd
0.6 0.7 0.8 0.9 1
−0.02
−0.01
0
0.01
0.02
0.03
Φ0.90
∆
(S
n
+
S r
)/
Sr
ef d
(b) ∆n(Sn + Sr)
0.6 0.7 0.8 0.9 1
−0.1
−0.05
0
0.05
Φ0.90
∆
S n
/S
re
f
d
(c) ∆nSn
0.6 0.7 0.8 0.9 1
−0.05
0
0.05
0.1
0.15
0.2
Φ0.90
∆
S r
/S
re
f
d
(d) ∆nSr
0.6 0.7 0.8 0.9 1
−0.01
0
0.01
0.02
0.03
0.04
Φ0.90
∆
S c
d/
Sr
ef d
(e) ∆nScd
0.6 0.7 0.8 0.9 1
−1
−0.5
0
0.5
1
1.5
2
2.5
x 10−4
Φ0.90
∆
(S
di
ff
,z
)/
Sr
ef d
 
 
R01
R01b
R01c
R01d
(f) ∆nSdiff,z
Figure 4.29: Displacement speed and its components for cases R00, R01, R02,
and R03
94
The drop in maximum of ∆n(Sn + Sr) for case R01c is attributed to the
reaction component Sr (refer to figure 4.29d). The normal diffusion component
maximum deviation for case R01c (refer to figure 4.29c), which occurs at about
Φ0.90 ≈ 0.77, is similar to the other R01 cases. However, the maximum deviation
for the case R01c reaction component max (∆nSr), which occurs at Φ0.90 ≈ 0.77,
is about 12% smaller compared to the other R01 cases.
95
4.9 Summary
One-dimensional simulations of unsteady and steady laminar premixed flames
have been performed. A library of steady laminar flame simulations with different
global equivalence ratio (0.5 ≥ ΦG ≥ 1.0) was obtained with the intention of
benchmarking the current reduced reaction mechanism as well as to create a
database of one-dimensional steady flames for unsteady flame comparison. The
unsteady laminar flame simulations were performed to investigate the effect of
mixture gradient that is parallel to the direction of flame propagation.
There is a local equivalence ratio dip within the flame for all the steady lam-
inar flames in the database. This equivalence ratio dip shows the presence of
preferential diffusion processes within the flame. The flame properties at the
c = 0.90 isosurface are used to compare the steady and unsteady laminar flames.
Using a well defined point within the flame for comparison alleviates the need to
interpolate across this equivalence ratio dip.
Two sets of unsteady laminar flame simulations were performed, one with
sinusoidal time varying inlet equivalence ratio (S cases), and the other with a
complementary error function time varying inlet equivalence ratio (R cases).
In the S cases, hysteresis of Sd is observed, where Sd becomes larger for a
flame that is propagating into a leaner mixture and becomes smaller for a flame
that is propagating into a richer mixture. The time averaged Sd for case S00
is 1% smaller than the value that corresponds to the time-averaged equivalence
ratio.
The main contributor to Sd is Sn + Sr, which represents 90% of the total.
The quantity Scd contributes about 9% while the rest is contributed by Sdiff,z.
However, the main contributors to the deviation ∆Sd are ∆(Sn + Sr) and ∆Scd.
The contribution of the extra diffusive term ∆Sdiff,z remains negligible.
In the R cases, the change of equivalence ratio magnitude ∆Φ influences the
trajectory of ∆Sd. The increase in ∆Φ produces a higher max(∆Sd). There is
also a maximum temporal gradient of Φ0.90 for a given change of equivalence ratio
magnitude ∆Φ. Hence, for a given ∆Φ, there is a maximum attainable ∆Sd.
The deviation of the normal component ∆Sn is controlled by ∆∂|∇c|/∂x,
the deviation of the reaction component ∆Sr is controlled by ∆ω˙CH4 , and the
96
deviation of the cross dissipation component ∆Scd is controlled by ∆∂z/∂x.
The surface density function |∇c| is the inverse of flame thickness, and be-
comes larger for flames propagating into a leaner mixture and smaller for flames
propagating into a richer mixture. This results in a thinning and thickening of
the flame during the unsteady transition. The effective reaction rate ω˙c and the
mixture gradient ∂z/∂x also increase for a rich to lean propagating flame and
vice versa.
The key species in the water-gas shift reactions and the system of oxygen-
hydrogen reactions are CO, CO2, H2, OH, and H. These radical species were
found to influence the trends of ∆∂z/∂x and ∆ρ. The mass diffusive fluxes of
these species also changed in the order of 10% with respect to the steady state
counterpart.
97
Chapter 5
Two-dimensional laminar flame
simulations
One-dimensional flame simulations described in the previous chapter have pro-
vided the benchmarks for steady and unsteady laminar premixed flames. In the
presence of temporally evolving mixture strength, the effects of normal stratifi-
cation on laminar premixed flames were discussed.
Extending from the one-dimensional flame cases, two-dimensional laminar
flame simulations were performed to incorporate the effects of flame geometry on
premixed flames with nonhomogeneous reactants. A gradient of mixture strength
perpendicular to the direction of flame propagation was imposed in the simula-
tions to study the effects of transverse mixture gradients on the flame.
This chapter will describe the configurations used for the two-dimensional
laminar flame simulations and the associated boundary conditions. The evolution
of the two-dimensional flame will be presented followed by the discussions of the
observed behaviour.
98
5.1 Two-dimensional inhomogeneous laminar flame
simulation configurations
The objective of the two-dimensional simulations is to study the evolution of the
initially planar, freely propagating premixed flame under the influence of trans-
versely stratified reactants. The computational domain consists of a rectangular
box measuring Lx =10 mm long by Ly =5 mm wide, discretised using a uniform
Cartesian mesh of 1000 by 500 points. The resulting square cell size is 10 µm by
10 µm.
The inhomogeneous reactants are introduced to the flame by imposing a
prescribed field of inhomogeneous reactants at the inlet boundary. The inflow
boundary conditions were set using the Navier-Stokes Characteristic Boundary
Condition (NSCBC) formulation. A constant mixture density was specified along
with a constant velocity component normal to the inlet boundary plane, while
the transverse velocity components were set to zero. The inlet mass fractions
were set to conform to a specified profile of reactant mixture strength which was
quantified using the equivalence ratio as defined in chapter 3.
The prescribed inlet equivalence ratio profile is shown in figure 5.1. The profile
consists of a peak region in the centre of the inlet corresponding to the maximum
mixture strength, regions of minimum mixture strength towards both outer edges
of the inlet, and regions of transition from the maximum to the minimum having
a prescribed gradient. The magnitude of the gradient is controlled by a gradient
coefficient α.
The equations used for prescribing the inlet profile of equivalence ratio are:
Φ(0 ≤ y < Ly/3) = M −A
Φ(Ly/3 ≤ y ≤ Ly/2) = β [ sin (2πθ/L) + α tanh (2παθ/L)] + M
Φ(Ly/2 < y ≤ Ly) = Φ(Ly − y)
(5.1)
99
0 0.2 0.4 0.6 0.8 1
0.5
0.6
0.7
0.8
0.9
1
non-dimensionalised cross-stream direction, yLy
In
le
tl
oc
al
eq
vu
ia
le
nc
e
ra
tio
Φ L
,i
n
 
 α = 0
α = 3
α = 5
α = 7
Figure 5.1: Equivalence ratio Φ specified at the inlet boundary.
where θ and β are given by:
θ = 3π/Ly(y − Ly/3)
β = 2A/(sin (π/2) + α tanh (απ/2)− sin (π/2)− α tanh (απ/2))
(5.2)
and where A and M are the specified amplitude and median of the inlet equiva-
lence ratio profile:
A = (Φmax − Φmin)/2
M = Φmed
(5.3)
The outflow boundary condition was set using the Navier-Stokes Characteris-
tic Boundary Condition (NSCBC) formulation, as discussed in chapter 3 to pro-
vide a non-reflecting condition for acoustic disturbances. The transverse bound-
aries were set to be periodic.
100
5.1.1 Test matrix of two-dimensional simulations
Six two-dimensional simulations were performed with the aim of assessing the
effects of the amplitude of the reactant inhomogeneity together with different
transverse mixture gradients. The relevant values used to prescribe the inlet
boundary conditions for these cases are given in Table 5.1.
Case Mean Median Max Min Φmax−Φmin
2
Gradient
Φ Φmed Φmax Φmin coefficient, α
A00 0.63 0.75 1.1 0.40 0.35 3
A01 0.63 0.75 1.1 0.40 0.35 5
A02 0.63 0.75 1.1 0.40 0.35 7
B00 0.67 0.75 1.0 0.50 0.25 3
B01 0.67 0.75 1.0 0.50 0.25 5
B02 0.67 0.75 1.0 0.50 0.25 7
Table 5.1: The inlet boundary conditions of the two-dimensional simulations.
There are two groups of simulations performed here; group A, where the max-
imum Φ is above unity and the minimum Φ is below the lean flammability limit,
and group B, where the maximum Φ is unity and the minimum Φ is approximately
equal to the flammability limit.
In each group, there are three simulations accounting for different gradient
coefficients. The higher gradient coefficient value will yield a steeper transverse
gradient between the maximum and minimum equivalence ratio. This makes the
inhomogeneous equivalence ratio profile more top-hat like.
The length scale of the distribution of inhomogeneous reactants is constant for
all cases, and is equal to one-third of the domain length, LΦ = Ly/3. This length
scale is similar to the length scale defined by Jimenez et. al. [67, 68]. The current
length scale of inhomogeneity is about three times larger than the flame length
scale associated with the mean equivalence ratio. This will ensure that the flame
will deform as it propagates towards the oncoming inhomogeneous reactants.
Note that diffusive effects will render the range of Φ experienced by the flame,
101
which is some distance downstream, to be lower than at the inlet.
5.1.2 Initial solution for the simulation
The initial conditions within the entire computational domain were set using
a precomputed planar laminar flame profile with global equivalence ratio of
ΦG = 0.80 for all cases tabulated in table 5.1. Additionally, the prescribed
equivalence ratio profile at the inlet is extended into the domain. The initial
equivalence ratio profile for case B00 is plotted in figure 5.2. The extension of
inhomogeneous reactants into the domain allows the reduction of required compu-
tational resources, which would otherwise be spent on calculating the convection
of inhomogeneous reactants towards the planar flame. The interface between the
initial inhomogeneous reactant distribution and the initial planar flame solution
is smoothed using an error function.
The prescribed inlet equivalence ratio profile can be seen in figure 5.2a, where
the red region denotes high equivalence ratio (Φ = 1.00) while the blue region
denotes low equivalence ratio (Φ = 0.50). The yellow region denotes the initial
flame equivalence ratio (Φ = 0.80). The extension of inhomogeneous reactants
into the domain is done with care so that it will not intrude the initial planar
flame solution. Note that the dip in the local equivalence ratio profile due to
preferential diffusion as discussed in the previous chapter 4 is clearly shown in
figure 5.2.
The initialised laminar flame has the following properties:
Φ0.90,init = ΦL = 7.575× 10−1
SL,init = SL = 1.954× 10−1 m s−1
δL,init = δL = 5.345× 10−4 m
tL,init = tL = 2.735× 10−3 s
(5.4)
These values will be used as the normalising factors in subsequent discussions.
102
(a) Case B00 initial equivalence ratio profile
(b) Plan view for case B00
Figure 5.2: The initial equivalence ratio profile for case B00. Φ = 0.50 in blue to
Φ = 1.00 in red.
103
5.2 The methane species mass fraction trans-
port budget
A closed budget provides strong evidence that the simulation and the post-
processing procedure is producing data with the required level of accuracy. The
budget of the methane mass fraction transport equation is plotted in figure 5.3
for case A00 at elapsed time t/tL = 2.2. The transport equation for methane
mass fraction is stated as equation 5.5 and the various terms are labelled:
unsteady
︷ ︸︸ ︷
∂ρYCH4
∂t +
convection
︷ ︸︸ ︷
∂ρukYCH4
∂xk
=
diffusion
︷ ︸︸ ︷
∂
∂xk
[
ρD∂YCH4∂xk
]
+
reaction
︷ ︸︸ ︷
ω˙CH4 (5.5)
From figure 5.3, it may be observed for the unsteady term (figure 5.3a), two
peaks were found near to the inlet boundary at a position close to that of the
maximum gradient of the inlet mixture strength, while a trough was found close
to the position of maximum mixture strength. It is interesting to note that the
convection term (figure 5.3b) exhibits a similar pattern to the unsteady term but
with the opposite sign.
The diffusion term (figure 5.3c) has its maximum value at the position where
the reaction term (figure 5.3d) also has its greatest magnitude. Note that the
reaction rate is negative since the methane fuel is being consumed within the
flame. The position of the maximum reaction rate magnitude can be used to
denote the location of the flame front. It is interesting to observe that the flame
front (figure 5.3d) is curved, with the middle of the flame front positioned nearer
to the inlet. This provides a strong indication that there are differences in the
propagation speed along the flame front.
The remaining plot (figure 5.3e) shows the residual error ǫ which results
from closing the budget of methane mass fraction numerically according to equa-
104
(a) unsteady term (b) convection term
(c) diffusion term (d) reaction term
(e) residual term, ǫ
Figure 5.3: Methane species mass fraction
105
tion 5.6.
ǫ = ∂ρYCH4∂t +
∂ρukYCH4
∂xk
− ∂∂xk
[
ρD∂YCH4∂xk
]
− ω˙CH4
(5.6)
The highest magnitude of the residual ǫ is found near the inlet boundary
where the inlet boundary condition on the species mass fraction distribution was
imposed. Moreover, the spatial differencing scheme near the boundary is of lower
accuracy than that in the domain. Note the scale on the z -axis, which indicates
that the numerical error present in the calculations amounts to about four orders
of magnitude less than that of any single term in the equation. This shows
clearly that indeed the simulation is accurate and that the budget of methane
mass fraction is closed.
5.2.1 Equivalence ratio and progress variable profile for
Case A00
The budget for Case A00 at elapsed time t/tL = 2.2 is closed and its progress
variable c and local equivalence ratio Φ at that instant is shown in figure 5.4.
Figure 5.4a is the progress variable plot for case A00 at t/tL = 2.2, where the
z -axis is the progress variable value ranging between 0.0 to 1.0, and the x and y-
axis show the spatial coordinates of the domain. The figure shown a curved flame
progress variable profile with flame elements at y/Ly = 0.50 reaching c = 1.0
before the flame elements at y/Ly = 0.00 and y/Ly = 1.00. The progress variable
is also well bounded from 0.0 to 1.0.
The stream-wise profile of c at four different transverse position (y/Ly =
0.5, 0.35, 0.25, and 0.00) is shown in figure 5.4b. Flame elements at y/Ly =
0.50 (blue line), which have the highest equivalence ratio, propagate the fastest
compared to the other three positions while flame elements at the transverse
periodic boundary (cyan line), which have the lowest reactant equivalence ratio,
propagate the slowest.
106
(a) Progress variable c field plot for case
A00 at te = 6.0 × 10−3s
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
x/Lx
c
 
 
y=0.50Ly
y=0.35Ly
y=0.25Ly
y=Ly
(b) Stream-wise profile of figure 5.4a at
various transverse positions.
(c) Local equivalence ratio Φ plot for
case A00 at t/tL = 2.2
0 0.2 0.4 0.6 0.8 1
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
x/Lx
Φ L
 
 
y=0.50Ly
y=0.35Ly
y=0.25Ly
y=Ly
(d) Stream-wise profile of figure 5.4c at
various transverse positions.
Figure 5.4: Field plots of Φ and c for case A00 at t/tL = 2.2.
107
The local equivalence ratio profile for case A00 at t/tL = 2.2 is plotted in
figure 5.4c. The z -axis is the local equivalence ratio Φ while the x and y-axes
represent the spatial coordinates of the domain. The prescribed inlet profile for
the equivalence ratio is clearly shown in the figure. The initially homogeneous
planar premixed flame is seen to have adjusted to the inhomogeneous reactants
prescribed at the inlet. As prescribed, the flame will experience an increase of
mixture strength in the middle of the transverse direction and will experience a
decrease (lowest local equivalence ratio prescribed) near to the transverse periodic
boundaries.
Similarly, the stream-wise profile of local equivalence ratio Φ at four different
transverse locations is plotted in figure 5.4c. Kinks in the profile are located at
about x/Lx = 0.20 for the stream-wise profile at y/Ly = 0.50, 0.35 and 0.25
in figure 5.4d. These kinks in the equivalence ratio corresponds to the position
of the flame. The location of the kink in local equivalence ratio agrees well
with the observation of the dip in equivalence ratio made in the one-dimensional
simulations.
Note that the flame near to the transverse boundaries, though it has a slight
kink in the equivalence ratio, is about to be quenched as the incoming equivalence
ratio is below the flammability limit.
5.3 Case B00 flame evolution
The evolution of the two-dimensional flame for case B00 is discussed in this
section. The position of the flame as defined by the iso-surface at c = 0.90 is
shown in figure 5.5 for case B00 at different elapsed times. It can be seen that
the curvature of the flame front is becoming more pronounced with increasing
elapsed time (blue solid line corresponds to the earliest time, blue dashed line
corresponds to the latest time). The increasing curvature of the c = 0.90 iso-
surface shows that all elements of the flame do not propagate at the same speed.
It is clear that the initially planar premixed laminar flame has been stretched out
by propagation at different speeds through the inhomogeneous field of reactants.
The effects of differential propagation as manifested in figure 5.5 provide evidence
confirming the hypothesis that there is a mechanism for the generation of flame
108
0.2 0.25 0.3 0.35 0.4
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x/Lx
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
Figure 5.5: Progress variable location (c = 0.90) for case B00
surface area in the presence of reactant mixture inhomogeneity.
The difference in propagation speed along the flame front increases the flame
surface area. The increase in flame surface area, or flame length in the present
two-dimensional context, has been quantified and is plotted against the elapsed
time in figure 5.6, where the total instantaneous flame surface (AΣ) is normalised
by its initial flame surface (AL) and plotted against the elapsed time for case
B00. The flame surface enhancement becomes significant (AΣ/AL > 1) at about
t/tL = 1.0. This is also observed in figure 5.5 where the c = 0.90 iso-surface
(purple solid line), becomes noticeably curved.
The flame surface area enchancement shows linear growth from t/tL ≥ 1.5,
which is consistent with the more curved position of the iso-surface as shown in
figure 5.5.
109
0 0.5 1 1.5 2 2.5 3
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
t/tL
A Σ
/A
L
Figure 5.6: Normalised total flame surface area (AΣ/AL) at c = 0.90 for case B00.
5.3.1 Time evolution of the equivalence ratio (Φ0.90) profile
The flame surface area for case B00 is shown to have increased with elapsed time
as shown on the c = 0.90 iso-surface plot in figure 5.5. With the flame curving
towards the reactants, the local equivalence ratio for the flame will inevitably
change from its initial value, hence displacement speed is expected to change
accordingly. The local equivalence ratio on the c = 0.90 iso-surface (Φ0.90) is
plotted in figure 5.7.
In this figure (5.7), the y-axis is the transverse direction of the simulation
domain while the x axis is the magnitude of local equivalence ratio on the c =
0.90 iso-surface (Φ0.90). Note that the Φ0.90 < ΦG, this is due to equivalence
ratio dip within the flame due to preferential diffusion (Refer to section 4.1 and
table 4.1). Flame elements near to the middle of the transverse direction (y/Ly =
0.50), which propagate into the richest part of the inhomogeneous mixture, are
shown to adapt to the oncoming changes and register the highest Φ0.90 magnitude.
Similarly, the flame elements closer towards the transverse periodic boundaries
y/Ly = 0.00, 1.00) are adapting to the oncoming reactants which were prescribed
to have the lowest value of Φ. This is consistent with the evolution of the c = 0.90
iso-surface where the iso-surface is curving towards the reactants.
110
0.55 0.6 0.65 0.7 0.75 0.8 0.85
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Φ0.90
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
Figure 5.7: Local equivalence ratio on the c = 0.90 iso-surface for case B00 at
different elapsed times.
5.3.2 Evolution of the surface-weighted averaged of equiv-
alence ratio (Φ0.90) profile
The evolving flame front observed in figure 5.5 has resulted in a non-constant
equivalence ratio on the isosurface (refer to figure 5.7). The effective local equiv-
alence ratio for the c = 0.90 isosurface at any given time can be represented using
the surface averaged quantity Φ0.90. The surface averaged local equivalence ratio
Φ0.90 against t/tL for case B00 is plotted in figure 5.8.
The time history started with Φ0.90 that corresponds to ΦG = 0.80. There is a
slight bump at t/tL ≈ 0.2 then the value of Φ0.90 rapidly decreases in magnitude.
The decrease happens in three stages; the first decreases almost linearly from
t/tL ≈ 0.25 to about 0.75, followed by another almost linear decrease but at a
slower rate for 1.00 < t/tL < 2.00 and the third for t/tL > 2.0 where the rate of
decrease is becoming more non-linear with time. These stages corresponds to the
stages of flame evolution which are discussed later in section 5.4.
111
0 0.5 1 1.5 2 2.5 3
0.66
0.68
0.7
0.72
0.74
0.76
t/tL
Φ 0
.9
0
Figure 5.8: The surface averaged local equivalence ratio Φ0.90 against t/tL for
case B00
5.3.3 Time evolution of the displacement speed profile
It is clearly shown that the distribution of Φ0.90 is changing with time, adapting
to the oncoming inhomogeneity in reactant mixture strength. The adaptation of
Φ0.90 gives rise to a non-constant Φ0.90 that will inevitably lead to a variation in
displacement speed Sd along the iso-surface. This variation of Sd along the flame is
a manisfestation of differential propagation speed. The displacement speed on the
c = 0.90 iso-surface is plotted in figure 5.9. Like figure 5.7, the y-axis represents
the transverse axis of the simulation domain while the x -axis represents Sd. Note
that the Sd plotted here is the magnitude of the displacement speed.
Figure 5.9 has shown that Sd along the iso-surface starts to vary as elapsed
time increases. The profile of Sd starts to resemble the equivalence ratio profile
as shown in figure 5.7. It is interesting to note that the Sd profile starts to
deviate from the Φ0.90 profile in figure 5.7 from t/tL > 1.5. Near to the periodic
boundaries, Sd has started to increase again while the value of Φ0.90 in that region
is showing a decreasing trend with time.
112
4 4.5 5 5.5 6 6.5 7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Sd/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
Figure 5.9: Displacement speed on the c = 0.90 iso-surface for case B00 at
different t/tL.
5.3.4 Time evolution of the displacement speed compo-
nents
The phenomenon of differential propagation along the iso-surface has been shown
figure 5.9. The displacement speed is expected to respond to the changes in
Φ0.90 monotonically. When there is a non-monotonic response of Sd towards the
evolution of Φ0.90, particularly for the region near to the periodic boundaries. It
is clear that additional factors are influencing the response of Sd. The underlying
factors for this behaviour are revealed by probing the different components of Sd.
The displacement speed components are plotted in figure 5.10. These com-
ponents are labelled on their corresponding sub-figures: Sd (figure 5.10a), Sn
(figure 5.10b), Sr (figure 5.10c), St (figure 5.10d), Scd (figure 5.10e), and Sdiff,z
(figure 5.10f). The y-axis is the transverse direction of the simulation domain
while the x -axis denotes the magnitude of the displacement speed components.
Note that the tangential component is non-trivial in the two-dimensional simu-
113
4 4.5 5 5.5 6 6.5 7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Sd/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(a) Sd
−20 −18 −16 −14 −12 −10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Sn/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(b) Sn
15 20 25
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Sr/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(c) Sr
−1 0 1 2 3
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
St/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(d) St
0.5 1 1.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Scd/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(e) Scd
−0.25 −0.2 −0.15 −0.1 −0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Sdiff,z/SL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
(f) Sdiff ,z
Figure 5.10: The snap-shots of displacement speed and its components for case
B00 at different elapsed times.
114
lation.
The normal diffusion component and reaction components, Sn and Sr, have
large magnitude but they are of opposite sign. The sum of these two components
has similar order of magnitude to that of Sd. The tangential component St and
the cross dissipation component Scd each contribute about a tenth of the Sd.
The diffusive components Sdiff,z have negligible contribution towards Sd. These
findings are similar to the one-dimensional simulations with the addition of the
St contribution due to curvature.
The sum of Sn+Sr is clearly contributing positively towards Sd on all elements
along the flame at all elapsed time. Similarly, Scd also contributes positively to
Sd in all regions at all the times since the start of the simulation. However,
the contribution from St is dependent on the curvature of the flame. The flame
has positive curvature near the middle of the transverse direction, therefore the
St contribution towards Sd is negative. A strong negative curvature is form-
ing towards the transverse periodic boundaries, and as expected, St contributes
positively in those regions.
It is observed that the main features of the Sd profile are predominantly
represented by Sn and Sr, especially for flame elements near y/Ly = 0.50. The
profile of Sn and Sr, both in aggregate and individually, behaves monotonically
towards the variation of Φ on the iso-surface. However, as the elapsed time
increases, St and Scd start to influence the shape of Sd, specifically on flame
elements near to the transverse periodic boundaries.
Figure 5.11 shows the value of ∆Sd normalised by Srefd for different time
instances. The value of ∆Sd is the difference between the two-dimensional flame
Sd values and the reference one-dimensional laminar unstrained premixed flame
Srefd . For example, if the two-dimensional flame elements are travelling at a faster
Sd than the reference counterpart, e.g. (Sd > Srefd ) the value of ∆Sd is positive
and vice versa. Details on this notation is given later in equation 5.8.
At the beginning of the simulation (t/tL < 0.50), the contribution towards
Sd comes from Sn, Sr, and Scd, with St ≈ 0. During this period, the initially
planar flame responds to the oncoming inhomogeneous reactants like the one-
dimensional simulations. From figure 5.11, the displacement speed for flame
elements at y/Ly = 0.50 is shown to register smaller values compared to flame
115
−0.2 0 0.2 0.4 0.6 0.8
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
∆Sd/Srefd
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
Figure 5.11: Normalised difference of displacement speed ∆Sd/Srefd on the c =
0.90 iso-surface for case B00 at different t/tL.
elements nearer to the transverse periodic boundaries. This is because the flame
near y/Ly = 0.50 is propagating into a richer mixture while the flame near to the
periodic boundaries is propagating into a leaner mixture.
After the initial response, the magnitude of St becomes non-trivial and starts
to contribute towards Sd as flame begins to have non-zero curvature. At t/te ≥
1.5, Sd for flame elements near to the transverse periodic boundaries is shown to
have increased in magnitude due to the contribution from St (figure 5.10d). This
profile is further accentuated by the non-monotonic Scd behaviour with respect
to Φ near that region (figure 5.10e). The combination of St and Scd in this
region (0.00 < y/Ly < 0.20, and 0.80 < y/Ly < 1.00) is the main reason for the
non-monotonic changes in Sd with Φ0.90 shown in figure 5.10a.
116
5.4 The two-dimensional flame evolution stages
The displacement speed Sd of flame elements does not behave monotonically with
time as observed in figure 5.9 where the Sd of flame elements located near to the
transverse periodic boundaries initially decreases in value but the trend is reversed
at t/tL = 1.50. The non-monotonicity between Sd and t/tL indicates that the
evolution of the two-dimensional flame happened in stages.
The values of Sd for flame elements at five different values of t/tL are plotted
against the local equivalence ratio Φ0.90 in figure 5.12. It is observed that the Sd
line encompasses a larger range of Φ0.90 with increasing t/tL; the flame is evolving
into a region of greater variation of Φ0.90. The dashed line in figure 5.12 represents
the steady homogeneous laminar flame Srefd . The two-dimensional flame Sd lines
(solid lines) are observed to deviate from the reference line (dashed line) in all
five instances. This shows that the two-dimensional flame is no longer a simple
ensemble of individually steady laminar flame elements.
0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
Φ
S d
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
reference
Figure 5.12: The displacement speed Sd against local equivalence ratio Φ for case
B00 at t/tL = 0.50 to 2.50 with an increment of 0.50. The black dashed lin is the
steady laminar flame reference displacement speed Srefd as a function of the local
equivalence ratio on c = 0.90 isosurface Φ0.90.
117
Although the two-dimensional flame Sd line deviates from the reference line
in figure 5.12, Sd still has a positive gradient with respect to the various Φ0.90
cases. The exception occurs at low values of Φ0.90 for higher values of t/tL, where
Sd shows negative gradient with respect to Φ0.90. From the observation of the
Sd − Φ0.90 relations, three stages of evolution can be identified. These stages are
1. Initial normal stratification stage
2. Flame stretching stage
3. Cusp formation stage
The initial normal stratification stage occurs when the initially planar flame is
required to adjust to the inhomogeneous reactant mixture it is propagating into.
At this stage, the flame elements of the two-dimensional flame resemble an ensem-
ble of unsteady one-dimensional flames that are propagating into a non-constant
mixture strength field.
At the end of the initial normal stratification stage, the flame will have dif-
ferential displacement speed along its surface. This marks the beginning of the
flame surface stretching stage, where the flame will have a non-trivial curvature.
As the flame stretches, cusps will form due to the occurrence of highly negative
curvature at positions where the mixture strength is the lowest (lowest value of
Φ0.90 along the surface). In the cusp formation stage, Sd near to the cusps will
have a strong relation to the curvature.
5.4.1 Normal stratification stage
The flame elements at the centre of the transverse domain propagate into a richer
mixture while flame elements near the transverse periodic boundaries propagate
into a leaner mixture. During the initial stage of the simulation, the flame ele-
ments will propagate into mixture with different strength (Φ), hence the flame
elements will have different Sd depending on the local value of Φ. Flame elements
with higher (lower) Φ will have a larger (smaller) value of Sd. During this phase,
the two-dimensional flame can be modelled as layers of one-dimensional flames,
each propagating into reactants with different mixture strength.
118
The values of Sd at t/tL = 0.50 and 1.00 for case B00 are plotted against Φ0.90
in figure 5.13.
0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95
0.8
1
1.2
1.4
1.6
1.8
2
Φ
S d
 
 
t/tL = 0.50
t/tL = 1.00
reference
Figure 5.13: The displacement speed Sd against local equivalence ratio Φ for case
B00 at t/tL = 0.50 and 1.00. The black solid line is the steady laminar flame
reference displacement speed Srefd as a function of the local equivalence ratio on
c = 0.90 iso-surface Φ0.90.
The displacement speed has a monotonic relation with Φ0.90 as expected.
The displacement speed for t/tL = 1.00 has lower values across the range of
Φ0.90 compared with t/tL = 0.50 and the gradients for both instances that lie
in the region of 0.70 < Φ0.90 < 0.80 are the same. The values of Sd differ from
the reference counterpart where a positive (negative) difference of Sd − Srefd is
observed for low (high) values of Φ0.90.
The overall larger Sd values across Φ for t/tL = 0.50 as shown in figure 5.13
occur because the flame has higher area-weighted average equivalence ratio Φ
than the flame at t/tL = 1.00 (refer to figure 5.8). The higher value of Φ yields
a greater equivalence ratio amplitude (Φ − Φin) across the flame, resulting in
a global positive shift of Sd. The inhomogeneous reactants average equivalence
ratio strength is set at Φin = 0.67.
These findings are consistent with the normal stratification effects observed
119
in the previous chapter. The deviation of flame-normal spatial gradient of equiv-
alence ratio from the reference counterpart is given as:
∆r
(∂Φ
∂n
)
= ∆r(∂nΦ) =
∂Φ
∂n −
(∂Φ
∂n
)ref
(5.7)
where ∆r denotes the deviation of flame properties from the reference counterpart.
A positive value of ∆r(∂nΦ) denotes that the flame has a less steep gradient when
compared with the reference counterpart. This is because ∂Φ/∂n has a negative
value for all flames at the c = 0.90 iso-surface. A less steep ∂Φ/∂n implies that
the flame is propagating into a mixture with higher equivalence ratio, and vice
versa.
Similarly, the deviation of Sd from the reference counterpart is given as:
∆r (Sd) = Sd − Srefd (5.8)
A positive ∆r (Sd) indicates that the flame on the c = 0.90 iso-surface is prop-
agating at a faster Sd than the reference steady counterpart that has the same
value of Φ0.90.
The unsteady flame under the normal stratification effect will register a neg-
ative (positive) value of ∆r(∂nΦ) when it is propagating into a mixture with
higher (lower) equivalence ratio. It is observed that a negative (positive) value of
∆r(∂nΦ) yields a positive (negative) value of ∆ (Sd). Mathematically:
∆r(∂nΦ) ∝ −∆r (Sd) (5.9)
The value of ∆r(Sd) is normalised by Srefd and ∆r(∂nΦ) is normalised by
(∂Φ/∂n)ref :
∆r,n(Sd) =
∆rSd
Srefd
= Sd/Srefd − 1
∆r,n(∂nΦ) =
∆r(∂nΦ)
(∂Φ
∂n
)ref = (∂nΦ)/(∂nΦ)
ref − 1
(5.10)
where ∆r,n denotes the normalised deviation of the flame property. The rela-
120
tion of the normalised Sd deviation (∆r,n(Sd)) and the normalised ∂nΦ deviation
(∆r,n(∂nΦ)) for case B00 at t/tL = 0.50 and 1.00 is shown in figure 5.14.
−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
∆r,n(∂nΦ)
∆
r,
n
(ρ
S d
)
 
 
t/tL = 0.50
t/tL = 1.00
Figure 5.14: The normalised Sd deviation (∆r,n(Sd)) against the normalised
∂Φ/∂n deviation (∆r,n(∂nΦ)) for Case B00 at t/tL = 0.5 and 1.0.
The x -axis for figure 5.14, ∆r,n(∂nΦ), is the percentage difference of ∂Φ/∂n for
case B00 from the laminar steady reference value. Flame elements that propagate
into a leaner (richer) mixture will have a positive (negative) ∆r,n(∂nΦ) value.
Therefore, the majority of the flame elements located near to the centre of the
domain will have ∆r,n(∂nΦ) < 0, whereas flame elements located near to the
transverse periodic boundaries will have ∆r,n(∂nΦ) > 0. The y-axis for figure 5.14,
∆r,n(Sd), measures the fractional difference of Sd for case B00 from the laminar
steady reference value. A negative ∆r,n(Sd) value indicates that the flame element
has a lower Sd value than the reference counterpart that has an equal value of
Φ0.90.
A linear relationship between ∆r,n(∂nΦ) and ∆r,n(Sd) is observed in figure 5.14
for the flame at t/tL = 0.50. This clearly shows that in this instance, the case B00
flame is predominantly influenced by the normal stratification effects. At t/tL =
1.00, ∆r,n(∂nΦ) and ∆r,n(Sd) start to deviate from the original linear relationship.
Moreover, non-monotonic behaviour of ∆r,n(Sd) is observed at high values of
121
∆r,n(∂nΦ). This indicates that there are additional factors which contribute to
the behaviour of ∆r,n(Sd) at high values of ∆r,n(∂nΦ). This also marks the end
of the normal stratification stage for case B00 flame.
The flame during this stage (0.00 ≥ t/tL < 1.00) is largely planar with cur-
vature not larger than |κδL| < 0.10. The small curvature confirms that normal
stratification effects are the dominant influence on the unsteady flame during this
period.
5.4.2 Flame stretching stage
The flame will undergo a stretching stage once the differential displacement speed
is established during the normal stratification stage. During the flame stretching
stage, which occurs for 1.0 < t/tL < 2.0, the adjacent flame elements will propa-
gate at a different Sd and the flame curvature will become non-trivial. The flame
surface area will also increase during this phase (refer to figure 5.6).
The normalised curvature κ× δL against the spanwise direction y/Ly for five
different instances are plotted in figure 5.15 for case B00.
−0.5 −0.4 −0.3 −0.2 −0.1 0 0.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
κ × δL
y/
L y
 
 
t/tL = 0.50
t/tL = 1.00
t/tL = 1.50
t/tL = 2.00
t/tL = 2.50
Figure 5.15: The normalised flame curvature κ× δL against transverse direction
y/Ly for case B00 at five different instances.
122
From figure 5.15, the curvature range is clearly observed to increase with
time. This is consistent with the observation that the flame surface increases
with time as shown in figure 5.6. This increasing κ trend indicates that the
curvature contribution becomes more prominent during the later stage of the
flame evolution. For example at t/tL = 2.50 (pink line), there are strong negative
values of κ towards the periodic boundaries (y/Ly =0 and 1), indicating the
formation of cusps.
At t/tL = 0.50, the curvature is slightly positive but the variation of curva-
ture across the flame is minimal. At this time, flame is undergoing the normal
stratification stage. The relatively small value of flame curvature confirms that
the flame is close to planar during the normal stratification stage. At the end of
the normal stratification stage, where t/tL = 1.0, a large variation in flame curva-
ture across the transverse direction is observed. The maxima and minima of the
flame curvature are each 10% of the value of the laminar flame thickness. The
non-trivial curvature indicates that the flame is starting to curve significantly.
During the flame stretching stage, where 1.0 < t/tL < 2.0, the variation of cur-
vature across the transverse direction increases. It can be observed that there
are very large negative value of curvature found towards the periodic boundaries
(y/Ly =0 and 1). The negative curvature is further enhanced in magnitude at
t/tL = 2.50, which culminates in the formation of cusp in this region. The cusp
formation stage is discussed in the next section.
The growth of curvature and the growth of flame surface area induce a finite
stretch rate on the flame surface. The flame surface stretch S˙ is given as:
S˙ = 1A
∂A
∂t = ∂(lnA)/∂t = aT + 2Sdκ (5.11)
where A is the instantaneous flame surface area. The surface-weighted averaged
of stretch is therefore the sum of aT and 2Sdκ. The surface-weighted averaged
stretch rate S˙ is plotted against t/tL in figure 5.16.
The surface averaged stretch rate S˙ is always positive, consistent with the
positive growth of the flame surface area shown in figure 5.6. The surface averaged
stretch rate S˙ is observed to hover about 45s−1 during the normal stratification
stage. The stretch is then observed to increase linearly with time during the flame
123
0 1 2 3
0.1
0.15
0.2
0.25
0.3
0.35
t/tL
a T
+
κS
d
×
δ L
/S
L
Figure 5.16: The surface averaged rate of flame surface stretch ∂(lnA)/∂t against
t/tL for case B00.
stretching stage (1.0 < t/tL < 2.0) before tapering off and decreasing with time
at about t/tL = 2.50. The behaviour of S˙ can be explained by the looking into
aT and 2Sdκ individually.
The surface averaged strain rate aT is plotted against t/tL in figure 5.17 where
aT has been smoothed to reduce high frequency wiggles that are caused by the use
of acoustically reflecting inlet boundary conditions in the simulation. The values
of aT shown in figure 5.17 hover about zero during the normal stratification stage
(0.0 < t/tL < 1.0), then increase linearly with time during the flame stretching
phase (1.0 < t/tL < 2.0). The growth of aT tapers off at about t/tL = 2.0 then
shows a decreasing trend for t/tL > 2.5.
It is observed in figure 5.17 that aT is small during the normal stratification
stage. This is consistent with the fact that the flame is largely planar during this
time, as observed from the small curvature range at t/tL = 0.5 and t/tL = 1.0
in figure 5.15. The surface averaged strain rate aT during the flame stretching
phase is observed to increase linearly with time. This observation is consistent
124
0 1 2 3
−0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
t/tL
a T
×
δ L
/S
L
Figure 5.17: The surface averaged strain rate aT against t/tL for case B00.
with the increasing curvature range observed in figure 5.15. Note that aT shows
a negative trend at the end of the flame stretching stage (t/tL = 2.5).
The surface average of the curvature component of the surface stretch 2κSd is
plotted against t/tL in figure 5.18 for case B00. The surface average of the product
of curvature and displacement speed 2κSd is generally positive throughout the
simulation, which indicates that this component contributes positively toward S˙.
Note that 2κSd is observed to be positive at the start of the simulation. This
is because the slight amount of curvature in this stage has been amplified by
displacement speed of different flame elements.
There is a dip in 2κSd during the normal stratification stage around t/tL =
0.50. This dip starts to recover going into the flame stretching stage. However,
2κSd is observed to decrease linearly with t/tL again at t/tL = 1.5, which is
during the flame stretching stage. The decrease becomes steeper after the flame
stretching stage.
The non-monotonic profile of 2κSd underlines that the κ and Sd relation
125
0 1 2 3
0.04
0.06
0.08
0.1
0.12
0.14
0.16
t/tL
κS
d
×
δ L
/S
L
Figure 5.18: The surface average of the curvature component of the surface stretch
2κSd against t/tL for case B00.
changes at different stages of the flame evolution. During the flame stretching
stage, 2κSd increases during the first half of the stage, then a decreasing trend is
observed at t/tL > 1.5.
5.4.3 Cusp formation stage
The flame stretching stage leads to the cusp formation stage as the curved flame
evolves. The onset of cusp formation is observed in figure 5.15 (pink line) at
t/tL ≥ 2.0. During this stage, the relation between Sd and Φ has become non-
monotonic (refer to figure 5.12), where there is a negative branch of Sd and Φ
at small value of Φ. This negative branch indicates that the cusp is starting
to form. Referring to figure 5.16, the mean stretch S˙ is starting to decrease in
value and both components of the stretch (aT and κSd) illustrate a decrease in
magnitude for t/tL ≥ 2.5. This is because the formation of cusp increases the
local displacement speed, which in turn reduces the stretch. Note that at this
126
stage, the curvature component of the stretch (refer to figure 5.18) decreases more
rapidly than the stretching stage (1.0< t/tL <2.0) and in a nonlinear fashion.
5.5 SDF scatter against progress variable across
the flame brush for case B00
The normalised surface density function |∇c| × δL scatter plot across the flame
brush for case B00 during the flame stretching stage (t/tL = 1.00) is presented in
figure 5.19.
0.0 0.2 0.4 0.6 0.8 1.0
c
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
|
∇
c
|
×
δ
L
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
Figure 5.19: Scatter plot of SDF σ = |∇c| across the flame brush coloured by
local equivalence ratio Φ. The mean of SDF conditioned on progress variable is
plotted in red asterisks.
The scatter plot shows that for a given isosurface, there are variations in
the SDF value. This indicates that the two-dimensional flame has varying flame
127
thickness. The inhomogeneous reactants are observed to have fully penetrated
the flame. Moreover, the local equivalence ratio has a clear correlation with the
SDF: where the upper part of the SDF scatter is associated with larger local
equivalence ratio while the lower part is associated with smaller local equivalence
ratio.
The SDF peak for flames with higher local equivalence ratio is located at about
c ≈ 0.7 while the SDF peak for flames with lower local equivalence ratio is located
slightly further towards the products relative to that of higher local equivalence
ratio counterpart. This is consistent with the one dimensional calculation (refer
to figure 4.7, where the peak of SDF is gradually shifted towards the products
with decreasing global equivalence ratio.
5.6 The budget of the progress variable trans-
port equations
The budget of the progress variable transport equation is plotted in figure 5.20
for case B00 at t/tL = 1.50, which is in the flame stretching stage of the two-
dimensional flame evolution. The progress variable transport equation (refer to
equation 3.34, can be recast to be the displacement speed equation (refer to
equation 3.40) as follows:
1
︷︸︸︷
Dρc
Dt =
2
︷ ︸︸ ︷
∂
∂xk
ρD ∂c∂xk
+
3
︷︸︸︷
ω˙c +
4
︷ ︸︸ ︷
2ρD∂ ln z∂xk
∂c
∂xk
+
5
︷ ︸︸ ︷
(1− c)
z
∂
∂xk
ρ(Dz −D)
∂z
∂xk
=
1
︷ ︸︸ ︷
ρ|∇c|Sd = ρ|∇c|(
2
︷ ︸︸ ︷
Sn + St +
3
︷︸︸︷
Sr +
4
︷︸︸︷
Scd +
5
︷ ︸︸ ︷
Sdiff,z)
(5.12)
The progress variable transport terms, which the displacement speed terms
were derived from, in equation 5.12 are labelled with the corresponding sub-figure
numbers below:
1. The combined unsteady and convection term, which is the displacement
speed term ρ|∇c|Sd, is plotted in figure 5.20a
128
2. The diffusion term, which is the sum of tangential and normal diffusion
term ∇(ρD∇c) = ρ|∇c|(St + Sn), is plotted in figure 5.20b
3. The reaction term ω˙c = ρ|∇c|Sr is plotted in figure 5.20c
4. The cross dissipation term 2ρD∇(ln z)∇c = ρ|∇c|Scd is plotted in fig-
ure 5.20d
5. The extra diffusive term (1− c)∇(ρ(Dz −D)∇z)/z = ρ|∇c|Sdiff,z is plotted
in figure 5.20e
The diffusion term is split into normal ρ|∇c|Sn and tangential ρ|∇c|St compo-
nents and they are plotted in figure 5.21a and 5.21b respectively. The combined
normal diffusion and reaction term ρ|∇c|(Sn + Sr) is plotted in figure 5.21c.
The material derivative of ρc, which is also the displacement speed term
ρ|∇c|Sd plotted in figure 5.20a shows a larger value at the middle (y/Ly = 0.5)
relative to the periodic boundaries (y/Ly =0.0, and 1.0). This is consistent with
the Sd against y/Ly plotted in figure 5.9 in which larger Sd on the c = 0.9 isosur-
face is found towards y/Ly = 0.5.
The diffusion term for reaction progress variable shown in figure 5.20b has an
inverted profile compared to the diffusion term for species mass fraction shown in
figure 5.3c. This diffusion term can be split into normal and tangential compo-
nents and they are plotted in figure 5.21a and 5.21b respectively. In the stream-
wise direction, the diffusion term is shown to have a positive contribution towards
Sd but changes sign and contributes negatively before returning to zero. Similar
to the displacement speed, the diffusion term registers its largest contribution in
the middle y/Ly = 0.5. This is also reciprocated by the reaction term.
The reaction term has a large positive contribution (almost twice the mag-
nitude of Sd) towards the displacement speed. It is positive due to the rela-
tionship ω˙c = −ω˙CH4/z. The reaction term registers its highest magnitude at
y/Ly = 0.50 and its lowest at y/Ly = 0.00 because at the middle of the domain,
the local equivalence ratio has the highest value as compared to that near to the
transverse periodic boundaries. This is consistent with the profile shown in the
diffusion term, where the diffusion term at the middle registers larger negative
magnitude than towards the periodic boundaries.
129
(a) Dρc/Dt = ρ|∇c|Sd
(b) ∇(ρD∇c) = ρ|∇c|(St + Sn) (c) ω˙c = ρ|∇c|Sr
(d) 2ρD∇(ln z)∇c = ρ|∇c|Scd
(e) (1 − c)∇(ρ(Dz − D)∇z)/z =
ρ|∇c|Sdiff,z
Figure 5.20: Budget for the r.h.s of progress variable transport equation for case
B00 in the flame stretching stage t/tL = 1.00 (1 of 2). Note that the SDF is
σ = |∇c| (refer to equation 3.46).
130
(a) −n.∇(ρD|∇c|) = ρ|∇c|Sn
(b) 2ρD|∇c|∇n = ρ|∇c|St
(c) −n.∇(ρD|∇c|)D+ ω˙c = ρ|∇c|(Sn+
Sr)
Figure 5.21: Budget for the r.h.s of progress variable transport equation for case
B00 in the flame stretching stage t/tL = 1.00 (2 of 2).
The cross dissipation term is plotted in figure 5.20d. There are two peaks
and a trough observed in the plot. The trough is found at the middle of domain
(y/Ly = 0.50) while the two peaks are found near the region of transition from
high to low local equivalence ratio (y/Ly = 0.33, and 0.66). Note that positive
magnitude of the cross dissipation term is found away from the peak and towards
the periodic boundaries. This behaviour is expected as the portion of the flame
in the middle of the domain is propagating into a stronger reactant mixture,
therefore registering a negative magnitude in the cross dissipation term. The
opposite is true for the portion of flame towards the sides of the domain. This
is consistent with the findings of the previous chapter, where an instantaneous
131
locally negative spatial gradient of equivalence ratio along the flame normal vector
will have a positive value of Scd, and vice versa.
The extra diffusion term represents the diffusion of mixture fraction with the
diffusion coefficient Dz−D. This diffusion coefficient is the difference between the
diffusion coefficient for the mixture fraction z and that for the progress variable c.
This term will disappear if the mixture fraction diffusion coefficient is the same as
the fuel mass fraction diffusion coefficient. The importance of this term relies on
the spatial gradient of the mixture fraction or the magnitude difference between
the diffusion coefficients. Also, this term has a 1−c factor (refer to equation 3.40).
Therefore the magnitude is relatively smaller towards the products.
The contribution of the extra diffusion term is about three times smaller than
the cross-dissipation term (term 5 in equation 5.12), and about two orders smaller
than the diffusion and reaction terms. However its contribution relative to other
terms is larger compared to the one-dimensional flame simulation. This is because
the gradient of the mixture fraction is larger in the two-dimensional cases.
The normal component of the diffusion term ρ|∇c|Sn plotted in figure 5.21a
is observed to be the major contributors towards the diffusion term. Its profile
closely resembles that of the diffusion term in figure 5.20b. The tangential com-
ponent of the diffusion ρ|∇c|St plotted in figure 5.21b is observed to have about
10% contribution towards the diffusion term. It has a positive value towards the
periodic boundaries y/Ly =0.0 and 1.0 because of negative curvature in these
regions. It has negative values in the middle y/Ly = 0.50 because of positive
curvature in the middle.
When the normal diffusion and reaction components are combined (refer to
figure 5.21c), they are observed to be the major contributors towards Sd in the
two-dimensional flame. This shows that the primary driver for the differential
propagation comes from Sn + Sr.
132
5.7 Scatter plot of displacement speed and its
components coloured by local equivalence
ratio
The scatter plot of displacement speed in terms of ρ|∇c|Sd and its components
against progress variable c across the flame brush for case B00 during the flame
stretching stage (t/tL = 1.50) is shown in figure 5.22. The mean behaviour
conditioned on progress variable is also plotted in the figure using red markers.
The scatter plot is coloured by its local equivalence ratio.
The scatter plot of ρ|∇c|Sd (refer to figure 5.22a) resembles that of the SDF
shown in figure 5.19. A strong preference of larger local equivalence ratio in the
upper part of the scatter is observed, similar to the SDF scatter. This indi-
cates that the displacement speed in the stretching phase (1.0< t/tL <2.0) has a
monotonic relation with the local equivalence ratio, as observed in figure 5.12.
Monotonicity of local equivalence ratio and the combined normal and reac-
tion component ρ|∇c|(Sn + Sr) is clearly observed in the scatter in figure 5.22b.
Compared to ρ|∇c|Sd, the peak of the mean of ρ|∇c|(Sn + Sr) is slightly shifted
towards the reactants. This can be attributed to the cross dissipation component
ρ|∇c|Scd in figure 5.22d, where the peak of the mean of ρ|∇c|Scd is found towards
the products. This underscores the importance of Scd as an additional source of
normal diffusion towards Sd
The mean behaviour of the cross dissipation component ρ|∇c|Scd has a mini-
mum at c ≈ 0.2 and maximum at c ≈ 0.8. From the scatter plot colour scale, the
cross dissipation component is observed not to have monotonic relation with local
equivalence ratio as there are abrupt changes of colour near the mean behaviour.
The non-monotonicity of Scd indicates that this term is influenced not only by
ΦL, but also the geometry of the flame, as observed in figure 5.10e. However,
the higher local equivalence ratio (Φ > ΦG) has a preference towards the lower
part of the scatter while the lower equivalence ratio (Φ < ΦG) has a preference
towards the upper part of the scatter. This is in agreement with the one dimen-
sion laminar flame where flames that propagate into leaner mixture have negative
Scd while flame that propagates into richer mixture has positive Scd (refer to fig-
133
0.0 0.2 0.4 0.6 0.8 1.0
c
0
100
200
300
400
500
600
700
800
ρ
|
∇
c
|
S
d
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
(a) ρ|∇c|Sd
0.0 0.2 0.4 0.6 0.8 1.0
c
0
200
400
600
800
ρ
|
∇
c
|
(
S
n
+
S
r
)
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
(b) ρ|∇c|(Sn + Sr)
0.0 0.2 0.4 0.6 0.8 1.0
c
−100
−50
0
50
100
150
ρ
|
∇
c
|
S
t
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
(c) ρ|∇c|St
0.0 0.2 0.4 0.6 0.8 1.0
c
−50
0
50
100
ρ
|
∇
c
|
S
c
d
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
(d) ρ|∇c|Scd
0.0 0.2 0.4 0.6 0.8 1.0
c
−60
−50
−40
−30
−20
−10
0
ρ
|
∇
c
|
S
d
i
f
f
,
z
Φ
0.50
0.55
0.60
0.65
0.70
0.75
0.80
0.85
0.90
0.95
1.00
(e) ρ|∇c|Sdiff,z
Figure 5.22: Scatter plot of r.h.s progress variable transport equation terms
against progress variable.
134
ure 4.11e on page 63). On the c = 0.90 isosurface, the contribution of Scd is about
20% of Sd.
The mean behaviour of the extra diffusion component ρ|∇c|Sdiff,z (refer to
figure 5.22e) is shown to be non-zero, and it has a non-negligible contribution
towards displacement speed for the progress variable range of 0.1 < c < 0.8.
However, the extra diffusion contribution on the c = 0.90 isosurface remains
negligible.
5.8 The surface density function transport equa-
tion
The balance equation for the surface density function (SDF) |∇c| is given as:
1
︷ ︸︸ ︷
∂|∇c|
∂t +
∂uk|∇c|
∂xk
=
2
︷ ︸︸ ︷
(δij − ninj)
∂ui
∂xj
|∇c|+
3
︷ ︸︸ ︷
Sd
∂nk
∂xk
|∇c|
4
︷ ︸︸ ︷
−∂Sdnk|∇c|∂xk
(5.13)
Like the budget of the progress variable balance equation shown in figure 5.20,
the surface density balance equation budget is plotted in figure 5.23 for case B00
during the flame stretching stage at t/tL = 1.50. These two budget plots are
of the same case and the same instant of time. The surface density balance
equation terms in equation 5.13 are labelled with the corresponding sub-figure
numbers below:
1. The unsteady and convection terms, figure 5.23a.
2. The curvature term, figure 5.23b.
3. The strain rate term, figure 5.23c.
4. The propagation term, figure 5.23d.
The combined unsteady and convective SDF term is shown in figure 5.23a.
A positive value of this term denotes the generation of SDF. It is observed that
generation of SDF is located towards the reactants while two troughs, which are
associated with the desctruction of SDF, are located towards the products.
135
(a) ∂|∇c|∂t +
∂uj |∇c|
∂xj (b) Sd ∂nk∂xk |∇c|
(c) (δjk − njnk) ∂uj∂xk |∇c| (d) −
∂Sdnk|∇c|
∂xk
(e) residual ǫ
Figure 5.23: Surface density function transport equations budget
136
The curvature term plot is shown in figure 5.23b. The curvature term is
positive for flame elements located in the middle of the transverse direction. It
shows that the flame is curved convex towards the reactants. Positive increase in
flame surface area in this region is observed. The curvature term shows strong
negative values towards the periodic boundaries, indicating the onset of cusp
formation and the destruction of SDF in that region. This is consistent with the
St plot shown in figure 5.10d. The total contribution from the curvature term is
in the same order as the peak of the strain rate term but an order lower than the
propagation term. Although small, the contribution from the curvature is finite,
given that the flame was initialised to be planar.
The strain rate term is plotted in figure 5.23c. In the current snapshot, there
is a high compressive strain rate towards the middle of the flame (y/Ly = 0.50)
followed by a strong extensive strain for flame elements adjacent to it. Aside from
these three peaks, the value of the strain rate term remains small, approximately
an order smaller than the propagation term.
The propagation term in the SDF balance equation is plotted in figure 5.23d.
It clearly shows that flame surface area is generated towards the reactants and
destroyed towards the products. The propagation term also reflects the behaviour
of Sd, where the peak found at the middle of the flame (y/Ly = 0.50) corresponds
to the highest Sd registered along a given iso-surface. Also, the high Sd value
towards the periodic boundaries due to the onset of the formation of cusps is also
shown, with a smaller peak in the plot. The propagation term contribution is an
order larger than the other two r.h.s terms of equation 5.13. This is because the
flame is laminar in nature and the contribution of SDF generation from strain
rate and curvature terms are small.
The residual of the budget is plotted in figure 5.23e. The residual is given by:
ǫ = ∂|∇c|∂t +
∂uk|∇c|
∂xk
− (δij − ninj)
∂ui
∂xj
|∇c| − Sd
∂nk
∂xk
|∇c|+ ∂Sdnk|∇c|∂xk
(5.14)
The residual is at least three orders smaller than the other terms in the SDF
balance equation. The highest residual is registered within the flame region, in
particular at y/Ly ≈ 0.33 and 0.66, where the maximum transverse gradient of
equivalence ratio is located. The small residual provides additional confidence in
137
the simulation results and the terms presented.
5.9 Scatter of the surface density function trans-
port equation across the flame brush.
The scatter of the terms in the surface density function transport equation (refer
to equation 5.13) is plotted against progress variable in figure 5.24 for case B00
during the flame stretching stage at t/tL = 1.50. The mean behaviour conditioned
on progress variable is also plotted in the figure in red markers. The scatter plot
is coloured by the local equivalence ratio.
The combined unsteady and convection term for the SDF transport equation
is plotted against progress variable in figure 5.24a. The mean behaviour shows
positive values towards the reactants (0.0< c <0.7) but negative values towards
the products (0.7< c <1.0). This is consistent with the budget in figure 5.23a.
Towards the products, the combined term is observed to have a monotonic re-
lation with the local equivalence ratio, whereby larger local equivalence ratio
corresponds to more negative values of the term for a given isosurface. How-
ever, this relation is not observed towards the reactants, where the larger local
equivalence ratio does not coincide with a larger contribution.
The strain rate term is plotted against progress variable in figure 5.24b. The
mean variation of the strain rate term conditioned on the progress variable across
the flame brush is observed to be negligible. The strain rate term scatter plot is
not showing a monotonic relation with the local equivalence ratio. This is consis-
tent with figure 5.23c where extensive strain rate is observed near the transition
region (y/Ly ≈ 0.33 and 0.66). Strong extensive strain rate is observed for flames
with larger local equivalence ratio (y/Ly = 0.50). Regions close to the periodic
boundaries, where the local equivalence ratio is smallest, have negligible strain
rate (denoted by blue coloured dots in figure 5.24b).
The curvature term contribution is plotted against progress variable in fig-
ure 5.24c. The mean variation of the curvature term conditioned on progress
variable across the flame is found to be about zero. In regions where the flame
elements are negatively curved (concave towards reactants), which corresponds to
138

5.10 Mixture fraction profile along the c=0.90
isosurface
The carbon