The bending effect in
turbulent flame propagation
Girish Venkata Nivarti
Darwin College
University of Cambridge
A thesis submitted for the degree of
Doctor of Philosophy
6 June 2017

Declaration
This dissertation is the result of my own work and contains nothing which
is the outcome of work done in collaboration except where specifically indi-
cated in the text. It has not been submitted for another qualification to this
or any other university. It contains nearly 33,300 words, 48 figures and 3 tables.
Girish Nivarti
1 June 2017

Abstract
In the present thesis, the sensitivity of flame propagation to the turbulent
motion of burning gases is investigated. The long-standing issue of the ‘bending
effect’ is focused upon, which refers to the experimentally-observed inhibition
of flame propagation velocity at high intensities of turbulence. Plausible mech-
anisms for the bending effect are investigated by isolating systematically the
effects of turbulence intensity. By providing a novel perspective on this topic,
the thesis addresses the fundamental limits of turbulent burning.
The investigation employs Direct Numerical Simulation (DNS), which en-
ables the basic conditions of burning to be controlled directly. A parametric
DNS dataset is designed and generated by increasing turbulence intensity over
five separate simulations. Effects of turbulent motion are isolated in this man-
ner, such that the bending effect is reproduced in the variation of flame propa-
gation velocity recorded. Subsequently, the validity of Damköhler’s hypotheses
is investigated to ascertain the mechanism of bending.
Analysis of the DNS dataset highlights the significance of kinematic flame
response in determining turbulent flame propagation. Damköhler’s first hy-
pothesis is found to be valid throughout the dataset, suggesting that the
bending effect may be a consequence of self-regulation of the flame surface.
This contradicts the dominant belief that bending occurs as a result of flame
surface disruption by the action of turbulence. Damköhler’s second hypothesis
is found to be valid in a relatively limited regime within the dataset, its va-
lidity governed by flame-induced effects on the prescribed turbulent flow field.
Therefore, this thesis presents turbulent flame propagation and the bending
effect as emergent from the dynamics of a flame surface that retains its internal
thermo-chemical structure.
Finally, experimental validation is sought for the proposed mechanisms of
bending. Comparisons have been initiated with measurements in the Leeds ex-
plosion vessel, based on which the widely accepted mechanism of bending was
hypothesized twenty-five years ago. Modifications to the DNS framework war-
ranted by this comparison have aided the development of novel computationally-
efficient algorithms. The ongoing work may yield insights into the key mech-
anism of the bending effect in turbulent flame propagation.
Acknowledgements
I value all the love I have received these four years in a land so far from home.
Foremost, I owe a deep debt of gratitude to my supervisor, Professor Stewart
Cant, who created with arms wide open a space for me to learn. It is in
this space that I was able to develop my ideas and find my personal brand
of research, influenced by his style and backed by his courage. Still Stewart’s
portrayal of science as a dull and mechanistic procedure remains simply ironic,
because science is just so much fun with him!
I am grateful to my teacher, Salil Srivastava, for imparting the foundations
and beauty of science, Professor Kendal Bushe for the guidance up to Cam-
bridge, and to Professor Epaminondas Mastorakos who conceived the basic
idea on our very first meeting that culminated eventually into this thesis. My
academic drive has benefited from the presence of Professor Simone Hochgreb,
and I thank her for all the conversations we have shared. I am also grateful to
Professor Andrei Lipatnikov for bridging the decades between us with passion
and wine ever since Birmingham. Lastly, I am inspired deeply by Professor
Derek Bradley and his incredible crew in Leeds. I am grateful to them, and to
him for accommodating us in his academic life.
At Darwin, I have learned much in the shadows of fellows and friends. I
thank sincerely Professor Andy Fabian for stargazing and conversations on
turbulence, and Professor Sir Harry Bhadeshia for his kind mentorship and
relentless humour. The content of this thesis has benefited significantly from
the conversations with Dr. Giles Shaw, Dr. Matteo Corso and all members of
Gazebo Society, and Dr. Richard Armstrong and the regulars at DarBar.
A special note of thanks goes to Eli-ran Bar-el and Stefanos Roimpas for the
life-changing conversations over fresh air and the foundation of Darwin Ring,
to Loren Held for sharing the love for turbulence, to Dr. Stoyana Alexandrova
– an ageless inspiration, to the residents of Wordsworth Grove, and to all the
Porters at Darwin. Elsewhere in Cambridge, I found, to my great fortune, Erik
Garrison over magic flights in chaos and aperiodic time crystals, Dr. András
Bárány over linguistic coffees, Claudia Bedin at the Scooter, Christian Owusu
whose love I will carry until the next flight to Cambodia, and Space Cadets
who have been an artistic inspiration.
Credit is due to my seniors from early days Ryan Griffiths, Dr. Pranay
Seshadri, and Dr. Ahmed Al-Shabab; I thank my colleagues and friends from
Ashby Laboratory: Andrea Masi for Brindisi, Ashley Scillitoe for Arsenal, and
Camille Bilger for Camille, Marseille, and the I♥Stewie group. Last but not
at all the least, I am grateful to Peter Benie for sharing with me a piece of the
wonder that he is, and to Katia Babayan for Pushkin.
The generous funding support of Cambridge International Trust has been
vital to my life in Cambridge. I thank also the Cambridge Philosophical So-
ciety and the Cambridge Society for Applied Research (CSAR) for their kind
donations and support.
My brothers in crime – Aamod Shanker and Ramkumar Ramachandra – are
always with me no matter how distant their habitats. I cherish the good times
with the core team members: Mikesh and Shubha. Friendships from Vancouver
have been inspirational throughout my time here; thank you Professor Artan
Seshmani, Dr. Dan Brox, and Sneha Sheth for nourishing these bonds beyond
Green College. My life with Annika Böddeling has been precious to me . . . she
is simply my favourite!
The love of my family empowers me; thank you mummy, papa and deedee.
to India:
the home of 0 and
the land of ∞

Contents
List of Figures xv
List of Tables xxi
Preface xxiii
0 Introduction 1
1 Theoretical Background 5
1.1 Concepts of turbulent burning velocity . . . . . . . . . . . . . . 5
1.1.1 Damköhler: the hypotheses of 1940 . . . . . . . . . . . . 7
1.1.2 Aerodynamic stretch and flame stretch . . . . . . . . . . 11
1.1.3 Regimes of premixed flame propagation . . . . . . . . . . 14
1.1.4 Bradley: the Leeds database . . . . . . . . . . . . . . . . 17
1.2 Simulations of turbulent burning velocity . . . . . . . . . . . . . 20
1.2.1 Model equations for flame dynamics . . . . . . . . . . . . 20
1.2.2 Peters: the Thin Reaction Zones regime . . . . . . . . . 25
1.2.3 What do Damköhler’s hypotheses mean today? . . . . . 27
1.2.4 Thesis objectives and approach . . . . . . . . . . . . . . 29
2 The Direct Numerical Simulation method 31
2.1 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.4 Benchmark solutions . . . . . . . . . . . . . . . . . . . . . . . . 47
xi
Contents
3 Parametric Simulations of the Bending Effect 55
3.1 Design conditions for bending . . . . . . . . . . . . . . . . . . . 55
3.1.1 Flame configuration and regime . . . . . . . . . . . . . . 58
3.2 Recorded effects of turbulence intensity . . . . . . . . . . . . . . 63
3.2.1 Images of the turbulent flame brush . . . . . . . . . . . . 63
3.2.2 Computation of turbulent burning velocity . . . . . . . . 73
3.3 Summary: on the bending effect recorded . . . . . . . . . . . . . 74
4 Damköhler’s First Hypothesis 77
4.1 Flame Surface Density approach . . . . . . . . . . . . . . . . . . 77
4.1.1 Computation of flame surface area . . . . . . . . . . . . 83
4.2 Stretch rate components . . . . . . . . . . . . . . . . . . . . . . 87
4.2.1 Statistics of strain and propagation terms . . . . . . . . 88
4.2.2 Balance of strain and propagation . . . . . . . . . . . . . 93
4.3 Summary: the role of DH1 in bending . . . . . . . . . . . . . . 96
5 Damköhler’s Second Hypothesis 97
5.1 Scalar flux transport in turbulence . . . . . . . . . . . . . . . . 98
5.1.1 Gradient transport hypothesis . . . . . . . . . . . . . . . 99
5.2 Counter-gradient transport . . . . . . . . . . . . . . . . . . . . . 101
5.2.1 Regime of counter-gradient transport . . . . . . . . . . . 106
5.2.2 Computation of the diffusion coefficient . . . . . . . . . . 109
5.3 Summary: validity regime of DH2 . . . . . . . . . . . . . . . . . 112
6 Towards Validation with the Leeds Experiment 115
6.1 Leeds: expanding spherical explosions . . . . . . . . . . . . . . . 115
6.1.1 Preliminary comparisons . . . . . . . . . . . . . . . . . . 117
6.2 Cambridge: freely propagating flames . . . . . . . . . . . . . . . 118
6.2.1 Multi-step chemical kinetics . . . . . . . . . . . . . . . . 118
6.2.2 Fast adaptive timestep control . . . . . . . . . . . . . . . 123
6.2.3 Linear turbulence forcing . . . . . . . . . . . . . . . . . . 127
6.2.4 Computational setup and regime . . . . . . . . . . . . . 130
6.3 Summary: the ongoing comparative study . . . . . . . . . . . . 134
xii
Contents
7 Conclusion 135
Appendix A 141
References 143
xiii

List of Figures
1.1 Photographs of turbulent Bunsen flames studied by Damköhler
(1940). Contrasting regimes of turbulent flame propagation are
observed upon increasing the Reynolds number (between top
and bottom rows). . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Illustration of a propagating flame surface interacting with a
field of vorticity. At large-scales (red vortical regions), the in-
teractions lead to wrinkling and aerodynamic stretch, whereas
at relatively small scales, internal gradients (zoomed inset) may
be modified leading to flame stretch. . . . . . . . . . . . . . . . 12
1.3 Regime diagram proposed by Borghi (1985) to categorise topo-
logical and structural effects of turbulence on flame propagation. 15
1.4 The explosion vessel investigated by Abdel-Gayed and Bradley
(1977) in Leeds to quantify effects of turbulence on flame propa-
gation. Four fans generate nearly homogeneous isotropic turbu-
lence in the reacting mixture filling the vessel volume. Turbu-
lent flame propagation subsequent to successful ignition in the
mixture can be observed through quartz windows. . . . . . . . . 17
1.5 Regime diagram as proposed by Peters (1999) demarcating the
extent of the Thin Reaction Zones (TRZ) regime using Ka = 1
and Kaδ = 1 isolines. . . . . . . . . . . . . . . . . . . . . . . . . 26
1.6 Regimes proposed by Damköhler (1940) illustrated on the regime
diagram proposed by Peters (1999). DH1 (`0  δL) and DH2
(`0  δL) correspond to the flamelet and Broken Reaction Zones
(BRZ) regimes respectively. . . . . . . . . . . . . . . . . . . . . 28
xv
List of Figures
2.1 Laminar propagating 1D flame solution modelled using single-
step kinetics and unity Lewis number. A spatially-thin reaction
zone (left) is simulated. The peak reaction rate occurs at c = 0.8
(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.2 Grid resolution study to determine the bounds for flame chem-
istry resolution. Convergence fails for grids with resolution
∆x & δL/5 where δL = 0.36 × 10−3m is obtained a priori on
the finest grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3 Simulation snapshots of a counter-clockwise rotating vortex ad-
vected through a 2D inflow-outflow domain with inlet velocity
u0 = 5ms
−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.4 Turbulent kinetic energy computed in a periodic 128×128×128
HIT simulation. Panels show a) snapshot of energy spectrum
in Fourier space and b) temporal decay of kinetic energy. . . . . 52
3.1 Inflow-outflow flame configuration used in the simulations. Left:
3D rendering of a propagating flame isosurface interacting with
upstream turbulent eddies in the domain. Right: 2D axial slice
showing the location of the flame and relative size of the domain. 58
3.2 Regime diagram showing the region of parameter space cov-
ered by the present simulations. The entire region between the
flamelet regime and broken reaction zones regime is spanned by
increasing u′ successively across five simulations labelled I – V. . 62
3.3 Contours of progress variable (colour) for each of the five simu-
lations (columns, with u′ value noted on top). Time evolution
of the respective flame occurs down each column, with elapsed
time noted on the left in each row. . . . . . . . . . . . . . . . . 64
3.4 Instantaneous snapshots of the turbulent flame leading edge c =
0.01 (coloured blue) and trailing edge c = 0.99 (coloured red)
in cases with the lowest and highest u′, respectively. The arrow
indicates the direction of flame propagation. . . . . . . . . . . . 66
xvi
List of Figures
3.5 Distributions of the instantaneous mean curvature hm and cur-
vature shape factor Sh for the leading edge (solid lines) as well
as trailing edge (dashed lines) surfaces in cases I and V (black
and red lines, respectively). . . . . . . . . . . . . . . . . . . . . 67
3.6 Magnified view of the instantaneous leading edge (blue) and
trailing edge (red) surfaces at t = 4τ0 as computed for case V
with u′ = 30sL. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.7 Reaction rate contours ω˙(x, z; y = 5L/6) in cropped axial slices
of the domain (above) and reaction rate scatter ω˙(c) across the
entire domain (below) at t = 4τ0 for cases I and V. . . . . . . . 71
3.8 Progress variable isosurface c∗ = 0.8, corresponding to the peak
reaction rate, in the lowest and highest intensity cases. Axial
sections of the domain show the variation of reaction rate ω˙
across the domain. . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.9 Computed variation of the normalized turbulent burning veloc-
ity sT/sL with normalized turbulence intensity u′/sL. The inset
shows a typical variation of sT (t) and the corresponding value
of sT chosen from this variation. . . . . . . . . . . . . . . . . . . 74
4.1 Expectation of FSD Σ as a function of the axial coordinate x
(left) and mean progress variable c (right) for three cases (I,
III, and V) spanning the dataset. The laminar flamelet solution
(red line) is shown for comparison. . . . . . . . . . . . . . . . . 79
4.2 Illustration of sampling methods (dashed lines). Left: sampling
over spanwise sections of the domain across the mean turbulent
flame brush (c¯∗0 < c¯ < c¯∗N). Right: sampling over instantaneous
progress variable c isosurfaces (solid lines) within the instanta-
neous turbulent flame brush (c∗0 < c < c∗N). . . . . . . . . . . . . 80
4.3 Instantaneous normalized Σ′(c) scatter (left column) at t = 4τ0
compared with the laminar flamelet profile (red line). Corre-
ponding contours of Σ′(x) (right column) are shown in an axial
section overlaid with progress variable contours (black and white
lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
xvii
List of Figures
4.4 Computed flame surface area ratio AT/AL compared with nor-
malised turbulent burning velocity sT/sL in the recorded bend-
ing effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5 Comparisons of turbulent flame structure and propagation char-
acteristics simulated using 1-step and 25-step reduced chemical
kinetic models for case III with u′ = 10sL. . . . . . . . . . . . . 86
4.6 Distribution of strain at and curvature hm over the peak reaction
rate isosurface c∗ = 0.8 for three cases (lowest u′: solid line,
highest u′: dotted line) spanning the entire dataset. . . . . . . . 88
4.7 Joint probability distribution of strain rate at (scaled with Tay-
lor scale strain rate u′/λ) and mean curvature hm (scaled with
Kolmogorov scale η) in the lowest and highest intensity cases. . 90
4.8 Distribution of displacement speed sd over the peak reaction
rate isosurface c∗ = 0.8 in three cases spanning the dataset
(solid line: lowest intensity and dotted line: high-intensity). . . . 91
4.9 Joint probability distribution of sd (scaled with laminar flame
speed sL) and hm (scaled with Kolmogorov scale η) in the lowest
and highest intensity cases. . . . . . . . . . . . . . . . . . . . . . 92
4.10 Variation of volume-integrated source-terms a and h of the trans-
port equation for Σ′ with increasing normalized turbulence in-
tensity u′/sL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.1 Schematic showing gradient (G) transport of the progress vari-
able flux. A disturbed flame surface (light dashed line) is wrin-
kled further (progressively darker solid curves) by vortical activ-
ity. The flame does not counteract this tendency, and is simply
convoluted by turbulent fluctuations. . . . . . . . . . . . . . . . 100
xviii
List of Figures
5.2 Schematic of the counter-gradient (CG) transport mechanism.
At sufficiently high heat release factors τ = ρu/ρb, the flame has
a tendency to act against gradient transport. An initially wrin-
kled flame surface (dashed convoluted curve) may un-wrinkle as
the lighter product peninsula is differentially accelerated in the
direction of products (along red arrows, away from the direction
of propagation). . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.3 Variation of mean scalar flux u˜′′c′′ with the mean scalar gradi-
ent ∂c˜/∂x. Colours denote corresponding values of the mean
progress variable c˜ (blue: leading edge c˜ = 0.01, red: trailing
edge c˜ = 0.99). . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.4 Instantaneous progress variable contours (black lines) and Favre
velocity fluctuations u′′ (background colours) in an axial cross-
section y = L
2
;L < x < 2L, z of the domain. . . . . . . . . . . . 105
5.5 Recorded correlations between u′′ and c′′ for case I (above) and
case V (below). Colours indicate the progress variable value
(blue: reactants, red: products) and size is proportional to the
reaction rate magnitude. The transition from positive to nega-
tive correlation corresponds to the G/CG transition. . . . . . . . 107
5.6 The Bray number NB for each case simulated in the present
dataset delineating the gradient/counter-gradient transport regimes
at unity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.7 Computed normalized turbulent burning velocity sT/sL com-
pared with predictions
√
DT/DL of DH2 in the gradient diffu-
sion regime NB < 1. DH2 is invalidated in the CG transport
regime denoted by NB > 1 (shaded red). . . . . . . . . . . . . . 111
6.1 Preliminary comparisons of methane-air flame propagation in
experiments (Leeds) with φ = 0.6, and DNS (Cambridge) with
φ = 1.0. The same turbulence intensity u′/sL is prescribed
(panel a) in both cases. Curvature distributions are normalized
using the Taylor length scale λ. . . . . . . . . . . . . . . . . . . 117
xix
List of Figures
6.2 Profiles of thermo-chemical variables: temperature and chemical
components (left) of the composition vector as well as radicals
(right) computed using a 25-step chemical kinetic model. . . . . 121
6.3 Grid resolution study to determine the spatial resolution re-
quirement for the 25-step methane-air chemical kinetic model.
Convergence fails for grids with ∆x & δL/22.3, where δL =
9.3mm denotes the value obtained a priori using a fine grid. . . 122
6.4 Timestep history (panels a,b, and c) as well as timestep ratios
(panel d) obtained using the PID controller αβγ and the mod-
ified PI4.2 controller α∗β∗ due to Söderlind (2002) for various
test problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
6.5 Turbulent forcing strategy. Left: time evolution of mean tur-
bulent kinetic energy k(t) for an HIT simulation with Q = 0
and with Q = ε/(2k0). Right: region of turbulent forcing in the
flame simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . 129
6.6 Resolution constraints of chemistry (black dashed lines) and tur-
bulence (red dashed line) as well as the corresponding compu-
tational expense (black and red solid lines, respectively) as a
function of DNS domain size. . . . . . . . . . . . . . . . . . . . 131
6.7 Combustion regimes corresponding to the Leeds experiment and
the Cambridge DNS in the regime diagram with a modified
definition of Ka. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
xx
List of Tables
3.1 Thermo-chemical input parameters. . . . . . . . . . . . . . . . . 57
3.2 Turbulence properties and dimensionless parameters. . . . . . . 60
6.1 Input parameters for the Leeds-Cambridge study. . . . . . . . . 133
xxi

Preface
Understanding combustion has been essential for sustaining progress this cen-
tury. Be it our long-anticipated trip to Mars, or tomorrow’s BA flight to
Mumbai, human exploration hinges upon controlled turbulent combustion. At
home, fire is necessary to cuisine, and far out in the cosmos, our very existence
appears to be tied to turbulent flames. Indeed, today we are re-discovering
our age-old relationship with combustion in its myriad forms. Using lasers,
supercomputers and advanced mathematics, the complexity of these forms is
being revealed once again. Unified conceptualisation of the turbulent flame is,
therefore, being pursued with great optimism and curiosity.
A mode of intense and efficient burning is found to occur when fuel and ox-
idiser are mixed prior to ignition. The drawback, however, is that subsequent
propagation of the ‘premixed’ flame is extremely sensitive to a variety of mix-
ture properties. Turbulence is one such property; the chaotic dynamics that
arise in a turbulent gaseous medium has manifold consequences in flame prop-
agation. Damköhler (1940) conjectured in his foundational hypotheses that
turbulence may only enhance burning. Yet over the 75 years since, evidence
for inhibitive effects of turbulence has grown. So, it is believed increasingly
that sufficiently-intense turbulence may extinguish a premixed flame. Bradley
(1992) summarized these possibilities in a review entitled: ‘How fast can we
burn?’
The present thesis addresses this very question using supercomputers. The
non-linear partial differential equations governing the motion of reacting gases
are solved on thousands of processing cores. As a result, the complex in-
teractions of the flame and turbulence are unveiled for scrutiny. By thereby
re-investigating the foundational hypotheses of flame propagation, the thesis
xxiii
Preface
contributes to a unified conception of the turbulent flame.
xxiv
Chapter 0
Introduction
Scientific conceptions of the ‘flame’ seek to reconcile two distinct spatio-
temporal perspectives. At large scales, the flame appears as a coherent surface
with well-defined physical properties. At relatively small scales, the flame is
described rather as a process that results in thermo-chemical transformation.
Dynamics of the flame at these distinct scales govern in unison how fast a
laminar gas mixture burns. In a turbulent burning mixture, the large- and
small-scale perspectives of the flame are interlaced by the motion of gases at
intermediate scales. As a result, difficulties arise in reconciling the ‘external’
surface and the ‘internal’ process perspectives of the flame. These difficulties
are addressed in the present thesis by isolating and studying systematically
the various scales of interactions between flame and turbulence.
Moments before the piston reaches top dead centre in an SI engine, a spark
ignites the compressed fuel-air mixture. Subsequently, a thin ‘premixed’ flame
propagates through the mixture, expanding it rapidly to restore the piston.
Propagation of the flame surface, therefore, results in producing mechanical
work. Mallard and le Chatelier (1883) pioneered a framework for calculating
this ‘burning velocity’ of the propagating flame by linking it with various mix-
ture properties. Interpreting this problem as being akin to heat conduction in
thermal waves, Michelson (1889) formulated the burning velocity as a function
of the temperature gradient across the flame surface. Hence, the perspective
of the premixed flame as a coherent surface wedged between unburned and
burned mixture states was born.
1
0. Introduction
Considered merely as a surface of thermal discontinuity, however, the flame
lacks the properties ascribed to it otherwise by virtue of combustion chemistry.
The activation energy required to spark-ignite a mixture, for instance, arises
from a chemical requirement. During the subsequent flame propagation, this
energy is sustained by exothermic chemical reactions. Indeed, under conditions
in which this chemical energy is lost, the flame may even extinguish. Evidently,
much of the sensitivity of the flame is linked to its chemical nature.
Lewis and von Elbe (1934) proposed, therefore, that the internal structure
of the flame is determined not only by thermal but also by molecular gradients.
Establishing this formally, Zel’dovich and Frank-Kamenetski˘i (1938) showed
that a relatively-thin ‘reaction zone’ depletes molecules within the internal
flame structure. This region releases heat whereas the preparatory ‘pre-heat
zone’ sustains the balance between thermal and molecular diffusion. These
studies provide a strong basis for the internal thermo-chemical structure that
arises from the flame considered as a process.
The narrative so far addresses the burning velocity of laminar flames, where
the motion of gases does not impart new scales to the problem of flame prop-
agation. Mixtures of gases, however, are often turbulent. In fact, the very
act of mixing itself is effected most commonly by the turbulent ‘cascade’ in
the first place (Kolmogorov, 1941). The various cascading scales of turbulent
motion blur the distinction between the large surface scales and the small pro-
cess scales of the flame. This poses a severe challenge to the formulation of a
consistent theory of burning velocity in turbulent flame propagation.
So in order to address turbulent flame propagation, Damköhler (1940) con-
ceived of two extreme cases. In the first of his seminal hypotheses, all scales
of turbulent motion are considered to be large, so that only the external sur-
face attributes of the flame are affected. In the second hypothesis, all scales
of turbulence are imagined to be small, so that only the internal process and
the associated structure is affected. Therefore, Damköhler’s conception divides
turbulent flame propagation into manageable regimes that provide a basis for
the prediction of turbulent burning velocity. In each of these regimes, turbu-
lence is hypothesized to ‘act’ on the flame to ‘enhance’ the burning velocity.
2
0. Introduction
This conception has underpinned the majority of investigations of turbulent
flame propagation undertaken to date.
With the surge of aircraft technology in the 1940s, afterburners were in-
vented. These help improve thrust, but may incur accidental flame blow-off at
sufficiently high speeds. Karlovitz et al. (1953) theorizes that the latter may
result from increased localized flame quenching due to steep velocity gradients
in high-speed gaseous flow. Hence, the notion of ‘stretch’ was conceptualised
so as to accommodate for the effects of gas motion inhibiting burning. This is
generalized by Klimov (1963) – turbulent flame propagation was envisaged to
be the consequence of interactions between ‘laminar flames in turbulence.’
Markstein (1964) and Williams (1985) improved upon this view, suggesting
that components of the burning mixture may mediate the effects of turbulence
on the flame. Furthermore, the possibilities that arise from increasing degrees
of turbulence-flame interaction are illustrated by Borghi (1985) in a ‘regime’
diagram. Meanwhile, the integrity in turbulent flow of the internal flame struc-
ture – the so-called ‘flamelet’ – was questioned increasingly. The possibility
of flame quenching in turbulence carries the implication that small-scale tur-
bulent motion disrupts the structure of the flamelet, thereby dismantling the
coherence of the flame surface. This challenges the entire conception, since
Damköhler (1940), of a propagating turbulent flame as a laminar flame af-
fected by turbulence. Hence, for decades to follow, the scientific community
became pre-occupied with the search for conditions which might limit this
hypothesis.
Abdel-Gayed and Bradley (1977) employed an explosion vessel in which the
intensity of turbulence can be controlled in a manner similar to Damköhler’s
theoretical conception. Correlating data from over a thousand experiments, it
was suggested that though large-scale turbulent motion enhances flame prop-
agation as per Damköhler’s first hypothesis, sufficiently small-scale turbulence
may, in fact, inhibit flame propagation. The resulting variation of burning
velocity over the entire range of turbulence regimes was termed ‘the bending
effect’ as the graph tapered away from linearity in the small-scale regime.
Bradley (1992) conjectures that ‘bending’ occurs because velocity gradients
3
0. Introduction
of small-scale turbulence cause disruption of the internal thermo-chemical pro-
cess of the flame. Yet Driscoll (2008) notes that in the twenty-five years since
Bradley’s conjecture, little direct evidence of total flame structure disruption
due to small-scale turbulence has emerged. So, the regime of ‘Thin Reaction
Zones’ has been proposed by Peters (1999), suggesting mitigated damage to
the flamelet hypothesis. This line of investigation continues to date, effects of
turbulence and the degree of loss of flamelet structure being seen as the central
problem of turbulent flame propagation.
The present investigation complements this dominant ongoing narrative.
In the systematic simulations conducted here, it is evident that as the variety
of turbulence scales become prominent, self-regulation of the flame follows
suit. Hence, the role of flame-response appears to be significant in determining
overall turbulent flame propagation. The bending effect emerges, therefore, as
a competition between turbulent enhancement of the flame surface and the
self-inhibition of the flame. In other words, it is not necessary for localized
flame quenching to occur to observe the bending effect.
Furthermore, Damköhler’s forgotten second hypothesis is examined. A re-
lated theme is observed here as a limiting factor of burning, namely the effect
of the flame on its surrounding turbulent flow. The resulting flame-induced ac-
celeration of burned gases limits the enhancement of turbulent burning velocity
in the small-scale turbulence regime. The present thesis highlights, therefore,
that turbulent flame propagation and the bending effect are governed not only
by the effects of turbulence on the flame, but equally by the effects of the flame
on itself as well as on turbulence.
In Chapter 1, the theoretical background is discussed beginning with the
foundational hypotheses laid out by Damköhler (1940). Furthermore, the
bending effect and Bradley’s work is elaborated upon. The direct simulation
approach employed to complement Bradley’s work is discussed in Chapter 2.
The simulation dataset is designed and generated in Chapter 3, and is anal-
ysed for the validity of Damköhler’s hypotheses in Chapter 4 and Chapter 5.
An outline of a collaborative investigation with Bradley and co-workers is pre-
sented in Chapter 6.
4
Chapter 1
Theoretical Background
1.1 Concepts of turbulent burning velocity
Beginning with Damköhler (1940), the turbulent burning velocity sT has been
seen as a conceptual extension of the laminar burning velocity, which may be
measured or computed for laminar flames. Indeed, a theoretical solution for the
burning velocity can be obtained analytically for only steady one-dimensional
laminar premixed flames (Zel’dovich and Frank-Kamenetski˘i, 1938). For this
idealised case, the steady-state equation for temperature can be solved to re-
late the laminar burning velocity with the thermal diffusivity α and the peak
reaction rate ω˙max of the burning mixture as
s0L ∼
√
α ω˙max/Ea, (1.1)
provided the mixture has a large activation energy Ea and a unity Lewis num-
ber Le ≡ α/D, i.e. equal thermal diffusivity α and molecular diffusivity D.
These turn out to be reasonable simplifications in most contexts related to
gaseous turbulent burning. Hence, the expression derived by Zel’dovich and
Frank-Kamenetski˘i (1938) provides a foundational building block towards the
formulation of expressions for the turbulent burning velocity sT .
Libby et al. (1979) have noted, however, that sT may have several concep-
tually different definitions to begin with. For instance, is sT is to reflect a local
or a global measure of burning rate? Is it to imply a displacement speed of the
5
1. Theoretical Background
flame surface, or the consumption speed of the flame brush? One may derive
a variety of formulations for sT along each of these definitions. Furthermore,
the sensitivities of the turbulent burning velocity to various factors, such as
mixture thermo-chemical properties, may be vastly different than the laminar
case. These sensitivities may vary considerably between different experimen-
tal configurations. As a result, difficulties arise in framing the problems of
turbulent flame propagation in a unified manner. The conceptually-distinct
burning velocity definitions have, in fact, been treated as the same in the past,
resulting in graphs that contain a large amount of scatter (Driscoll, 2008).
To aggravate the problem, certain features of turbulent flows that ren-
der challenging the computation or measurement of turbulent burning veloc-
ity. Typical configurations for turbulent flames do not lend themselves in a
straightforward manner to the computation of burning velocity as an eigen-
value, as is possible in the case of one-dimensional steady laminar flames (Pe-
ters, 2000). Furthermore, a variety of time and length scales are involved in
turbulent flames, implying that the assumptions inherent to Equation 1.1 may
not always be applicable to turbulent flames (Kerstein, 2002). Finally, the
equations for turbulent flame propagation are evidently more complex to the
extent that a direct solution approach is tedious (Poinsot, 1996). Due to these
reasons, the problem of turbulent burning velocity remains to date without a
consistent overarching formulation despite its significance to the field of com-
bustion (Bray, 1990; Lipatnikov and Chomiak, 2002; Driscoll, 2008; Bradley
et al., 2011).
Explanations of many phenomena related to turbulent flame propagation
remain contested today because of the above-mentioned challenges. The fo-
cus of this thesis – the bending effect – presents such an example, bearing
fundamental importance to the theory of turbulent burning. The bending ef-
fect refers to the variation of sT with increasing turbulence intensity u′ of the
unburned mixture. Although a ‘bending’ of sT (u′) has been witnessed and
quantified in many investigations (Lipatnikov and Chomiak, 2002; Driscoll,
2008), it continues to be debated by several researchers to date. In the follow-
ing, the conceptual background of this problem is presented, beginning and
6
1.1. Concepts of turbulent burning velocity
ending with an evaluation of the seminal work of Damköhler (1940).
1.1.1 Damköhler: the hypotheses of 1940
A basic apparatus for investigating the turbulent burning velocity is the Bun-
sen burner (Lockermann, 1956). In this apparatus (shown in Figure 1.1), pre-
mixed gaseous mixtures are issued from a vertical pipe such that a flame may
be stabilised at the burner rim. As a result, the flame propagates constantly
against the direction of reactants flowing out of the pipe. The turbulent burn-
ing velocity is, therefore, related directly to the global mass transformation
rate of the unburned gas m˙ in this context, so that
sL,T =
m˙
ρuYuA0
, (1.2)
Figure 1.1 Photographs of turbulent Bunsen flames studied by Damköhler
(1940). Contrasting regimes of turbulent flame propagation are observed upon
increasing the Reynolds number (between top and bottom rows).
7
1. Theoretical Background
for a known inlet mass flow rate ρuYuA. Here, ρu is the unburned gas density,
Yu is the mixture fuel mass fraction, and A0 is the cross-sectional area of the
pipe. Equation 1.2 holds for laminar L as well as turbulent T flames (note
the subscripts in Equation 1.2). Damköhler (1940) exploited this equivalence,
comparing the value obtained for turbulent flames with the laminar case.
In Figure 1.1, the laminar flame (top row) is seen as having a distinct
surface whose area AL can be measured to a fair degree of accuracy. If every
point in the laminar Bunsen flame satisfies Equation 1.1, the overall burning
velocity may be computed as
sL ∼ AL
√
D ω˙max, (1.3)
upon substituting α with the molecular diffusivity D. Turbulent flames, char-
acterised by a high Reynolds number Re ≡ UL/ν, have a relatively ambiguous
structure and shape in comparison (as seen in the bottom row of Figure 1.1).
Damköhler (1940) proceeded, therefore, to relate these differences in terms
of length scales expected in a turbulent flow relying on concepts proposed by
Prandtl (1925). Measured values of sT were compared with the computation of
sL based on Equation 1.3 in order to highlight differences which were then used
to elaborate physical mechanisms of turbulent flame propagation. It was noted
that sT > sL in all cases, leading to the conclusion that ‘turbulence always
increases the flame speed referred to the average flow cross section.’ Further,
the effects of turbulence were separated into two regimes, hypothesizing the
enhancement due to turbulence to occur ‘in a two-fold manner.’
Coarse-body turbulence
The first of these two ‘manners’ concerns a regime where the turbulence length
scales ` are large compared to the laminar flame thickness δ0L = ν/s0L, i.e. for
` δ0L. Addressing this regime, Damköhler (1940) proposed that
Hypothesis 1 ‘for coarse-turbulence bodies exceeding the laminar flame thick-
ness, the flame surface is effectively enlarged, which is, by roughening [the
flame surface.] Thus the coarse-body turbulence effects a rise of the flame
8
1.1. Concepts of turbulent burning velocity
velocity, . . . , even if the flame velocity in the individual surface elements is
entirely the same as in laminar flow.’
wherein the flame velocity of individual surface elements refers to the local
burning velocity s0L or, equivalently, the local laminar flame speed. Damköhler
(1940) argued further that the local velocity fluctuations u′ which give rise
to the ‘roughening’ balance the enhanced effective local burning velocity sT ,
stabilizing the flame at the burner rim. This led to a linear relation sT ∼ u′.
Hence, Damköhler’s first hypothesis (hereby, DH1) may be summarised as
follows
sT
sL
=
AT
AL
=⇒ sT ∼ u′ (1.4)
where AT is the enlarged turbulent surface area and AL is the appropriate flow
cross-section area. In this regime, it is the enhanced area ratio AT/AL that
acts to magnify the laminar flame speed sL to its equivalent sT in turbulence.
The statement of DH1 assumes, therefore, that the internal flame structure
is preserved in a manner that leaves unchanged the thermo-chemical terms α
and ω˙max in Equation 1.3. Framed in this manner, DH1 forms the basis of
nearly all turbulent flame modelling to date.
Fine-body turbulence
The second ‘manner’ of sT enhancement – Damköhler’s second hypothesis
(here called DH2) – concerns a limiting regime of high u′ turbulence wherein
the turbulence length scales `  δ0L. Addressing this regime, which contrasts
with DH1, Damköhler (1940) proposed that
Hypothesis 2 ‘for fine-body turbulence, . . . , an amplification of every micro-
scopic transport process in the flame surface must occur, especially between the
reaction zone and the pre-heat zone. These transport processes [then] . . . govern
the flame speed decisively.’
wherein the microscopic transport processes referred to may be molecular dif-
fusion or thermal conduction. Reasoning along these lines, Damköhler (1940)
described the enhancement of sT in ‘fine-body’ turbulence using the formula-
9
1. Theoretical Background
tion
sT
sL
=
√
DT
DL
=⇒ sT ∼
√
u′`0 ∼
√
Re, (1.5)
which is obtained by assuming that small-scales do not affect the area or reac-
tion rate in Equation 1.3. The area AL is not affected because the turbulence
length scales are considered to be too small, whereas the reaction rate ω˙ does
not change because the time scale of reaction is considered to be too small
following Zel’dovich and Frank-Kamenetski˘i (1938). The statement of DH2
further applies the mixing length hypothesis of Prandtl (1925), by setting
DT ∼ u′`0 so that sT in this regime may also be computed for a specified u′.
It is implicit in this statement that the ratio DT/DL remains positive upon
specifying the model for DT . This ensures that the mathematical validity of
DH2 is preserved.
The discussion above highlights that Damköhler (1940) not only conjec-
tured the regimes of turbulent flame propagation, but also provided formula-
tions for the burning velocity variation sT (u′). By relating these regimes and
the respective formulations to a few parameters such as u′ and `0, a platform
for modelling flame propagation was pioneered. Soon after, Schelkin (1943)
developed further the model of turbulent burning velocity sT as a function of
turbulence intensity u′. In an attempt to unify the two regimes, an expression
was provided as
sT
sL
=
√
1 +
(
u′
sL
)2
, (1.6)
which heralded subsequently an era of seeking a unified model of sT (u′) appli-
cable to all regimes of combustion (Bray, 1990). These attempts continue to
date as discussed by Lipatnikov and Chomiak (2002). The expectation is that
scaling relationships satisfied by sT (u′, `0) reflect intimately the underpinning
mechanisms that influence turbulent flame propagation. A primary focus has,
therefore, been to find better fits to the experimental sT variation so that the
related propagation mechanisms may be revealed.
Damköhler’s conception may be summarised as follows. Firstly, the tur-
bulent burning velocity differs from the laminar burning velocity due to the
effect of turbulence on the flame. Secondly, the burning velocity attained in
10
1.1. Concepts of turbulent burning velocity
turbulence is always enhanced in comparison to the laminar burning veloc-
ity. Lastly, the level of enhancement is governed by flame surface area and
turbulent diffusivity in two contrasting regimes of turbulence.
1.1.2 Aerodynamic stretch and flame stretch
A decade after Schelkin (1943), inhibitive effects of velocity gradients in the
burning mixture were considered by Karlovitz et al. (1953). The notion of
‘stretch’ was developed and quantified in terms of
K =
δ0L
U
∂U
∂y
, (1.7)
a function of velocity gradient ∂U/∂y modulated further by the flame thickness
δ0L and mean inflow velocity U . Attributing to stretch a negative effect on the
laminar burning velocity, Karlovitz et al. (1953) were able to predict laminar
flame extinction in a Bunsen burner. Klimov (1963) generalised this idea to
obtain the stretch rate
s˙ =
1
A
dA
dt
, (1.8)
as the fractional rate of change of flame surface area A over time. Clearly, this
expression relies on the notion of the flame surface area, so that the concept of
stretch necessarily assumes DH1 to be valid in some sense. A persistent issue
with stretch-based turbulence modelling has, therefore, been the difficulty in
being able to address the invalidation of DH1 due to flame quenching from
within the framework of DH1. Nevertheless, the work of Williams (1985) has
provided a strong mathematical basis for stretch.
The stretch rate s˙ defined in Equation 1.8 can be decomposed into the
component prescribed by the flow and the component experienced as a conse-
quence of the flame response (Williams, 1985). This decomposition is written
as
s˙ = at + sdhm, (1.9)
where the ‘prescribed’ component takes the form of tangential surface strain
at = ∇ · u− n̂ n̂ : ∇u, (1.10)
11
1. Theoretical Background
and the component ‘experienced’ due to the flame response occurs as a result
of normal propagation of the surface elements with mean curvature
hm =
1
2
∇ · n̂, (1.11)
where n̂ is the local flame surface normal vector (pointing in the direction of
unburned reactants) and sd is the surface displacement speed in this direction.
Stretch rate thus obtained has been useful in describing a range of flame phe-
nomena (Law, 1989). However, the effects of stretch, as well as the relationship
between its components, are subtle and at times difficult to comprehend.
Extensive strain rate at > 0, for instance, increases the flame surface area
by definition, but it may also reduce internal gradients causing the local dis-
placement speed sd to diminish. Strain has thereby proved vital in character-
ising extinction probabilities in flame propagation (Abdel-Gayed et al., 1984;
Bradley et al., 1992). Figure 1.2 shows how both large-scale surface deforma-
tion (due to aerodynamic stretch) as well as small-scale process modification
(due to flame stretch) may operate at the same time in a turbulent flame.
Mean Flow
Flame Propagation
Reactants Products
Figure 1.2 Illustration of a propagating flame surface interacting with a field
of vorticity. At large-scales (red vortical regions), the interactions lead to
wrinkling and aerodynamic stretch, whereas at relatively small scales, internal
gradients (zoomed inset) may be modified leading to flame stretch.
12
1.1. Concepts of turbulent burning velocity
Further, strain may also effect part of the local flame surface curvature hm
(Pope, 1988), so that the two stretch components are not entirely independent
of each other.
In a similar manner, the curvature component hm can also play a sup-
portive or inhibitive role on the local flame propagation velocity. Markstein
(1964) characterised the role of mixture diffusive properties (through the Lewis
number Le) in regulating flame surface curvature effects on flame propagation.
Unifying these concepts, it was proposed (Clavin, 1985) that, for flame thick-
ness δ0L much smaller than the characteristic wrinkle size, the local normal
displacement speed sd varies according to
s0L − sd
s0L
=
1
A
dA
dt
δ0L
s0L
Ma (1.12)
where the Markstein number Ma is related to the mixture Lewis number Le.
In turbulent flames, Clavin and Williams (1979) derived the large-scale
behaviour of turbulent premixed flame propagation for mildly stretched flames.
As a more general formulation, Bray and Cant (1991) proposed that flame
stretch effects may be absorbed simply into DH1 using a factor I0 . 1 under
low levels of stretch, such that
sT
sL
= I0
AT
AL
, (1.13)
thereby modifying DH1 without invalidating its basic premise. Further, in
order to address the difficulties in computing the exact area AT of convoluted
flame surfaces, a statistical approach was adopted. The mean Flame Surface
Density (FSD, denoted by Σ) is defined, for example, as the mean flame surface
area per unit volume allowing DH1 to be written (Bray and Cant, 1991) as
sT
sL
= I0
1
AL
∫
V
Σ dv, (1.14)
which is the same as Equation 1.4 with AT replaced by the integral of FSD
over the flame brush volume V . This description has helped account to some
degree for the inhibitive turbulent stretch effects on the local chemical structure
13
1. Theoretical Background
(Driscoll, 2008). It remains unclear, however, what effects precisely lead to this
variation of I0. The predominant role is believed to be the action of strain,
which may quench flamelets locally (Abdel-Gayed et al., 1984). Curvature
effects have been largely neglected, in comparison, and considered to ‘even
out’ in the context of turbulence (Bradley et al., 1992).
Accounting for the increasingly disruptive effects of turbulence, a more
general class of sT (u′) models as that developed by Schelkin (1943) emerged.
In an atmospheric pressure environment, a non-linear relationship
sT
sL
∼
(
1 +
(
u′
sL
)n)(1/n)
(1.15)
has been suggested to hold, proposing the exponent n as a parameter that
may fit the scattered data. These models attempted to address the various
bending curves of sT (u′) observed in experimental measurements by Karpov
and Severin (1980); Kido et al. (1989) as well as by Gülder (1991); Kobayashi
et al. (1996); Shy et al. (2000b). The precise variations captured in these
datasets have been reviewed in detail by Lipatnikov and Chomiak (2002). One
such work is focused upon here: that of Bradley et al. (1992) and co-workers.
A phenomenological conception of turbulent flame propagation developed at
nearly the same time helped visualise subsequently the physical mechanisms
responsible for such variations in sT (u′).
1.1.3 Regimes of premixed flame propagation
Following Damköhler (1940), the classification of flame propagation regimes
has been pursued based on the comparison of length scales involved in turbu-
lent reacting flows (Borghi, 1985). Effects of turbulent stretch are considered
in this framework as the result of an interaction between the various length
scales ` of turbulence and those of flame chemistry (such as δ0L). The turbu-
lence scales can be computed based on the k-ε model of turbulence (Jones and
Launder, 1972). The flame speed sL (referred to in DH1) is then expected to
be ‘the same as in laminar flow’, i.e. correspond to s0L, only if the Kolmogorov
14
1.1. Concepts of turbulent burning velocity
Figure 1.3 Regime diagram proposed by Borghi (1985) to categorise topological
and structural effects of turbulence on flame propagation.
scale
η =
ν3
ε
, (1.16)
and all scales larger than it, is larger than δL (where ε is the turbulent dissi-
pation rate). The Karlovitz number
Ka ≡ δ2L/η2 (1.17)
is a dimensionless marker that delineates the regimes of turbulence-flame in-
teraction by comparing the smallest turbulent microscale η with the flame
length scale δL. In essence, the size of the Kolmogorov scales of turbulence is
determined by the Reynolds number
Re0 = u
′`0/ν (1.18)
for a given integral length scale `0 and kinematic viscosity ν. Similarly, large-
scale convective time scale may be compared with the flame time scale, using
the Damköhler number
Da ≡ τf
τc
=
sL`0
u′δL
(1.19)
15
1. Theoretical Background
where τf = `0/u′ and τc = δL/sL. Many different diagrams have been proposed
all loosely based on the variety of turbulence flame interactions possible for dif-
ferent ranges of these dimensionless groups (Borghi, 1985; Abdel-Gayed et al.,
1989). A classical diagram is shown in Figure 1.3 with various visualizations
of the flame structure associated with the respective regimes.
In the Borghi diagram as Figure 1.3 is called, isolines of unity Ka, Da and
Re are used to demarcate the regimes of turbulence-flame interaction (Borghi,
1985). The ‘flamelet’ regime satisfies the Klimov-Williams criterion Ka < 1
and, hence, DH1 is expected to be valid in this regime based on the phe-
nomenological description. With increasing disruption of the local internal
flame structure, different regimes of flame propagation are evident until at the
highest intensities of turbulence, ‘distributed’ combustion prevails. Many ex-
periments have reported a significant departure from the linear relation sT ∼ u′
implied by DH1 in turbulent premixed flames with Ka > 1. This so-called
‘bending effect’ of the sT (u′) curve has been isolated to occur when all param-
eters but u′ are maintained constant within the parameter space consisting
of u′, `0, sL, and δL (Bradley, 1992; Lipatnikov and Chomiak, 2002; Driscoll,
2008). It was inferred that stretch effects begin to play an important role at
high intensities, where the deviation from linearity is observed in sT (u′).
Renormalisation group theories attempted to use such disparities in time
scales to find the scaling behaviour of sT in various regimes of combustion.
A rigorous approach pursued by Yakhot (1988) led to a relatively complex
relationship of the form
(
sT
sL
)2
log
(
sT
sL
)
=
(
u′
sL
)2
, (1.20)
that has also been derived heuristically (Kerstein et al., 1988) without con-
sidering the detailed internal structure. However, with the development of
detailed internal chemical structure (Williams, 2000) of the premixed flame,
the renormalisation group theories are faced with the challenge of accounting
for various other scales in the problem, so that these have not been able to
account successfully for the bending effect.
16
1.1. Concepts of turbulent burning velocity
Figure 1.4 The explosion vessel investigated by Abdel-Gayed and Bradley
(1977) in Leeds to quantify effects of turbulence on flame propagation. Four
fans generate nearly homogeneous isotropic turbulence in the reacting mixture
filling the vessel volume. Turbulent flame propagation subsequent to successful
ignition in the mixture can be observed through quartz windows.
1.1.4 Bradley: the Leeds database
In England, motivated by the need to verify the inhibitive effects of small-scale
turbulence, Andrews et al. (1975) devised an explosion vessel. The principal
advantage of the vessel is the ability to generate homogenous isotropic turbu-
lence using four fans located along the vertices of a tetrahedron as shown in
Figure 1.4. Flame kernels initiated using a spark in the centre would prop-
agate spherically in the field of turbulence, and by varying the speed of the
fans, turbulence intensity u′ may be regulated. By avoiding direct contact
with walls and indirect interaction through wall shear generated turbulence,
the experiment was able to reduce anisotropies in turbulence. Abdel-Gayed
and Bradley (1977) followed with experiments of turbulent flame propagation
under increasing turbulence intensity in this vessel.
Summarising the sensitivity of turbulent burning velocity sT to various
17
1. Theoretical Background
mixture properties, Abdel-Gayed et al. (1985) presented the relation
sT
sL
= f
(
u′
sL
,Re0,Le, β
)
(1.21)
where u′/sL is the normalised turbulence intensity used in the Borghi diagram
in Figure 1.3, Re0 is the Reynolds number defined in Equation 1.18, Le is the
mixture Lewis number and β, the Zeldovich number, is a function of the acti-
vation energy Ea. Abdel-Gayed et al. (1987) conducted the same experiments
under a wide range of equivalence ratios and turbulence intensities; a bending
effect in sT (u′/sL) was evident. Further still, at high enough intensities, the
premixed flame would fail to grow. This supported the hypothesis that the
bending effect occurs due to local flame quenching.
Bradley et al. (1992) sought to characterise this flame quenching due to
the contribution of tangential strain rate to the overall stretch rate. Hence, a
strain rate scale was assumed to be related to the Taylor scale λ as provided
by Batchelor (1952). The Karlovitz number based on this scale Ka = δ0L/λ is
derived as
Ka = 0.157 Re
− 1
2
0
(
u′
sL
)2
. (1.22)
Further, along the early work of Wohl et al. (1953), Abdel-Gayed et al. (1989)
pointed out the role of Lewis number in accentuating the effect of hydrody-
namic strain on turbulent premixed flames. Hence, quenching criteria were
developed based on correlations with 1650 experiments as follows
KaLe ≥ 6.0 (1.23)
for which flames quench globally. Here, the Karlovitz number Ka is meant to
denote the role of strain and the Lewis number Le that of the thermo-chemical
mixture properties. Curvature and flame propagation effects are neglected.
Indeed, in the expression of Equation 1.9 in terms of flame kernel radius R
s˙ =
(
U
R
+∇tU
)
+
sd
R
, (1.24)
where the first term includes the strain rate term computed using the tan-
18
1.1. Concepts of turbulent burning velocity
gential derivative operator ∇t and the second is the curvature term. Bradley
et al. (1992) hypothesised that ‘as flow and turbulence velocities increase, the
bracketed terms will increase relative to [the second term].’ The bracketed
terms representing prescribed flow field effects such as tangential strain were
hypothesised to dominate in comparison to the flame curvature term. Based
on this, it suggested that bending marks the onset of a departure from the lo-
cally flamelet structure – the foundational basis of DH1 – leading to quenching
of the flamelet due to excessive straining (Abdel-Gayed et al., 1985).
Two aspects of the work in explosion vessels have been explored since.
First, a relationship of turbulent burning velocity with the flame configura-
tion appears to be inevitable (Bedat and Cheng, 1995; Filatyev et al., 2005).
Hence, it is believed increasingly that a general explanation of the bending ef-
fect without considering turbulence generation mechanisms may be infeasible
(Driscoll, 2008). Further, the relationships between strain and flame prop-
agation have been investigated, primarily to determine the degree to which
transient strain determines the overall flame propagation. Recent work of
Steinberg and Driscoll (2009) has revealed through cinema stereography, for
instance, aspects of the non-linear relationship between strain and curvature,
which are found to be tightly coupled.
Further, while the work of Bradley and co-workers has had significant im-
pact in terms of the explanation of bending, a consistent general correlation of
sT (u
′) was never provided. Recently, Shy et al. (2000b) considered radiation
losses in methane-air premixed flames when mixed with diluents and found
that a correlation using Damköhler number better approximates methane-air
sT/sL data across different values of equivalence ratio (φ). The relation is as
follows
sT − sL
u′
∼ Da0.59 (1.25)
where the constants were obtained through correlations across many experi-
ments. Shy et al. (2000a) observe a similar bending effect in sT (u′) with two
distinct slopes in low and high u′ regimes. The mechanism and physical ram-
ifications of this scaling remain to be understood. In this regard, simulations
provide a window to access detailed information of the flow field and flame
19
1. Theoretical Background
process.
1.2 Simulations of turbulent burning velocity
1.2.1 Model equations for flame dynamics
Direct Numerical Simulations
With advances in computing, an important approach has emerged that allows
turbulent premixed flames to be simulated exactly under prescribed conditions
(McMurtry et al., 1986). Emerging dynamics of flame propagation can then be
inspected to ascribe physical mechanisms to the phenomena of interest. This
provides a significant amount of information accessible for analysis as opposed
to experimental techniques.
In Direct Numerical Simulation, or DNS, the governing equations for tur-
bulent gaseous reactive flow are prescribed first. Under certain initial and
boundary conditions, these equations may be solved to a high degree of nu-
merical accuracy. The equations include the Navier-Stokes equations as well as
the transport equations for the various chemical species involved in the com-
bustion reactions. Local mass-fractions Yα of an arbitrary number of active
chemical species may be computed according to the transport equation
∂ρYα
∂t
+∇ · (ρuYα) = −∇ ·J α + ω˙α (1.26)
where α = 1, . . . , N is the species index, J α is the molecular diffusion flux of
the species α. Hydrocarbon combustion may involve N > 50 such chemical
species. In order to accurately model the transport of these species, a corre-
sponding reaction mechanism must be evaluated that provides the local reac-
tion rates ω˙α. The solution of these additional chemical transport equations
poses a severe computational challenge over and above the expensive direct
numerical integration of the Navier-Stokes equations. Nevertheless, simplified
2D simulations have been feasible when employed in conjunction with reduced
chemical kinetic models (Poinsot et al., 1995). These reduced kinetic models
consider the transport and chemistry of N ≤ 25 chemical species thereby re-
20
1.2. Simulations of turbulent burning velocity
ducing computational expense. Such simulations have provided insight into
the mechanisms that underpin the observed dynamics of flame propagation.
In premixed flames, a further simplification may be adopted assuming that
the overall reaction relies on just a few key chemical species (Peters, 1988). It
is assumed often that the reaction progress variable
c =
Yα − Yα,b
Yα,u − Yα,b (1.27)
may denote sufficiently the progress of the overall combustion reaction locally
from an unburned u to burned b state provided α corresponds to a chemical
species that ensures monotonicity of c. Hence, the problem of simulating
premixed turbulent combustion may be reduced faithfully (Bray and Cant,
1991) to a single tractable equation
∂ρc
∂t
+∇ · (ρuc) = −∇ · (ρD∇c) + ω˙c (1.28)
alongside equations that describe the flow. Here, Fick’s law has been used to
model the diffusive mass flux in terms of the gradient ∇c of reaction progress
variable. Such a framework forms the basis for turbulent premixed combustion
modelling and, despite its simplicity, has provided much insight into the dy-
namics of premixed flame propagation as discussed by Bray (1990); Lipatnikov
and Chomiak (2002).
For instance, 2D flame propagation in a laminar flow with a single vortex
has been simulated by Poinsot et al. (1991a). Dynamics of vortex-flame in-
teractions were investigated for a range of vortex intensities in order to gain
insight into the regime of validity for DH1. Hence, as Poinsot (1996) notes, the
effects of excessive strain on flame quenching may be studied without having
to consider a turbulent flow. Interestingly, the area enhancement mechanism
of DH1 was shown to remain valid despite local flamelet quenching because
quenching was found to be intermittent in nature. This reinforced the conclu-
sions of Meneveau and Poinsot (1991) who, through DNS of vortex-pair pre-
mixed flame interactions, also observed that Kolmogorov scales did not have
any effect on the flame front barring short-lived high strain rates imposed by
21
1. Theoretical Background
these scales.
The flamelet regime Ka ≤ 1 corresponds to η ≥ δ0L, so it is argued that
turbulence can not disrupt the local laminar flamelet structure in this regime
Peters (1999). Experimental evidence by Bedat and Cheng (1995); Shepherd
et al. (2002) as well as DNS evidence from Poinsot et al. (1995) not only
supports this argument but emphasises further that local flamelet structure
persists even for Ka > 1. At the same time, some experimental works (Gülder
and Smallwood, 2007) question the persistence of flamelet structure in similar
or even milder conditions. Models have been developed (Cant et al., 1991;
Duclos et al., 1993) on the assumption that DH1 persists (in fact, up to Ka ≈
16) so long as quenching ‘does not inhibit the growth of the active flame surface’
(Meneveau and Poinsot, 1991). Experimental studies have since correlated sT
data within the extended DH1 framework (1.13 and 1.14), reformulating the
stretch factor I0 as a ‘burning factor’ Pb ∼ I20 instead (Bradley et al., 2011).
Further to computing exactly the solutions under prescribed hypothetical
conditions, as noted by Kerstein (2002), DNS also allows for investigation of
exact but unclosed equations of flamelet evolution. Diagnostic tests of these
equations which are arrived at using various assumptions provide further in-
sights into phenomena such as flame quenching (Poinsot et al., 1991b). More-
over, such unclosed equations have been developed to be applied for turbulent
combustion modelling in other frameworks (Veynante and Vervisch, 2002).
Three common examples (Bilger et al., 2005) are provided below.
Bray-Moss-Libby model
Proposed over several joint investigations (Bray and Moss, 1977; Bray, 1980),
the BML model presents turbulent flame propagation in terms of the Favre-
averaged form of Equation 1.28 in steady state
∇ · ρ¯u˜c˜ = −∇ · ρ¯u˜′′c′′ + ω˙, (1.29)
which is valid for Re  Da  1 so that the diffusive term may be ignored.
Favre mean quantities q˜ = ρq/ρ¯ represent density-weighted averages whereas
q′′ denotes the respective fluctuations, i.e. q = q˜ + q′′. In this framework,
22
1.2. Simulations of turbulent burning velocity
propagation of the mean flame surface c˜ emerges as a consequence of convection
(l.h.s.) and the two terms on the r.h.s. which denote respectively the turbulent
scalar flux −∇ · ρ¯u˜′′c and the chemical source-term ω˙.
Analysis of DNS datasets using the BML framework has led to insights
related to the various terms that play a role in flame propagation. One exam-
ple is provided by the manifestation of counter-gradient turbulent transport
(discussed in Chapter 5) proposed originally by Bray (1980). Evidence for
counter-gradient transport has since been found in various DNS datasets (Rut-
land and Cant, 1994; Veynante et al., 1997; Chakraborty and Cant, 2009), and
this has allowed further analysis of the mechanisms responsible for it (Bray,
1995; Veynante and Poinsot, 1997; Veynante and Vervisch, 2002).
Counter-gradient transport relates to the effect of the flame motion and
heat release on the turbulence. Still further, this relates to the concept of
flame generated turbulence initiated long ago by Karlovitz et al. (1951). Again,
DNS provides a framework for investigating the underlying mechanisms of
such effects (Louch and Bray, 2001). The BML model therefore provides a
valuable tool for investigating key mechanisms involved in turbulent flame
propagation. In the present work, this framework is employed to understand
the ramifications on DH2 by elaborating on its validity limit in Chapter 5.
Flame Surface Density
Another exact but unclosed formulation has been provided in terms of the
transport of Flame Surface Density (Pope, 1988; Cant et al., 1991; Candel and
Poinsot, 1990). This formulation provides an evolution equation of the average
local flame surface-to-volume ratio
∂Σ
∂t
+∇ ·
(
〈X˙〉sΣ
)
= 〈s˙〉sΣ, (1.30)
in terms of the mean stretch rate 〈s˙〉s, which acts as the source-term. The
quantity X˙ = u + sdn̂ is the local velocity of the surface under consideration
such that sd is the displacement speed along the normal direction n̂. In contrast
to the BML formulation, the FSD framework provides a perspective in terms of
the evolution of an isosurface. DNS datasets may be investigated for budgets
23
1. Theoretical Background
of the FSD equation in a Favre-averaged formulation (Trouvé and Poinsot,
1994). As a result, the FSD framework allows for investigation of stretch rate
influences on particular isosurfaces. Unsteady investigations have also been
conducted allowing the possibility of understanding concepts such as negative
flame speed (Gran et al., 1996; Peters et al., 1998).
Furthermore relationships between individual stretch rate components (in-
troduced in Section 1.1.2) can be investigated based on this framework. For
instance, Echekki and Chen (1996) concluded that negative curvature reduces
the rate of stretch, and its contribution is further enhanced by the nonlinear
behaviour of sd as a function of curvature. Further, the structure and the
effects on flame topology and resulting flame propagation rates were inferred
(Echekki and Chen, 1996).
While Lewis number effects have been observed in turbulent flame curva-
tures (Haworth and Poinsot, 1992), significant evidence has also been gathered
supporting purely kinematic effects of curvature (Cant, 1999). These possibil-
ities have been investigated further in 3D by Chakraborty and Cant (2004,
2009) within the same framework. Thus, this allows verification of earlier
claims made by Bradley et al. (1992) regarding flame propagations in regions
where stretch rate may dominate a predictable quenching stretch rate s˙q. These
concepts are, therefore, further investigated in the present work by linking the
concept of stretch rate source-terms with the bending effect behaviour.
Level set approach
Instead of computing locally the mean progress of reaction, the flame surface
itself may be tracked using a level-set approach (Kerstein et al., 1988). As-
suming the flame surface to be thin compared to the flow scales, the convective
transport of a passive scalar G provides
∂G
∂t
+ X˙ ·∇G = 0 (1.31)
where the local velocity X˙ = u + sdn̂ as before is the vector sum of local
flow velocity u and the local displacement velocity sdn̂ of the flame surface
along its normal n̂ = ∇G/|∇G|. By analysing the dynamics of this interface it
24
1.2. Simulations of turbulent burning velocity
is, as with FSD, possible to inspect local surface behaviour, but here changes
to the local velocity of the isosurface. Indeed, the FSD can be shown to be
proportional to the mean of the gradient modulus |∇G| (Cant and Mastorakos,
2008).
Indeed, it is interesting to note that the experimentally observed bending
effect can be reproduced by tuning long-establised models such as the Coher-
ent Flame Model (Duclos et al., 1993) and the G-equation model (Wenzel and
Peters, 2000). Despite these apparent successes, a mechanism of bending ex-
pressed using physically-intuitive terms remains elusive. DNS offers a useful
tool to investigate such mechanisms.
1.2.2 Peters: the Thin Reaction Zones regime
The level-set approach formed a key part of an investigation conducted by
Peters (1999). A scaling analysis was adopted to differentiate the Thin Reac-
tion Zones (TRZ) regime from the corrugated and distributed regimes identi-
fied earlier (see Section 1.1.4). The overall flame time scale was identified as
τF ∼ D/s2L in relation to diffusion, whereas the inner reaction time scale was
associated with the reaction zone thickness `δ ∼ 0.1δL. This divide between
the flame time scale and the reaction time scale was exploited to formulate two
closely related Karlovitz number expressions. A Karlovitz number defined as
Ka ≡ τF
τη
=
δ2L
η2
, (1.32)
the ratio of a chemical scale (τF or δL) to the Kolmogorov scale (τη or η)
determines the overall flame related stretch. In contrast, the ratio
Kaδ ≡ τδ
τη
∼ `
2
δ
η2
∼ 0.01δ
2
L
η2
, (1.33)
denotes a comparison of the Kolmogorov time scale with the reaction zone time
scale `δ. According to Peters (1999), the TRZ regime lies between Ka = 1 and
Kaδ = 1. Figure 1.5 shows the Borghi diagram of Figure 1.3 as conceptualized
in terms of the regimes proposed by Peters (1999).
The Broken Reaction Zones (BRZ) regime Ka ≥ 100 corresponds to η ≤
25
1. Theoretical Background
Figure 1.5 Regime diagram as proposed by Peters (1999) demarcating the
extent of the Thin Reaction Zones (TRZ) regime using Ka = 1 and Kaδ = 1
isolines.
δL/10 and, hence, turbulence microscales may be expected to affect the local
internal flame structure. Recent DNS evidence (Poludnenko and Oran, 2011)
shows that close to Ka ≈ 100, the effects of turbulence are limited to disruption
of the preheat zone. The reaction zone remains largely unaffected and retains
flamelet structure locally.
Investigating the scaling proposed in Damköhler’s hypotheses, Peters (1999)
further employed the notion of TRZ to understand the bending effect in terms
of the level-set approach. This was related, thereby, to the difference in effects
of kinematic restoration and scalar dissipation in the corrugated and TRZ
regimes, respectively. There is reason however to re-investigate these results
because the level-set approach limits the description of the simulated flames
to a surface with no internal structure. In another study, Wenzel and Pe-
ters (2000) arrived at a dataset that captured bending across the same regime
boundary. However, instead of using the G equation as a diagnostic tool, it
was used as the equation for simulation, thereby requiring a model input of
local burning velocity. Hence, it remains to be settled that bending indeed
occurs across this boundary of regimes. Other diagnostic frameworks outlined
above may be expected to shed light on this problem. Hence in Chapter 4, the
26
1.2. Simulations of turbulent burning velocity
bending effect is investigated in the framework of the FSD approach.
Several recent DNS studies have resolved to investigate flames at high tur-
bulence intensity. These have included DNS of laboratory-scale flames (Bell
et al., 2005, 2007; Day et al., 2012) as well as canonical configurations in-
cluding spherical flame kernels (Jenkins and Cant, 1999; Fru et al., 2011) and
statistically planar flames (Poludnenko and Oran, 2011; Aspden et al., 2011;
Savre et al., 2013; Lapointe et al., 2015). More realistic configurations are
also being addressed under similar levels of turbulence (Sankaran et al., 2007;
Hawkes et al., 2012). These efforts have paved a way for systematic studies
which attempt to isolate the parameters critical to the bending effect (Ham-
lington et al., 2011; Poludnenko, 2015; Lapointe and Blanquart, 2016). As yet,
however, a consolidated framework for addressing the bending effect is lacking.
Furthermore, several computational investigations have corroborated the
resilience of the flamelet structure (especially, the reaction zone) to intense
turbulence (Wenzel and Peters, 2000; Lipatnikov and Chomiak, 2002; Bell
et al., 2007; Poludnenko and Oran, 2011; Wang et al., 2016a). As a result,
a description of bending in the absence of flamelet quenching may be sought
(Filatyev et al., 2005). More recent DNS simulations and their comparison
with experiments at high u′ (Wang et al., 2016b,a) support the plausibility of
such an explanation. Some experiments, however, have provided contrary evi-
dence and disputed the persistence of flamelet structure as well as the validity
of DH1 in moderate to intense turbulence (Gülder and Smallwood, 2007; Yuen
and Gülder, 2013). Subsequent DNS studies have, therefore, focused on illu-
minating further the role of quenching at high u′ (Aspden et al., 2011; Savard
and Blanquart, 2015; Lapointe et al., 2015). To shed light on this problem,
Damköhler’s foundational hypotheses (Damköhler, 1940) are reframed with
regard to these investigations.
1.2.3 What do Damköhler’s hypotheses mean today?
Damköhler’s regimes, DH1 as well as DH2, correspond to regions in the Peters
regime diagram with contrasting values of length scales (`0  δL and `0  δL)
along a constant Re isoline (as discussed in §5 of Damköhler (1940)). As seen
27
1. Theoretical Background
Figure 1.6 Regimes proposed by Damköhler (1940) illustrated on the regime
diagram proposed by Peters (1999). DH1 (`0  δL) and DH2 (`0  δL) corre-
spond to the flamelet and Broken Reaction Zones (BRZ) regimes respectively.
in Figure 1.6, these belong to two extremities of the regime diagram, such that
DH1 corresponds to the flamelet regime and DH2 to the BRZ regime.
The size of these turbulent microscales is inversely related to the intensity u′
of a homogeneous isotropic turbulent flow field through the relation η ∝ Re−3/4
for a given integral scale `0 (Frisch, 2004). Hence, DH1 is expected to be strictly
valid for a turbulent flow with a specified `0 if u′ is low enough that η > δL. It
is unclear, however, how far into the TRZ regime the first hypothesis continues
to be valid. Furthermore, Damköhler proposed each of the two hypotheses to
elucidate the primary mechanism of sT enhancement in their respective regimes
of turbulence. While a primary mechanism operates, the secondary mechanism
was expected to ‘recede in the background’ contributing alongside albeit to a
lesser degree. Over the years, the range and upper u′ limit of validity for
DH1 has received significant attention. The lower u′ limit of validity for DH2,
however, remains unknown still.
Recent investigations in intense (Poludnenko and Oran, 2011; Savard and
Blanquart, 2015) and extreme u′ turbulence (Wabel et al., 2016) have high-
lighted the significance of diffusive transport processes in determining sT under
such conditions. Appropriate modifications of DH2 are being proposed so as to
28
1.2. Simulations of turbulent burning velocity
fit the sT data acquired. One may ask, if the transition between Damköhler’s
two proposed regimes lies within the TRZ regime, then what is the appropriate
marker for recording such a transition? Third, does such a transition mark
an onset of validity for DH2? How far down into the TRZ does DH2 remain
valid? The aim of the present work is to gain some insight by studying the
relationship between sT and underlying mechanisms.
1.2.4 Thesis objectives and approach
Evidently, several important questions regarding turbulent burning velocity
remain open. The present thesis investigates specifically the bending effect in
turbulent burning velocity. The following questions are addressed
1. Is bending real? Firstly, a DNS dataset is designed carefully to capture
systematically the variation of turbulent burning velocity sT with tur-
bulence intensity u′. As a consequence, the very possibility of recording
a bending effect in sT (u′) is explored. The range of u′ considered span
the entire Thin Reaction Zones regime conceptualised by Peters (1999).
The DNS method employed is discussed in the following chapter, and
the simulation dataset is designed and generated in Chapter 3.
2. Is bending kinematic or chemical? Secondly, upon registering a bending
effect, the mechanism governing this variation of sT (u′) is investigated.
The validity of Damöhler’s first hypothesis is verified within the dataset,
providing an insight into whether the bending phenomenon retains the
thermo-chemical structure of the flame. The plausibility of a kinematic
mechanism of bending is explored within the framework of the Flame
Surface Density approach in Chapter 4. Due consideration is given to
the possible role of flame surface curvature.
3. Is Damköhler’s second hypothesis valid? Finally, the validity of the
neglected hypothesis of Damköhler (1940) is investigated in the DNS
dataset. The objective is to shed further light on the mechanisms un-
derpinning the recorded bending effect. The associated role of flame-
induced self-acceleration in limiting the validity of DH2 is studied within
29
1. Theoretical Background
the framework of the Bray-Moss-Libby model in Chapter 5.
The thesis concludes with an ongoing work that aims to validate the pro-
posed explanations with experiment. The Leeds explosion vessel is considered
for this comparison in Chapter 6. By probing the possible mechanisms of
the bending effect, this thesis therefore aims to provide an insight into the
long-standing problem of turbulent flame propagation.
30
Chapter 2
The Direct Numerical Simulation
method
As discussed in the previous chapter, Direct Numerical Simulation (DNS) pro-
vides the ability to capture and analyse flame dynamics in detail by resolving
the various scales of gas motion in space and time. Until recently, DNS has
yielded significant insights into turbulent flame propagation (Chakraborty and
Cant, 2004; Poludnenko and Oran, 2011). Yet its application to the bending
effect remains to be addressed (Lipatnikov and Chomiak, 2002; Driscoll, 2008).
Hence, in the present work, the bending effect is investigated using the method-
ological framework of DNS as follows.
The basic equations governing combustion chemistry and species transport
are written in a 3D form. In order to close the basic equations, a set of
well-posed boundary conditions (BCs) is supplied as appropriate for the flame
configuration. Further, an initial condition (IC) is applied to the variables of
the basic equations. This is done in such a manner that the computational
effort required to reach a statistically significant solution is minimised. Upon
the imposition of an IC and appropriate BCs, the system of equations are
converted to a high-order discretised form in space and time. This enables
a numerical approach to obtaining the solution. The numerical algorithm is
implemented in a FORTRAN77 computer code called SENGA2 developed by Cant
(2013). Time-dependent solutions of the desired problem are obtained by
running the code (also called solver) on supercomputing clusters.
31
2. The Direct Numerical Simulation method
Each of these aspects specific to the presently used DNS framework are
elaborated below. The presentation of the basic equations follows the style of
Buckmaster and Ludford (1982). Following the presentation of the framework,
a set of benchmark solutions conducted to verify the computational method is
discussed. These benchmarks are presented in a manner leading progressively
up to the primary DNS investigations of the current work.
2.1 Basic equations
Consider first the temporal evolution of a homogeneous gaseous mixture un-
dergoing chemical reaction. At any time t during the reaction process, the
mixture is described by its composition vector Y (t) = {Yα(t) : α = 1, . . . , N}T
which comprises mass-fractions of each component (chemical species) denoted
by the index α. At any such time,
N∑
α=1
Yα(t) = 1, (2.1)
providing a constraint to the temporal evolution of the N mixture compo-
nents. Hence, any N − 1 components of the gaseous mixture determine the
composition vector Y (t) completely. The density ρ(t) of the gaseous mixture
is computed as
ρ(t) =
N∑
α=1
ρα(t)Yα(t), (2.2)
the mass-weighted-average of component densities ρα(t). The temperature
T (t) may be assumed to be the same for all components and is obtained using
the caloric equation of state. Together, the local fields Y (t), ρ(t), and T (t) form
a basis for the time-evolving thermo-chemical state-space of a homogeneous
reacting gaseous fuel-air mixture.
Kinetic models pertinent to hydrocarbon combustion chemistry are pre-
sented in the sections below. By providing a model for the reaction rate
ω˙, these descriptions address the time-evolution of the thermo-chemical state
space.
32
2.1. Basic equations
Combustion Chemistry
The overall chemical reaction for the combustion of a stoichiometric mixture,
i.e. with equivalence ratio φ = 1, of methane and air may be expressed as
CH4 + 2(O2 + 3.76N2)
−∆H◦c−−−−→ 2H2O + CO2 + 2 (3.76N2), (2.3)
where it is assumed that air is composed of O2 and N2 in a molar proportion
of 1:3.76, and that N2 remains inert throughout the reaction.
The overall reaction given by Equation 2.3 proceeds exothermically, re-
leasing ∆H◦c as the heat of combustion for methane. The rate of the overall
reaction ω˙ denotes the rate of consumption of reactants or, equivalently, the
rate of formation of products. However, the overall reaction is composed of
many elementary reactions where intermediate chemical species are formed and
consumed. The ‘full’ chemical mechanism has been elaborated for methane-air
combustion, for example in the GRI-MECH 3.0 (Smith et al., 2017) and it con-
sists of M > 300 elementary reactions and N > 50 species. A reduced kinetic
model (with M < 50) is computationally tractable while at the same time
it can describe some of the primary underlying physics of flame propagation.
Two such reduced chemical kinetic models are considered in the present thesis.
One-step Kinetics
The most basic description of the overall combustion reaction is
R
−∆H◦c−−−−→ P, (2.4)
which models the chemical mechanism of combustion as a one-step irreversible
transformation from unburned reactants R to burned products P . The rate ω˙
of a one-step reaction is expressed simply in the Arrenhius form
ω˙α = AρYαe
−Ea/(R◦T ), (2.5)
based on the overall activation energy Ea and overall collision factor A. The
index α may correspond to R or P , which can be assigned the thermo-chemical
33
2. The Direct Numerical Simulation method
properties (such as enthalpy of formation ∆H◦f ) of the reactants and products
of the overall reaction given by Equation 2.3. This allows the overall heat
release ∆H◦c and its rate to be accounted for by this kinetic model.
Upon the specification of an appropriate Lewis number Le and activation
energy Ea, the flamelet structure and propagation speed may be also com-
puted with reasonable accuracy (this is discussed in Section 2.4). By account-
ing for these salient aspects of combustion processes (Le, Ea, and ∆H◦c ), this
seemingly simple kinetic model has provided satisfactory explanations of many
combustion phenomena related to laminar as well as turbulent flame propaga-
tion. These are presented in the review by Lipatnikov and Chomiak (2002). As
a further bonus, computational expense is minimal because the model requires
only one species to be transported, with the other determined as YP = 1−YR.
Due to these merits, the one-step kinetic model is chosen in the paramet-
ric simulations of the bending effect in Chapter 3. This has allowed further
investigation of the Damköhler’s hypotheses in the recorded bending effect
(Chapter 4 and Chapter 5). Nevertheless, the role of intermediate chemical
species in the bending effect can be elaborated only by considering multi-step
kinetics.
Multi-step Reduced Kinetics
An improvement upon the one-step kinetic model is to use multi-step kinetic
models. These are obtained by reducing systematically the ‘full’ chemical
mechanism of combustion using steady-state assumptions for elementary re-
actions and fast chemistry approximations for intermediate species. In this
manner it may sometimes, but not always, provide a more realistic descrip-
tion.
A 25-step kinetic model due to Smooke and Giovangigli (1991) is used
in Chapter 6. The consideration of various reaction pathways helps model
the role of intermediate species transport in the bending effect. Furthermore,
two other kinetic models were investigated: the 4-step kinetic model due to
Peters and Williams (1987), and the 68-step model due to Sandia National
Laboratories (Kee et al., 1986). These models were not adopted because the
34
2.1. Basic equations
former was found to be computationally stiff and the latter computationally
prohibitive. The single-step model is the one used primarily in the present
investigation. Initial comparisons using such the 25-step kinetic model have
been considered in Chapter 4, but the description of multi-step kinetics is
deferred until Chapter 6 where it is used for experimental validation.
In a volume that exhibits flame propagation, the state {ρ,Y , T} changes
spatially as well, so that the local state is affected by species transport. This
spatio-temporal evolution is addressed in the following sections. Adiabatic
systems are considered so that volumetric heat transfer may be neglected.
Continuum Description of Species Transport
In turbulent flame propagation, spatial variations of the composition vector
Y (x, t) may be described using its transport equation. This is determined by
the equation of mass-balance for each species α
∂ρYα
∂t
+∇ · (ρvαYα) = ω˙α, (2.6)
where vα is the local species velocity and ω˙α is its rate of formation source-term.
The variable ω˙α is modelled using one of the chemical kinetic models discussed
above. A further relation is needed, however, for the species velocity vα. The
precise molecular contribution in vα can be made explicit by defining the local
bulk velocity as the mass-averaged velocity of the mixture components
v ≡
N∑
α
Yαvα, (2.7)
in similar fashion to Equation 2.2 for the local bulk gas density. The molecular
diffusion velocity Vα = vα − v is introduced subsequently so that the species
mass-transport equation may be re-written as
∂ρYα
∂t
+∇ · (ρvYα) = ω˙α −∇ · (ρVαYα), (2.8)
where the unclosed term (the second term on the r.h.s.) is explicitly a molec-
ular diffusion term. In non-reacting mixtures, the diffusive flux ρVαYα is ex-
35
2. The Direct Numerical Simulation method
pressed as a consequence of the local concentration gradient of the component
α. This formulation (Bird et al., 1960) is due to Fick, and is applied to reacting
mixtures as is for the sake of simplicity. Hence, in the present work
Vα,jYα = −Dα∂Yα
∂xj
, (2.9)
neglecting the Soret effect (molecular diffusion due to temperature gradient).
The species diffusion coefficient
Dα =
λ
ρCpLeα
(2.10)
is related to its thermal conductivity λ and heat capacity Cp via the Lewis
number Leα. The ratio λ/Cp expressed as a function of temperature
λ
Cp
= Aλ
(
T
T ◦
)r
, (2.11)
holds for the local mixture with constants Aλ and r specified a priori. Due
to its relative simplicity, this approach is adopted in the present work from
Chapter 3 through Chapter 5 assuming a unity mixture Lewis number. In
Chapter 6, this framework is changed. The Fick’s law of Equation 2.9 is
modified so that a dependence on the local composition vector is reflected.
Similarly, the values of Dα and λ are obtained based on the local properties of
the mixture.
Mass, Momentum and Energy Transport
The sum of mass balances for each species in Equation 2.8 yields the continuity
equation for the bulk fluid written as
∂ρ
∂t
+∇ · (ρv) = 0, (2.12)
under the compatibility conditions
∑N
α=1 ω˙α = 0 and
∑N
α=1 VαYα = 0. Simi-
larly, the local momentum balance equation for the mixture can be obtained by
adding individual momentum-balance equations for each component species.
36
2.1. Basic equations
This is given by the Navier-Stokes equations in compressible form
∂ρv
∂t
+∇ · (ρv ⊗ v) = −∇p+∇ · τ , (2.13)
upon neglecting gravitational and other body forces. The variable p introduced
here is the local hydrodynamic pressure, whereas τ is the viscous stress-tensor
of the bulk fluid. Both variables must be modelled in order to close the system.
The local pressure p emerges, as a consequence of Dalton’s Law for single-phase
systems (Buckmaster and Ludford, 1982), as the sum
p =
N∑
α=1
pα, (2.14)
of partial pressures given by the ideal gas law for a perfect gas
pα = ρR
◦TYα/Wα, (2.15)
for each individual component α with molar mass Wα. Assuming the mixture
to be a Newtonian fluid, the viscous stress tensor τ is computed as
τji = µ
(
∂ui
∂xj
+
∂uj
∂xi
− 2
3
∂uk
∂xk
δji
)
, (2.16)
where the viscosity
µ =
λ
Cp
Pr, (2.17)
is expressed as a linear function of λ/Cp using the Prandtl number Pr, which
is assumed constant.
Lastly, the equation of energy-balance is derived for the chemically reacting
flow system. The local internal energy for a compressible flow is defined as
e ≡ 1
2
ujuj − p
ρ
+
N∑
α=1
Yαhα, (2.18)
where hα is the local enthalpy of species α given by the integral
hα =
∫ T
T0
CP,αdT + h
0
α (2.19)
37
2. The Direct Numerical Simulation method
for a specified value of formation enthalpy h◦α. The mass-specific heat capacity
Cp,α is approximated using a polynomial in temperature tabulated a priori for
each species (McBride and Gordon, 1967). The transport equation for e is also
obtained in a similar manner as for the mass and momentum balance, and can
be written as
∂ρe
∂t
+∇ · (ρve) = −∇ · pv +∇q + τ : ∇v (2.20)
assuming the work done by molecular interaction forces is negligible (Buck-
master and Ludford, 1982). The heat flux vector q is introduced here as new
unclosed variable. Neglecting radiative transfer and the Dufour effect (heat
flux due to concentration gradient), the components of q are expressed in
terms of the variables ρ, T , Vα, Yα and hα as follows
qj = −λ ∂T
∂xj
+
N∑
α=1
ρVα,jYαhα, (2.21)
using Fourier’s conduction law (first term on r.h.s.) and the sum of diffusion
enthalpies (second term). Hence, a closed system of equations is obtained.
Summary of equations in closed form
Conservative variables ρ, ρui, ρe, and ρYα are described as
∂ρ
∂t
+
∂
∂xi
(ρui) = 0, (2.22)
∂ρui
∂t
+
∂
∂xj
(ρujui) = − ∂
∂xi
(p) +
∂
∂xj
(τji), (2.23)
∂ρe
∂t
+
∂
∂xj
(ρuje) = − ∂
∂xj
(puj) +
∂
∂xj
(qj) +
∂
∂xj
(τjiui), and (2.24)
∂ρYα
∂t
+
∂
∂xj
(ρujYα) = ω˙α +
∂
∂xj
(ρVα,jYα), (2.25)
38
2.2. Boundary Conditions
where e ≡ 1
2
ujuj − p/ρ+
N∑
α=1
Yαhα and the variables p, τji and qj are modelled
as
p = ρR◦T
N∑
α=1
Yα/Wα (2.26)
τji = µ
(
∂ui
∂xj
+
∂uj
∂xi
− 2
3
∂uk
∂xk
δji
)
, and (2.27)
qj = −λ ∂T
∂xj
+
N∑
α=1
ρVα,jYαhα (2.28)
and, finally, the variables Vα,j and ω˙α are expressed as
Vα,jYα = −Dα∂Yα
∂xj
, and (2.29)
ω˙α = AρYαe
−Ea/(R◦T ), (2.30)
where λ, Cp, Dα and µ are functions of temperature. All other constants are
known a priori from the mixture components and chemical kinetic models.
2.2 Boundary Conditions
Equations of transport – Equation 2.22 through Equation 2.25 – follow the
general conservative form
∂U
∂t
+∇ ·C = D+ S, (2.31)
in which U = {ρ, ρui, ρe, ρYα}T are the conserved variables, C and D are
the convective and diffusive fluxes, and S is the vector of source-terms. These
equations do not require any further information in an unbounded domain with
infinite extents in all dimensions; however, such a problem is not amenable to
numerical solution in a finite time. Hence, unbounded domains of infinite
extent must be truncated to lend computationally-tractable domains. This is
done by exploiting various symmetries in the solution expected. For instance,
the periodic BC exploits translational symmetry in the specified direction.
Such BCs are used routinely to render computationally feasible 3D simulations
39
2. The Direct Numerical Simulation method
of homogeneous isotropic turbulence (Pope, 2001).
Basic Concepts and Flame Configuration
In directional flows such as vortex advection and flame propagation, the trans-
lational symmetries are broken in at least one dimension. Hence, periodic BCs
cannot be applied for the simulation of such problems. In such circumstances,
because the domain must nevertheless be truncated, a different kind of BC
must be specified at the artificial boundaries.
The most basic configuration for modelling a propagating flame is a 1D
domain with inflow and outflow boundaries on either end (discussed in Sec-
tion 2.4). In this configuration, laminar inflow balances the propagating flamelet
between the boundaries. This problem can be extended to two dimensions by
considering a planar flame where the added span-wise dimension is assigned pe-
riodic BCs. Extending the 2D configuration to three dimensions, a 3D inflow-
outflow configuration is obtained. Such a domain is necessary for modelling
realistically the turbulent energy spectrum. A statistically planar flame may
propagate in the longitudinal or x direction (inflow-outflow) and the remaining
directions are assigned periodic BCs. This configuration is used for the flame
simulations conducted in the present work (see Section 3.1.1).
Primitive Variables
At the inflow and outflow boundaries, the solution is specified in terms of
primitive variables U = {ρ, ui, p, Yα}T rather than conservative variables U.
This is because the imposed values are generally known in primitive rather
than conservative form. The primitive variables U follow a transport equation
that may be obtained by transforming the conservative equations 2.31 using
the Jacobian matrix P = ∂U/∂U
Conservative Form
P−1−−→←−−
P
Primitive Form, (2.32)
to arrive at an equivalent system of equations
∂U
∂t
+ A · (∇U) = D + S, (2.33)
40
2.2. Boundary Conditions
written using primitive variables, where A = P−1∂C/∂U is the corresponding
flux-transformation matrix. Since the two forms are equivalent, closure of one
system of equations is informed by the other.
The primitive system consists of hyperbolic (unsteady diffusive terms) as
well as non-linear terms (reaction rate source-terms). The corresponding equa-
tions require an appropriate number of constraints for closure (Poinsot and
Lele, 1992). The imposition of these constraints at the boundaries, however,
introduces unwanted waves in the domain of solution that may contaminate
the numerical solution over time. A challenge in flame simulations, therefore,
is to suppress all such unwanted effects of the artificial boundary (Sutherland
and Kennedy, 2003).
Navier-Stokes Characteristic Boundary Conditions
One approach to reduce spurious boundary-effects is to use the characteristic
wave form of the primitive equations. The idea is to be able to specify the
amplitudes of the known waves while diminishing the amplitudes of the un-
known or unwanted waves. In an inflow-outflow configuration, waves travel in
the longitudinal x direction, so that the primitive form of equations can be
written in the boundary-normal form as
∂U
∂t
+ A(n) ·
(∇(n)U)+ A(t) · (∇(t)U) = S, (2.34)
upon neglecting the diffusive terms D in Equation 2.33. Subscripts n and t
indicate normal and transverse components to the boundaries. The character-
istic decomposition of this equation is written as
∂U
∂t
+

1
ρa
(
L(x)5 − L(x)1
)
L(x)3
L(x)4
L(x)2 +
1
a2
(
L(x)5 + L(x)1
)
L(x)5 + L(x)1
L(x)5+α

= S (2.35)
41
2. The Direct Numerical Simulation method
in the locally one-dimensional inviscid (LODI) formulation. Each of the terms
L(x)i denote here the characteristic wave-amplitude variations in the axial di-
rection; transverse terms as well as diffusive terms are neglected in the LODI
formulation. The terms L(x)i
L(n) =

L(x)1
L(x)2
L(x)3
L(x)4
L(x)5
L(x)5+α

=

λ
(x)
1 ·
1
2
(
∂p
∂x
− ρa∂u
∂x
)
λ
(x)
2 ·
1
2
(
∂ρ
∂x
− 1
a2
∂p
∂x
)
λ
(x)
3
2
·
(
∂v
∂x
)
λ
(x)
4
2
·
(
∂w
∂x
)
λ
(x)
5 ·
1
2
(
∂p
∂x
+ ρa
∂u
∂x
)
λ
(x)
5+α
2
·
(
∂Yα
∂x
)

(2.36)
are characteristic wave-amplitudes with characteristic velocities λ(x)k
λ
(x)
1 = u− a, λ(x)2 = λ(x)3 = λ(x)4 = λ(x)5+α = u, λ(x)5 = u+ a. (2.37)
in the boundary-normal direction n, which refers to the longitudinal/x direc-
tion here. Given the sound speed a, the incoming and outgoing waves at a
given boundary may be distinguished based on the sign of the characteristic
velocities λ(x)i . In subsonic flow u < a, all waves except L(x)5 are incoming at
the left inflow boundary, whereas all waves except L(x)1 are outgoing at the
right outflow boundary. In order to specify the BCs for Equation 2.35, only
the incoming characteristic wave-amplitudes need to be specified, while the
outgoing waves can be computed using the solution in the domain. These
characteristic wave-amplitudes arise in the characteristic form
∂U
∂t
+ L(n) = S, (2.38)
42
2.2. Boundary Conditions
of Equation 2.33 obtained by transforming the inviscid primitive form as
Primitive Form
(Sn)−1−−−−→←−−−−
Sn
Characteristic Form (2.39)
using the matrix of right eigenvectors Sn of the flux transformation matrix A.
The characteristic variables U = Sn−1U have source-terms S (Sutherland and
Kennedy, 2003). Knowing the characteristic variables and the source-terms,
the wave-amplitudes Li can be estimated. These may then be specified at the
appropriate domain boundaries. Once the characteristic wave-amplitudes are
computed, these are supplied in the primitive equation which is transformed
back to the conservative form and solved numerically at the boundary.
Subsonic constant-density inflow
At the inlet, the characteristic wave-amplitudes L1, L2, L3, L4 and L5+α must
be specified whereas the outgoing wave-amplitude L5 is computed from the
characteristic equation. In a constant-density inflow BC, the density ρ to
be prescribed is known at the inlet as well as the species mass-fractions Yα
and the velocity components ui. These are used in the computation of the
characteristic wave-amplitudes Li along with the Taylor’s hypothesis which
allows calculation of the velocity components and the velocity gradients.
With the longitudinal component u1 set to a specified mean velocity 〈u1〉
in flame simulations, the laminar flame speed sL, the unsteady velocity field
in the other two directions is obtained by scanning a frozen field of homoge-
neous isotropic turbulence. Relating the unsteady local velocity field to spatial
traversal in a frozen velocity field is valid strictly when
u′  〈u1〉 (2.40)
in which case Taylor’s hypothesis provides the unsteady velocity
∂ui
∂t
= 〈u1〉∂ui
∂x
(2.41)
in terms of the spatial velocity gradient required for closing the inflow boundary
43
2. The Direct Numerical Simulation method
condition. The temperature at the inlet is calculated using the LODI form of
the energy equation.
Subsonic non-reflecting outflow
At the outlet, all waves are outgoing except L1. The value of this amplitude
can be obtained using a steady-state approximation of the corresponding char-
acteristic equation. The value L1 = S1 is then obtained upon setting all other
wave-amplitudes to 0. However, it is observed that with such a prescription,
the pressure in the domain tends to drift (Poinsot and Lele, 1992). Hence,
in order to relax the boundary condition towards the steady-state value, a
pressure-correction term is added as
L(x)1 = S1 +
σ
2L
a
(
1−M2) (p− p∞) (2.42)
which accounts for the initial differences in pressure that may not satisfy the
steady state equation of the characteristic. Rudy and Strikwerda (1980) pro-
vided an approximate constant value for σ = 0.287. The ambient pressure is
set to p∞ = 1 bar , and the Mach number is specified asM = 5× 10−3 (Cant,
2013).
2.3 Numerical method
Upon prescribing the boundary conditions, the closed system of boundary-
value-problem equations can be solved numerically as a time-varying ODE
∂U
∂t
= R(U, t) (2.43)
such that the r.h.s. is composed of source-terms and first and second order
spatial derivatives fˆ ′ ≡ ∂/∂xi and fˆ ′′ ≡ ∂/∂xi∂xj of primitive variables U.
The exact re-arranged form of the basic equations consists of a skew-symmetric
decomposition of the terms so as to reduce numerical error. Further details
are discussed in the SENGA2 User Manual (Cant, 2013).
44
2.3. Numerical method
Spatial Discretisation
The derivatives f ′ and f ′′ are estimated at coordinate points obtained by dis-
cretising the domain uniformly along each direction. At any interior mesh
point, the derivatives are computed to 10th order accuracy using a centred
finite difference scheme. The first order approximation
f ′ =
5∑
m=1
a′m
fi+m − fi−m
2mhi
, (2.44)
as well as the second approximation
f ′′i =
5∑
s=1
a′′m
fi+m − 2fi + fi−m
(mhi)
2 , (2.45)
of the derivatives require 5 points on each side of the chosen point. Similar
expressions are formed for mixed second order derivatives ∂/∂xi∂xj, i 6= j.
The respective coefficients a′m and a′′m are obtained by solving the a system
of equations formed by the Taylor series expansions at each of the 11 points.
Points on a periodic boundary are treated the same as general interior points.
For points on a non-periodic boundary, if the chosen point is less than 5 points
away from a non-periodic boundary (inlet and outlet), the order of the method
is sacrificed. At points less than 3 points away from the boundary, the 4th order
centred scheme (used 3 points away) is skewed towards the interior so as to
maintain the order of the method. The coefficients a′m and a′′m are altered to
accommodate for these changes. Values for these coefficients are provided in
the SENGA2 User Manual (Cant, 2013).
Time Marching
Upon computing the spatial differences, the r.h.s. of Equation 2.43 is known
completely, and the ODE
dU
dt
= R(U, t), (2.46)
can be solved using a 4th-order explicit 5 stage explicit Runge-Kutta scheme
(classified as RK4(3)5[2R+]C) due to Kennedy and Carpenter (2003). In this
45
2. The Direct Numerical Simulation method
method, the conservative solution U is advanced at each time step tn by a
discrete interval ∆t over successive stages. At any intermediate stage i, the
solution is calculated as
U(i) = U(n) + ∆t(n)
i−1∑
j=1
a
(i)
j R
(j) (2.47)
and in the final stage, the time-advanced solution is obtained using all the
solutions at previous stages as
U(n+1) = U(n) + ∆t(n)
5∑
j=1
bjU
(j) (2.48)
where the stage coefficients a(i)j and bj are pre-computed to match the order
of the method. The Butcher array is used to obtain the coefficients in a
straightforward manner as described in the SENGA2 User Manual (Cant, 2013).
Hence, in principle, this can be performed using two registers – one for the
solution vector U(i) and another for the r.h.s. vector R(j) required at each
stage.
Adaptive Time-step Control
The error incurred over the numerical time marching can be estimated using
the truncation error with respect to a lower-order scheme. As the solution is
advanced using a 4th-order scheme as shown above, a 3rd order scheme is used
in parallel. This approach is referred to as ‘embedded method.’ The difference
between the two solutions in an embedded method provides an estimate of the
error. This error may be used to compute adaptively the timestep for the next
iteration of the algorithm. A PID controller due to Kennedy and Carpenter
(2003) takes the form
∆t(n+1) = TOL ·
(
ˆ
rˆn+1
)α(
ˆ
rˆn
)β(
ˆ
rˆn−1
)γ
∆t(n) (2.49)
where α = 0.49/p, β = 0.34/p and γ = 0.1/p given p = 3 is the order of
the embedded method. The error-bound ̂ = 10−3 and the tolerance variable
46
2.4. Benchmark solutions
TOL = 0.9 allows for a safety-factor in the timestep estimation. The use of
such controllers in the time-marching schemes for PDEs has a rich background
in dynamical control theory. This PID controller is used as is for Chapter 3
through Chapter 5. An improved form of the controller (Söderlind, 2002) is
developed in the final chapter.
Computational Algorithm
Numerical solutions of the governing equations may be computed using the
schemes described above. In order to exploit parallel computing resources, a
parallel decomposition of the domain is performed. The numerical solution is
thus computed simultaneously over various sub-domains mapped to individ-
ual cores on a supercomputing cluster. In order to compute gradients using
finite difference approximations described above, points close to the edges of
the sub-domain require information from neighbouring sub-domains. For this
purpose, a HALO array is constructed where the necessary information from
the neighbouring sub-domain is stored during the computation procedure. For
optimised fast I/O operations and optimised overall performance of the CFD
solver, the HDF5 library is used. This allows collective MPI I/O operations using
an API. The DARWIN (Intel Sandy Bridge E5-2670) cluster and the massively
parallel ARCHER (Cray XC30) supercomputer were used to conduct all simula-
tions in the present work. Post-processing of data is conducted using various
tools including Python, Jupyter and an in-house FORTRAN95 code called SATN.
2.4 Benchmark solutions
Standard test cases have been performed using SENGA2 to validate the frame-
work. These include, in order of increasing flow complexity, a) 1D laminar
propagating flame solutions, b) 2D vortex advection study, and c) homoge-
neous isotropic decaying turbulence. These serve to verify respectively: a)
combustion chemistry calculations, b) the inflow and outflow boundary condi-
tions, c) computation of 3D Navier-Stokes equations with turbulent spectra.
Hence, all necessary components of a 3D turbulent flame propagation problem
47
2. The Direct Numerical Simulation method
are benchmarked in this manner.
1D Laminar Flame
A benchmark for flame chemistry models is the 1D laminar flame propagation
solution. The 1D domain with inflow and outflow boundary conditions is set
up. A steady inlet velocity is set approximately equal to the unstretched lam-
inar flame speed sL expected a priori under the specified ambient conditions
p◦ = 1 bar and T ◦ = 300K. The initial thermochemical field is specified as
YR(x; t = 0) =
1
2
[
1 + erf
(
x− x0
δ
)]
(2.50)
for the reacting species R. The density and temperature are set in a similar
manner using error function profiles. Using a single-step description of chem-
istry, a converged solution for 1D laminar flame propagation is obtained for
these conditions after approximately 1600000 time steps. Figure 2.1 shows the
key profiles obtained in terms of flame coordinates x (Figure 2.1a) as well as in
the progress variable space c (Figure 2.1b). A thin reaction zone is evident (red
line in Figure 2.1a). A peak reaction rate of ω˙max ≈ 2000 s−1 is computed for
the stoichiometric mixture. In terms of the progress variable, this peak occurs
(a) Axial Coordinates (b) Progress variable space
Figure 2.1 Laminar propagating 1D flame solution modelled using single-step
kinetics and unity Lewis number. A spatially-thin reaction zone (left) is sim-
ulated. The peak reaction rate occurs at c = 0.8 (right).
48
2.4. Benchmark solutions
Figure 2.2 Grid resolution study to determine the bounds for flame chemistry
resolution. Convergence fails for grids with resolution ∆x & δL/5 where δL =
0.36× 10−3m is obtained a priori on the finest grid.
at c = 0.8; this is an important aspect that will be used later on in the sub-
sequent analysis. A laminar flame speed sL ≈ 0.39ms−1 and flame thickness
δL ≈ 0.36mm are computed. All other relevant thermo-chemical parameter
values will be discussed in Chapter 3. The discussion here is intended solely
to verify the framework.
In order to verify further the dependence of the solution on the mesh prop-
erties, a mesh convergence study was conducted. This helps also find the
resolution constraint imposed by the chemical flame structure of a 1D laminar
flame. The laminar flame thickness δL was computed first using a fine mesh
with nearly 40 points across this flame width δL. This corresponds to the
right most point δfineL /∆x ≈ 40 in Figure 2.2. Subsequently, the mesh was
coarsened by increasing ∆x so that δfineL /∆x decreases (from right to left in
Figure 2.2). In the coarsest mesh, (left most point in Figure 2.2), the flame was
resolved using nearly 5 points. Inclusive of the fine and coarse extremes, all
intermediate mesh configurations provided the same converged laminar flame
speed and similar laminar flame structure as the fine-mesh result. The added
benefit of using coarser meshes is that the number of steps Nc (ordinate in
Figure 2.2) required to attain convergence is reduced. Meshes with coarser
49
2. The Direct Numerical Simulation method
resolution failed to provide a converged solution (marked using red text).
2D Vortex Advection
A benchmark case for boundary conditions is performed in a 2D domain with
the left and right boundaries specified with inflow-outflow NSCBCs along the
longitudinal direction and periodic BCs along the span-wise directions. The
initial condition for a 2D vortex is given by the velocity field
u
v
 =
u0
0
+ 1
ρ

∂ψ
∂y
−
∂ψ
∂x
 , (2.51)
and is determined using the potential flow coordinate
ψ = C exp
−x2 + y2
2R2c
 , (2.52)
where C determines the vortex strength and Rc the core radius. The pressure
p for this flow field is initialised as
p− p∞ = ρC
2
R2c
exp
(
−x
2 + y2
2R2c
)
(2.53)
where p∞ = 1 bar is the ambient pressure. The mean axial inflow u0 = 5ms−1
is subsonic unlike in the study by Poinsot and Lele (1992).
Figure 2.3 shows various snapshots of vorticity during the time-evolution of
the snapshot. As the vortex approaches the outflow boundary, it is distorted
and squeezed out of the domain through its lower half (x velocities are higher
in this region of the vortex). Acoustic waves (small vertical lines at t = 1ms)
are traced travelling away from the inflow boundary in one of the snapshots.
A region of weak vorticity is generated near the inflow from the interaction of
this wave with the boundary. These vorticity regions may contaminate the so-
lution and are the subject of ongoing investigation. In particular, an improved
formulation is sought for the outflow boundary condition prescribed in the
50
2.4. Benchmark solutions
(a) t = 0 (b) t = 0.6ms
(c) t = 1.0ms (d) t = 1.2ms
Figure 2.3 Simulation snapshots of a counter-clockwise rotating vortex ad-
vected through a 2D inflow-outflow domain with inlet velocity u0 = 5ms−1.
current work. At present, however, this serves to verify the general framework
of NSCBCs for computing solutions in an inflow-outflow configuration.
3D Turbulent Non-reacting Flows
Finally, a benchmark is performed to verify the 3D computation of turbulent
energy spectra. In turbulent flows, the initial turbulent velocity field v(x, t =
0) for a cubic domain of size L is computed by specifying an energy spectrum
due to Batchelor and Townsend (1948) as follows
E(κ) =
c0L
κ0
(
κ
κ0
)4
exp
[
−2
(
κ
κ0
)2]
, (2.54)
51
2. The Direct Numerical Simulation method
in the wavenumber κ space such that κ0 denotes the wavenumber with the
maximum turbulent kinetic energy. The constant c0 is corresponds to the
initial mean turbulent kinetic energy k0 ≡ 1/2〈v ·v〉 at t = 0 using the relation
c0 =
32
3
√
2
pi
k0. (2.55)
For a specified initial energy spectrum E(κ), one instantiation v(κ) of
the possible velocity fields is obtained in wavenumber space using the method
proposed by Orzsag discussed in Cant (2013). A Fourier transform allows
the initial velocity field v(κ, t = 0) to be transformed into the velocity field
v(x, t = 0) in Cartesian space. Hence, a homogeneous isotropic turbulence
(HIT) field is obtained as the initial condition. With the periodic BCs specified
in all three directions, the turbulent decay of this HIT field can be simulated.
The various regions of the energy spectrum are illustrated in Figure 2.4a for
a time-evolved solution of an HIT field. The central region (red) denotes the
energetic modes κE whereas the inertial subrange (orange/yellow) denotes the
modes where dissipation is independent of the wavenumber. The spectrum
extends up to κmax where the dissipation range terminates.
The instantaneous turbulent kinetic energy k(t) is expected to follow the
(a) Energy spectrum (b) Temporal energy decay
Figure 2.4 Turbulent kinetic energy computed in a periodic 128 × 128 × 128
HIT simulation. Panels show a) snapshot of energy spectrum in Fourier space
and b) temporal decay of kinetic energy.
52
2.4. Benchmark solutions
temporal relationship
k = k0
(
t
t0
)−n
, (2.56)
with the exponent n constant for a developed flow signifying a constant dissi-
pation rate  (Pope, 2001). Figure 2.4b shows the temporal evolution of mean
kinetic energy k(t) obtained using DNS with grid size 128× 128× 128. Initial
kinetic energy k0 represents a flow with u′ = 10sL. The value of the decay
exponent n ≈ 1.85 achieved after a time t ≈ t0 is found to be close to that
expected in the final period of decay (Batchelor and Townsend, 1948; Pope,
2001). The exponent stays constant for the remaining time at least until the
kinetic energy attains a small fraction k ≈ 0.05k0 of the initial value. This
serves to verify the 3D computations of the Navier-Stokes solutions in the
present framework.
Summary: computational framework
The computational framework used in the subsequent investigations is laid
out. The Navier-Stokes equations are framed in a 3D compressible form along
with thermo-chemical transport (Cant, 2013). A 10th-order centred finite dif-
ference operator evaluates spatial derivatives and reduces to 4th-order at the
inlet and outlet faces. In the stream-wise direction, inflow-outflow bound-
aries is specified using the Navier-Stokes characteristic boundary conditions
(Sutherland and Kennedy, 2003). At the inlet, a reflecting constant-density
condition with specified velocity components is imposed, whereas at the out-
let, a non-reflecting outflow condition is imposed. The transverse directions
are specified with periodic boundary conditions. A 4th-order explicit Runge-
Kutta method is used for time-marching along with adaptive time-stepping
implemented using an error controller (Kennedy and Carpenter, 2003).
53

Chapter 3
Parametric Simulations of the
Bending Effect
The goal of the present chapter is to use the computational framework de-
scribed in the previous chapter to simulate turbulent premixed flame propa-
gation. Selection of relevant simulation parameters has been informed by the
conditions of bending recalled below. Recorded images of the turbulent flames,
the pre-heat and reaction zone structures, as well as the associated turbulent
burning velocities are discussed.
3.1 Design conditions for bending
Generalising the expression discussed in Chapter 2, Driscoll (2008) has pro-
vided a general dependence of sT on physico-chemical parameters as
sT = f(u
′, `0, U¯in; sL, δL,Le;L) (3.1)
where u′, `0, U¯in are physical properties, namely the turbulence intensity, in-
tegral length scale and mean inlet velocity of the premixture, whereas sL, δL
and Le are thermo-chemical properties of the mixture, namely the unstrained
laminar flame speed, laminar flame thickness, and the Lewis number. Further,
the dependence on flame configuration established experimentally by Shepherd
et al. (2002) is considered within a characteristic length scale L.
55
3. Parametric Simulations of the Bending Effect
The bending effect refers (Bradley, 1992; Lipatnikov and Chomiak, 2002;
Driscoll, 2008) to the non-linear variation of the curve sT (u′). Hence, in or-
der to capture the bending effect, it is essential that `0, U¯in, sL, δL,Le and
the flame configuration parameter L remain unchanged between any two mea-
surements. This constraint has been emphasised by Driscoll (2008) in light
of recent experiments as well as simulations wherein other parameters vary
simultaneously with u′. Such simulations may result in capturing a different
parametric variation of sT and, therefore, may not capture the bending effect
or lead to different conclusions on the same. Therefore, the focus here is to
design a set of simulations wherein u′ is the only parameter varied systemat-
ically. Furthermore, values of u′ are varied in a regime where the emergent
curve sT (u′) bends, i.e. sT records a bending effect. The selection of each
parameter maintained across the simulations is discussed as follows.
Thermo-chemical parameters
The thermochemical parameters: Lewis number Le, laminar flame speed sL
and laminar flame thickness δL are related to the choice of fuel-air mixture.
Combustion of stoichiometric methane-air mixtures is considered here for which
the Lewis number Le ∼ 1. Effects of pyrolysis as well as radiative heat transfer
are negligible for methane-air combustion. Hence, a relatively simple thermo-
chemical model can describe flame propagation accurately for these mixtures.
This has motivated the investigation of methane-air flames in many experi-
mental and simulation databases as listed by Lipatnikov and Chomiak (2002).
Subsequently, methane-air combustion chemistry is described using a single-
step kinetic model with unity Lewis number (see Chapter 2). The assumptions
of single-step chemistry together with a unity Lewis number drastically sim-
plify the modelling framework. Indeed, this provides a basic set of conditions
in which to investigate the salient mechanisms of a particular combustion phe-
nomenon. So this was thought to be a reasonable starting point for the in-
vestigation of the bending effect. It has been pointed out since that stretch
effects, especially the contributions of local flame curvature, may be signif-
icantly different for non-unity Lewis number flames. Hence, in Chapter 6,
56
3.1. Design conditions for bending
Table 3.1 Thermo-chemical input parameters.
Ambient pressure, Pa 1.01× 105 Pa
Inlet temperature, Tu 300 K
Pre-exponential constant, A 2.75× 107 m3kg−1s−1
Activation energy, Ea 108 Jkg−1
Gas constant, R 287.1Jkg−1K−1
Lewis number, Le 1.0
Prandtl Number, Pr 0.7
Schmidt Number, Sc 0.7
Heat release parameter, τ 6.7
Equivalence ratio, φ 1.0
Laminar flame speed, sL 0.39 ms−1
Laminar flame thickness, δL 3.6× 10−4 m
Reaction zone FWHM, δ 1.65× 10−4 m
these assumptions are lifted. Nevertheless, the relative simplicity of the sim-
ulations discussed here and until Chapter 6 has provided several insights in a
relatively inexpensive manner.
An unstrained 1D laminar flame solution profile calculated using the frame-
work discussed above is used as the initial thermo-chemical field: density ρ,
temperature T , and species mass-fractions Yα. The laminar flame speed sL
is calculated by integrating the unburned mixture reaction rate ω˙R across an
unstrained 1D laminar premixed flame
sL = − 1
ρuYR,u
∫ ∞
−∞
ω˙R(x)dx (3.2)
where ρu is the density of the unburned mixture and YR,u is the reactant mass
fraction in this mixture. The unstrained laminar flame thickness δL = δ0L is
calculated using the maximum temperature gradient method
δL =
Tb − Tu
(|dT/dx|)max
(3.3)
where Tu = 300K is the unburned gas temperature and Tb = 2300K is the
adiabatic burned gas temperature. The values of sL and δL are provided in
Table 3.1 along with other relevant thermo-chemical parameters.
57
3. Parametric Simulations of the Bending Effect
3.1.1 Flame configuration and regime
Flame configuration plays an important role in determining the burning veloc-
ity sT of the flame. Several realistic flame configurations have been considered
by previous researchers as discussed in Section 1.2. A freely propagating flame
is considered here instead for the following reasons
1. Jet flames and bluff-body stabilised flames require a huge domain to
account for the turbulent spread. Freely propagating flames require only
that the axial extent of the domain be several times the turbulent flame
brush thickness. Undue computational expense is, therefore, mitigated
in the latter.
2. Turbulent flow properties vary considerably over the flame surface in
realistic flame configurations. Freely propagating flames offer a means
to specify a single turbulence intensity u′ and quantify its effects on sT .
3. Effects of mean curvature are absent in statistically-planar propagating
flames. Therefore, the problem of turbulent burning velocity in this
configuration is the true analogue of the 1D propagating laminar flame.
Figure 3.1a illustrates the flame configuration using a typical 3D snapshot
obtained from one of the simulations. The flame surface (denoted by a single
(a) 3D snapshot (b) Axial slice
Figure 3.1 Inflow-outflow flame configuration used in the simulations. Left:
3D rendering of a propagating flame isosurface interacting with upstream tur-
bulent eddies in the domain. Right: 2D axial slice showing the location of the
flame and relative size of the domain.
58
3.1. Design conditions for bending
progress variable isosurface: red contour) propagates into upstream turbulence
(denoted by high-enstrophy contours: grey eddies) with a speed sT (red arrow).
In Figure 3.1b, a 2D axial slice is shown. The axial length of the domain is 3
times its width in spanwise directions. This helps accommodate flame prop-
agation in the axial direction and, at the same time, avoids flame-boundary
interactions. As discussed in Chapter 2, the transverse directions are specified
with periodic boundary conditions, whereas the upstream and downstream
faces (marked using blue arrows in Figure 3.1a) are specified with inflow and
outflow NSCBCs respectively.
Turbulence initialisation
Turbulence is initialised using a procedure outlined by Rogallo (1981). First,
a velocity field of desired turbulence intensity u′ is generated using an en-
ergy spectrum for homogeneous isotropic turbulence proposed by Batchelor
and Townsend (1948). As discussed in Chapter 2, the spectrum is specified
in Fourier space using two parameters: 1) u′, which specifies the turbulent
kinetic energy and 2) k0, which specifies the wavenumber of peak turbulent
kinetic energy. Next, the corresponding Cartesian field is computed using a
Fourier transform and applied to a cubic periodic domain of width L. The
Cartesian periodic velocity field is allowed to evolve for a brief period. No
forcing has been employed in the present simulations so that the turbulence
decays ever so slightly during this period. This time-evolved homogeneous
isotropic turbulence (HIT) solution is used as the initial velocity field in the
flame domain.
Since the stream-wise extent of the flame domain is 3 times its width L,
three identical HIT solutions are copied adjacent to each other for initialisation.
Further, an identical HIT solution is convected through the inlet. Due to the
time-evolution of the initial HIT solution, fewer transients are induced in the
corresponding flame simulation.
The parameters k0 and L together determine the integral length scale `0 and
the Taylor length scale λ for the initialised spectrum (Batchelor and Townsend,
59
3. Parametric Simulations of the Bending Effect
Ta
bl
e
3.
2
Tu
rb
ul
en
ce
pr
op
er
ti
es
an
d
di
m
en
si
on
le
ss
pa
ra
m
et
er
s.
C
as
e
−→
I
II
II
I
IV
V
In
le
t/
In
it
ia
lu
′
u
′ i/
s L
1.
50
7.
50
10
.0
0
20
.0
0
30
.0
0
le
ad
in
g
ed
ge
u
′
u
′ le
/s
L
0.
94
4.
50
6.
00
12
.0
0
18
.0
0
Tu
rb
ul
en
t
K
E
k
m
2
/s
2
0.
6
11
.7
20
.6
82
.5
18
5.
6
P
ea
k
K
E
w
av
en
um
be
r
κ
0
2.
0
In
te
gr
al
le
ng
th
sc
al
e
` 0
m
m
←−
0.
99
−→
K
ol
m
og
or
ov
sc
al
e
η
µ
m
74
.5
9
33
.3
5
28
.8
9
20
.4
2
16
.6
7
D
om
ai
n
le
ng
th
L
x
=
3
cm
1.
5
D
om
ai
n
w
id
th
L
y
,L
z
cm
←−
0.
5
−→
G
ri
d
si
ze
N
x
×
N
y
×
N
z
←−
28
8
×
96
×
96
−→
A
xi
al
m
es
h
re
so
lu
ti
on
∆
x
µ
m
52
.2
6
Sp
an
-w
is
e
re
so
lu
ti
on
∆
y
,∆
z
µ
m
52
.6
3
LE
T
O
T
τ 0
m
s
1.
59
0.
34
0.
25
0.
12
0.
08
C
he
m
ic
al
ti
m
e
τ c
m
s
0.
93
0.
93
0.
93
0.
93
0.
93
To
ta
lr
un
ti
m
e
T
=
4τ
0
m
s
6.
36
1.
28
1.
0
0.
48
0.
32
D
am
kö
hl
er
nu
m
be
r
D
a
1.
71
0.
36
0.
27
0.
13
0.
09
K
ar
lo
vi
tz
nu
m
be
r
K
a
1.
20
12
.3
6
19
.0
8
53
.9
6
99
.8
2
†
A
rr
ow
s
in
di
ca
te
th
at
em
pt
y
fie
ld
s
ha
ve
th
e
sa
m
e
va
lu
e
as
th
e
ce
nt
re
co
lu
m
n
fie
ld
of
th
e
co
rr
es
po
nd
in
g
ro
w
.
*
A
cr
on
ym
s:
L
E
T
O
T
(L
ar
ge
-E
dd
y
T
ur
n-
O
ve
r
T
im
e)
an
d
K
E
(K
in
et
ic
E
ne
rg
y)
.
60
3.1. Design conditions for bending
1948; Cant, 2013) as follows
`0 =
L√
2pik0
; λ =
L
2pik0
. (3.4)
Therefore, these parameters have been fixed throughout the dataset with the
intention of controlling `0 and λ of the prescribed turbulent flow field across
the entire dataset (refer to Table 3.2 for values). In fact, L is chosen so that
the integral length scale is several times the laminar flame thickness: `0 ≈ 3δL.
As the flame simulation evolves in time, the inflow establishes a statistically-
stationary spatial profile of turbulence from the inlet to the leading edge of the
flame after a few eddy turn-over times. The intensity u′ of turbulence varies to
noticeable degree over this region. Hence, a leading edge turbulence intensity
u′le calculated at a planar section corresponding to c = 0.01 has been supplied
alongside the inlet turbulence intensity u′i in Table 3.2. The length scales `0
and λ vary to a relatively minor extent in comparison. The leading edge values
of these scales have, therefore, not been reported. Conforming with standard
practice, all following values of u′ correspond to the inlet or initial intensity.
Since the integral length scale `0 is fixed across all simulations, the Kol-
mogorov scale η becomes very small in high u′ turbulence. Due to this, com-
plete resolution of high u′ turbulence up to the Nyquist frequency imposes
a severe computational challenge. Hence, a fixed mesh spacing of the order
∆x ∼ O(η) of Kolmogorov scales is specified throughout the dataset. In the
high u′ cases, ∆x is up to a few times the Kolmogorov length scale η calculated
based on the Batchelor-Townsend spectrum as
η =
[
ν2
15u′2pi2k
2
0
]1/4
(3.5)
where k0 is the wavenumber of peak turbulent kinetic energy normalised by
the domain width L. Coarse-grained simulations such as these have been
conducted previously (Poludnenko and Oran, 2011). The relatively larger grid
sizes contribute small-scale dissipation by numerical means, so that this frame-
work is numerically valid. Nevertheless, fine-grained simulations are being
61
3. Parametric Simulations of the Bending Effect
conducted (see Chapter 6) to ascertain the ramifications of the present study.
Combustion regime
The selection of all parameters on the r.h.s. of Equation 3.1 has been dis-
cussed above. The integral length scale `0 together with the chosen turbulence
intensity u′ specifies the regime of the simulation. The combustion regime of
simulations is estimated here using a Karlovitz number
Ka ≡
(
u′
s0L
) 3
2
(
δ0L
`0
) 1
2
, (3.6)
based on the Taylor scale strain rate (Bradley, 1992). The Borghi diagram in
Figure 3.2 shows the regime of each of the simulations marked as a single point
(labelled I - V). The choice of u′ values for the simulations has been such that
by increasing u′ successively between each simulated case, the dataset bridges
the flamelet regime and the Broken Reaction Zones (BRZ) regime. These
values of the Karlovitz number were calculated using inlet/initial turbulence
properties. These values have been reported along with the Damköhler number
Figure 3.2 Regime diagram showing the region of parameter space covered by
the present simulations. The entire region between the flamelet regime and
broken reaction zones regime is spanned by increasing u′ successively across
five simulations labelled I – V.
62
3.2. Recorded effects of turbulence intensity
Da for the entire dataset in Table 3.2.
The dataset as prepared above consists of 5 different simulations each with a
specified inlet u′. The DARWIN (Intel Sandy Bridge E5-2670) and ARCHER (Cray
XC30) supercomputers were used to compute the solutions. Each simulation
used approximately Tc ∼ 2.5× 104 core-hours. The simulations were each run
for t = 4τ0 where τ0 = `0/u′ is the integral time scale. During this period, an
initial phase of monotonic growth in sT is observed up to t . 2τ0, followed by
a phase of oscillations in sT about an average value specific to each case. All
statistics are collected during this phase at the instant t = 4τ0 as it provided
representative results for the period of statistical-stationarity. Averaging, for
instance, is conducted at the snapshot corresponding to this instant over the
spanwise directions.
3.2 Recorded effects of turbulence intensity
3.2.1 Images of the turbulent flame brush
Figure 3.3 shows cropped slices of the domain corresponding to each of the sim-
ulations along the columns with the turbulence intensity noted on top. The
time evolution of each of the flames is shown up to 4 eddy turn-over times τ0
across the different rows. Each panel shows a 2D slice of the instantaneous
progress variable contours at the corresponding time noted on the left. It is ev-
ident that the turbulent flame brush thickness increases across the simulations
as u′ is increased from 1.5sL to 30sL (left to right along any row in Figure 3.3).
This thickness also increases over time in each simulation up to a certain max-
imum thickness as presented in (Nivarti and Cant, 2015a). Whereas wrinkling
governs the flame brush thickness in the low to moderate intensity cases, con-
siderable pocket-formation takes place in the high-intensity cases u′ > 10sL
determining the overall flame brush thickness.
A closer look at the time-evolution of flame wrinkling behaviour shows
that wrinkling takes place at a much smaller time scale and much smaller
length scales in the high-intensity cases. Low-intensity cases with u′ < 10sL
are dominated by slowly evolving, large-scale wrinkling behaviour; i.e. each
63
3. Parametric Simulations of the Bending Effect
Figure 3.3 Contours of progress variable (colour) for each of the five simulations
(columns, with u′ value noted on top). Time evolution of the respective flame
occurs down each column, with elapsed time noted on the left in each row.
wrinkle evolves over a significantly greater number of rows in the columns to
the left. This highlights an interesting divide between low- and high-intensity
cases. Given the constraints of the physical geometry, large-scale wrinkles of
low-intensity turbulence cases keep the flame wrinkled for longer durations of
time. In contrast, the wrinkling is dominant at smaller scales in the high-
64
3.2. Recorded effects of turbulence intensity
intensity cases so that the flames are flatter on average than the low-intensity
cases.
In order to inspect further the physical differences between the low- and
high-intensity cases, three dimensional instantaneous snapshots at t = 4τ0
of the lowest and highest u′ flame brushes are presented in Figure 3.4. The
differences are striking. In the low u′ case (Figure 3.4a, above), the leading edge
c = 0.01 isosurface (coloured blue) and the trailing edge c = 0.99 isosurface
(coloured red) are close to and roughly parallel to each other: the flame brush
is thin. Moreover, the wrinkles scale with the domain size `w ∼ O(L). Similar
large-scale wrinkling was also observed previously in laminar flame simulations
of low-Karlovitz flame-vortex interaction (Poinsot et al., 1991b).
In contrast, the high u′ flame brush (Figure 3.4b, below) has an extremely
irregular and fragmented leading edge. Only a few large-scale wrinkles are
noticeable (analysed quantitatively later). The instantaneous flame surface
consists instead of many small-scale wrinkles `w  O(L). The trailing edge
is also highly wrinkled, but the size of these wrinkles is larger than those on
the leading edge. The two edges of the flame brush are much farther apart
than in the low u′ case. These images resemble DNS results obtained for
premixed single-step flames with u′ ≈ 35sL in the BRZ regime by Poludnenko
and Oran (2011). This sharp contrast between the two instantaneous flame
brushes illustrates visually the relationship between u′ and the dominant scale
of flame wrinkling `w.
Flame brush edges
In Figure 3.4, the effects on the leading edge appear to be different from those
on the trailing edge. A quantitative analysis of these surfaces is presented in
Figure 3.5. Distributions of mean curvature hm and curvature shape factor Sh
are shown for low-intensity (black lines) and high-intensity (red lines) instan-
taneous flame brush edges. The leading edge (solid lines) has been taken as
the c = 0.01 contour and the trailing edge (dashed lines) has been taken as the
c = 0.99 contour. Previous DNS studies have adopted related definitions with
slightly different choices of progress variable contour values (Lipatnikov and
65
3. Parametric Simulations of the Bending Effect
(a) Lowest intensity case I: u′ = 1.5sL
(b) Highest intensity case V: u′ = 30sL
Figure 3.4 Instantaneous snapshots of the turbulent flame leading edge c = 0.01
(coloured blue) and trailing edge c = 0.99 (coloured red) in cases with the
lowest and highest u′, respectively. The arrow indicates the direction of flame
propagation.
66
3.2. Recorded effects of turbulence intensity
Chomiak, 2002; Driscoll, 2008). It has been verified in the present cases that
the precise choice does not have a significant effect on the statistics obtained.
Mean curvature statistics
Figure 3.5a shows a comparison of distribution of mean curvature given by
hm = (hl + hu)/2, the mean of local principal curvatures hl and hu, for the
flame surfaces shown in Figure 3.4. The low-u′ flame brush edges (black lines)
exhibit little variation in hm about the mean in comparison to the broad range
of mean curvatures observed in the PDF for the high u′ case (red lines). This
shows that the flame surface is wrinkled at smaller scales at higher intensities;
i.e. the scales of flame surface wrinkling are inversely related to the intensity
of turbulence. In both cases, however, the mean value of mean curvature is
close to zero for both leading (solid lines) and trailing edges (dashed lines).
In the higher u′ case, the peak probability of hm occurs closer to zero than
in low u′ for both edges. This represents quantitatively the reduced large-
scale curvature in the higher u′ case noticed qualitatively in Figure 3.3 and
Figure 3.4.
Furthermore, in the low-u′ case (black lines in Figure 3.5a), the leading
edge (solid black line) has a high likelihood of positive mean curvature – the
(a) Mean Curvature hm (b) Shape Factor Sh
Figure 3.5 Distributions of the instantaneous mean curvature hm and curvature
shape factor Sh for the leading edge (solid lines) as well as trailing edge (dashed
lines) surfaces in cases I and V (black and red lines, respectively).
67
3. Parametric Simulations of the Bending Effect
trailing edge (dashed black line) follows this distribution except for a slightly
higher probability of negative curvature. In the high u′ case (red lines), both
leading as well as trailing edges have significant probability of negative mean
curvature. The distribution of hm on the trailing edge in the high u′ case
(dashed red line in Figure 3.5a) resembles the distribution on the leading edge
of the low u′ case (solid black line in Figure 3.5a) but with the mean value
negative 〈hm〉 < 0.
Curvature shape factor statistics
Figure 3.5b presents statistics of the curvature shape factor Sh = hl/hu defined
as the ratio of the minimum (hl) and maximum (hu) local principal curvatures
by magnitude. Broadly, three shapes are captured: a) Sh = 1 represents
locally spherical topology, b) Sh = 0 represents cylindrical topology, and c)
Sh = −1 represents saddle shapes. These topologies are evident upon closer
visual examination of the flame brush surfaces in Figure 3.4.
The wrinkles in the low u′ flame brush (pictured in Figure 3.4a) form
smooth and rounded ridges. Notably, the wrinkles do not protrude into the re-
actants or products but instead are gently warped in between. No “finger”-like
structures have been seen here such as those observed more recently elsewhere
(Lipatnikov et al., 2015). Since the trailing edge is parallel to the leading edge
in this case, both surfaces have a similar distribution P(Sh) (black lines in
Figure 3.5b). Except for a small probability of saddle shapes (Sh = −1), the
surfaces tend to have cylindrical topology (|hl|  |hu|) as is apparent in the
large-scale thin rounded ‘folds’ in Figure 3.4a described earlier. This distri-
bution of Sh has been noted previously in the simulation of material surfaces
evolving in turbulence (Pope et al., 1989). It has also been observed more
recently in the DNS of a slot-jet premixed flame in shear turbulence (Griffiths
et al., 2015). A non-negligible probability of spherical shapes Sh = 0 is noted
here in addition.
Spherical and ellipsoidal shapes with hl, hu > 0 are less common in high u′
flames (red lines in Figure 3.5b) as opposed to the low u′ flames (black lines).
Instead, the high u′ flame brush shows a greater probability of saddle-like
68
3.2. Recorded effects of turbulence intensity
(a) Leading edge, c = 0.01
(b) Trailing edge, c = 0.99
Figure 3.6 Magnified view of the instantaneous leading edge (blue) and trailing
edge (red) surfaces at t = 4τ0 as computed for case V with u′ = 30sL.
(given by Sh ≈ −1) as well as of cylindrical (Sh ≈ 0) shapes. Unlike the mean
curvature PDFs, the distributions of P(Sh) on the trailing edge of the high u′
case (dashed red line in Figure 3.5b) is not similar to that on the leading edge
of the low u′ case (solid black line in Figure 3.5b).
69
3. Parametric Simulations of the Bending Effect
The topological features of the high u′ case are inspected further in Fig-
ure 3.6 using a magnified view of the surfaces. The leading edge (coloured blue
in Figure 3.6a) is highly convoluted and consists mainly of twisted ‘folds’ and
the edges of these folds have a locally cylindrical topology. The sides of these
folds are twisted at scales larger than the maximum principal curvature at the
edge. Evidently, these folds contribute to the increased probability of Sh = 0
(cylindrical shapes) and Sh < 0 in the high u′ distribution in Figure 3.5b.
The relatively smooth trailing edge consists of long extended ‘fingers’ instead.
The fingers have spherical tips which contribute to Sh ≈ 1; the sides of these
fingers are cylindrical, and the base forms a saddle shape. Some of the fingers
merge into each other and form large-scale toroidal structures. Ongoing work
is focused on investigating further the preference for these topological shapes.
Reaction zone contours
Evidently, the 20-fold increase in u′ from case I to case V has a marked influence
on the turbulent flame brush. The topological observations noted above may
suggest that the internal turbulent flame structure might also differ markedly.
In Figure 3.7, reaction rate contours are shown (panels above) using black and
white shades in axial sections of the domain. Progress variable contours (black
lines) shown alongside indicate that reaction ω˙ 6= 0 is confined to thin regions
between c = 0.5 and c = 0.99. At higher u′ (Figure 3.7b), this zone of reaction
appears broadened yet it remains confined to within the same contour lines.
The reaction rate scatter ω˙(c) shown below in Figure 3.7 indicates only slight
deviations from the laminar flamelet profile (red line). The reaction rate peak
occurs at c ≈ 0.8, exactly as for a flamelet. Even for the highest u′ flame
(Figure 3.7d), a magnified view of the scatter reveals very minor differences
from the flamelet profile. In other words, the reaction zones are not affected
even though the surface topology changes drastically.
In Figure 3.8, 3D contours of the progress variable isosurface c = 0.8 are
shown to emphasise this point further. In the backdrop, slices of the reaction
rate ω˙(x) have been shown to reveal the reaction profile alongside the convo-
luted surface. In the low u′ flame (above), the c = 0.8 isosurface is attached
70
3.2. Recorded effects of turbulence intensity
(a) Contours for case I: u′ = 1.5sL (b) Contours for case V: u′ = 30.0sL
(c) Scatter for case I (d) Scatter for case V
Figure 3.7 Reaction rate contours ω˙(x, z; y = 5L/6) in cropped axial slices of
the domain (above) and reaction rate scatter ω˙(c) across the entire domain
(below) at t = 4τ0 for cases I and V.
to the reaction zone as expected. This continues to be the case even in the
highest u′ case (below) where the c = 0.8 isosurface undergoes considerable
wrinkling. The reaction rate contours (visible as a thin red band in the slices)
continue to remain smooth and connected with the exception of minor local
variations due to occasional flame self-intersections.
Both visual and statistical assessments presented here indicate clearly that
even though the turbulent flame brush is affected severely under increasing u′,
the reaction zone (i.e. the region corresponding to 0.7 < c < 0.9) retains a
laminar flamelet structure even up to u′ = 30sL. These observations, therefore,
71
3. Parametric Simulations of the Bending Effect
(a) Case I: u′ = 1.5sL.
(b) Case V: u′ = 30.0sL.
Figure 3.8 Progress variable isosurface c∗ = 0.8, corresponding to the peak
reaction rate, in the lowest and highest intensity cases. Axial sections of the
domain show the variation of reaction rate ω˙ across the domain.
72
3.2. Recorded effects of turbulence intensity
conform with the scaling arguments for the TRZ regime (Peters, 1999) as
discussed in Section 1.2.2. This provides an indication that flame propagation
and the associated turbulent burning velocity sT is not affected by quenching
since that would necessarily cause the reaction zone to change.
3.2.2 Computation of turbulent burning velocity
The turbulent burning velocity sT is calculated as the integral
sT ≡ m˙F
ρuYu,FA0
= − 1
ρuYu,FA0
∫
V
ω˙F dV, (3.7)
thereby defining sT as the global fuel consumption speed expressed using the
mass transformation rate m˙F or, equivalently, the integrated reaction rate
ω˙F of fuel in an unburned premixture of density ρu and fuel mass fraction
Yu,F ; finally, A0 is the flow cross-section area. Here, F denotes the fuel which
corresponds to the unburned reactant mixture.
In each of the simulations, sT is computed using the expression above. It
is observed that sT grows with time up to a period of t ∼ 2τ0 after which the
value of sT oscillates about a value specific to the specific case. In Figure 3.9,
the computed values for turbulent burning velocity at t = 4τ0 are shown for
each of the cases. The inset shows the temporal variation for a typical case;
oscillations after a certain period of time are within a small error margin of
the value at t = 4τ0 (grey line and circle) used to plot the variation of sT
across the dataset. Since all parameters other than u′ are fixed across the set
of simulations, the graph in Figure 3.9 shows the curve sT (u′). Clearly, the
bending effect has been captured – the deviation from linear enhancement of
sT is prominent for cases IV and V with u′ > 10sL.
The normalized turbulent burning velocity sT/sL values shown in Figure 3.9
are similar in order of magnitude O(sT/sL) to the values obtained in previous
DNS studies (Aspden et al., 2008; Poludnenko and Oran, 2011; Savard and
Blanquart, 2015). The order of magnitude O(sT/sL) obtained using DNS is
typically smaller than that obtained experimentally (Bradley, 1992; Shy et al.,
2000b; Gülder and Smallwood, 2007). For example, stoichiometric methane-air
73
3. Parametric Simulations of the Bending Effect
Figure 3.9 Computed variation of the normalized turbulent burning velocity
sT/sL with normalized turbulence intensity u′/sL. The inset shows a typical
variation of sT (t) and the corresponding value of sT chosen from this variation.
flames have been considered experimentally in downward propagating flames
(Shy et al., 2000b) up to u′/sL ≈ 15. The corresponding measurements indicate
sT/sL ≈ 10 − 12; i.e. a few times the predictions obtained using DNS. These
differences may arise due to the differences in Lewis numbers of chemically
reacting species between experimental and simulated flames. Perhaps more
importantly, these differences may be due to the dependence of sT on the
geometry. In other words, the relatively low values of sT/sL observed here in
comparison to experiments could be related to the small geometrical extent of
the DNS. The relatively larger integral length scales `0 ≈ 10 − 20 cm in the
experiments of Shy et al. (2000b) are as yet unable to be considered by DNS.
Nevertheless, the present DNS dataset provides the ability to investigate the
mechanism of the observed bending effect.
3.3 Summary: on the bending effect recorded
A parametric DNS dataset has been designed with due considerations of the
notion of the bending effect. The most basic conditions that can reproduce the
bending effect are prescribed upon an initially-laminar stoichiometric methane-
74
3.3. Summary: on the bending effect recorded
air flame. All controllable thermo-chemical and turbulence parameters are
maintained across the set of simulations; these include the flame configuration
and domain size as well as numerical parameters such as the mesh spacing.
The turbulent intensity u′ is the only parameter varied across the simulations.
The curve sT (u′) thus obtained records a bending effect which is, therefore, a
consequence of increasing u′ only.
The present thermo-chemical framework relies on the twin assumptions
of unity Lewis number and single-step chemistry. This has helped simplify
the description and analysis and narrow the potential causes of bending to
kinematic rather than chemical effects such as local quenching. Since the
framework indeed records a bending effect (Figure 3.9), it follows that the
mechanism responsible for the bending effect recorded must be kinematic and
not chemical in nature. It is plausible that the underlying mechanism in the
present dataset may also operate in flames with more complicated thermo-
chemical properties although the significance of such a mechanism remains to
be evaluated by using precisely such a framework (see Chapter 6). In the
following two chapters, the validity of Damköhler’s hypotheses is investigated
in the present DNS dataset with the objective of highlighting the mechanism
of the bending effect recorded.
75

Chapter 4
Damköhler’s First Hypothesis
In this chapter, the validity of Damköhler’s first hypothesis (DH1) is investi-
gated within the bending effect recorded in the previous chapter. Recall the
modified form of DH1
sT
sL
= I0
AT
AL
(4.1)
where I0 the so-called stretch factor is a function of Le. The equivalent laminar
flame surface area AL is the square cross-section area of the domain. Hence,
given a fuel-air mixture, the primary task that remains to validate DH1 is to
compute the turbulent flame surface area AT . This computation is pursued
within the framework of the Flame Surface Density (FSD) approach. The
FSD approach was discussed in Chapter 1; it provides a measure of the flame
surface area AT . Moreover, it provides the factors influencing this measure;
namely, the stretch rate components. Hence this allows for a conceptualisation
of the bending effect in terms of the stretch rate components: strain rate and
curvature propagation. In the passages below, relevant concepts are discussed
followed by its application to the bending effect recorded.
4.1 Flame Surface Density approach
Any progress variable isosurface c = c∗ within the turbulent flame brush may
be considered as a constant-property material surface (Pope, 1988) and may
hence be used for the computation of AT . In the framework of the FSD ap-
77
4. Damköhler’s First Hypothesis
proach, each point on such an instantaneous surface is specified by an in-
finitesimal area δA within a fixed infinitesimal volume δV . The local surface-
to-volume ratio may then be defined as
Σ′(x, t; c = c∗) ≡ δA(x, t; c = c
∗)
δV
, (4.2)
at each point x|c=c∗ on the surface. The local surface area δA may be cal-
culated using components of the local gradient ∇c(x, t) of the isosurface. It
has been found (Boger et al., 1998), however, that in the limit of LES filter
width lim∆→0 Σ′ = |∇c|. This approximation holds in DNS because ∆ → 0
necessarily. Hence, Σ′ = |∇c| in the present study.
The quantity Σ′ is a local measure of the degree to which a progress variable
isosurface c = c∗ is wrinkled because points on the surface that have a larger
surface area within a local infinitesimal volume necessarily have higher local
mean curvatures. At the same time, it is noteworthy that the progress variable
gradient ∇c is high only where reaction rate is non-zero. Hence, the flame
surface-to-volume ratio when defined using the progress variable gradient Σ′ ≡
|∇c| combines two aspects – the physical: surface wrinkling and the chemical:
reaction rate.
Generalised FSD profiles
The expectation of Σ′ is denoted by Σ which represents the mean flame surface
to volume ratio. It is calculated here as
Σ(x, t) ≡ 〈Σ′(x, t)〉 = 〈|∇c(x, t)|〉, (4.3)
and it is, therefore, a statistical quantity – also known as Flame Surface Density
(FSD) – that is often measured in experiments (Driscoll, 2008), and is central
to combustion modelling (Trouvé and Poinsot, 1994; Veynante and Vervisch,
2002). Indeed, the transport of this quantity was discussed earlier in Section 1.2
as a modelling framework. The averaging in Equation 4.3 may be performed
temporally over different ensembles or spatially over statistically homogeneous
regions in the domain. In the present study, the averaging is performed at a
78
4.1. Flame Surface Density approach
specific instant (t = 4τ0) over span-wise planes of the domain which provide a
statistically homogeneous sample. In such statistically-stationary flames, the
FSD is expected to be a function Σ(x) of axial distance only.
The profile Σ(x) is compared with the laminar profile (red line) in Fig-
ure 4.1a for three cases of the present dataset (I: u′ = 1.5sL, III: u′ = 10sL
and V: u′ = 30sL). The broadening of the Σ(x) profile with increasing u′ indi-
cates that the flame surface extends over a larger portion of the domain, i.e.
the mean turbulent flame brush thickens with increasing u′. This is in agree-
ment with the images shown earlier in Figure 3.4 as also with experimental
measurements of FSD discussed by Driscoll (2008).
Conventionally, the variation of FSD over mean progress variable c values
is compared with the laminar flamelet FSD to assess the effects of turbulence
on the flame brush. Figure 4.1b shows Σ as a function of c for three repre-
sentative cases (I, III and V) exhibits a peak in Σ close to c ≈ 0.6 similar
to the laminar flamelet profile (blue line). In comparison, the Σ peak shifts
towards a lower value of c at u′ = 10sL (dashed line). At higher intensities
still (dotted line), the peak shifts further towards the leading edge. While the
figure includes only three of the cases for clarity, the remaining cases in the
(a) Expected FSD, Σ along x. (b) Expected FSD, Σ along c.
Figure 4.1 Expectation of FSD Σ as a function of the axial coordinate x (left)
and mean progress variable c (right) for three cases (I, III, and V) spanning
the dataset. The laminar flamelet solution (red line) is shown for comparison.
79
4. Damköhler’s First Hypothesis
dataset also follow this trend. Such a shift in the peak of Σ with increasing
u′ has been reported in a recent experimental study (Kheirkhah and Gülder,
2015) where it is linked with counter-gradient and gradient transport regimes.
Further inspection of Figure 4.1b shows that the FSD values for turbulent
flames are lower than the laminar flamelet values. This again conforms with
the measurements of the aforementioned experimental study (Kheirkhah and
Gülder, 2015). However, the trend observed here is that the FSD values for
turbulent flames increase with increasing u′ which stands in contrast to the
previous measurements (Kheirkhah and Gülder, 2015).
Scatter in progress variable space
A key challenge in the analysis of turbulent premixed flames is to compare
the local internal flame structure with a 1D laminar flamelet profile. Due to
wrinkling of the flame surface, the surface normal vector is not always aligned
with the streamwise direction. It is common practice to sample over span-wise
sections instead of progress variable isosurfaces. The section-sampling method
is often applied for evaluating mean profiles of Σ as a function of c¯∗ obtained by
averaging Σ′(x, y = y∗, z = z∗) throughout a slice which corresponds to an av-
(a) Section-sampling (b) Surface-sampling
Figure 4.2 Illustration of sampling methods (dashed lines). Left: sampling over
spanwise sections of the domain across the mean turbulent flame brush (c¯∗0 <
c¯ < c¯∗N). Right: sampling over instantaneous progress variable c isosurfaces
(solid lines) within the instantaneous turbulent flame brush (c∗0 < c < c∗N).
80
4.1. Flame Surface Density approach
erage progress variable c¯∗ (illustrated in Figure 4.2a). Since wrinkling changes
the direction of the local structure from that of the slices, the slice-sampling
method can provide a misleading comparison between the local turbulent flame
structure and the laminar flamelet structure. A similar concern was expressed
previously (Shepherd et al., 2002) in order to avoid incorrect assessments of
“flame broadening” relying on slice-sampling. Hence, Σ′(x; y∗, z∗) is not an
appropriate benchmark of comparison with the laminar flamelet profile Σ′L(x)
even in the statistically 1D case considered here.
Instead, the scatter of points over the progress variable space Σ′(c) makes
for a more direct comparison with the laminar flamelet profile. Geometrically,
this method is analogous to sampling over progress variable isosurfaces c =
c∗ as shown in Figure 4.2b and is, therefore, termed surface-sampling. The
comparison with the laminar flamelet profile is direct and any deviations from
the flamelet profile can be noticed straightforwardly. Hence, even though the
exact spatial point corresponding to a given value of FSD may not be clear in
the scatter, this is the apt and preferred method of comparison of turbulent
local internal structure with laminar flamelets.
The scatter Σ′(c) obtained using the surface-sampling method is presented
in the left column in Figure 4.3 for the three cases I, III and V. For the low
u′ case I, the overall Σ′ profile is not disturbed noticeably in comparison to
the laminar flamelet profile Σ′L = |∇c|L (red line). As u′ is increased, the
deviation of Σ′ from the flamelet profile at any given value of c (isosurface)
increases. For u′ ≥ 10sL, large deviations from the flamelet profile are evident
close to the leading edge c ≤ 0.4. The reaction zone c ≥ 0.7, however, is much
less affected under increasing u′. Even for the highest u′ flames, only minor
deviations are visible. These correspond to the minor deviation in reaction
rate profile observed in Figure 3.7. This implies that the flame brush tends
to be wrinkled by length scales that do not disturb significantly the internal
progress variable gradient ∇c in the reaction zone.
Instantaneous slices of Σ′(x, y = 5/6Ly, z = 5/6Lz) at t = 4τ0 are shown
alongside on the right column in Fig. 4.3 for the same three cases. The extent
of the turbulent flame brush is shown by overlaying progress variable con-
81
4. Damköhler’s First Hypothesis
(a) Case I: u′ = 1.5sL
(b) Case III: u′ = 10.0sL
(c) Case V: u′ = 30.0sL
Figure 4.3 Instantaneous normalized Σ′(c) scatter (left column) at t = 4τ0
compared with the laminar flamelet profile (red line). Correponding contours
of Σ′(x) (right column) are shown in an axial section overlaid with progress
variable contours (black and white lines).
82
4.1. Flame Surface Density approach
tours in white. The contours in black represent the interior flame structure
0.4 < c < 0.8. The colour bars show clearly that the peak magnitude Σ′max
(red colour) in a given slice increases with u′. Since this magnitude is higher
than the that in laminar flames, these regions indicate steepening of the gra-
dients and, thereby, thinning of flame thickness by turbulence. In addition,
the spatial region where the peak is observed appears to be increasingly frag-
mented, occurring selectively in regions of negative curvature (note: the flame
is propagating to the left). This localisation of Σ′max = |∇c|max has been at-
tributed to a preferential “tight packing” of progress variable contours due to
turbulence (Poludnenko and Oran, 2011). Elsewhere (Steinberg et al., 2012),
the same behaviour is referred to as the preferential alignment of the compres-
sive eigenvector of strain rate with the c gradient. Further investigation of
flame surface topology will help establish reasons for this observation.
4.1.1 Computation of flame surface area
As discussed in the previous chapter as well as in the passages above, the
reaction zone appears to retain integrity across the entire dataset. Hence, the
turbulent flame surface area AT may be estimated to within a scaled constant
using an isosurface within the reaction zone (Filatyev et al., 2005; Driscoll,
2008). In the FSD approach, the area AT can be expressed as the volume
integral
AT =
∫ ∞
−∞
Σ′δ(c− c∗)dV, (4.4)
of the surface-to-volume ratio Σ′ conditioned on the progress variable value
c = c∗ corresponding to the reaction zone isosurface considered. Presently,
this is chosen as the peak reaction rate isosurface, i.e. c∗ = 0.8, but it was
found that any value 0.6 < c∗ < 1.0 corresponding to the reaction zone is
appropriate. The Dirac delta function is evaluated by conditioning the integral
to the spatial region corresponding to c = c∗±0.01. Hence, the DH1 expression
sT
sL
= I0
AT
AL
= I0
1
AL
∫ ∞
−∞
Σ′δ(c− c∗)dV (4.5)
83
4. Damköhler’s First Hypothesis
Figure 4.4 Computed flame surface area ratio AT/AL compared with nor-
malised turbulent burning velocity sT/sL in the recorded bending effect.
may be validated using the conditioned Σ′ volume integral. This approach may
only be followed in DNS studies where the local progress variable corresponding
to the local FSD Σ′ is available for computation. Hence, this equation is used
here and is consistent with with other DNS studies (Aspden et al., 2011; La-
pointe et al., 2015; Poludnenko, 2015). In experimental studies (Bagdanavicius
et al., 2015), however, this approach using Equation 4.5 is often not feasible.
Therefore, in such studies, a mean flame surface area to volume ratio of the
entire flame surface is computed (also, confusingly, denoted Σ in such studies)
and integrated over the entire volume to calculate AT . These differences in
the computation procedure are considered in Chapter 6 where comparisons are
being made with experiments in the explosion vessel references above.
In Fig. 4.4, it is seen that the ratio of the reaction zone surface area
AT (c
∗ = 0.8) to that of the laminar surface area AL indeed matches closely
the burning velocity ratio sT/sL. Hence, Damköhler’s first hypothesis (Equa-
tion 1.4, Equation 4.1) is valid for the entire range of u′ simulated with the
stretch factor I0 ≈ 1 throughout. Any deviations in sT/sL from AT/AL are
within a few percent of the sT/sL value as was also noted by previous DNS
works (Sankaran et al., 2007; Poludnenko and Oran, 2011). Even in the highest
u′ case, which lies at the TRZ/BRZ regime boundary, AT predicts accurately
84
4.1. Flame Surface Density approach
the burning velocity sT . Such a close correspondence of sT/sL and AT/AL
stands in contrast to some recent experiments in “V” and Bunsen configura-
tions (Filatyev et al., 2005; Gülder and Smallwood, 2007; Bradley et al., 2013).
Ongoing work (Ahmed et al., 2017) is aimed at understanding the role of tur-
bulence length scale `0 and flame configuration in effecting these differences.
Preliminary validation with detailed chemistry
The validity of DH1 throughout the bending effect carries significant impli-
cations. The present framework considers a simplistic 1-step chemical kinetic
model elaborated in Chapter 2. In order to verify the scope and validity of fur-
ther analysis conducted on this dataset, comparisons were made (Nivarti and
Cant, 2015b) using a detailed 25-step chemical kinetic model due to Smooke
and Giovangigli (1991) as discussed briefly in Chapter 2 and elaborated in
Chapter 6. Parametric simulations of flame propagation using this detailed
model have not been feasible due to high computational expense. Therefore,
the case III with u′ = 10sL was chosen as a point of verification. Using the
identical flame configuration, ICs and BCs a case was set up with the only dif-
ference that the initial laminar flame was described using the detailed chemical
kinetic model.
Results of the chemical mechanism comparisons are shown in Figure 4.5.
A snapshot at t = 4τ0 of the propagating flame as recorded using 1-step chem-
istry in the present dataset is shown in Figure 4.5a. In line with the discussions
above, the pre-heat zone appears thickened (between the black progress vari-
able contour lines 0.01 and 0.5), whereas the reaction zone (between contour
lines 0.5 and 0.99) remains thin and the respective contour lines remain paral-
lel to each other. The corresponding snapshot for detailed chemistry is shown
in Figure 4.5b. The similarity with the 1-step chemistry simulation is striking.
Most importantly, the reaction zone remains thin throughout, even though a
similar degree of pocket formation is evident in both cases.
Investigating further, the time-evolving turbulent burning velocity sT (t) is
plotted against the turbulent flame surface area AT (t) for both cases in Fig-
ure 4.5c. As expected from the discussions so far, the 1-step chemistry flame
85
4. Damköhler’s First Hypothesis
(a) 1-step kinetics (b) 25-step kinetics
(c) Comparison of sT (AT ) (d) YCH2O × YOH
Figure 4.5 Comparisons of turbulent flame structure and propagation charac-
teristics simulated using 1-step and 25-step reduced chemical kinetic models
for case III with u′ = 10sL.
(black circles) displays conformance with Damköhler’s first hypothesis (DH1)
at all times. This corresponds to a stretch factor value I0 ≈ 1. The con-
formance with DH1 remains true, however, for the 25-step chemistry flame as
well, with the only difference that the stretch factor I0 ≈ 0.8 is slightly lower in
comparison. Qualitatively, the conformance with DH1 is also noted in the con-
tours of heat release indicated using the product of mass fractions YCH2O×YOH.
No signs of localised flame extinctions have been evident. These comparisons
seem to suggest that differences in sT recorded by adopting detailed chemistry
models may be subject to thermo-diffusive effects which influence I0 rather
86
4.2. Stretch rate components
than localised flame extinctions which influence the validity of DH1. Still, the
study pursued in Chapter 6 aims to address comprehensively the flame prop-
agation behaviour using detailed chemical kinetics. At present, however, this
comparative study serves to give additional weight to the results obtained in
the present dataset using single-step chemistry.
Based on the discussions above, it is apparent that the reaction zone retains
globally its internal structure. DH1 is satisfied for all turbulence intensities
within the dataset. Hence, the wrinkled turbulent reaction zone can be treated
as an ensemble of locally-laminar reaction zones without serious loss of accu-
racy. This implies that the bending effect recorded within the dataset must
also be explained in terms of a flamelet description of the turbulent flame.
Further investigation of the FSD equation and its source-terms is carried out
in order to find possible underlying mechanisms of the bending effect.
4.2 Stretch rate components
An equation for the turbulent transport of the local surface-to-volume ratio
Σ′(x) of a surface propagating in turbulence has been derived previously (Pope,
1988; Candel and Poinsot, 1990) as
∂Σ′
∂t
+∇ ·
(
X˙Σ′
)
= s˙Σ′ (4.6)
where X˙(x) = u+sdn̂ is the total velocity of the point which propagates with
displacement speed sd along its normal n̂ given by n̂(x)|c=c∗ = −(∇c/|∇c|)|c=c∗ .
Equation 4.6 shows that stretch rate s˙ regulates the local production and de-
struction of surface-to-volume ratio Σ′(x, t > t◦) in proportion to the available
local surface-to-volume ratio at a given instant t◦. The role of the stretch rate
is, therefore, investigated across the dataset. First, it is recalled that stretch
rate s˙ can be decomposed (Pope, 1988) into
s˙ = at + sdhm, (4.7)
where at and sdhm correspond to local tangential strain rate and local flame
87
4. Damköhler’s First Hypothesis
propagation rate respectively. These are given by
at = ∇ · u− n̂ n̂ : ∇u, and (4.8)
hm =
1
2
∇ · n̂, (4.9)
and the local surface normal n̂ is defined (as above) to be positive in the
direction of the unburned mixture. These individual components of the stretch
rate source term are analysed in the following passages.
4.2.1 Statistics of strain and propagation terms
Tangential strain
The distribution of tangential strain rate at on the c∗ = 0.8 isosurface is shown
in Figure 4.6a. The strain rate has been normalized using the Taylor scale
strain rate u′/λ. It is evident that turbulent straining of the flame surface oc-
curs close to the Taylor length scale λ as suggested previously (Bradley, 1992;
Bray and Cant, 1991). With increasing u′ the Taylor scale strain rate u′/λ
increases, so that the similarity of distributions implies that the prescribed
tangential strain rate at increases in accordance. All the cases exhibit an ex-
(a) Normalized tangential strain rate at (b) Normalized mean curvature hm
Figure 4.6 Distribution of strain at and curvature hm over the peak reaction
rate isosurface c∗ = 0.8 for three cases (lowest u′: solid line, highest u′: dotted
line) spanning the entire dataset.
88
4.2. Stretch rate components
tensional mean surface strain 〈at〉c=c∗ > 0 consistent with theory on turbulent
transport of material surfaces (Pope, 1988) and with previous results on flame
surfaces (Echekki and Chen, 1998; Chakraborty and Cant, 2007). Further, the
distribution shapes are similar to those observed previously in lower intensities
of turbulence (Bray and Cant, 1991) throughout the dataset.
Mean curvature
The variation of local mean curvature hm of the c∗ = 0.8 isosurface is shown in
Figure 4.6b for the same three cases. For consistency with strain rate scaling,
the Taylor length scale λ is chosen here as the scaling parameter. For all cases,
the surface average of mean curvature 〈hmλ〉c=c∗ ≈ 0. With increasing u′, the
likelihood of negative mean curvatures increases (indicated by red arrow). At
the same time, the spread of the curvature distribution also increases. The
shapes of these distributions also conform with previous observations (Bray
and Cant, 1991). The variation of strain (as shown in Figure 4.6a) is contrary
to the variation of curvature (as shown in Figure 4.6b) with increasing u′.
Configurations of strain and curvature
In order to further clarify the relationship between strain at and curvature hm
in the discussion above, it is important to note the configuration of strain and
curvature components observed in turbulent premixed flames. The effect of
strain rate at on surface stretch is direct: extensive strain generates surface area
and compressive strain destroys it. The mechanism by which propagation rate
sdhm contributes to surface stretch, however, is more subtle. The propagation
rate term only contributes if a wrinkle is formed – surface elements with no
mean curvature do not generate surface area due to propagation.
The relationship between curvature and strain is investigated further in Fig-
ure 4.7. The correlation between and hm and at scaled using the Kolmogorov
scale η and Taylor scale λ respectively are shown for low- and high-intensity
flames. The negative correlation at low intensities (Figure 4.7a) is striking.
Previous studies conform with this observation, indicating similarly that at
and hm are negatively correlated up to moderate values of turbulence (Peters
89
4. Damköhler’s First Hypothesis
(a) Case I: u′ = 1.5sL (b) Case V: u′ = 30sL
Figure 4.7 Joint probability distribution of strain rate at (scaled with Taylor
scale strain rate u′/λ) and mean curvature hm (scaled with Kolmogorov scale
η) in the lowest and highest intensity cases.
et al., 1998; Chakraborty and Cant, 2004; Steinberg and Driscoll, 2009). The
high-intensity case (Fig. 4.7b) continues to show a significant negative corre-
lation although the correlation coefficient is lower. Hence, the effects of at and
hm act to oppose each other locally: for hm > 0, strain rate is compressive
(at < 0) and for hm < 0, strain rate is extensive (at > 0). The underlying rea-
sons for the preference of such configurations are being investigated presently.
It is emphasised that the mean curvature hm is not a mean in the sense
of Reynolds averaging. While the mean hm = 1/2(hl + hu) may be zero, the
principal curvatures hl and hu that constitute the mean may not be zero. In
such cases, the flame surface may have a net propagation effect which may
contribute to a net stretch rate. Going back to Figure 3.6a, it is seen that a
finite amount of the surface undergoes these forms of curvature. At present,
we ignore this possible complication in flame surface propagation but it may
be worthy of further investigation.
Displacement speed
In the overall stretch rate s˙, the displacement speed sd acts to regulate the
propagation stretch effects of curvature hm. According to Peters (2000), the
expression for the displacement speed sd of a point x on an isosurface c = c∗
90
4.2. Stretch rate components
is given by
sd(x; c = c
∗) =
[∇ · (ρD∇c) + ω˙
ρ|∇c|
]
c=c∗
, (4.10)
where D is the local molecular diffusivity. The first term ∇ · (ρD∇c) in the
numerator is molecular diffusion and the second term ω˙ is the reaction rate.
Figure 4.8 shows profiles of displacement speed sd at the isosurface c∗ = 0.8
for the same three cases I, III and V spanning the lowest to highest turbu-
lence intensity. In the lowest intensity case (solid line), a relatively narrow
distribution of sd/sL about the expected value corresponding to the density
ratio at this surface is evident. Under higher levels (dashed line) of turbulence
intensity u′ = 10sL, the variance of sd|c∗=0.8 increases due to increased tur-
bulent fluctuations experienced by the reaction layer. In the highest intensity
case (dotted line), the variance spans a significantly larger range denoting the
increased scatter due to turbulent fluctuations. In accordance, negative sd be-
comes increasingly common in higher intensities of turbulence. This is in line
with previous 2D DNS studies (Echekki and Chen, 1996; Peters et al., 1998).
The surface mean 〈sd〉c=c∗ is close to the value expected for a laminar
flamelet at the isosurface c∗ = 0.8 for case I (solid line in Figure 4.8). A
Figure 4.8 Distribution of displacement speed sd over the peak reaction rate
isosurface c∗ = 0.8 in three cases spanning the dataset (solid line: lowest
intensity and dotted line: high-intensity).
91
4. Damköhler’s First Hypothesis
minor reduction of the peak value is noted for case III (dashed line). In the
highest intensity case V, however, a further reduction is observed indicating a
greater role played by turbulent fluctuations. The increasing scatter in sd(c)
in high u′ may be attributed to the diffusive term in Equation 4.10 because
the reaction rate term ω˙ retains to a great extent the laminar flamelet profile
(as discussed in Section 1.2.2). Ongoing investigation is aimed at addressing
these observations.
Configuration of curvature and displacement speed
The net stretch rate contribution due to the propagation term sdhm relies
on both aspects: the creation of mean curvature hm and the propagation of
these curved elements sd. The consequent curvature propagation behaviour
may be revealed in the correlations between quantities hm and sd appearing
in the propagation term sdhm shown in Figure 4.9. It is clear that at low
intensities (Figure 4.9a), sd and hm are negatively correlated; i.e. sd is higher
when mean curvatures are negative hm < 0, and vice versa. In high intensities
(Figure 4.9b), the correlation continues to remain negative, but a considerable
amount of negative displacement speed sd < 0 is also observed. Nevertheless, of
(a) Case I: u′ = 1.5sL (b) Case V: u′ = 30sL
Figure 4.9 Joint probability distribution of sd (scaled with laminar flame speed
sL) and hm (scaled with Kolmogorov scale η) in the lowest and highest intensity
cases.
92
4.2. Stretch rate components
the four possible configurations between ±sd and ±hm, the two configurations
with opposite signs (sd < 0, hm > 0 and vice versa) are most common. Some
sections of the surface with sd > 0 and hm > 0 are observed, but nearly no
sections with sd < 0 and hm < 0 are observed.
Assuming sd has a laminar flamelet behaviour in the mean, the local mean
curvature hm determines the sign of the propagation rate term sdhm: a posi-
tively curved surface element hm > 0 contributes productively Σ (given sd > 0)
and a negatively curved surface element hm < 0 contributes destructively to
Σ.
4.2.2 Balance of strain and propagation
It has been argued by Bradley et al. (1992) that the increase in tangential strain
rate at with u′ causes localised flame extinctions at sufficiently high u′. This
has been attributed historically as the primary cause of the bending effect in
sT (u
′). In particular, increased at reduces sd locally causing the flame surface
area AT to decrease by virtue of holes. The role of curvature hm has been
considered minimal thus far. Informed by the evolution of non-reacting scalar
surfaces, turbulence has been expected to even out flame surface curvatures.
This picture conforms with the observation that 〈hm〉 ≈ 0 for all u′ values
considered in the present dataset.
However, the true effect of curvature on the flame surface area AT (as
computed in Equation 4.5) is manifested through propagation rate sdhm and
not merely curvature hm. Hence, the actual individual contributions of the
stretch rate components at and sdhm are identified (Nivarti and Cant, 2017)
here by calculating volume-averaged quantities
a =
1
V
∫
V
at dV, and h =
1
V
∫
V
sdhm dV, (4.11)
where V is the domain volume. These volume integrals are graphed together
in Figure 4.10 over the entire range of u′ simulated presently.
On the one hand, it is observed that the volume-averaged tangential strain
rate a is positive for all the turbulence intensities considered, and its mag-
93
4. Damköhler’s First Hypothesis
Figure 4.10 Variation of volume-integrated source-terms a and h of the trans-
port equation for Σ′ with increasing normalized turbulence intensity u′/sL.
nitude |a| increases monotonically with u′. Hence, net turbulent strain rate
contributes productively to Σ′ over the entire range of simulations. This con-
forms indeed with the previous studies as well as the figures discussed above.
On the other hand, the volume-averaged propagation rate h contributes de-
structively to Σ′ and its magnitude increases by several orders over the same
range of u′. Therefore, the two source-terms appear to regulate together the
total flame surface area AT . A balance between the two terms determines the
value of this area at a given u′. Since DH1 is satisfied as per Figure 4.4, the
surface area AT governs the turbulent burning velocity sT . Hence, it may be
inferred that the bending effect recorded in sT (u′) is related to this balance of
strain rate at and propagation rate sdhm.
A phenomenological explanation of bending
The results presented above outline a self-regulating mechanism which main-
tains on average the value of Σ in a prescribed turbulent flow field. For a given
u′, production of Σ due to extensive turbulent strain is opposed by destruction
of Σ due to flame propagation in negative curvature regions. The overall effect
of strain and curvature is to regulate Σ. As u′ increases, the magnitude of
94
4.2. Stretch rate components
extensive strain at increases resulting in increased production of Σ. This is
balanced by increased destruction of Σ because larger portions of the brush
have sdhm < 0. The steady-state Σ attained for any given u′ is directly related
to the flame surface area AT and, consequently, determines sT . Hence, the
variation of sT (u′) relies on a competition between strain and curvature – the
source-terms of Σ transport – that determines the steady-state Σ attained.
The conditions of bending may therefore be sought by posing the question:
what are the relationships that determine the Σ attained for a given u′?
The key mechanism of bending seen in Figure 3.9 is highlighted in unison
by Figure 4.4 and Figure 4.10. These figures indicate that turbulent burning
velocity sT is governed by the flame surface area AT (Figure 4.4) but the mono-
tonic increase in extensive strain rate at with u′ does not result in a monotonic
increase in AT due to the destructive role of negative mean curvature (Fig-
ure 4.10). This non-monotonic variation of AT results in the bending of sT (u′)
(seen earlier in Figure 3.9).
For a given geometry L and integral length scale `0, the most significant
change with increasing u′ is the production of an increasing number of small-
scale eddies. It may be concluded, therefore, that these small-scale eddies are
less efficient at generating flame surface area. This reduced efficiency may be
due to the fact that these eddies have shorter lifetimes in comparison to the
overall time scale of flame wrinkling as proposed by Meneveau and Poinsot
(1991). The reduced efficiency may also result from a shift in balance between
the relevant time scales of strain rate as well as of propagation rate. For ex-
ample, the images in Section 3.2.1 show that the shapes and sizes of wrinkles
change with increasing u′. The ‘folds’ and ‘fingers’ seen in high u′ turbulence
act to increase the significance of propagation rate effects over strain rate
effects. Further, small-scale wrinkling causes higher propagation speeds (Fig-
ure 4.8) due to enhanced molecular diffusion between flames in close proximity.
This combination of increased small-scale wrinkling and increased propagation
speed may, therefore, also reduce the time scale of self-intersections.
95
4. Damköhler’s First Hypothesis
4.3 Summary: the role of DH1 in bending
Turbulent surface strain rate at and flame surface propagation rate sdhm have
opposing effects on the surface area AT . While strain tends to produce AT ,
flame propagation tends to destroy AT . With increasing u′, propagation rate
effects dominate and tend to reduce the overall flame surface density Σ at-
tained. An explanation for this imbalance and the resulting decreased Σ in
high u′ turbulence is sought. Two distinct signatures of the governing mecha-
nism are
1. increased negative mean propagation (Figure 4.10) in high u′ turbulence,
2. and increased displacement speed sd in regions of negative curvature
(Figure 4.8 and Figure 4.9) which exaggerates the previous effect.
Hence, a possible mechanism of the bending effect is framed in terms of kine-
matic effects as opposed to chemical effects such as localised extinctions. Since
all parameters but u′ are held fixed, an increased number of small-scale ed-
dies are expected with increasing u′. Further investigation will be focused
on investigating the role of small-scale eddies on flame surface topology and
propagation.
96
Chapter 5
Damköhler’s Second Hypothesis
Damköhler proposed both hypotheses (DH1 and DH2) as primary mechanisms
for sT enhancement in the respective regimes of operation. As the primary
mechanism operates in a given regime, the secondary mechanism was expected
to ‘recede in the background’ contributing alongside albeit to a lesser degree.
This raises an important question: when does one primary mechanism tran-
sition to being a secondary mechanism? Alternatively, is there a limit to the
applicability of a single primary mechanism? Over the years, the range and
upper u′ limit of validity for DH1 has received significant attention, and this
has been the focus of the previous chapters. The validity limit for DH2 remains
unknown.
As discussed in Chapter 1, Damköhler (1940) reasoned that the enhance-
ment of sT in ‘fine-body’ turbulence arises as a consequence of the enhancement
of the local internal transport process. Therefore, DH2 concerns the effects of
turbulence on the diffusive component of the expression sL ∼ A
√
Dω˙. The
local reaction rate ω˙ as well as the global flame surface area AT were assumed
to be unaffected. Reasoning along these lines, the ratio of flame speed in ‘fine-
body’ turbulence to the unstrained laminar flame speeds may be expressed
as
sT
sL
=
√
DT
DL
, (5.1)
where DT and DL represent the transport coefficients in turbulent and laminar
flames respectively. Recent investigations in intense (Poludnenko and Oran,
97
5. Damköhler’s Second Hypothesis
2011; Savard and Blanquart, 2015) and extreme (Wabel et al., 2016) intensities
of turbulence have highlighted the possible significance of diffusive transport
processes in determining sT in such conditions. Appropriate modifications of
DH2 are, therefore, being proposed so as to fit the sT data acquired from
experiments and DNS. Motivated by these investigations, the present work
focuses on outlining the regimes where DH2 is applicable. In this chapter,
the scalar flux transport recorded using the present DNS is analysed and its
implication for the minimal conditions of validity for DH2 is discussed.
5.1 Scalar flux transport in turbulence
In cold non-reacting flows, the ratio DT/DL > 1 because turbulent mixing and
transport of the scalar flux enhances effectively the diffusive transport coeffi-
cient. This enhanced turbulent scalar flux transport leads to increased propa-
gation rates of scalar isosurfaces. Most importantly, the enhancement abides
by the mathematical validity of Equation 5.1 which requires that DT/DL > 0.
Bray-Moss-Libby formulation
In laminar premixed flames, the analogous scalar which designates propagation
of flame isosurfaces is the reaction progress variable c introduced earlier in
Chapter 1. For turbulent premixed flames, the Favre mean progress variable c˜
provides an appropriate framework (Veynante and Vervisch, 2002) to quantify
the mean transport of the propagating isosurfaces. Favre mean quantities
q˜ = ρq/ρ¯ represent density-weighted averages whereas q′′ denotes the respective
fluctuations so that q = q˜+q′′. As introduced in Section 1.2, the Favre average
framework provides an equation for the transport of the scalar c˜ as
∇ · ρ¯u˜c˜ = −∇ · ρ¯u˜′′c′′ + ω˙, (5.2)
according to the Bray-Moss-Libby (BML) formulation due to Bray (1980). In
this framework, the viscous term is neglected under the assumption Re 
Da 1, and the unsteady term is neglected considering a statistically-steady
state. Although these assumptions may not hold strictly in the present dataset,
98
5.1. Scalar flux transport in turbulence
this framework serves to highlight the key physical mechanisms governing the
propagation of the mean flame surface. In the present study, c denotes the
progress variable based on reactant mass-fraction c = (YR,u−YR)/(YR,u−YR,b)
where the subscripts u and b stand for unburned and burned components
respectively. Alternatively, the temperature T may also be used in which case
the diffusivity will correspond to thermal diffusivity.
In the BML framework, flame propagation is linked with the convection
term on the l.h.s. in Equation 5.2 as well as the two terms on the r.h.s. which
denote
1. −∇·ρ¯u˜′′c′′ is the scalar flux term which describes turbulent mixing across
the Favre mean c˜ isosurface, and
2. ω˙ is the reaction source-term which describes the mean chemical trans-
formation of the reactive scalar c˜ from 0→ 1.
Hence, Equation 5.2 extends the mechanics for laminar flame propagation to
the context of mean flame propagation in turbulence. Indeed, this framework
forms the basis for several modelling approaches for turbulent flame propaga-
tion (Veynante et al., 1997; Lipatnikov and Chomiak, 2002).
The discussions so far suggest that the second term ω˙ is not influenced
significantly in the present dataset (see Figure 3.7). The first term, however,
may be expected to change in turbulence if the effects on transport of the scalar
flux u˜′′c′′ are considered. Supporting evidence is provided by the previous
discussions which indicate that the pre-heat zone thickens considerably (see
Figure 4.1). Diffusive transport processes occur in this very region and, hence,
the role of turbulence must be examined in this context.
5.1.1 Gradient transport hypothesis
The enhancement of molecular diffusivity DT due to turbulence may be seen
in the same vein as the diffusion of passive scalars due to molecular transport.
The application of Fick’s law to turbulent transport of passive scalars suggests
that scalar fluxes must be transported favourably down the scalar gradient
99
5. Damköhler’s Second Hypothesis
(Pope, 2001); i.e. from higher to lower concentrations. Such a model of tur-
bulent transport conforms with the qualitative picture where turbulence aids
mixing.
A ‘gradient transport’ expression for the turbulent reactive scalar flux u˜′′c′′
has been hypothesised as
u˜′′c′′ = −DT ∂c˜
∂x
, (5.3)
based on an analogy with Fick’s law. This expression models the transport
of the scalar flux u˜′′c′′ down the mean scalar gradient ∂c˜/∂x for which the
diffusion coefficient DT > DL > 0; i.e. turbulence enhances mixing of different
concentrations of the reactive scalar according to gradient (G) transport.
Effect on flame propagation
In premixed flames, the propagation of the mean flame brush is enhanced by
the mixing of products and reactants . Hence, gradient transport of reactive
scalar is favourable to flame propagation. This is visualized in Figure 5.1.
Since u˜′′c′′ < 0 in gradient transport, it follows that turbulent fluctuations
u′′ that arise from turbulent vortical motion (red arrows in Figure 5.1) sup-
port progress variable fluctuations c′′ convoluting further an initially disturbed
c˜ = 0
Mean Flow
Flame Propagation
c˜ = 1
c00 > 0
c00 < 0
u00 < 0
u00 > 0
Reactants Products
c00 < 0
u00 > 0
Figure 5.1 Schematic showing gradient (G) transport of the progress variable
flux. A disturbed flame surface (light dashed line) is wrinkled further (progres-
sively darker solid curves) by vortical activity. The flame does not counteract
this tendency, and is simply convoluted by turbulent fluctuations.
100
5.2. Counter-gradient transport
flame surface. Hence, the turbulent diffusivity DT > 0 based on gradient hy-
pothesis in Equation 5.3. As a consequence, the ratio DT/DL is applicable
under the square-root operation in the DH2 expression (Equation 5.1).
5.2 Counter-gradient transport
Numerical (Libby and Bray, 1981; Rutland and Cant, 1994) and experimental
studies (Shepherd et al., 1982) suggest that gradient transport is not the only
mechanism of scalar flux transport in premixed flames. In other words, the
turbulent scalar flux u˜′′c′′ (l.h.s. in Equation 5.3) does not always have the
opposite sign (direction) to the scalar gradient ∂c˜/∂x (r.h.s. in Equation 5.3).
In such a case, Equation 5.3 provides DT < 0 posing a threat to the validity
of the DH2 expression (Equation 5.1). Such a mode of scalar flux transport is
referred to as counter-gradient (CG) transport (Libby and Bray, 1981).
Mean Flow
Flame Propagation
Reactants Products
⇢ = ⇢u/⌧⇢ = ⇢u
Figure 5.2 Schematic of the counter-gradient (CG) transport mechanism. At
sufficiently high heat release factors τ = ρu/ρb, the flame has a tendency to
act against gradient transport. An initially wrinkled flame surface (dashed
convoluted curve) may un-wrinkle as the lighter product peninsula is differ-
entially accelerated in the direction of products (along red arrows, away from
the direction of propagation).
101
5. Damköhler’s Second Hypothesis
Mechanism of counter-gradient transport
CG transport has been attributed (Veynante et al., 1997) to a counter-acting
mechanism in premixed flames which acts to inhibit flame surface wrinkling.
This mechanism arises due to gas expansion across the flame, and it competes
with the wrinkling action of turbulent vortical motion (as seen in G transport).
DNS analysis of mean scalar flux transport has highlighted the role of pres-
sure gradients (Rutland and Cant, 1994; Veynante et al., 1997) in the CG/G
transition. The concept that emerges from these studies has been illustrated
in Figure 5.2. Thermal expansion gives rise to a self-induced pressure gradi-
ent which scales with the heat release factor τ ∼ ρu/ρb > 1. Wrinkled flame
surfaces may be characterised by a peninsula of lighter products (as shown in
Figure 5.2); these regions are accelerated differentially in comparison to neigh-
bouring reactant gulfs under the same pressure gradient. Hence, an initially
convoluted flame surface (dark contour) may un-wrinkle if the effect of heat
release is large enough to compensate for the wrinkling effect of turbulence.
Hence, CG transport results from a shift in balance between flame-induced
gas expansion (represented by τ) and flow-imposed turbulent fluctuations (rep-
resented by u′). An appropriate model that captures the accompanying change
in sign of u˜′′c′′ is written (Veynante et al., 1997) as follows
u˜′′c′′ = c˜(1− c˜) (τsL − 2αu′) , (5.4)
wherein the ratio NB = τsL/2αu′ (called the Bray number) dictates the CG/G
transition in scalar flux transport (since c˜(1− c˜) > 0 ∀c˜ ∈ [0, 1]).
In more recent works, the influence of flame-generated vorticity has been
investigated (Louch and Bray, 2001). Wenzel and Peters (2000) have verified
the role of gas expansion using G-equation model simulations of flame prop-
agation with added heat release (discussed in Section 1.2). The simulations
recorded a bending effect in sT (u′) akin to the present dataset. The simulations
of Wenzel and Peters (2000) also recorded a CG/G transition in the scalar flux
transport, but the implication of this transition on Damköhler’s formulations
for sT enhancement was not addressed. In fact, a connection of the CG/G
102
5.2. Counter-gradient transport
transition with the DH2 expression for sT has not been considered to date.
Scalar flux and gradient computations
In Figure 5.3, each panel shows the variation of Favre averaged turbulent
scalar flux u˜′′c′′ with the Favre averaged scalar gradient ∂c˜/∂x for the specified
case with corresponding u′. Each point (∂c˜/∂x, u˜′′c′′) is coloured according
to the corresponding value of c˜, indicating its position in the mean flame
brush 0.01 < c˜ < 0.99. The leading edge c˜ ≈ 0 (cold unburned reactants)
is coloured blue, whereas the trailing edge c˜ ≈ 1 (hot burned products) is
coloured red. In all of the cases, ∂c˜/∂x > 0 throughout the mean turbulent
flame brush, i.e. points lie to the right of the dotted vertical line at ∂c˜/∂x =
0. An exception is the small region close to c˜ ≈ 0.6 (marked pink points
in Figure 5.3c) in case IV. This exception may be attributed to excessive
entrainment of unburned reactants within the preheat zone of the mean flame
brush (visualised below in Figure 5.4). In contrast to the scalar gradient, the
scalar flux u˜′′c′′ transitions from being entirely positive to entirely negative
across the different cases. These four representative cases (I, III, IV and V)
are discussed below to highlight the salient features of this transition.
In case I (Figure 5.3a), u˜′′c′′ > 0 throughout the Favre mean flame brush.
In addition, u˜′′c′′ appears to have a direct correspondence with mean scalar
gradient ∂c˜/∂x, both terms attaining a maximum in the preheat zone between
c˜ ≈ 0.4 and c˜ ≈ 0.6. Similarly, case II (not shown) as well as case III (Fig-
ure 5.3b) exhibit u˜′′c′′ > 0 across most of the brush; i.e. all points lie above
the dotted horizontal line. The correspondence of u˜′′c′′ with ∂c˜/∂x follows
in similar fashion, but with the exception that u˜′′c′′ < 0 close to the leading
edge, exhibiting a local minimum close to c˜ ≈ 0.01 (marked dark blue dots
in Figure 5.3b). Clearly, the upstream turbulence intensity u′ in case IV be-
gins to overturn the counter-gradient (CG) transport u˜′′c′′ > 0 that existed
throughout the flame brush in the low u′ case I. This shift in the behaviour of
u˜′′c′′ in the preheat zone is, therefore, a signature of the imminent transition
to gradient (G) transport at higher intensities.
Upon increasing u′, the magnitudes of u˜′′c′′ and ∂c˜/∂x do not change sig-
103
5. Damköhler’s Second Hypothesis
(a) Case I: u′ = 1.5sL (b) Case III: u′ = 10sL
(c) Case IV: u′ = 20sL (d) Case V: u′ = 30sL
Figure 5.3 Variation of mean scalar flux u˜′′c′′ with the mean scalar gradient
∂c˜/∂x. Colours denote corresponding values of the mean progress variable c˜
(blue: leading edge c˜ = 0.01, red: trailing edge c˜ = 0.99).
nificantly (as seen in the axes of panels in Figure 5.3). Instead, the sign of
u˜′′c′′ and the correspondence between u˜′′c′′ and ∂c˜/∂x is affected considerably.
In case IV (Figure 5.3c), u˜′′c′′ < 0 through most of the flame brush (barring
a small region between c˜ ≈ 0.2 and c˜ ≈ 0.6) whereas in case V (Figure 5.3d),
u˜′′c′′ < 0 across the entire flame brush. The high degree of scatter that ac-
companies this transition marks a departure from a smooth wrinkled flame to
a highly convoluted flame. This follows from the discussion above: in gradient
transport wrinkling is enhanced by turbulence resulting in convoluted flames,
whereas in CG transport wrinkling is diminished resulting in relatively smooth
flames.
104
5.2. Counter-gradient transport
(a) Case I: u′ = 1.5sL (b) Case III: u′ = 10sL
(c) Case IV: u′ = 20sL (d) Case V: u′ = 30sL
Figure 5.4 Instantaneous progress variable contours (black lines) and Favre
velocity fluctuations u′′ (background colours) in an axial cross-section y =
L
2
;L < x < 2L, z of the domain.
Figure 5.4 shows flame surface contours corresponding to the cases I, III, IV
and V within streamwise slices of the domain. The flame is propagating to the
left in all panels. The departure from a smooth wrinkled flame in case I to the
highly convoluted flame in cases IV and V is seen readily. This corresponds to
the contrast in scatter of points in the respective cases noticeable in Figure 5.3.
Such a change in the pattern of wrinkling has also been noted in the CG/G
transition captured in other simulations conducted previously (Veynante et al.,
1997; Chakraborty and Cant, 2009).
Each panel of Figure 5.4 also shows the velocity fluctuations u′′ using
colours (blue: negative, and red: positive) in the backdrop of the contour
105
5. Damköhler’s Second Hypothesis
lines. The decreasing size of turbulence small-scales with increasing u′ may
be observed in the gradually changing distribution of coherently coloured flow
field regions. On closer inspection, it is seen additionally that while positive
velocity fluctuations (red colour) appear predominantly behind the flame in
case I (Figure 5.4a), such fluctuations appear ahead of the flame to a signifi-
cantly higher extent in case V (Figure 5.4d). This transition provides visual
evidence of a shift in balance within the dataset between velocity gradients of
two varieties: 1) those generated behind the flame due to gas expansion, and
2) those present ahead of the flame in the turbulent flow field.
5.2.1 Regime of counter-gradient transport
The role of the balance between the flame imposed gas expansion and the flow
imposed turbulence in the CG/G transition was investigated by Wenzel and
Peters (2000). By using a G-equation model of the flame surface with added
heat release, they were able to compute the two corresponding components
of turbulent velocity fluctuations u′′ = u′′f + u′′g . The shift in favour of the
gas expansion induced fluctuations u′′g as opposed to the flow field fluctuations
u′′f explained the CG/G transition recorded between two separate simulation
cases. In the simulations of Wenzel and Peters (2000), the gas expansion factor
τ was varied. Here the turbulence intensity u′ is varied instead with τ held
constant across the dataset. Hence, u′′g = g(τsL) and its corresponding flux
component u˜′′gc′′ is expected to be constant, while u′′f = f(u′) and u˜′′fc′′ are
expected to change as u′ increases from case I to case V. The shift in balance
between these quantities accompanies the CG/G transition recorded in the
present dataset.
In case I (Figure 5.5a, low u′), u′′ and c′′ are positively correlated through
most of the flame brush. Either u′′ < 0 and c′′ < 0 (seen in the long trail of
blue dots from the leading edge) or u′′ > 0 and c′′ > 0 (seen predominantly
in the big red dots of the reaction zone region). A small number of points
(in the preheat and reaction zones) display a negative correlation, but their
contribution to the mean scalar flux u˜′′c′′ is insignificant as the corresponding
correlations do not carry over to the averaged values (shown in Figure 5.3a).
106
5.2. Counter-gradient transport
(a) Case I: u′ = 1.5sL
(b) Case V: u′ = 30sL
Figure 5.5 Recorded correlations between u′′ and c′′ for case I (above) and case
V (below). Colours indicate the progress variable value (blue: reactants, red:
products) and size is proportional to the reaction rate magnitude. The transi-
tion from positive to negative correlation corresponds to the G/CG transition.
107
5. Damköhler’s Second Hypothesis
Figure 5.6 The Bray number NB for each case simulated in the present dataset
delineating the gradient/counter-gradient transport regimes at unity.
A more complicated correlation is evident in case V (Figure 5.5b, high u′)
which exhibits G transport. Regions with positive correlations continue to be
present, indicating the constant contribution of gas expansion in both cases.
On the other hand, a large portion of the preheat and reaction zone shifts to a
negative correlation, representing the overturning of this contribution by the
turbulent flow field. In fact, the negative correlations are more significant in
the mean (as seen in Figure 5.3d). From the preceding discussion, it is clear
that this qualitative change in the different behaviour of the scalar flux u˜′′c′
occurs at a certain threshold u′ between 10sL and 20sL within the dataset.
Veynante et al. (1997) provided a simple criterion to delineate the CG/G
regimes with contrasting u˜′′c′ behaviour. According to Equation 5.4, u˜′′c′ has
opposite signs across a Bray number of unity, NB ≡ τsL/2αu′ = 1. In the
present context, u′ is the only parameter relevant to the criterion that is being
varied between any two cases. Moreover, the integral length scale `0 is also
held constant (`0/δL = 2.75) across all the cases. The corresponding graph
NB(u
′/sL) is shown in Figure 5.6. With α = 0.25±0.1, the predicted transition
matches (to within ±5sL between 10sL < u′ < 20sL) the transition observed
between cases III (Figure 5.3b)and IV (Figure 5.3c).
108
5.2. Counter-gradient transport
Veynante et al. (1997) hypothesised α to be a function of `0/δL. Ignoring
further the possible dependence on u′/sL, they found that NB = 1 delineates
the CG/G transition in their dataset for α(`0/δL = 11) = 0.5. The present
analysis, with α(`0/δL ≈ 3) = 0.25, indicates that α is a gradually increasing
function of `0/δL. A recent experimental investigation (Kheirkhah and Gülder,
2015) of counter-gradient transport in V-shaped flames suggests a logarithmic
variation of α with `0/δL supporting the assumption of a monotonic variation
considered here. Kheirkhah and Gülder (2015) further report α(`0/δL = 4) ≈
0.3 which is very close to the values obtained here. This value of α, therefore,
provides a valuable datapoint pertinent to subsequent investigations on the
general behaviour of the efficiency function α. It is noted that the CG/G
transition is reported here for a planar flame, and such transitions may occur
locally and under slightly different conditions for more complicated geometries.
For example, Kheirkhah and Gülder (2015) observe a CG/G transition at
u′/sL ≈ 3 which is lower than the conditions u′/sL ≈ 10 observed here. The
shift in balance between τ and u′ in the CG/G transition at NB(α = 0.25) = 1
occurs as u′ is increased with `0 held constant in this dataset. Hence, the
recorded CG/G transition is linked to a transition in the role of turbulence
small-scales.
5.2.2 Computation of the diffusion coefficient
Knowing the scalar flux u˜′′c′′ and the Favre mean progress variable gradient
∂c˜/∂x, the corresponding turbulent diffusivity DT may be computed. This can
be compared with the laminar diffusivity DL to provide the overall enhance-
ment of diffusivity due to turbulence.
The molecular diffusivity DL of a species in laminar flow is calculated as
DL =
λ
ρCpLe
, (5.5)
given its Lewis number (here, Le = 1.0). It is recalled that DH2 refers to the
turbulent amplification of molecular diffusion ‘between the preheat zone and
the reaction zone’. In order to produce the diffusion coefficients for evaluating
109
5. Damköhler’s Second Hypothesis
DH2, values must then be chosen within this (0.6 ≤ c ≤ 0.9) region. In
the present work, DL is calculated as DL = (λ/ρCp)|c=0.8 at the c = 0.8
axial location where the reaction rate attains a maximum in a laminar flame.
Similarly, the turbulent diffusivity DT is evaluated at c˜ ≈ 0.8 in the turbulent
premixed flame brush.
In Damköhler’s original work (Damköhler, 1940), DT was estimated for tur-
bulent flames using Prandtl’s mixing length hypothesis (Prandtl, 1925). Ac-
cording to this hypothesis, the diffusivity may be approximated as the product
of a length scale and a turbulent velocity scale. These scales are chosen respec-
tively as the integral length scale `0 and the turbulence intensity u′ to provide
the expression
DT ∼ u′`0 (5.6)
which is expected to be positive for all turbulent flows because u′, `0 > 0
necessarily. For laminar flames, the relevant scales may be chosen as the flame
thickness δL and flame speed sL respectively to provide DL ∼ δLsL ∼ ν. These
expressions may then be used to compute sT as
sT
sL
∼
√
DT
DL
∼
√
u′`0
sLδL
∼
√
u′`0
ν
(5.7)
where ν ∼ sLδL is based on a constant Prandtl number assumption. In fact,
Damköhler used this very line of reasoning to arrive at sT ∼
√
Re0, the final
term in the above equation, for ‘fine-body turbulence.’ More recently, im-
proved formulations of DT have been used instead of Prandtl’s method using
parameters of the k- model of turbulence (Lipatnikov and Chomiak, 2002;
Savard and Blanquart, 2015).
While these approaches to calculateDT may be adopted here, the focus here
is to understand what it means to use the gradient hypothesis instead. Since
the gradient hypothesis is closely linked with flame propagation behaviour, it
has the ability to shed light on the underlying physical mechanisms and their
limitations. The counter-gradient transport mechanism was one such instance
which exposed a limitation of the gradient hypothesis and consequently sug-
gested the existence of a different mode of flame propagation. Therefore, DT
110
5.2. Counter-gradient transport
Figure 5.7 Computed normalized turbulent burning velocity sT/sL compared
with predictions
√
DT/DL of DH2 in the gradient diffusion regime NB < 1.
DH2 is invalidated in the CG transport regime denoted by NB > 1 (shaded
red).
is computed here using the gradient hypothesis expression of Equation 5.3 as
DT = − u˜
′′c′′
∂c˜/∂x
∣∣∣∣
c˜=0.8
. (5.8)
at the Favre mean reaction zone c˜ = 0.8. In the presence of a CG/G transition,
DT is expected to change in sign. Given the formulation of DH2 as presented
in Equation 5.1, the applicability of a negative DT is prohibited.
DH2 states that the role of ‘fine-body turbulence of size `  δL’ is to en-
hance sT ‘by amplifying the microscopic (diffusive) transport processes’ within
the flame brush. As discussed above, u˜′′c′′ < 0 only in a small portion of the
dataset corresponding to the cases IV and V (as shown in Figure 5.6). Numer-
ical evaluation of Equation 5.1 for these cases is shown in Figure 5.7. For the
remaining cases, NB > 1 and the diffusion coefficient DT < 0. Hence, Equa-
tion 5.1 cannot be evaluated. In other words, Damköhler’s second hypothesis
(DH2) applies strictly to the gradient transport regime where NB < 1. This
is illustrated in Figure 5.7 by marking the region of inapplicability NB > 1 in
red.
111
5. Damköhler’s Second Hypothesis
Consequently, it appears that the general claim that ‘turbulence always
acts to enhance the flame speed’ might only apply in certain special limiting
conditions. In fact, for NB > 1, the turbulence scales may have an inhibitive
effect on sT because u˜′′c′′ > 0 is transported against the gradient −∂c˜/∂x
favourable to flame propagation. The difference between values of
√
DT/DL
and sT/sL in Figure 5.7 opens a further question: to what extent do DH1 and
DH2 influence sT individually in this region of overlapping validity? As yet, the
evidence is inconclusive and further investigation is necessary to understand
fully the significance of DH1 and DH2 in the thin reaction zones regime. An
investigation of higher u′ in a regime corresponding precisely to DH2, i.e. with
`  δL, will help understand when it is that the diffusive processes begin to
govern decisively the enhancement in sT . The present analysis elucidates rather
the minimal criteria for the applicability of Damköhler’s second hypothesis.
5.3 Summary: validity regime of DH2
In low u′ turbulence, DH1 applies and turbulent motion both roughens and
stabilizes the surface. These simultaneous effects result in a positive correlation
between u′′ and c′′, i.e. CG transport prevails. Referring to Equation 5.3, this
implies that the l.h.s. u˜′′c′′ > 0 whereas the r.h.s. ∂c˜/∂x > 0. Hence, since
DT < 0, it is clear that DH2 cannot operate.
In high u′ turbulence, DH1 continues to apply as the large-scales serve to
deform the surface. In this regime, however, G transport (u˜′′c′′ < 0) prevails.
Hence, DH2 is applicable and operates alongside DH1. The transition between
CG/G transport, therefore, delineates the regimes of applicability of DH2.
Hence, contrary to Damköhler’s original conclusion, it appears that ‘fine-
body turbulence with ` δL’ may not always act to enhance the propagation
speed sT by ‘amplifying microscopic transport processes.’ In fact, this is plau-
sible solely beyond the threshold u′ (i.e. for NB < 1) in a regime where G
transport prevails.
Two aspects solicit further investigation: 1) the nature (positive or nega-
tive) of effect of ‘fine-body’ turbulence ` δL on sT below the threshold u′ (i.e.
112
5.3. Summary: validity regime of DH2
for NB > 1) and its associated expression, and 2) the underlying causes that
govern the shift between different effects of turbulence on flame propagation.
113

Chapter 6
Towards Validation with the Leeds
Experiment
In the previous chapters, a mechanism of the bending effect has been proposed
based on the validity of DH1. Further, a recorded transition in the scalar
flux transport has revealed a limitation of DH2. These investigations suggest
that flame propagation can be limited by physical mechanisms inherent to the
kinematics of propagation. Experimental investigations of the Leeds explosion
vessel have suggested (Bradley et al., 1992), however, that thermo-chemical
processes such as flame quenching constitute the primary mechanism of bend-
ing. Hence, in order to validate these two contrasting explanations of bending,
a comparative study has been initiated with Leeds group.
In the following, the Leeds experiment and its recent 3D reconstruction ca-
pability are described along with preliminary comparisons with the parametric
simulations conducted so far. Subsequently, the host of modifications made
to the DNS method for carrying out the present work are elaborated. These
modifications have resulted in the design of a computationally-efficient DNS.
The simulation results will soon be compared with experimental results.
6.1 Leeds: expanding spherical explosions
As described in Section 1.1.4, the Leeds explosion vessel provided 2D flame
surface measurements of an outward expanding spherical explosion previously
115
6. Towards Validation with the Leeds Experiment
(Bradley et al., 2011). Ongoing PhD work conducted by Ben Thorne has fo-
cused on conducting experimental measurements leading to 3D reconstructions
of the flame surface (Nivarti et al., 2017). The key features of this new setup
are described in the following passages.
Experimental setup
The Leeds explosion vessel provides a constant pressure environment for the
initiation and outward propagation of spherical flame kernel in a turbulent
premixed charge. With the overall construction of the vessel largely unchanged
(Gillespie et al., 2000), the present setup allows a spark to be generated at the
centre using a high-intensity laser issue through one of the quartz faces (see
Figure 1.4). The fans rotate at a constant frequency during the ensuing flame
propagation, so that the turbulence intensity u′ and the integral length scale
`0 are maintained throughout the experiment.
A laser sheet can be issued through one of the quartz windows. This
sheet excites unburned seed particles in the mixture so that an instantaneous
approximate flame surface contour may be located in the vessel (Gillespie et al.,
2000). The laser sheet can also be made to swing across the explosion vessel
by reflecting it on a rotating mirror. As the mirror rotates, the pulsating
laser sheet captures several discrete slices across the vessel volume. A set of
2D flame surface contours is obtained as a result. Provided the laser sheet
sweeps the vessel fast enough in comparison to the propagation speed of the
flame, these contours may be assumed to correspond to an instantaneous flame
surface; subsequently, the individual contours are spliced together to form a
3D reconstruction of the flame surface.
Given the constraints of the laser pulse frequency and the angular velocity
of the rotating mirror, it is noted that stoichiometric flames propagate too
fast for this assumption to remain valid. Hence, lean methane-air flames with
φ = 0.6 are considered which have a relatively lower propagation velocity even
at high intensities of turbulence. This chemical composition was, therefore,
agreed upon by the two groups as the basis of the comparative study. Some
comparisons were made, nevertheless, using the previous stoichiometric flame
116
6.1. Leeds: expanding spherical explosions
dataset obtained in Chapter 3. These preliminary comparisons have been
discussed below.
6.1.1 Preliminary comparisons
Snapshots of the reconstructed flame surface are compared with the recorded
data obtained in the parametric DNS. The case u′ = 10sL is chosen as the
point of comparison as shown in Figure 6.1a – surfaces obtained from the
Leeds experiment and Cambridge DNS (at c = 0.5 at t = 4τ0) are illustrated
(a) Chosen intensity of comparison (b) Mean curvature distributions
(c) Reconstructed spherical flame surface (d) Isosurface obtained using DNS
Figure 6.1 Preliminary comparisons of methane-air flame propagation in exper-
iments (Leeds) with φ = 0.6, and DNS (Cambridge) with φ = 1.0. The same
turbulence intensity u′/sL is prescribed (panel a) in both cases. Curvature
distributions are normalized using the Taylor length scale λ.
117
6. Towards Validation with the Leeds Experiment
alongside. In Figure 6.1b, these surfaces are compared quantitatively using the
curvature distribution hm obtained upon removing the average mean curvature
component 〈hm〉. Indeed hm is found to have a very similar distribution on both
flame surfaces. This indicates that the effect of various scales of turbulence on
a given flame surface is similar for identical intensities. Yet it is evident that
the spherical reconstruction has a wider spread of curvatures around the mean
in comparison to the planar DNS isosurface. Such differences may be affected
once the equivalence ratios of the two flames considered are matched.
In Figure 6.1c and Figure 6.1d the flame isosurfaces are visualized in order
to aid a qualitative comparison. Local mean curvature hm normalised with
the Taylor length scale is shown using colours (red: positive, and blue: nega-
tive) on the isosurface. Significant similarities are observed in the two surfaces.
Indeed, the scale of wrinkling appears to be the same for both surfaces. How-
ever, the planar isosurface tends to have locally cylindrical shapes whereas
the spherical kernel exhibits more spherical shapes. Improved comparisons are
sought between the experiment and DNS to elaborate on these differences. In
the following, a DNS case is prepared so as to make suitable comparisons with
experimental measurements in Leeds.
6.2 Cambridge: freely propagating flames
The following DNS setup is based on the method discussed in Chapter 2 with
a range of modifications made to address better the experimental results.
6.2.1 Multi-step chemical kinetics
As discussed above, lean premixed turbulent flames display several phenomena
that cannot be modelled adequately by the single-step chemical kinetic model
used in the previous chapters. For instance, effects of non-unity species Lewis
numbers are expected to be prominent in lean premixed flames. These effects
may result in the modification of local mean curvatures on the leading and
trailing edges. Moreover, effects of differential diffusion of various chemical
species may alter flame propagation. Finally, the different reaction pathways
118
6.2. Cambridge: freely propagating flames
of these species may influence further the overall turbulent flame propagation.
In order to record the experimentally observed dynamics due to such ef-
fects in lean premixed flames, a 25-step reduced chemistry model by Smooke
and Giovangigli (1991) has been used here in place of the single-step method
adopted previously. The key implication of this change is that two terms
of the species transport Equation 2.25 are now modelled differently. These
are the reaction rate ω˙α and the diffusive mass flux ρVα,kYα. In the previ-
ous chapters, ω˙α was modelled using a single-step Arrenhius formulation (see
Equation 2.5). The latter term ρVα,kYα was modelled using a constant Lewis
number assumption within the framework of Fick’s law (see Equation 2.9).
The new formulations for these terms are described below.
Reaction Rate
In an M -step reduced chemical kinetic description the overall combustion re-
action (Equation 2.3) is modelled using a series of elementary reactions. A
single elementary reaction may be written
N∑
α=1
ν ′α,mMα −−→
N∑
α=1
ν ′′α,mMα, (6.1)
where m = 1, 2, . . .M the index of the reaction. Stoichiometric coefficients of
the various species α are denoted ν ′α,m in the reactant side and ν ′′α,m in the
product side. The rate of this mth elementary reaction is given then by a
generalised Arrenhius expression
ω˙m = km
N∏
α=1
[
ρYα
Wα
]ν′α,m
, (6.2)
where
km = AmT
nme−
Ea,m
R0T (6.3)
and the constants Am, nm, and Ea,m prescribed a priori for each elementary
reaction m. It follows that the rate of formation of each chemical species α is
119
6. Towards Validation with the Leeds Experiment
expressed as
ω˙α = −Wα
M∑
m=1
(
ν ′′α,m − ν ′α,m
)
ω˙m, (6.4)
which is the species-weighted sum of the rates of reaction ω˙m over all M steps.
It is seen therefore, that the fuel consumption rate ω˙F for instance would de-
pend on each elementary reaction step it participates in within this description.
As a consequence, the overall rate ω˙F , which dictates the value of sT , will de-
pend on the availability of reacting species participating in those elementary
steps. Hence, the M -step reduced chemical mechanism provides a more de-
tailed consideration of the dependence of sT on local mixture composition.
In this manner it may provide a more realistic description of the combustion
process and its burning velocity.
Diffusive Mass Fluxes
In order to compute the diffusive mass flux ρVα,kYα for species α, a mixture-
averaged transport formulation in the absence of thermal diffusion flux (Soret
effect) and diffusion thermal flux (Dufour effect) is considered here. The molec-
ular transport coefficients for the mixture are assembled from molecular trans-
port data for each individual species. This approach is more accurate, and
also allows for greater generality in the treatment of molecular transport phe-
nomena. Hence, for the diffusive mass flux, Fick’s Law is replaced by
ρVα,kYα = −ρDα∂Yα
∂xk
+ ρDαYα
1
Wm
∂Wm
∂xk
+ ρDαYα
(
1− Wα
Wm
)
1
P
∂P
∂xk
(6.5)
where Dα is now the mixture-averaged diffusion coefficient for species α and
Wm is the effective molar mass for the mixture given by
1
Wm
=
N∑
α=1
Yα
Wα
. (6.6)
For completeness it should be noted that the viscosity is also now a mixture-
averaged quantity. Values of viscosity µα and thermal conductivity λα for each
species, together with values of the binary diffusion coefficients Dαβ for each
pair of species, are obtained using a set of polynomial fits (Cant, 2013).
120
6.2. Cambridge: freely propagating flames
(a) Temperature, major components (b) Oxygen, radicals
Figure 6.2 Profiles of thermo-chemical variables: temperature and chemical
components (left) of the composition vector as well as radicals (right) computed
using a 25-step chemical kinetic model.
In the 25-step chemical kinetic model chosen here (Smooke and Giovangigli,
1991), the composition vector consists of 16 chemical species. This is given by
[CH4,O2,CO2,H2O,H2,O,OH,H,HO2,H2O2,CO,CH2O,HCO,CH3,CH3O,N2].
The computation of these species including formaldehyde CH2O and hydroxyl
radical OH allows for a head-to-head comparison with experimental measure-
ments. The detailed chemical reaction rate as well as diffusion flux description
in this framework allows for a more accurate simulation of the pre-heat as well
as reaction zone structure. The chemical mechanism used here (Smooke and
Giovangigli, 1991) has, in fact, been tailored for lean flames and excludes C2
chemistry. Figure 6.2 shows the profiles of thermo-chemical variables computed
based on this chemical mechanism.
Chemical Resolution Constraints
Initial thermo-chemical field (density ρ, temperature T , and species mass-
fractions Yα) in the turbulent flame simulation is specified with an unstrained
1D laminar flame solution profile calculated a priori. In the present case, the
equivalence ratio φ = 0.6 and inlet unburned mixture conditions p = 1 bar
and T = 300K. The laminar flame speed sL = 0.122ms−1 and thickness
δL = 0.928mm are obtained a priori from a 1D simulation with a fine mesh
121
6. Towards Validation with the Leeds Experiment
N = 1024 within a domain of length L = 1 cm. Such a fine resolution is in
principle unnecessary for the turbulent flame simulation and a coarser grid
may be used.
In order to find the coarsest grid that can stably simulate the chemical
structure, a convergence study is conducted using 1D laminar flame simulations
in this domain of length L = 1 cm. Beginning with an extremely fine grid N =
1024 that corresponds to a thermal flame thickness resolution δfineL /∆x ≈ 100,
the grid is coarsened by factors of 2 until the solution becomes unstable due to
lack of adequate resolution. Close to the point of stability-loss, grids are refined
gradually until the resolution requirement is narrowed to within δfineL /∆xc±1.
Figure 6.3 shows the number of steps required Nc to obtain the converged
value of laminar flame speed sfineL = 0.122ms
−1 for each resolution level which
is quantified by δfineL /∆x for a fixed domain length. The resolution is found to
be inadequate at δfineL /∆x ≈ 22.3 where no convergence is achieved. Hence,
this value is taken to be a resolution constraint of the present chemical kinetic
model. This constraint is also used later for estimating the optimal domain
size for the DNS.
Figure 6.3 Grid resolution study to determine the spatial resolution require-
ment for the 25-step methane-air chemical kinetic model. Convergence fails
for grids with ∆x & δL/22.3, where δL = 9.3mm denotes the value obtained a
priori using a fine grid.
122
6.2. Cambridge: freely propagating flames
6.2.2 Fast adaptive timestep control
Numerical solution of 16 additional transport equations imposes a severe com-
putational expense on the simulation. It is critical, therefore, that the desired
statistically-stationary solution state be achieved in as few timesteps as possi-
ble.
A direct approach to reducing the number of timesteps would be to use large
timesteps. The magnitude of timestep ∆t is restricted, however, by the sta-
bility range of the numerical method. The domain of stability is characterised
by a range of the Courant-Friedrichs-Lewy (CFL) number CFL ≡ U∆t/∆x
specific to the numerical integration scheme, where U is a characteristic ve-
locity such as the sound speed a. In implicit integration methods, the range
of CFL stability is much higher in comparison to the explicit method used
presently. However, the additional stability comes with added computational
expense and complexity involved in the matrix inversion that is required in
implicit integration.
In explicit methods, the timestep size ∆t(n) may be calculated adaptively
based on a pre-specified stable CFL number for a uniform grid ∆x. The adap-
tive timestep control strategy introduced in Section 2.3 presents yet another
approach to computing timestep sizes adaptively. Rather than using a fixed
CFL number, a fixed error tolerance ̂ is specified. Thus, appropriate timestep
sizes ∆t(n) are computed based on PID control that maintains this error toler-
ance ̂. This precludes a priori knowledge of stability regions of the numerical
scheme.
Variation of the timestep size ∆t(n) computed using PID control (Equa-
tion 2.49) is shown in Figure 6.4 using black lines. The profiles of timestep
history in different simulations (shown in Figure 6.4a, Figure 6.4b and Fig-
ure 6.4c) are denoted by the PID exponents α, β and γ. These profiles il-
lustrate a general trend; after an initial period of small step sizes the PID
controller reaches a steady-state where the step size fluctuates about a fixed
value. In a 1D laminar flame solution (black line in Figure 6.4a), this cor-
responds to CFL ≈ 1.2, whereas in a periodic HIT simulation (black line in
Figure 6.4b) this timestep corresponds to CFL ≈ 0.2. Furthermore, a much
123
6. Towards Validation with the Leeds Experiment
lower CFL ≈ 0.01 is achieved in a fully-resolved turbulent flame propagation
case using PID control.
This highlights that the numerical solutions so far have been performing
inadequately due to an inefficient adaptive strategy of selecting timestep sizes.
The timestep sizes ∆t ∼ 10−10 s are rather small and render the overall solution
to be computationally infeasible. Therefore, the control strategy is modified
in order to improve timestep sizes. This is done by investigating three aspects:
limiters, cutoffs and exponents.
Limiters and cutoffs
In the PID controller used so far, timestep selection has been constrained
by imposing a limiter that restricts the maximum ratio of step size σn =
∆t(n+1)/∆t(n) provided by the controller to a small value σ̂n. In other words,
the selected timestep is given by ∆t(n+1) = min(σ̂n, σn)∆t(n). Limiters of this
kind avoid the selection of relatively large timsteps so that large errors may
be avoided at subsequent steps of the numerical solution. However, this turns
out to be a rather severe restriction. The timestep profiles discussed above
each indicate a period of slow growth to a steady timestep size (black lines
in Figure 6.4). Such a phase of relatively-small timesteps is neither necessary
numerically.
The reason for the gradual growth of ∆t(n) is evident in Figure 6.4a.
Timestep size ratios σn ≡ ∆t(n+1)/∆t(n) provided by the controller (black
dots) are restricted to σ̂ = 1.1 even though the corresponding error ratios
r̂n/̂ are small (x axis in Figure 6.4d). The error ratio shown here denotes the
ratio of the maximum L2 norm of the solution r̂n to the pre-scpecified error
tolerance ̂. In contrast, the unlimited profile with σ̂ = ∞ shown alongside
exhibits the same errors but with much larger step size ratios (red dots in Fig-
ure 6.4d). This suggests that the use of limiters may be abandoned entirely,
saving essential computational time.
Despite the performance of the unlimited scheme (red dots in Figure 6.4d),
a smooth limiter proposed by Söderlind and Wang (2006) has been imple-
mented presently. Limiting of this sort allows large timestep ratios to be
124
6.2. Cambridge: freely propagating flames
(a) 1D laminar flame (b) Periodic HIT
(c) 3D propagating turbulent flame (d) Timestep ratio and error
Figure 6.4 Timestep history (panels a,b, and c) as well as timestep ratios (panel
d) obtained using the PID controller αβγ and the modified PI4.2 controller
α∗β∗ due to Söderlind (2002) for various test problems.
accepted in a safe manner, avoiding disastrous consequences at later steps in
the numerical integration. In this approach, the timestep ratio obtained from
the controller is limited as
σ̂n = 1 + LIM · tan−1
(
σn − 1
LIM
)
, (6.7)
where LIM is a pre-specified constant value. The formulation above restricts
timestep ratios σn to within a factor of the coefficient LIM. This means that
the timestep ratio provided by the controller is bounded within 1/LIM < σn <
LIM, but in a smooth manner.
125
6. Towards Validation with the Leeds Experiment
Similarly, cutoffs are used normally to impose a restriction on the timestep
selection whenever the error ratio r̂n/̂ > 1 (as indicated by the vertical line
in Figure 6.4d). Since large errors are not usually incurred in the present
framework, cutoffs have not been investigated further.
Tuning exponents
Choice of controller exponents is a central issue in dynamic control of timestep
selection. As part of the investigation of small step sizes provided by the
previous PID controller, it was concluded that the tuning exponents were in-
sufficiently tuned. For instance, the variation of a single exponent, say γ,
would result in significant improvements of the timestep the controller con-
verges upon. Hence, exponents suggested by Söderlind and Wang (2006) have
been used in the present simulation. In fact, these correspond to the PI4.2
controller where the exponent γ = 0 (Söderlind, 2002), i.e.
∆t(n+1) = TOL ·
(
ˆ
rˆn
)α∗(
ˆ
rˆn−1
)β∗
∆t(n) (6.8)
where α∗ = 0.6 and β∗ = −0.2 as discussed by Söderlind and Wang (2006),
The error tolerance is shifted to a smaller value than before ̂ = 10−5 and the
safety factor is also constrained further TOL = 0.8. So as to make obvious the
change between the previous and present controller, this is referred to as PI
controller as opposed to the PID controller.
The timestep profiles obtained in this manner are presented alongside the
profiles shown earlier in Figure 6.4 using red lines. In all profiles – the 1D
laminar solution (Figure 6.4a), the 3D HIT solution (Figure 6.4b) and the
3D flame propagation solution (Figure 6.4c) – the timesteps rise immediately
to the value converged upon due to the absence of limiting. In other words,
there is no period of gradual growth beginning with the initially prescribed
value as seen in the black lines. Moreover, the profiles are all much smoother
than the corresponding profiles obtained using the PID controller with strict
limiting (black lines). The smooth profile may also be attributed to the loss
of instabilities associated with the γ term.
126
6.2. Cambridge: freely propagating flames
In Figure 6.4, it is also evident that the actual magnitude of the timestep
converged upon by the PI controller is higher than the PID controller for
most cases. For instance, apart from the 1D laminar flame solution case (Fig-
ure 6.4a), the timestep size improves nearly threefold for the HIT simulation
(Figure 6.4b) and almost by two orders of magnitude in the turbulent flame
propagation case (Figure 6.4c). Due to these merits of the PI controller, it
is implemented in the present simulations. Further improvements are sought,
however, including the ability to reject solutions that incur too much error.
These are related to the cutoffs neglected above.
6.2.3 Linear turbulence forcing
In the Leeds setup, the intensity u′0 ≡ u′|t=0 of homogeneous isotropic turbu-
lence (HIT) is maintained throughout flame propagation by keeping the fans
on. In the DNS simulations conducted so far, the initial u′0 was subject to de-
cay. Furthermore, the rate of decay is a function of velocity gradients, so that
the high-intensity cases are prone to faster turbulent decay. This is undesir-
able; the altered intensity u′ may change the expected regime of combustion.
Investigating this effect further, the source of turbulent decay is located in
the viscous term in the momentum equation
∂ρui
∂t
+
∂
∂xj
(ρujui) = − ∂
∂xi
(p) +
∂
∂xj
(τji) + fb, (6.9)
recalled here from Chapter 2, with the exception that a body-force fb is also
considered here. In the absence of body-forces fb = 0, as in the simulations so
far, the HIT simulation exhibits an exponential temporal decay of turbulent
kinetic energy ∂k/∂t at the rate of turbulent dissipation
ε = ν 〈∇u : ∇u〉 , (6.10)
as shown in Figure 2.4. So as to mitigate the extent of turbulent decay, the
flame simulations so far employed an inflow BC which convected downstream
a frozen turbulence state with intensity u′0. Nevertheless, the convected tur-
bulence undergoes decay from the inlet to the flame leading edge (as reported
127
6. Towards Validation with the Leeds Experiment
in Table 3.2). In other words, despite the inflow BC, the leading edge inten-
sity u′le 6= u′0. The relatively small domain reduces the scope of this decay,
but the domain size warranted by the present comparisons do not allow this
simplification. Hence, a strategy of turbulent forcing is adopted here.
In an HIT simulation based on the solution of Equation 6.9, the average
rate of change of turbulent kinetic energy is given by
∂ρ〈u · u〉
∂t
= −ρε+ 〈ρu · f b〉, (6.11)
in the presence of a body-force fb. Evidently, the body-work term 〈ρu ·fb〉 acts
in addition to dissipation ρε. The strategy for turbulence forcing proposed by
Lundgren (2003) exploits this aspect. An ‘artificial’ body-force is added in such
a manner that the decay due to turbulent dissipation is balanced. Evidently,
a body-force fb = Q(u − 〈u〉) proportional to the local turbulent fluctuation
results in the corresponding work term of form Q〈u · u〉 provided 〈u〉 = 0.
Hence, the coefficient Q may then be chosen such that turbulent decay due to
dissipation is balanced exactly by the work done due to forcing leading to
0 = −ρε+Q〈ρu · u〉, (6.12)
which provides Q = ε/〈u ·u〉 = ε/2k for incompressible flows where the turbu-
lent kinetic energy k ≡ 〈u · u〉/2. In this manner, Lungren’s forcing strategy
(Lundgren, 2003) may be used to maintain the initial turbulence intensity u′0
throughout the simulation.
This approach to turbulent forcing has been found satisfactory in many
simulations (Rosales and Meneveau, 2005; Ishihara et al., 2009). In order to
extend its applicability further, Petersen and Livescu (2010) have proposed a
modification to address compressible statistically-stationary turbulent flows.
Such a modification has been found challenging to implement as well as un-
necessary here due to the low Mach numbers involved. Carroll and Blanquart
(2013), however, investigated strategies to improve the rate of convergence to-
wards statistical stationarity. A modification of the linear forcing term was
proposed as fb = Qk0/ku where k0 is the initial kinetic energy. The outcome
128
6.2. Cambridge: freely propagating flames
(a) Kinetic energy in HIT simulation (b) Forcing in 3D flame simulation
Figure 6.5 Turbulent forcing strategy. Left: time evolution of mean turbulent
kinetic energy k(t) for an HIT simulation with Q = 0 and with Q = ε/(2k0).
Right: region of turbulent forcing in the flame simulation.
of this modification is that the value Q = ε/2k0 thus obtained from Equa-
tion 6.11 is constant. In Figure 6.5a, results from using this forcing term are
compared with the the decaying HIT simulations conducted previously. It is
clear that the forced simulation (red dots) undergoes little to no fluctuation,
whereas the decaying simulation (black dots) loses 90% of the energy in the
same time.
Lundgren’s strategy (Lundgren, 2003) along with the modification proposed
by Carroll and Blanquart (2013) is followed in the present simulations. This
contrasts with some other forcing strategies developed previously. For instance,
Eswaran and Pope (1988) follow a ‘naturally expected’ forcing condition in tur-
bulent flows. Energy is injected only at low wavenumbers in spectral space.
The Navier-Stokes equations are allowed to re-distribute subsequently the en-
ergy over smaller scales. Such a strategy, however accurate, is not adopted here
for the sake of convenience. Nevertheless, it would be interesting to explore
these alternatives in future.
Due to the merits of Lundgren’s forcing strategy (Lundgren, 2003), it has
been employed recently in several DNS studies. For instance, Savard and
Blanquart (2015) uses this strategy to force the turbulent flow in their sim-
ulation of high-intensity turbulent flame propagation. The forcing term is
129
6. Towards Validation with the Leeds Experiment
effected throughout the domain of flame propagation. Given that the forcing
strategy is valid only for incompressible HIT simulations, the validity of such
approaches may be questioned. Furthermore, concerns (Quadrio et al., 2016)
have arisen recently regarding the effect turbulence forcing may have on the
statistics collected. The present simulations, therefore, take due consideration
of this, restricting the forcing to a small region of the domain upstream of the
flame. The location of forcing is shown in Figure 6.5b. Both the inlet as well as
the flame leading edge are avoided by using a multiplier coefficient ζ∗ ∈ [0, 1]
in fb = ζ∗Qk0/ku such that ζ∗ 6= 0 only in the specified region. With the
exception of a mean axial velocity, this region of forcing indeed resembles an
HIT simulation case. This approach is expected to mimic the effect of fans in
the Leeds experiment.
6.2.4 Computational setup and regime
Domain size optimisation
Small-scale turbulence-flame interactions are critical to the present compari-
son. Yet the role of large-scales must be simulated at the same time for at
least two reasons. Firstly, the simulated energy spectrum represents the cas-
cade process realistically if a large separation between integral and Kolmogorov
length scales is considered (Pope, 2001). Secondly, the problem of flame prop-
agation in the Leeds vessel corresponds to a regime between coarse-body and
fine-body turbulence. Hence, it is required that `0/δL be as large as possible
while η/δL < 1. Therefore, the quality of comparison with experiment is a
function of the computational domain size.
It is computationally not feasible, however, to match the experiment ex-
actly in size. Hence, an attempt is made to maximise the possible size of the
domain given the constraints of small-scale resolution. This may be framed as
an optimisation problem with the following constraints
1. Resolution of chemistry: the maximum spatial step size that can stably
simulate the 25-step reduced chemical kinetic model is ∆xc ≈ δL/23 (as
obtained in Figure 6.3).
130
6.2. Cambridge: freely propagating flames
2. Resolution of turbulence: the maximum spatial step size that can faith-
fully resolve small-scale turbulent motion is ∆xf = η/2.0 where η =
f(L, u′) (as provided in Equation 3.5). This is directly proportional to
L and inversely proportional to u′.
3. Expense of computation: the number of points N used to resolve the
above criteria along the smallest dimension L determines the overall
computational expense. For an inflow-outflow domain with a relatively
longer axial length sL, the total computational expense varies as sN3.
The variation of these constraints with domain size is shown in Figure 6.6.
Adequate resolution of chemistry imposes a constant constraint on the mesh
size (dashed blue line) regardless of the domain size. Resolution of turbulence
Kolmogorov scale η, however, imposes a monotonically diminishing constraint
(dashed red line). In other words, a larger mesh size ∆x may be used to resolve
η expected in a domain of larger width L.
The variation of computational expense N3 corresponding to each of these
constraints N = L/∆x is shown along side in Figure 6.6 using solid lines. The
step size imposed by turbulence ∆xf = η/2 increases as a power α = 1/4 of
Figure 6.6 Resolution constraints of chemistry (black dashed lines) and tur-
bulence (red dashed line) as well as the corresponding computational expense
(black and red solid lines, respectively) as a function of DNS domain size.
131
6. Towards Validation with the Leeds Experiment
Lα. Hence, the related computational expense (solid red line) varies as L1−α.
The step size imposed by chemistry ∆xc = δL/23, however, remains fixed
for all L. The number of points required due to this constraint (solid black
line) increases, therefore, in direct proportion to L. As a consequence of these
varying rates of growth in computational expense, the chemistry constraint
becomes more expensive at a critical L (red square in Figure 6.6). In fact, this
critical value L ≈ 0.015m corresponds to the largest configuration for which
N is determined by the flow resolution constraint ∆xf , and at the same time
satisfies the chemical constraint. Hence, this is the domain size chosen in the
present simulation.
In principle, this algorithm can be generalised, allowing the a priori choice
of a few essential physical parameters most relevant to the simulation, such
as the regime of combustion. The task of the algorithm would then be to
provide the computational parameters such as domain size and integral length
scale `0. Such algorithms are being investigated presently in order to provide
computationally efficient design of DNS studies.
Combustion regime
In order to set abreast the Leeds estimates (Bradley et al., 2013; Bagdanavicius
et al., 2015) of the combustion regime, a Karlovitz number is used as follows
K =
1
4
(
u′
sL
)2
Re−
1
2 (6.13)
where the Reynolds number Re ≡ `0u′/ν. The use of this form of Ka is mo-
tivated further by its historical use in previous experimental studies (Bradley
et al., 2011). The formulation of Equation 6.13 would predict localized flame
extinctions to expected to occur for K ≥ 1.0 in unity Lewis number mixtures.
In order to compute the normalised Reynolds number, the kinematic viscosity
ν is computed as ν ∼ δLsL, maintaining consistency with experiments (Bag-
danavicius et al., 2015).
A regime diagram is presented in Figure 6.7 using the new Karlovitz number
definition. The value K = 1.0 is taken to separate the flamelet and distributed
132
6.2. Cambridge: freely propagating flames
Figure 6.7 Combustion regimes corresponding to the Leeds experiment and
the Cambridge DNS in the regime diagram with a modified definition of Ka.
Table 6.1 Input parameters for the Leeds-Cambridge study.
Parameter Cambridge Leeds
Ambient pressure, Pa 1 bar 1bar
Inlet temperature, Tu 300K 300K
Kinematic viscosity, ν 1.585× 10−5 m2/s 1.62× 10−5 m2/s
Equivalence ratio, φ 0.6 0.6
Laminar flame speed, sL 0.122 ms−1 0.119 ms−1
Laminar flame thickness, δL 1.29× 10−3 m 1.36× 10−3 m
Integral length scale, `0 0.002m 0.02m
Turbulence intensity, u′ 1.25m/s 1.25m/s
Karlovitz number, K 2.1 0.7
regimes of combustion. The normalised turbulence intensity u′/sL = 10.0 for
both studies. On the other hand, the Karlovitz numbers are slightly different
due to the disparity in integral length scales. Values of these parameters as
relevant to the present study have been tabulated in Table 6.1.
133
6. Towards Validation with the Leeds Experiment
6.3 Summary: the ongoing comparative study
Initial comparisons with the Leeds experiment demonstrate a similarity in
curvature distributions between two flame surfaces. The ongoing investigation
has warranted several modifications of the DNS solver.
Firstly, a 25-step reduced chemical kinetic model is studied and imple-
mented for simulating lean methane-air combustion (Smooke and Giovangigli,
1991). The detailed description of chemical pathways and Lewis number effects
in this framework allows for realistic simulations comparable to experiments.
Secondly, a modified form (Carroll and Blanquart, 2013) of Lundgren’s forcing
strategy is implemented to maintain the intensity of upstream turbulence.
Furthermore, the computational expense of the necessary large-scale de-
tailed simulations are reduced by an improved parameter selection in the
adaptive timestep controller (Söderlind and Wang, 2006). Finally, in order
to better address the experimental setup, the largest possible domain size is
designed. Resolution constraints of turbulence and chemistry are satisfied in
a computationally feasible manner using a novel strategy.
It remains to be seen if fully-developed flames of the same turbulence inten-
sity will exhibit similar propagation behaviours as recorded in the preliminary
comparisons. As pursued with the previous dataset, Damköhler’s hypotheses
will be verified as a starting point.
134
Chapter 7
Conclusion
This thesis investigates the effects of turbulence intensity u′ on flame prop-
agation in gas mixtures. It focuses upon ‘the bending effect’ of turbulent
burning velocity sT (u′). A DNS approach is adopted in light of the recent
computational advances. In the present DNS approach, governing equations
of turbulent reactive flow are solved in compressible form. Using this method,
two lines of investigation are pursued.
The first comprises the design, generation and analysis of a parametric
DNS dataset that simulates flame propagation in the TRZ regime. The bend-
ing effect of sT (u′) is captured upon increasing u′ systematically across this
regime (see Chapter 3). Subsequently, the validity of Damköhler’s first and
second hypotheses is investigated in the dataset generated. Damköhler’s first
hypothesis (DH1) is found to be valid throughout the dataset (see Chapter 4).
Motivated by this finding, the framework of the FSD approach is used to study
the underlying mechanisms governing flame propagation. The bending effect
is thereby explained in terms of a competition between components of the FSD
source-term. Damköhler’s second hypothesis (DH2) is found to have limited
validity in the dataset (see Chapter 5). Minimal conditions for the validity of
DH2 are expressed in terms of scalar flux transport regimes.
The second line of investigation aims to verify the explanation proposed for
the bending effect so far. Experiments are conducted in the Leeds explosion
vessel under conditions of bending informed by the parametric study. In order
to facilitate comparisons, a new DNS is designed to simulate flame propagation
135
7. Conclusion
under appropriate conditions (see Chapter 6). Several modifications are made
to the DNS approach adopted previously so as to address comprehensively the
goals of this investigation. Upcoming comparative results may help contrast
the proposed explanation with the mechanism of bending proposed historically
by the Leeds group.
Conclusions drawn from each line of investigation are summarized below.
Parametric study of the bending effect
Firstly, the bending effect of sT (u′) does not depend necessarily upon localized
flame quenching as is widely believed (Bradley, 1992). Indeed, the reaction
zone may remain intact (see Figure 3.7 and Figure 3.8) even in turbulence with
u′ > 10sL. The emergent flame propagation may exhibit the bending effect
nevertheless (see Figure 3.9).
Secondly, turbulent flame surface area AT and Flame Surface Density Σ
continue to govern sT (according to DH1) in the bending region (see Fig-
ure 4.4). This contrasts with recent experimental measurements in u′ > 10sL
(Gülder and Smallwood, 2007; Wabel et al., 2016). The present computa-
tions make it plausible that bending occurs as a consequence of the kinematic
response of the flame rather than its structural disruption.
Thirdly, the contributions to Σ of local flame surface curvatures hm do not
necessarily ‘even out’ as presumed in previous explanations of the bending
effect (Bradley et al., 1992). Indeed, local displacement speed sd regulates the
overall curvature contribution sdhm to the source-term of Σ. This contribution
is negative and increases in magnitude with u′, whereas the contribution of
strain rate at remains positive. The overall Σ for any given u′ as well as
the corresponding sT emerges from the balance between destructive curvature
propagation sdhm < 0 and extensive tangential strain at > 0 (see Figure 4.10).
Lastly, the enhancement of sT due to turbulent diffusion as dictated by DH2
is limited to high-intensity turbulence where gradient transport operates (see
Figure 5.7). In low-intensity turbulence corresponding to the Bray number
NB > 1, counter-gradient transport prevails, so that DT < 0, and the DH2
expression is invalidated. This suggests that Damköhler’s mechanisms may
136
7. Conclusion
neither act ‘always’ nor ‘enhance’ sT as hypothesized (Damköhler, 1940).
Ongoing comparisons with the experimental measurements in Leeds might
help validate these conclusions. Insights gained so far are summarized below.
Comparisons with the Leeds experiment
Initial comparison between the experimentally reconstructed flame surface and
the corresponding parametric DNS isosurface show striking similarities (see
Figure 6.1 and Figure 6.1b). This suggests that curvature effects may indeed
play an important role in overall flame propagation characteristics. Conclu-
sions are reserved at present pending the completion of the ongoing simulation.
Design of the corresponding DNS has provided the following conclusions.
Firstly, the detailed 25-step chemical kinetic model (Smooke and Gio-
vangigli, 1991) has a resolution constraint of ∆xc . δfineL /22.3 (see Figure 6.3).
The lean methane-air flame simulated using this model addresses non-unity
Lewis number effects thereby complementing the parametric study.
Secondly, the modified form (Carroll and Blanquart, 2013) of the forcing
scheme introduced by Lundgren (2003) maintains turbulent kinetic energy in
an HIT simulation (see Figure 6.5a). The upstream region in an inflow-outflow
flame configuration may be forced selectively to maintain the regime of com-
bustion (see Figure 6.5b).
Thirdly, adaptive timestep selection based on the PID controller (Kennedy
and Carpenter, 2003) used in the parametric investigation is improved using a
controller developed by Söderlind and Wang (2006). Not only is the timestep
history more stable, but the magnitude of timesteps also increases by nearly
two orders of magnitude in the flame simulation (see Figure 6.4).
Lastly, an optimal size of the DNS domain can be obtained by considering
various constraints due to resolution requirements on the one hand and com-
putational requirements on the other. On these lines, a basic investigation has
provided a relatively large domain configuration with size L = 0.015m three
times the previous parametric study (see Figure 6.6). This optimized domain
helps enhance the length scale ratio `0/δL.
Several possibilities for further investigation have emerged. These comprise
137
7. Conclusion
ongoing work being pursued in collaboration with researchers in the interna-
tional community as well as proposals for future work.
Ongoing and future work
Firstly, the effect of flame configuration on efficiency factors is being investi-
gated. These include the stretch factor I0 (seen in Figure 4.5) and efficiency
function α (seen in Figure 5.6). DNS datasets of a V-flame and a stagnation
flame are being analysed for this purpose in collaboration with Dr. Umair
Ahmed (University of Manchester). The objective is to provide insights for
flamelet modelling of turbulent premixed flames.
Secondly, the flame response to transient strain is being investigated under
the guidance of Dr. Alan Kerstein (Consultant, USA). Correlations between
the normal displacement speed sn and the local strain at are being studied
(on the lines of Figure 4.6 and Figure 4.9). This may help address the effect
of the strain component at on the propagation component sdhm of the FSD
source-term.
Thirdly, a comparison of the present DNS dataset with a similarly obtained
dataset (at the University of Lund, Sweden) using an incompressible formula-
tion (Yu et al., 2015) is being pursued. The comparisons are expected to shed
light on the role of gas expansion on computed sT values in the bending effect.
At the same time, recent discussions with Prof. Ö. Gülder and co-workers
(University of Toronto) have highlighted similarities recorded in the computed
Σ profile (see Figure 4.1b) with experimentally measured profiles. The recent
study conducted by Kheirkhah and Gülder (2015) links the increasingly up-
stream peaking of Σ with a transition to counter-gradient transport (as also
presently recorded). This line of investigation is being pursued further.
Furthermore, improvements in adaptive timestep size (see Section 6.2.2)
have led to a detailed discussion with Dr. Marc Day (Lawrence Berkeley
National Laboratory, USA). Assessment of the solution error and convergence
properties is being pursued to formalize the new control procedure.
Finally, the success of the optimization study (see Section 6.2.4) shows
138
7. Conclusion
that it is possible to automate the design of a DNS configuration. Algorithms
developed for this purpose may also be used in the context of experiments in
which an optimized search for design criteria is usually the centre point of all
measurements. Indeed, similar algorithms may also be applied to experimental
and DNS datasets to reveal critical relationships in the data.
To c or not to c?
A central issue of turbulent flame propagation modelling is the applicability of
the flamelet description under increasingly-intense turbulence. In this regard,
the scientific community has been bent on demarcating the regime where flame
velocity is inhibited by turbulence due to flame structure disruption, which
would thereby invalidate the framework of the progress variable c. The present
thesis informs us that inhibition of flame propagation velocity, such as that
observed in the bending effect, may occur without flame structure disruption.
As opposed to the effect of turbulence on the flame, the effect of the flame
on turbulence as well as the self-inhibition of the flame is highlighted. This
suggests that we may need to consider all possible interactions between the
flame and turbulence in order to address the question: ‘How fast we can burn?’
139

Appendix A
List of Publications
The contents of this thesis have been published the following communications:
1. Nivarti, G. V. and Cant, R. S. (2015) Aerodynamic quenching and burn-
ing velocity of turbulent premixed methane-air flames, Proceedings of the
ASME Turbo Expo 2015: Turbomachinery Technical Conference and Ex-
position, GT2015-43416 (Chapter 3).
2. Nivarti, G. V. and Cant, R. S. (2015) Premixed flame propagation in
high-intensity turbulence: investigating the role of detailed chemistry,
Proceedings of the 25th International Colloquium on Dynamics and Ex-
plosions of Reactive Systems, paper 028 (Chapter 3)
3. Nivarti, G. V. and Cant, R. S. (2017) Direct Numerical Simulation of the
bending effect in turbulent premixed flames, Proceedings of the Combus-
tion Institute 36: 1903 – 1910 (Chapter 3 and Chapter 4)
4. Nivarti, G. V. and Cant, R. S. (2017) Analysis of the bending effect
in turbulent burning velocity of premixed flames using Direct Numerical
Simulation, Combustion Theory and Modelling (under review, Chap-
ter 4)
5. Nivarti, G. V. and Cant, R. S. (2017) Scalar transport and the validity
of Damköhler’s second hypothesis in turbulent premixed flames, Physics
of Fluids 29: 085107 (Chapter 5)
141
Appendix A
6. Ahmed, U., Nivarti, G. V., Prosser, R. and Cant, R. S. (2017) Flame con-
figuration effects on flamelet modelling in premixed turbulent combustion,
Flow, Turbulence and Combustion (in preparation, Chapter 1)
The content has also been presented at the following international conferences:
1. Nivarti, G. V. and Cant, R. S. (2015) ASME Turbo Expo 2015: Turbo-
machinery Technical Conference and Exposition, Montréal
2. Nivarti, G. V. and Cant, R. S. (2015) 25th International Colloquium on
Dynamics and Explosions of Reactive Systems, Leeds
3. Nivarti, G. V. and Cant, R. S. (2015) 16th International Conference on
Numerical Combustion, Avignon
4. Nivarti, G. V. and Cant, R. S. (2016) Joint British, Spanish and Por-
tuguese Section Combustion Meeting, Cambridge
5. Nivarti, G. V., and Cant, R. S. (2016) 36th International Symposium on
Combustion, Seoul
6. Nivarti, G. V., and Cant, R. S. (2017) 17th International Conference on
Numerical Combustion, Orlando
7. Nivarti, G. V., Thorne, B., Cant, R. S., Lawes, M., Bradley, D. (2017)
17th International Conference on Numerical Combustion, Orlando
The dataset generated in this work is available for investigation upon request.
142
References
Abdel-Gayed, R., Bradley, D., and Lawes, M. (1987). Turbulent burning ve-
locities: a general correlation in terms of straining rates. In Proceedings
of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences, 414, 389–413.
Abdel-Gayed, R. G., Al-Khishali, K. J., and Bradley, D. (1984). Turbulent
burning velocities and flame straining in explosions. In Proceedings of the
Royal Society of London A: Mathematical, Physical and Engineering Sci-
ences, 391, 393–414.
Abdel-Gayed, R. G. and Bradley, D. (1977). Dependence of turbulent burning
velocity on turbulent Reynolds number and ratio of laminar burning velocity
to RMS turbulent velocity. 16th Symposium (International) on Combustion,
1725–1735.
Abdel-Gayed, R. G., Bradley, D., Hamid, M. N., and Lawes, M. (1985). Lewis
number effects on turbulent burning velocity. 20th Symposium (Interna-
tional) on Combustion, 505–512.
Abdel-Gayed, R. G., Bradley, D., and Lung, F.-K. (1989). Combustion regimes
and the straining of turbulent premixed flames. Combustion and Flame, 76,
213–218.
Ahmed, U., Nivarti, G. V., Prosser, R., and Cant, R. S. (2017). Effects of flame
configuration on flamelet modelling of turbulent premixed flames. Flow,
Turbulence and Combustion. (in preparation).
Andrews, G., Bradley, D., and Lwakabamba, S. (1975). Measurement of tur-
bulent burning velocity for large turbulent Reynolds numbers. 655–664.
143
References
Aspden, A. J., Bell, J. B., Day, M. S., Woosley, S. E., and Zingale, M. (2008).
Turbulence-flame interactions in type Ia supernovae. The Astrophysical
Journal, 689, 1173–1185.
Aspden, A. J., Day, M. S., and Bell, J. B. (2011). Turbulence-flame interactions
in lean premixed hydrogen: transition to the distributed burning regime.
Journal of Fluid Mechanics, 680, 287–320.
Bagdanavicius, A., Bowen, P. J., Bradley, D., Lawes, M., and Mansour, M. S.
(2015). Stretch rate effects and Flame Surface Densities in premixed turbu-
lent combustion up to 1.25 MPa. Combustion and Flame, 162, 4158–4166.
Batchelor, G. K. (1952). The effect of homogeneous turbulence on material
lines and surfaces. In Proceedings of the Royal Society of London A: Math-
ematical, Physical and Engineering Sciences, 213, 349–366.
Batchelor, G. K. and Townsend, A. A. (1948). Decay of Turbulence in the
Final Period. Proceedings of the Royal Society A: Mathematical, Physical
and Engineering Sciences, 194, 527–543.
Bedat, B. and Cheng, R. (1995). Experimental study of premixed flames in
intense isotropic turbulence. Combustion and Flame, 100, 485–494.
Bell, J. B., Day, M. S., Grcar, J. F., Lijewski, M. J., Driscoll, J. F., and
Filatyev, S. A. (2007). Numerical simulation of a laboratory-scale turbulent
slot flame. Proceedings of the Combustion Institute, 31, 1299–1307.
Bell, J. B., Day, M. S., Shepherd, I. G., Johnson, M. R., Cheng, R. K., Grcar,
J. F., Beckner, V. E., and Lijewski, M. J. (2005). Numerical simulation of
a laboratory-scale turbulent V-flame. Proceedings of the National Academy
of Sciences of the United States of America, 102, 10006–10011.
Bilger, R. W., Pope, S. B., Bray, K. N. C., and Driscoll, J. F. (2005). Paradigms
in turbulent combustion research. Proceedings of the Combustion Institute,
30, 21–42.
Bird, R. B., Stewart, W. E., and Lightfoot, E. N. (1960). Transport Phenom-
ena. John Wiley and Sons, NY.
144
References
Boger, M., Veynante, D., Boughanem, H., and Trouvé, A. (1998). Direct
Numerical Simulation analysis of Flame Surface Density concept for Large
Eddy Simulation of turbulent premixed combustion. 27th Symposium (In-
ternational) on Combustion, 917–925.
Borghi, R. (1985). On the structure and morphology of turbulent premixed
flames. In Recent Advances in the Aerospace Sciences, Casci, C. and Bruno,
C., editors, 117–138. Springer US.
Bradley, D. (1992). How fast can we burn? 24th Symposium (International)
on Combustion, 247–262.
Bradley, D., Lau, A. K. C., and Lawes, M. (1992). Flame stretch rate as
a determinant of turbulent burning velocity. Philosophical Transactions of
the Royal Society of London A: Mathematical, Physical and Engineering
Sciences, 338, 359–387.
Bradley, D., Lawes, M., Liu, K., and Mansour, M. S. (2013). Measurements
and correlations of turbulent burning velocities over wide ranges of fuels and
elevated pressures. Proceedings of the Combustion Institute, 34, 1519–1526.
Bradley, D., Lawes, M., and Mansour, M. S. (2011). The problems of the
turbulent burning velocity. Flow, Turbulence and Combustion, 87, 191–204.
Bray, K. N. C. (1980). Turbulent flows with premixed reactants. In Turbulent
Reacting Flows, Libby, P. A. and Williams, F. A., editors, 44, 115–183.
Springer Berlin Heidelberg.
Bray, K. N. C. (1990). Studies of the turbulent burning velocity. Proceedings
of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences, 431, 315–335.
Bray, K. N. C. (1995). Turbulent transport in flames. Proceedings of the Royal
Society of London A: Mathematical, Physical and Engineering Sciences, 451,
231–256.
145
References
Bray, K. N. C. and Cant, R. S. (1991). Some applications of Kolmogorov’s
turbulence research in the field of combustion. Proceedings of the Royal
Society A: Mathematical, Physical and Engineering Sciences, 434, 217–240.
Bray, K. N. C. and Moss, J. B. (1977). A unified statistical model of the
premixed turbulent flame. Acta Astronautica, 4, 291–319.
Buckmaster, J. D. and Ludford, G. S. S. (1982). Theory of Laminar Flames.
Cambridge University Press.
Candel, S. M. and Poinsot, T. J. (1990). Flame stretch and the balance equa-
tion for the flame area. Combustion Science and Technology, 70, 1–15.
Cant, R. S. (2013). SENGA2: User guide. Technical Report CUED-THERMO-
2012/04 - 2, Cambridge University Engineering Department.
Cant, R. S. and Mastorakos, E. (2008). An introduction to turbulent reacting
flows. Imperial College Press.
Cant, R. S., Pope, S. B., and Bray, K. N. C. (1991). Modelling of flamelet
surface-to-volume ratio in turbulent premixed combustion. 23rd Symposium
(International) on Combustion, 809–815.
Cant, S. (1999). Direct Numerical Simulation of premixed turbulent flames.
Philosophical Transactions of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences, 357, 3583–3604.
Carroll, P. L. and Blanquart, G. (2013). A proposed modification to Lundgren’s
physical space velocity forcing method for isotropic turbulence. Physics of
Fluids, 25, 105114.
Chakraborty, N. and Cant, R. S. (2004). Unsteady effects of strain rate and
curvature on turbulent premixed flames in an inflow-outflow configuration.
Combustion and Flame, 137, 129–147.
Chakraborty, N. and Cant, R. S. (2007). A priori analysis of the curvature
and propagation terms of the Flame Surface Density transport equation for
large eddy simulation. Physics of Fluids, 19, 105101.
146
References
Chakraborty, N. and Cant, R. S. (2009). Effects of Lewis number on scalar
transport in turbulent premixed flames. Physics of Fluids, 21, 035110.
Clavin, P. (1985). Dynamic behaviour of premixed flame fronts in laminar and
turbulent flows. Progress in Energy and Combustion Science, 11, 1–59.
Clavin, P. and Williams, F. A. (1979). Theory of premixed-flame propagation
in large-scale turbulence. Journal of Fluid Mechanics, 90, 589–604.
Damköhler, G. (1940). Der einfluss der turbulenz auf die flam-
mengeschwindigkeit in gasgemschen. Zeitschrift für Elektrochemie und ange-
wandte physikalische Chemie, 46, 601–652. (English translation NACA TM
1112, 1947).
Day, M., Tachibana, S., Bell, J., Lijewski, M., Beckner, V., and Cheng, R. K.
(2012). A combined computational and experimental characterisation of
lean premixed turbulent low swirl laboratory flames: I. Methane flames.
Combustion and Flame, 159, 275–290.
Driscoll, J. F. (2008). Turbulent premixed combustion: Flamelet structure and
its effect on turbulent burning velocities. Progress in Energy and Combustion
Science, 34, 91–134.
Duclos, J., Veynante, D., and Poinsot, T. (1993). A comparison of flamelet
models for premixed turbulent combustion. Combustion and Flame, 95,
101–117.
Echekki, T. and Chen, J. H. (1996). Unsteady strain rate and curvature effects
in turbulent premixed methane-air flames. Combustion and Flame, 106, 184–
190.
Echekki, T. and Chen, J. H. (1998). Structure and propagation of methanol–air
triple flames. Combustion and Flame, 114, 231–245.
Eswaran, V. and Pope, S. B. (1988). An examination of forcing in Direct
Numerical Simulations of turbulence. Computers and Fluids, 16, 257–278.
147
References
Filatyev, S. A., Driscoll, J. F., Carter, C. D., and Donbar, J. M. (2005).
Measured properties of turbulent premixed flames for model assessment,
including burning velocities, stretch rates, and surface densities. Combustion
and Flame, 141, 1–21.
Frisch, U. (2004). Turbulence. Cambridge University Press.
Fru, G., Thévenin, D., and Janiga, G. (2011). Impact of turbulence intensity
and equivalence ratio on the burning rate of premixed methane–air flames.
Energies, 4, 878–893.
Gillespie, L., Lawes, M., Sheppard, C. G. W., and Woolley, R. (2000). As-
pects of laminar and turbulent burning velocity relevant to SI engines. SAE
Technical Paper, 2000-01-0192.
Gran, I. R., Echekki, T., and Chen, J. H. (1996). Negative flame speed in
an unsteady 2-d premixed flame: A computational study. 26th Symposium
(International) on Combustion, 323–329.
Griffiths, R. A. C., Chen, J. H., Kolla, H., Cant, R. S., and Kollmann, W.
(2015). Three-dimensional topology of turbulent premixed flame interaction.
Proceedings of the Combustion Institute, 35, 1341–1348.
Gülder, Ö. L. (1991). Turbulent premixed flame propagation models for dif-
ferent combustion regimes. 23rd Symposium (International) on Combustion,
743–750.
Gülder, Ö. L. and Smallwood, G. J. (2007). Flame surface densities in premixed
combustion at medium to high turbulence intensities. Combustion Science
and Technology, 179, 191–206.
Hamlington, P. E., Poludnenko, A. Y., and Oran, E. S. (2011). Interactions
between turbulence and flames in premixed reacting flows. Physics of Fluids,
23, 125111.
Hawkes, E. R., Chatakonda, O., Kolla, H., Kerstein, A. R., and Chen, J. H.
(2012). A petascale Direct Numerical Simulation study of the modelling of
148
References
flame wrinkling for Large-eddy Simulations in intense turbulence. Combus-
tion and Flame, 159, 2690–2703.
Haworth, D. C. and Poinsot, T. J. (1992). Numerical simulations of Lewis
number effects in turbulent premixed flames. Journal of Fluid Mechanics,
244, 405–436.
Ishihara, T., Gotoh, T., and Kaneda, Y. (2009). Study of high-Reynolds
number isotropic turbulence by Direct Numerical Simulation. Annual Review
of Fluid Mechanics, 41, 165–180.
Jenkins, K. W. and Cant, R. S. (1999). Direct Numerical Simulation of tur-
bulent flame kernels. In Recent Advances in DNS and LES, Kinight, D. and
Sakell, L., editors, 54, 191–202.
Jones, W. P. and Launder, B. (1972). The prediction of laminarization with a
two-equation model of turbulence. International Journal of Heat and Mass
Transfer, 15, 301–314.
Karlovitz, B., Denniston, D. W., Knapschaefer, D., and Wells, F. E. (1953).
Studies on turbulent flames: A. Flame propagation across velocity gradients
B. Turbulence measurement in flames. 4th Symposium (International) on
Combustion, 613–620.
Karlovitz, B., Denniston Jr, D. W., and Wells, F. E. (1951). Investigation of
turbulent flames. Journal of Chemical Physics, 19, 541–547.
Karpov, V. and Severin, E. (1980). Effects of molecular-transport coefficients
on the rate of turbulent combustion. Combustion, Explosion, and Shock
Waves, 16, 41–46.
Kee, R. J., Dixon-lewis, G., Warnatz, J., Coltrin, M. E., and Miller, J. A.
(1986). A Fortran computer code package for the evaluation of gas-phase,
multicomponent transport properties. Technical Report SAND86-8246, San-
dia National Laboratories.
149
References
Kennedy, C. A. and Carpenter, M. H. (2003). Additive Runge-Kutta schemes
for convection-diffusion-reaction equations. Applied Numerical Mathematics,
44, 139–181.
Kerstein, A. R. (2002). Turbulence in combustion processes: Modelling chal-
lenges. Proceedings of the Combustion Institute, 29, 1763–1773.
Kerstein, A. R., Ashurst, W. T., and Williams, F. A. (1988). Field equation
for interface propagation in an unsteady homogeneous flow field. Physical
Review A, 37, 2728.
Kheirkhah, S. and Gülder, Ö. L. (2015). Consumption speed and burning
velocity in counter-gradient and gradient diffusion regimes of turbulent pre-
mixed combustion. Combustion and Flame, 162, 1422–1439.
Kido, H., Kitagawa, T., Nakashima, K., and Kato, K. (1989). An improved
model of turbulent mass burning velocity. In Memoirs of the Faculty of
Engineering, Kyushu University, 49, 229–247.
Klimov, A. M. (1963). Laminar flame in turbulent flow. Prikladnoy Mekhaniki
i Tekhnicheskoy Fiziki Zhurnal, 3, 49 – 58. (English translation AD-A200
241 Foreign Technology Division, Air Force Systems Command, 1988).
Kobayashi, H., Tamura, T., Maruta, K., Niioka, T., andWilliams, F. A. (1996).
Burning velocity of turbulent premixed flames in a high-pressure environ-
ment. 26th Symposium (International) on Combustion, 389–396.
Kolmogorov, A. N. (1941). The local structure of turbulence in incompressible
viscous fluid for very large Reynolds numbers. Doklady Akadademii Nauk
SSSR, 30, 301–305.
Lapointe, S. and Blanquart, G. (2016). Fuel and chemistry effects in high
Karlovitz premixed turbulent flames. Combustion and Flame, 167, 294–307.
Lapointe, S., Savard, B., and Blanquart, G. (2015). Differential diffusion ef-
fects, distributed burning, and local extinctions in high Karlovitz premixed
flames. Combustion and Flame, 162, 3341–3355.
150
References
Law, C. K. (1989). Dynamics of stretched flames. 22nd Symposium (Interna-
tional) on Combustion, 1381–1402.
Lewis, B. and von Elbe, G. (1934). On the theory of flame propagation. The
Journal of Chemical Physics, 2, 537–546.
Libby, P. A. and Bray, K. N. C. (1981). Countergradient diffusion in premixed
turbulent flames. AIAA Journal, 19, 205–213.
Libby, P. A., Bray, K. N. C., and Moss, J. B. (1979). Effects of finite reaction
rate and molecular transport in premixed turbulent combustion. Combustion
and Flame, 34, 285–301.
Lipatnikov, A. and Chomiak, J. (2002). Turbulent flame speed and thickness:
phenomenology, evaluation, and application in multi-dimensional simula-
tions. Progress in Energy and Combustion Science, 28, 1–74.
Lipatnikov, A. N., Chomiak, J., Sabelnikov, V. A., Nishiki, S., and Hasegawa,
T. (2015). Unburned mixture fingers in premixed turbulent flames. Proceed-
ings of the Combustion Institute, 35, 1401–1408.
Lockermann, G. (1956). The centenary of the bunsen burner. Journal of
Chemical Education, 33, 20–22.
Louch, D. S. and Bray, K. N. C. (2001). Vorticity in unsteady premixed
flames: vortex pair-premixed flame interactions under imposed body forces
and various degrees of heat release and laminar flame thickness. Combustion
and Flame, 125, 1279–1309.
Lundgren, T. S. (2003). Linearly forced isotropic turbulence. In Annual Re-
search Briefs, CTR Stanford, 461–473.
Mallard, F. E. and le Chatelier, H. (1883). Recherches expérimenatales et
théoriques sur le combustion des mélanges gazeux explosifs. Annales des
Mines, 4, 274–378.
Markstein, G. H. (1964). Nonsteady flame propagation. Pergamon Press.
151
References
McBride, B. J. and Gordon, S. (1967). FORTRAN IV program for calculation
of thermodynamic data. Technical Report NASA TN D-4097, NASA.
McMurtry, P. A., Jou, W.-H., Riley, J., and Metcalfe, R. (1986). Direct Nu-
merical Simulations of a reacting mixing layer with chemical heat release.
AIAA Journal, 24, 962–970.
Meneveau, C. and Poinsot, T. (1991). Stretching and quenching of flamelets
in premixed turbulent combustion. Combustion and Flame, 332, 311–332.
Michelson, V. A. (1889). Ueber die normale entzündungsgeschwindigkeit ex-
plosiver gasgemischen. Annalen der Physik, 273, 1–24.
Nivarti, G. V. and Cant, R. S. (2015a). Aerodynamic quenching and burn-
ing velocity of turbulent premixed methane-air flames. Proceedings of the
ASME Turbo Expo: Turbomachinery Technical Conference and Exposition,
GT2015–43416.
Nivarti, G. V. and Cant, R. S. (2015b). Premixed flame propagation in high-
intensity turbulence: Investigating the role of detailed chemistry. In 25th
International Conference for Dynamics of Explosions and Reactive Systems.
Nivarti, G. V. and Cant, R. S. (2017). Direct numerical simulation of the
bending effect in turbulent premixed flames. Proceedings of the Combustion
Institute, 36, 1903–1910.
Nivarti, G. V., Thorne, B., Cant, R. S., Lawes, M., and Bradley, D. (2017).
Comparisons of Flame Surface Density measurements with Direct Numerical
Simulations of a lean methane-air flame in high-intensity turbulence. In 16th
International Conference on Numerical Combustion, Orlando.
Peters, N. (1988). Laminar flamelet concepts in turbulent combustion. 21st
Symposium (International) on Combustion, 1231–1250.
Peters, N. (1999). The turbulent burning velocity for large-scale and small-
scale turbulence. Journal of Fluid Mechanics, 384, 107–132.
Peters, N. (2000). Turbulent combustion. Cambridge University Press.
152
References
Peters, N., Terhoeven, P., Chen, J. H., and Echekki, T. (1998). Statistics of
flame displacement speeds from computations of 2-D unsteady methane-air
flames. 27th Symposium (International) on Combustion, 833–839.
Peters, N. and Williams, F. A. (1987). The asymptotic structure of stoichio-
metric methane-air flames. Combustion and Flame, 68, 185–207.
Petersen, M. R. and Livescu, D. (2010). Forcing for statistically stationary
compressible isotropic turbulence. Physics of Fluids, 22, 116101.
Poinsot, T. (1996). Using Direct Numerical Simulations to understand pre-
mixed turbulent combustion. 26th Symposium (International) on Combus-
tion, 219–232.
Poinsot, T., Candel, S., and Trouvé, A. (1995). Applications of Direct Numer-
ical Simulation to premixed turbulent combustion. Progress in Energy and
Combustion Science, 21, 531–576.
Poinsot, T. and Lele, S. K. (1992). Boundary conditions for direct simulations
of compressible viscous flows. Journal of Computational Physics, 101, 104–
129.
Poinsot, T., Veynante, D., and Candel, S. (1991a). Diagrams of premixed
turbulent combustion based on direct simulation. 23rd Symposium (Inter-
national) on Combustion, 613–619.
Poinsot, T., Veynante, D., and Candel, S. (1991b). Quenching processes and
premixed turbulent combustion diagrams. Journal of Fluid Mechanics, 228,
561–606.
Poludnenko, A. Y. (2015). Pulsating instability and self-acceleration of fast
turbulent flames. Physics of Fluids, 27, 014106.
Poludnenko, A. Y. and Oran, E. (2011). The interaction of high-speed tur-
bulence with flames: Turbulent flame speed. Combustion and Flame, 158,
301–326.
153
References
Pope, S. B. (1988). The evolution of surfaces in turbulence. International
Journal of Engineering Science, 26, 445–469.
Pope, S. B. (2001). Turbulent flows. Cambridge University Press.
Pope, S. B., Yeung, P. K., and Girimaji, S. S. (1989). The curvature of material
surfaces in isotropic turbulence. Physics of Fluids A, 1, 2010–2018.
Prandtl, L. (1925). Bericht über untersuchungen zur ausgebildeten turbulenz.
Zeitschrift für Angewandte Mathematik und Mechanik, 5, 136–139. English
translation NACA Technical Memorandum 1231, 1949.
Quadrio, M., Frohnapfel, B., and Hasegawa, Y. (2016). Does the choice of
the forcing term affect flow statistics in DNS of turbulent channel flow?
European Journal of Mechanics-B/Fluids, 55, 286–293.
Rogallo, R. S. (1981). Numerical experiments in homogeneous turbulence.
Technical Report NASA TM 81315, NASA Ames Research Center.
Rosales, C. and Meneveau, C. (2005). Linear forcing in numerical simulations
of isotropic turbulence: Physical space implementations and convergence
properties. Physics of Fluids, 17, 095106.
Rudy, D. H. and Strikwerda, J. C. (1980). A nonreflecting outflow boundary
condition for subsonic Navier-Stokes calculations. Journal of Computational
Physics, 36, 55–70.
Rutland, C. J. and Cant, R. S. (1994). Turbulent transport in premixed flames.
In Proceedings of the Summer Program at the Centre for Turbulence Re-
search, Stanford, 75–94. NASA Ames/Stanford University.
Sankaran, R., Hawkes, E. R., Chen, J. H., Lu, T., and Law, C. K. (2007).
Structure of a spatially developing turbulent lean methane-air Bunsen flame.
Proceedings of the Combustion Institute, 31, 1291–1298.
Savard, B. and Blanquart, G. (2015). Broken reaction zone and differential
diffusion effects in high Karlovitz n-C7H16 premixed turbulent flames. Com-
bustion and Flame, 162, 2020–2033.
154
References
Savre, J., Carlsson, H., and Bai, X.-S. (2013). Tubulent methane/air pre-
mixed flame structure at high Karlovitz numbers. Flow, Turbulence and
Combustion, 90, 325–341.
Schelkin, K. I. (1943). On combustion in a turbulent flow. Zhurnal Tekhnich-
eskoi Fiziki, 13, 520–530. (English translation NACA TM 1112, 1947).
Shepherd, I., Cheng, R., Plessing, T., Kortschik, C., and Peters, N. (2002).
Premixed flame front structure in intense turbulence. Proceedings of the
Combustion Institute, 29, 1833–1840.
Shepherd, I. G., Moss, J. B., and Bray, K. N. C. (1982). Turbulent transport in
a confined premixed flame. 19th Symposium (International) on Combustion,
423–431.
Shy, S. S., I, W. K., and Lin, M. L. (2000a). A new cruciform burner and its
turbulence measurements for premixed turbulent combustion study. Exper-
imental Thermal and Fluid Science, 20, 105–114.
Shy, S. S., Lin, W. J., and Wei, J. C. (2000b). An experimental correlation of
turbulent burning velocities for premixed turbulent methane-air combustion.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering
Sciences, 456, 1997–2019.
Smith, G. P., Golden, D. M., Frenklach, M., Moriarty, N. W., Eiteneer, B.,
Goldenberg, M., Bowman, C. T., Hanson, R. K., Song, S., William C. Gar-
diner, J., Lissianski, V. V., and Qin, Z. (2017). GRI-MECH 3.0. Available
at http://www.me.berkeley.edu/gri_mech/.
Smooke, M. D. and Giovangigli, V. (1991). Formulation of the premixed and
nonpremixed test problems. In Reduced Kinetic Mechanisms and Asymptotic
Approximations for Methane-Air Flames, Smooke, M. D., editor, 384, 1–28.
Springer Berlin Heidelberg.
Söderlind, G. (2002). Automatic control and adaptive time-stepping. Numer-
ical Algorithms, 31, 281–310.
155
References
Söderlind, G. and Wang, L. (2006). Adaptive time-stepping and computational
stability. Journal of Computational and Applied Mathematics, 185, 225–243.
Steinberg, A. M. and Driscoll, J. F. (2009). Straining and wrinkling processes
during turbulence–premixed flame interaction measured using temporally-
resolved diagnostics. Combustion and Flame, 156, 2285–2306.
Steinberg, A. M., Driscoll, J. F., and Swaminathan, N. (2012). Statistics and
dynamics of turbulence–flame alignment in premixed combustion. Combus-
tion and Flame, 159, 2576–2588.
Sutherland, J. C. and Kennedy, C. A. (2003). Improved boundary conditions
for viscous, reacting, compressible flows. Journal of Computational Physics,
191, 502–524.
Trouvé, A. and Poinsot, T. (1994). The evolution equation for the Flame Sur-
face Density in turbulent premixed combustion. Journal of Fluid Mechanics,
278, 1–31.
Veynante, D. and Poinsot, T. (1997). Effects of pressure gradients on turbulent
premixed flames. Journal of Fluid Mechanics, 353, 83–114.
Veynante, D., Trouvé, A., Bray, K. N. C., and Mantel, T. (1997). Gradient
and counter-gradient scalar transport in turbulent premixed flames. Journal
of Fluid Mechanics, 332, 263–293.
Veynante, D. and Vervisch, L. (2002). Turbulent combustion modeling.
Progress in Energy and Combustion Science, 28, 193–266.
Wabel, T. M., Skiba, A. W., and Driscoll, J. F. (2016). Turbulent burning ve-
locity measurements: extended to extreme levels of turbulence. Proceedings
of the Combustion Institute, 36, 1801–1808.
Wang, H., Hawkes, E. R., and Chen, J. H. (2016a). Turbulence-flame interac-
tions in DNS of a laboratory high Karlovitz premixed turbulent jet flame.
Physics of Fluids, 28, 095107.
156
References
Wang, H., Hawkes, E. R., Zhou, B., Chen, J. H., Li, Z., and Aldén, M. (2016b).
A comparison between Direct Numerical Simulation and experiment of the
turbulent burning velocity-related statistics in a turbulent methane-air pre-
mixed jet flame at high Karlovitz number. Proceedings of the Combustion
Institute, 36, 2045–2053.
Wenzel, H. and Peters, N. (2000). Direct numerical simulation and modeling
of kinematic restoration, dissipation and gas expansion effects of premixed
flames in homogeneous turbulence. Combustion Science and Technology,
158, 273–297.
Williams, F. A. (1985). Combustion Theory. Westview Press.
Williams, F. A. (2000). Progress in knowledge of flamelet structure and ex-
tinction. Progress in Energy and Combustion Science, 26, 657–682.
Wohl, K., Shore, L., Von Rosenberg, H., and Weil, C. (1953). The burning
velocity of turbulent flames. 4th Symposium (International) on Combustion,
620–635.
Yakhot, V. (1988). Propagation velocity of premixed turbulent flames. Com-
bustion Science and Technology, 60, 191–214.
Yu, R., Bai, X.-S., and Lipatnikov, A. N. (2015). A Direct Numerical Simula-
tion study of interface propagation in homogeneous turbulence. Journal of
Fluid Mechanics, 772, 127–164.
Yuen, F. T. and Gülder, O. L. (2013). Turbulent premixed flame front dy-
namics and implications for limits of flamelet hypothesis. Proceedings of the
Combustion Institute, 34, 1393–1400.
Zel’dovich, Y. B. and Frank-Kamenetski˘i, D. A. (1938). Teori˘ia teplovogo
rasprostraneni˘ia plameni. Zhurnal Fizicheskoi Khimii, 12, 100–105. (En-
glish translation Selected Works of Yakov Borisovich Zeldovich, Volume I,
Princeton University Press, 1992.).
157