University of Southern California Dissertations and Theses

Accuracy and feasibility of combustion studies under engine relevant conditions

(USC Thesis Other)

Accuracy and feasibility of combustion studies under engine relevant conditions

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

ACCURACY AND FEASIBILITY OF COMBUSTION STUDIES UNDER ENGINE RELEVANT
CONDITIONS

by
Abtin Ansari

A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements of the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)

AUGUST 2018

Copyright 2018 Abtin Ansari

ii
Acknowledgments
First and foremost, I would like to thank my advisor Prof. Fokion Egolfopoulos for his constant support
and guidance throughout my doctoral work. I was granted the opportunity to work under his supervision
in an exciting and productive research environment. Without his guidance, this work would not have been
possible.
I would like to extend a special thanks to Prof. Larry Redekopp for his kindness, concern, and for
sharing his insight into science and giving me additional guidance. I also want to thank my dissertation
committee members Profs. Paul Ronney, Ivan Bermejo-Moreno, and Jiin-Jen Lee.
Special thanks to Dr. Vyaas Gururajan who was a mentor, providing me with valuable lessons in
coding, computational fluid dynamics, numerical simulations, and research during my doctoral studies.
He happily shared his research ideas with me and reviewed all my papers as well as this dissertation. I
also want to thank Dr. Jagannath Jayachandran for teaching me how to use his reacting flow code and Dr.
Kevin Chen for giving me valuable lessons in collaboration software programming and sharing his
stability analysis libraries.
Many thanks to my colleagues at University of Sothern California Dr. Roe Burrell, Dr. Christodoulos
Xiouris, and Dr. Francesco Carbone who helped me facilitating the experiments.
Special acknowledgements go to Drs. Gareth Oskam, Pryank Saxena, and Chiping Li for financially
supporting my research.
There are no words to express my gratitude to my parents Mohammadreza Ansari and Nay Bayani
and my nanny Parvin Ramezani for their endless love, constant support, and encouragement. They always
encouraged me to believe that I could do anything that I put my mind to and provided me with the best
support no matter what I have chosen to do. Finally, I want to thank to Sarah Mirzaie and Tannaz Telikany
for their patience and for walking this path with me.

iii
Table of Contents
1. Introduction
1.1. Background
1.2. Validation of chemical mechanisms and mass transport rates
1.3. Counterflow flame ignition
1.4. End gas ignition by spherically expanding flames in a constant volume chamber
1.5. References
2. Flame ignition in the counterflow configuration: Reassessing the experimental assumptions
2.1. Introduction
2.2. Experimental approach
2.2.1. Method and apparatus
2.2.2. Flow and flame visualization
2.3. Modeling approach
2.4. Results and discussion
2.4.1. Bifurcation effects
2.4.2. Simulation results for geometry (a)
2.4.3. Effects of jet exit parabolic temperature profile
2.4.4. Simulation results for geometry (b)
2.4.5. Simulation results for geometry (c)
2.5. Concluding remarks
2.6. Appendix A: Momentum balance
2.7. Appendix B: Details of modeling approach
2.7.1. Geometry (a) Specifications
2.7.2. Geometry (b) Specifications
2.7.3. Geometry (c) Specifications
2.8. Appendix C: Effect of intrusive temperature measurements
2.9. Appendix D: Further results and discussion
2.10. Acknowledgements
2.11. References
3. Effects of confinement, geometry, inlet velocity profile, and Reynolds number on the asymmetry of
opposed-jet flows
3.1. Introduction
3.2. Experimental approach
3.3. Modeling approach
3.3.1. 1-D simulations
1
1
2
2
5
11
14
14
18
18
20
22
25
25
29
31
34
38
41
41
44
46
49
51
52
55
59
60

66
66
72
75
75

iv
3.3.2. 2-D simulations
3.3.2.1. Transient flow simulations
3.3.2.2. Bifurcation analysis
3.3.2.3. Geometry, initial and boundary conditions and cases investigated
3.4. Results and discussion
3.4.1. Experimental results
3.4.2. Numerical results
3.4.2.1. Momentum balanced conditions
3.4.2.1.1. Confinement effects
3.4.2.1.2. Aspect ratio effects
3.4.2.1.3. Inlet velocity profile effects
3.4.2.1.4. Geometrical asymmetry and imperfection effects
3.4.2.2. Sensitivity to momentum imbalance
3.5. Concluding remarks
3.6. Acknowledgements
3.7. References
4. Parameters influencing the burning rate of laminar flames propagating into a reacting mixture
4.1. Introduction
4.2. Modeling approach
4.3. Results and discussion
4.4. Concluding remarks
4.5. Appendix A: Details of transient numerical simulations
4.6. Appendix B: Details of the developed sensitivity analysis approach
4.7. Appendix C: Transient effects of stopping and restarting the simulation
4.8. Appendix D: Sensitivity analysis results for 𝜏
4.9. Appendix E: Sensitivity analysis results for 𝑚 ̇
4.10. Appendix F: Predicting the 𝑚 ̇ behavior using LSC values
4.11. Acknowledgements
4.12. References
5. Concluding Remarks and Recommendations
5.1. Counterflow configuration
5.2. End gas ignition by spherically expanding flames in a constant volume chamber
5.3. References
75
75
75
77
82
82
87
87
91
92
92
94
95
98
99
100
104
104
106
107
115
115
117
123
125
126
127
127
128
131
131
132
133

1
Chapter 1
Introduction
1.1. Background
Worldwide dependence on energy supply and its ever increasing demand is a known and over-
discussed topic (e.g., [1]). Fossil fuels are particularly of interest due to their high energy density. The
U.S. Energy Information Administration reports a 50% increase of worldwide energy consumption in the
span of 30 years from 2010 [2]. Although nuclear and renewable sources of energy are growing alongside,
fossil fuels will still need to satisfy around 80% of the energy needs by 2040.
Energy production via fossil fuels are viable through combustion processes in a controllable
environment. Ignition is one of the fundamental combustion processes. Simplest example of ignition event
can be a mundane lighting a gas stove. On the other side of the spectrum, ignition is fundamental in the
complex process of mid-air re-ignition of a gas turbine engine of an aircraft.
High dependency on fossil fuels however has brought about problems of global warming and air
pollution which in turn has imposed severe restrictions on emissions. Such requirements have forced
combustors designers not to rely on semi-empirical or simplified combustion models. Therefore, the motto
in numerical modeling is shifting from being predictive to becoming an optimization tool to control the
desirable variables in combustion processes in real time. Achieving such comprehensive models is not a
trivial task since combustion processes involve coupling of chemistry, mass, heat and momentum transfer
and scales as small as those dealt with in quantum mechanics realm all the way to large scales of fluid
mechanics [3].
In principle, ignition is a limit phenomenon of transition of non-burning state of a combustible mixture
to a burning state which occurs instantaneously due to exponential dependence of chemistry on
temperature. Generally, ignition can be divided into “auto-ignition” and ‘forced ignition”. Auto-ignition
occurs when the temperature of a combustible mixture is raised by some adiabatic means such as isentropic
compression. When heat or radical production exceeds their loss, the whole mixture explodes
instantaneously. Forced ignition takes place when energy is deposited into the mixture locally by some
external source such as heated surface, laser beam or electrical discharge. Subsequently a flame front will
be initiated locally and it propagates through the rest of the mixture. There is an extensive literature on
ignition that an interested reader can consult (e.g., [4-10]).

2
1.2. Validation of chemical mechanisms and mass transport rates
One of the main goals of combustion research is to provide combustor designers with reliable
chemistry models (reaction mechanisms) and mass transport models which can comprehensively describe
the combustion of various fuel systems. Combustion models should be able to predict the observables
measured in well-controlled environments such as canonical experimental setups. Only in these setups,
the dimensionality of the physics involved can be reduced and effect of each variable can be studied
independently. The experimental results can give feedback to kinetic and transport models to further
reduce the discrepancies.
Real engines and combustors have highly turbulent flows. Hence, they are impractical to be
systematically studied due to highly coupled physics. For certain setups and experimental conditions, by
comparison between relevant time scales such as fluid mechanics and radical pool production rate,
unimportant disciplines can be neglected from the analysis. Such simplification allows the use of simple
models to describe the phenomena observed experimentally. In the middle of the spectrum between
simplified homogenous systems and complex real-life applications, there stands laminar combustion field
which can be achieved using canonical configurations with very well-defined boundary conditions (BC).
In this regime, reaction mechanisms and transport models can be targeted for validations and pertinent
scaling parameters can be derived to aid the turbulent combustion community.
The main validation target and scaling parameter in ignition studies performed in laminar combustion
field is “ignition delay time” for homogenous systems (e.g., [4]) and “ignition temperature” for
convective-diffusive systems (e.g., [9]). In a homogenous and adiabatic system, if the temperature of a
reacting mixture is instantaneously raised to a certain temperature, the chemical reactions start to take
place and the mixture will ignite after some time which is defined as ignition delay time. When there exist
losses due to transport phenomena, the system is no longer homogenous and the ignition occurrence
depends on the competition between transport losses and chemical gains. For a certain time scale
determined by the rate of losses, there exist a minimum temperature required at which chemical reactions
overcome the losses and the mixture ignites. This temperature is referred to as ignition temperature.
1.3. Counterflow flame ignition
In the realistic environments, such as internal combustion engines, ignition occurs in the presence of
significant spatial inhomogeneities namely temperature and concentration gradients. Therefore, the effect
of transport processes on ignition cannot be neglected. The counterflow configuration (e.g., [11]) provides

3
the freedom of evaluating the effect of transport processes on ignition while maintaining the low-
dimensionality.
In this configuration, two opposed gaseous jets of fuel/inert mixture at room temperature and
electrically-heated air impinge on each other to create a stagnation-type flow field (e.g., [12]) as
schematically shown in Fig. 1.1.

Fig. 1.1. Schematic of a typical axisymmetric counterflow.

Ignition is achieved by gradually increasing the temperature of the air jet until a flat flame appears
close to the stagnation plane. Theoretically, the response of counterflow flames can be presented typically
in the form of S-curve [9] by plotting either the maximum reaction temperature versus a system Damkohler
number, Da, a non-dimensional ratio between a characteristic flow time and a characteristic reaction time.
Figure 1.2 is a schematic of a typical S-curve.

Fig. 1.2. Schematic of counterflow response (maximum temperature) as a function of the Damkohler
number.
Air
Fuel/inert
Stagnation
Plane
Symmetry
Axis

4
At low temperatures, the system exhibits no reactivity (LA branch in Fig. 1.2). By increasing the
temperature, chemical reactions take place until ignition is achieved and the system jumps to a vigorously
burning state (IF path in Fig. 1.2). This sudden transition is an indication of “saddle point bifurcation” and
has been introduced in the analysis of dynamical systems (e.g., [13]); this is a manifestation of the inherent
non-linear dependency of chemistry to temperature.
This configuration is, in principle, a well-controlled environment in which global flame properties
such as laminar flame speeds, extinction strain rates, and ignition temperatures can be defined and
measured (e.g., [9,14-16]). Moreover, the aforementioned quantities can be used as validation targets for
kinetic models (e.g., [9,14-16]). The counterflow configuration is usually modeled by assuming that the
flow-field is steady and (quasi) one-dimensional (1-D) [17]. However, these assumptions must be first
verified before any comparison can be made between the experimental and simulation results. Therefore,
one of the main objective of this thesis was to assess the validity of 1-D modeling under wide range of
conditions, and to provide guidelines to experimentalist for designing counterflow setups that conform
best with the underlying assumptions.
Chapter 2 contains the results of experiments and direct numerical simulations that assess the validity
of the assumptions of 1-D formulation for counterflow configuration. Experimentally, particle image
velocimetry, shadowgraph, and a high-speed camera were employed to characterize the flow field before
ignition, and to capture the ignition position and further evolution of the flame. The modeling involved
axisymmetric numerical simulations using detailed molecular transport and chemical kinetic models.
Both experiments and simulations revealed that if solid surfaces are present in the vicinity of the jets exit,
the flow separates generating recirculation zones that are unstable and result in the bifurcation of the flow
field. As a result, for a given set of boundary conditions at the burners’ exits, there exists two possible
stable states of the flow field which have different velocity and scalars distribution, and the fuel
concentration at which ignition occurs was determined to differ for these two states. A novel approach is
proposed to correct for the unavoidable radial non-uniformity of the temperature profile at the exit of the
heated jet and the conditions that do not result in bifurcation are outlined, so that the results from 1-D
codes can be compared to the data with confidence.
Chapter 3 contains the results of experiments and numerical simulations that investigate the behavior
and controlling parameters of counterflowing isothermal air jets for various nozzle designs, Reynolds
numbers, and surrounding geometries. The flow field in the jets’ impingement region was analyzed in
search of instabilities, asymmetries, and two-dimensional effects that can introduce errors when the data
are compared with results of 1-D simulations. The modeling involved transient axisymmetric numerical
simulations along with bifurcation analysis, which revealed that when the flow field is confined between

5
walls, local bifurcation occurs, which in turn results in asymmetry, deviation from the one-dimensional
assumption, and sensitivity of the flow field structure to boundary conditions and surrounding geometry.
Particle image velocimetry was utilized and results revealed that for jets of equal momenta at low
Reynolds numbers of the order of 300, the flow field is asymmetric with respect to the middle plane
between the nozzles even in the absence of confining walls. The asymmetry was traced to the asymmetric
nozzle exit velocity profiles caused by unavoidable imperfections in the nozzle assembly. The asymmetry
was not detectable at high Reynolds numbers of the order of 1000 due to the reduced sensitivity of the
flow field to boundary conditions. The cases investigated computationally covered a wide range of
Reynolds numbers to identify designs that are minimally affected by errors in the experimental procedures
or manufacturing imperfections and the simulations results were used to identify conditions that best
conform to the assumptions of quasi one-dimensional modeling.

1.4. End gas ignition by spherically expanding flames in a constant volume chamber
Most of the understanding gained on ignition phenomena is based on assumption of homogeneity in
the systems. Shock tubes and rapid compression machines (RCM) are usually assumed to be homogenous
(e.g., [18-20]). In these facilities, ideally there are no spatial gradients therefore the differential equations
describing the state of the gas could be marched only in the temporal domain. The main motivation for
reducing the dimensionality of the problem is the large size of reaction mechanisms and species involved
which requires solving simultaneously thousands of coupled, nonlinear differential equations, at numerous
grid points in space and time in which obtaining a realizable solution is a daunting task due to
nonlinearities [13] and their accompanying wide spectrum of time-scales [3].
It should be noted that homogenous ignition refers to simultaneous reaction of a mixture that is
spatially under uniform thermodynamic conditions. However, ignition in realistic systems usually occurs
in the presence of inhomogeneities even in well-established homogenous experiments like the shock tubes
and RCMs. In most cases, heat release is localized in space followed by the development of different
modes of propagating reacting fronts (e.g., [21-32]). Griffiths et al. [32] discovered that for fuel/air
mixtures that exhibit the negative temperature coefficient (NTC) behavior, ignition can occur at a location
that has a lower temperature relative to the surrounding. Numerical simulation performed utilizing a
turbulence model and a reduced chemical scheme to model ignition in a RCM demonstrated that the
reactivity originated close to the wall where the temperature was lower compared to the adiabatic core
[33,34]. Subsequent studies were aimed at characterizing temperature fields in RCMs and demonstrated
that temperature inhomogeneities that originated from BL-vortex interaction were leveled out during the
first stage heat release for fuel/air mixtures that displayed NTC behavior [29,34]. More recently it has

6
been shown that for conditions that can be encountered in these “legacy experiments”, the position of
ignition highly depends on the fuel type and thermodynamic conditions [35].
Since shock tubes and RCMs are key experimental facilities to obtain data for kinetic model
development and validation and considering the recent findings that ignition in these facilities can be
inhomogeneous as well as in experiments designed to provide insight into the physics of the knock (e.g.,
[24,25]), a novel approach is needed to perform homogenous ignition studies that is devoid of mentioned
complexities.
The explosion of adiabatically compressed combustible gas mixtures is one of the simplest and most
important combustion processes that can be utilized to study homogenous ignition (e.g., [36]). During the
flame propagation in a constant volume spherical chamber, the unburned gas is compressed resulting in
an increase in its temperature and pressure as shown schematically in Fig. 1.3.

7
(a)

(b)

Fig. 1.3. Schematic of the compression of the unburned end-gas by flame propagation in a spherical
chamber at different times after initiation of the flame, Physics are represented by the temperature maps
in the chamber; (a) 10 ms; (b) 25 ms.

If this increase in pressure and temperature in the unburned gas is sufficient to generate conditions that
the time scale of chemical reaction in the unburned gas becomes comparable to that of the flame
propagation, auto-ignition of the unburned gas can occur. The viability of this approach as a mean to study
ignition propensity of combustible mixtures and to validate chemical kinetic models has been
demonstrated first by Keck et al. [36] and has been used in subsequent studies using spherical chamber
(e.g., [37-39]), constant volume pancake type chamber [40] and knocking SI engines (e.g., [41-43]).
C o mp r e ss i o n

8
Hu et al. [37] utilized the pressure time history measured in the chamber to solve the species and
temperature evolution equations of the end gas numerically with a realistic chemical mechanism.
Although numerical determination of the state of the explosion was based on only an order of magnitude
analysis between the sensible energy release and isentropic compression rate, the agreement between the
zero-dimensional model and measurements was remarkable. However due to poor pressure readout
accuracy and sensitivity of data to the absolute temperature determination, the uncertainty of the reported
data was significant.
Although the readout accuracy from oscilloscopes window was later eliminated by incorporating
more advanced signal processing techniques [39], modeling using full set of chemistry is still to come. SI
engines community has successfully used realistic chemical mechanism to determine the conditions of the
end gas during knock but the interpretation of the results is not straightforward due to stochastic events,
turbulence uncertainty, and cycle to cycle variations (e.g., [41-43]).
The adiabatically compressed mixture by a spherically expanding flame (SEF) in a constant volume
spherical chamber emerges as a preferred candidate to study ignition owing to the simplicity of the
experiment. Experiments are performed in a totally spherical chamber. The chamber was successfully
used to measure the flame speeds in the compression region (e.g., [44]). The typical schematic of the
experiment is shown in Fig. 1.4.

Fig. 1.4. Typical Schematic of the spherical chamber experiments.

9
The entire chamber is located in an oven which is heated by resistance heaters mounted on the oven
wall. The oven temperature is controlled by closed loop temperature controller system capable of keeping
the oven at the desired temperature. The chamber is first filled with a combustible mixture using the partial
pressure method (e.g., [39,44]). The mixture is then ignited at the center by a spark between two extended
radial electrodes. Since there is no optical access, the only observable is the pressure time history. At low
initial temperatures and pressures, the mixture burns smoothly to completion as shown schematically by
the pressure vs. time trace in Fig. 1.5a. At high initial pressures and temperatures, the situation is markedly
different and an abrupt increase in pressure followed by large oscillations occurs prior to completion of
burning as shown in Fig. 1.5b. The pressure oscillation is an indication of auto-ignition of the end gas as
demonstrated in constant volume chambers [36-39] and SI engines operating under knocking conditions
(e.g., [41-43]).

(a)

(b)
Fig. 1.5. Typical pressure trace in the constant volume chamber; (a) normal flame propagation; (b)
auto-ignition.

This system unlike others is devoid of fluid mechanic interferences like vortices introduced by a
moving piston as observed in RCMs [45] and shock bifurcations [33] as in shock tube experiments thus
creating a “clean” environment to study ignition. Thanks to new technologies for measuring partial
pressures accurately and the availability of sophisticated uncertainty analysis [44], the initial and final
state of the unburned gas before ignition can be determined accurately to distinguish between the ignition
propensity of various fuels. Moreover, due to improved computational power, full set of reaction
mechanisms with variable rates as a function of pressure can be used to exactly model the evolution of the
end gas during compression.
Furthermore, since the flame is compressing the reacting end gas, the flame is essentially propagating
into a reacting mixture. Hence, the physics of flame propagation can be studied in the reactive

10
environment. This is relevant to spark ignition engines conditions in which considerable reactivity may
be present in the unburned mixture, introducing thus challenges due to couplings of auto-ignition and
flame propagation phenomena. Therefore, another objective of this thesis was to study SEF propagation
in a constant volume spherical chamber both in terms of feasibility and accuracy of performing end gas
ignition experiments.
Due to rather simple physics in constant volume chamber experiment, Direct Numerical Simulation
(DNS) can be performed to study the evolution of the whole mixture as the flame propagates through the
spherical chamber [44,46]. This is particularly useful to study the BL effects due to relatively colder
chamber wall [35], considering the mentioned effects of inhomogeneous ignition for fuels having NTC
behavior. Moreover, due to presence of the steep temperature gradients near the wall, the concentrations
of mixture components can exhibit stratifications near the wall for reacting mixtures that have contrasting
molecular weights due to Ludwig-Soret (L-S) diffusion [47]. When such effects are significant, the
explosion limits determined by DNS can be different from the zero-dimensional models predictions [35].
However, when the thermal BL effects are not present, the pressure-time history generated by DNS is
expected to be recovered by the rather simpler any relatively low cost models which only use
thermodynamics to determine the frozen end gas pressure and temperature (e.g., [48]). As a result, the
predicted pressure trace obtained from these simplified models can be successfully used to predict the
explosion limits using zero-dimensional models that predict reactivity evolution in the end-gas.
Chapter 3 contains the results of numerical simulations for the propagation of transient, 1-D laminar
flames into a reacting unburned mixture. The purpose of the investigation was to identify the parameters
influencing the flame burning rate in the conduction-reaction controlled regime at constant pressure. It
was found that the fuel chemical classification significantly influences the burning rate. More specifically,
for hydrogen flames, the “evolution” of the burning rate does not depend on the initial unburned mixture
temperature. On the other hand, for n-heptane flames that exhibit low temperature chemistry, the burning
rate depends on the instantaneous temperature and composition of the unburned mixture in a coupled way.
A new approach was developed allowing for the decoupling the flame chemistry from the ignition
dynamics as well as for the decoupling of parameters influencing the burning rate, so that meaningful
sensitivity analysis could be performed. It was determined that the burning rate is not directly affected by
fuel specific reactions even in the presence of low temperature chemistry whose effect is indirect through
the modification of the reactants composition entering the flame. The controlling parameters include but
not limited to mixture conductivity, enthalpy, and the species composition evolution in the unburned
mixture.

11
1.5. References
[1] K. Wu, I. Storey, Energy security in China's capitalist transition: Import dependence, oil diplomacy,
and security imperatives, London, 2008.
[2] U.S Energy Information Administration International Energy Outlook 2014 DOE/EIA-0484 (2014).
[3] C.F. Curtiss, J. Hirschfelder, Integration of stiff equations, Proc. Nat. Acad. Sci 38 (1951) 235-243.
[4] N.N. Semenov, Chemical kinetics and chain reactions, The Clarendon Press, Oxford, 1935.
[5] J.H. Van't Hoff, Etudes de Dynamique Chimique, 1884, p. 16l.
[6] J. Taffanel, C. Le Floch, Compt. Rend., (1913) 156-1544.
[7] D.A. Frank-Kamenetskii, Diffusion and heat transfer in chemical kinetics, Plenum, New York, 1969.
[8] B. Lewis, G. von Elbe, Combustion, flames, and explosion of gases, Academic Press, New York,
1987.
[9] C. G. Fotache, Ignition in convective-diffusive systems, Ph.D. thesis, 1997
[10] C.K. Law, Combustion physics, Cambridge University Press, 2008.
[11] U. Niemann, K. Seshadri, F.A. Williams, Accuracies of laminar counterflow flame experiments,
Combust. Flame 162 (2015) 1540-1549.
[12] J. Carpioa, A. Liñán, A.L. Sánchez, F.A. Williams, Aerodynamics of Axisymmetric Counterflowing
Jets, Combust. Flame (177) (2017) 137-143.
[13] S.T. Strogatz, Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry, and
engineering, Westview Press, 2000, pp. 44-93.
[14] O. Park, P.S. Veloo, H. Burbano, F.N. Egolfopoulos, Studies of premixed and non-premixed
hydrogen flames, Combust. Flame 162 (2015) 1078-1094.
[15] S.H. Won, B. Jiang, P. Diévart, C.H. Sohn, Y. Ju, Self-sustaining n-heptane cool diffusion flames
activated by ozone, Proc. Combust. Inst. 35 (2015) 881-888.
[16] R.R Burrell, Studies of Methane counterflow flames at low pressures, Ph.D. thesis, Mechanical
Engineering, University of Southern California, Los Angeles, California 2017.
[17] R.J. Kee, J.A. Miller, G.H. Evans, G. Dixon-Lewis, A computational model of the structure and
extinction of trained, opposed flow, premixed methane-air flames, Symp. (Int.) Combust. 22 (1988)
1479–1494.
[18] J.F. Griffiths, S.K. Scott, Thermokinetic interactions: Fundamentals of spontaneous ignition and
cool flames, Prog. Energy Combust. Sci., 13 (1987) 161-197.
[19] W. Jost, Explosion and combustion processes in gases, McGraw-Hill, New York, 1946.
[20] J.F. Griffiths, Reduced kinetic models and their applications to practical combustion systems, Prog.
Energy Combust. Sci. 21 (1995) 25-107.

12
[21] V.V. Voevodsky, R.I. Soloukhin, On the mechanism and explosion limits of hydrogen-oxygen chain
self-ignition in shock waves, Symp. (Int.) Combust. 10 (1965) 279-283.
[22] J.W. Meyer, A.K. Oppenheim, On the shuck-induced ignition of explosive gases, Symp. (Int.)
Combust. 13 (1971) 1153-1164.
[23] R. Blumenthal, K. Fieweger, K.H. Komp, G. Adomeit, Gas dynamic features of self ignition of non
diluted fuel/air mixtures at high pressures, Combust. Sci. Tech. 113 (1996) 137-166.
[24] R. Schieβl, U. Maas, Analysis of endgas temperature fluctuations in an SI engine by laser-induced
fluorescence, Combust. Flame 133 (2003) 19-27.
[25] N. Kawahara, E. Tomita, Y. Sakata, Auto-ignited kernels during knoking combustion in a spark-
ignition engine, Proc. Combust. Inst. 31 (2007) 2999-3006.
[26] S.M. Walton, X. He, B.T. Zigler, M.S. Wooldridge, A. Atreya, An experimental investigation of iso-
octane phenomena, Combust. Flame 150 (2007) 246-262.
[27] X. Gu, D.R. Emerson, D. Bradley, Modes of reaction front propagation from hot spots, Combust.
Flame 133 (2003) 63-74.
[28] P. Dai, Z. Chen, Supersonic reaction front propagation initiated by a hot spot in n-heptane/air mixture
with multistage ignition, Combust. Flame 162 (2015) 4183-4193.
[29] J. Clarkson, J.F. Griffiths, J.P. MacNamara, B.J. Whitaker, Temperature fields during the
development of combustion in a rapid compression machine, Combust. Flame 125 (2001) 1162-
1175.
[30] G. Mittal, M.P. Raju, C.J. Sung, CFD modeling of two-stage ignition in a rapid compression
machine: Assessment of zero-dimensional approach, Combust. Flame 157 (2010) 1316-1324.
[31] R. A. Strehlow, A. Cohen, Limitations of the reflected shock technique for studying fast chemical
reactions and its application to the observation of relaxation in nitrogen and oxygen, J. Chem. Phys.
30 (1959) 257-265.
[32] J.F. Griffiths, D.J. Rose, M. Screiber, J. Meyer, K.F. Knoche, Novel features of end-gas autoignition
revealed by computational fluid dynamics, Combust. Flame 91 (1992) 209-212.
[33] K.P. Grogan, M. Ihme, Weak and strong ignition of hydrogen/oxygen mixtures in shock-tube
systems, Proc. Combust. Inst. 35 (2015) 2181-2189.
[34] K. Nishiwaki, Y. Yoshihara, K. Saijyo, Numerical analysis of the location of knock initiation in SI
engines, SAE Technical Paper, (2000), paper 2000-01-1897.
[35] J. Jayachandran, F.N. Egolfopoulos, Thermal and Ludwig–Soret diffusion effects on near-boundary
ignition behavior of reacting mixtures, Proc. Combust. Inst. 36 (2017) 1505-1511.

13
[36] J.C. Keck, H. Hu, Explosions of adiabatically compressed gases in a constant volume bomb, Symp
(Intl.) Combust. 21 (1986) 521-529.
[37] H. Hu, J.C. Keck, Autoignition of adiabatically compressed combustible gas mixtures. Society of
Automotive Engineers (1987).
[38] H. Hu, The autoignition of hydrocarbon-oxygen-diluent mixtures in a constant volume bomb, Ph.D
Thesis, 1986.
[39] A. Moghaddas, C. Bennett, K. Eisazadeh-Far, H. Metghalchi, Measurement of Laminar Burning
Speeds and Determination of Onset of Auto-Ignition of Jet-A/Air and Jet Propellant-8/Air Mixtures
in a Constant Volume Spherical chamber, J. Energy Resour. Technol 134 (2012) 022205.
[40] Y. Moriyoshi, S. Kobayashi, Y. Enomoto, Analysis of knock phenomenon induced in a constant
volume chamber by local gas temperature measurement and visualization, JSME International
Journal Series B 48 (2005) 695-700.
[41] F. Bozza, D. Siano, M. Costa, Heat Transfer, Knock Modeling and Cyclic Variability in a
“Downsized” Spark-Ignition Turbocharged Engine, Advances in Applied Mathematics and
Mechanics 3 (2011) 310-326.
[42] F. Bozza, G. Fontana, E. Galloni, E. Torella, 3D-1D analyses of the turbulent flow field, burning
speed and knock occurrence in a turbocharged SI engine, SAE Technical Paper (2007), paper 2007-
24-0029.
[43] H. Schapertons, W. Lee. Multidimensional modelling of knocking combustion in SI engines. No.
Society of Automotive Engineers (1985) paper TP-850502.
[44] C. Xiouris, T. Ye, J. Jayachandran, F.N. Egolfopoulos, Laminar flame speeds under engine-relevant
conditions: Uncertainty quantification and minimization in spherically expanding flame
experiments, Combust. Flame 163 (2016) 270-283.
[45] K. B. Mansfield, M. S. Wooldridge, H. Di, X. He, Low-temperature ignition behavior of iso-octane,
Fuel 139 (2015) 79-86.
[46] J. Jayachandran, R. Zhao, Determination of laminar flame speeds using stagnation and spherically
expanding flames: molecular transport and radiation effects, F.N. Egolfopoulos, Combust. Flame
(2014) 2305-2316.
[47] D.E. Rosner, Transport Processes in Chemically Reacting Flow Systems, Dover Publications, 2000.
[48] D. Bradley, A. Mitcheson, Mathematical solutions for explosions in spherical vessels, Combust.
Flame 26 (1976) 201-217.

14
Chapter 2
Flame ignition in the counterflow configuration: Reassessing the
experimental assumptions
2.1. Introduction
The counterflow configuration is in principle a well-controlled environment in which global flame
properties such as ignition temperatures, laminar flame speeds, and extinction strain rates can be measured
and used as validation targets for kinetic models under the assumption that the flow-field is (quasi) one-
dimensional (1-D). If the axial velocity and scalars exhibit zero radial gradients along the axis of
symmetry, then a similarity solution exist and as a result, only the axis of symmetry needs to be modeled
[1] reducing thus the computational cost and removing ambiguities associated with boundary conditions
(BC) that need to be characterized in multi-dimensional configurations (e.g., [2]).
The validity of the 1-D assumption has been assessed in several studies both in the counterflow
configuration (e.g., [3-11]) and stagnation flows (e.g., [12-17]) for conditions that could be encountered
in experiments. Chelliah et al. [3] demonstrated that the requirement of zero radial gradient of the axial
velocity in the counterflow configuration can be relaxed by implementing the measured value of this
quantity as a BC. It has been shown that parabolic burner exit velocity profiles could result in notable
radial non-uniformities of velocities and scalars in the flow-field [4,5,7]. For non-premixed flames, it has
been shown [8] that for an M-shape burner exit velocity profile, the stagnation plane and the stoichiometric
surface exhibit notable curvature at large radii. As a result, heat release, scalar dissipation rates, and
temperature peak away from the centerline, which apparently violates the 1-D assumptions. Mittal et al.
[9] demonstrated that a non-uniform M-shape burner exit velocity profile, caused by a dip in the centerline,
induces a temperature curvature that has a significant effect on laminar flame speeds obtained through
extrapolations.
Niemann et al. [10] demonstrated that careful design of large diameter straight tube burners with flow
straightening screens, can produce velocity profiles in counterflowing non-reacting cold air jets, which
satisfy the 1-D assumptions within an accuracy of 5%. Johnson et al. [11] showed computationally that
the burner-design proposed by Niemann et al. [10] exhibit negligible two-dimensionality for non-
premixed H2 flames by thoroughly comparing the terms in momentum and energy equations generated by
1-D and 2-D approaches. The required design considerations for appropriately reproducing one-
dimensionality using contoured-nozzles have also been identified in the same study [11].

15
Bergthorson et al. [12] reported that 1-D formulation assumptions would be guaranteed in premixed
flames in the stagnation flow configuration for properly designed contoured-nozzles [13] if the measured
velocity BC’s at suitably chosen point are implemented in the simulations. Success of 1-D formulation in
such designs was further verified by comparing to the measured velocity, temperature and NO
concentration profiles [13,17]. Furthermore, Bergthorson et al. [16] verified the validity of the 1-D
formulation in non-reacting stagnation flows by applying the approach of Chelliah et al [3]. On the other
hand, Bouvet et al. [15] reported that due to variations of pressure curvature along the centerline, 1-D
models could not reproduce the measured velocity profiles in stagnation flames. However, the
conclusions for the non-reacting cases were flawed since the 1-D model was used outside its domain of
applicability [11,16,18-20].
While the aforementioned studies on potential multi-dimensional effects are meritorious, they have
been carried out in the presence of a flame with some exceptions ([7,10,16]). On the other hand, when
the counterflow configuration is used to study flame ignition, the physical and chemical processes that
lead to ignition evolve in the absence of a flame and are controlled by fluid mechanics, heat and mass
transfer to and from the ignition kernel whose temperature is modest compared to a flame, and the
relatively slow evolution of the radical pools. Additionally, the presence of one heated jet could add to
the complexity of the flow field. Studies of multi-dimensional effects on the ignition of counterflow
flames have yet to appear. Furthermore, it is generally accepted that stretch dampens any disturbances
resulting thus in a steady flow field (e.g., [10]).
However, what may be less known is that while the stagnation flow is a stable configuration [21], the
counterflow configuration in the absence of a flame, is in principle hydrodynamically unstable resulting
in two stable states for the same BC’s (e.g., [22-25]). More specifically, the symmetry of the flow field
could break, and for identical flow BC’s, the stagnation plane could be established at two different
locations away from the symmetry plane between the nozzles. This asymmetric behavior is observed
when a critical Reynolds number based on nozzle diameter (Recrit) is reached (e.g., [23,25]). More
specifically, experiments and simulations suggest that even when the momenta of two opposing jets are
equal, the resulting flow field may be asymmetric and exhibit multiple steady states. The transition from
a single steady state to “bi-stability” is an indication of “pitchfork bifurcation” and has been introduced in
the analysis of dynamical systems (e.g., [26]); this is a manifestation of the inherent non-linearity of the
Navier-Stokes equations. Figure 2.1 depicts the schematic of a typical stability diagram for counterflow
jets indicating that beyond the Recrit value the single stable solution bifurcates (e.g., [23,25]).

16

Fig. 2.1. Schematic of a typical stability diagram for counterflowing jets [23,25]. The two stable solutions
are distinct after the bifurcation.
In bifurcation studies, the state at which the stagnation point is closer to the lower burner is denoted
as the “lower branch” while the other state is called the “upper branch” (e.g., [23,25]). Since 1-D codes
simulate only the symmetric branch, after the Recrit value has been reached, the fidelity of computed results
based on 1-D assumptions must be assessed.
Physically, this bifurcation is partially attributed to the Coanda effect, which is the tendency of the
flow stream to attach to adjacent surfaces (e.g., [27]). At high enough velocities, the flow exiting the
nozzles can eventually separate and create recirculation zone(s) next to any solid surface. These
recirculation zones create radial velocity profiles that satisfy the necessary condition to develop local
instability according to Rayleigh and Fjortoft criterion as shown in classical shear flow studies (e.g.,
[28,29]).
On the other hand, Rolon et al. [22] showed experimentally that for counterflowing ambient-
temperature air jets, the stagnation plane is not located in the middle plane between the nozzles even when
all surfaces surrounding the nozzle exits and resulting recirculation zones are removed. It was further
shown, that the offset of the stagnation plane location increases with the nozzle separation distance (L).
In counterflow flame experiments, flow-surface interactions are frequently unavoidable. First, the two
impinging jets leave the nozzles that have rims with a finite thickness. Second, cooling jackets and coils
are used to protect the nozzles from thermal stresses (e.g., [8,30]), and insulation jackets are an essential
part of the counterflow ignition experiments in which high temperatures need to be reached at the exit of
one nozzle (e.g., [31]). Finally, the experiments are performed typically in closed chambers in order to
minimize ambient disturbances, establish inert environment for safety reasons, and reach pressures of
relevance to applications (e.g., [32-34]).
Upper Branch
Lower Branch
Re
crit
Stagnation Point Position

17
Ciani et al. [25] reported a significant hysteresis between the two stable solutions in counterflowing
air jets. It was shown that once the flow establishes itself at one state, it could not be switched to the other
unless the momentum at the exit of the nozzle that is closer to the stagnation plane increases by a certain
amount. The excess amount of momentum required for the hysteretic transition was found to increase
with Re.
Bifurcation of non-isothermal counterflowing jets has been reported also. Conditions that result in
steady [35] and unsteady behavior of the flow-field have been identified [24,36,37]. Those conditions
were shown to depend on the velocity and temperature at the nozzle exit as well as L.
Stability analysis of counterflowing non-reacting jets has been performed for a wide range of L values.
The appropriate geometric parameter commonly used is the aspect ratio   L/D, where 𝐷 is the jet
diameter. Nevertheless, the range where flame studies are being performed is 0.3 <  < 2, i.e. outside the
“free-floating” regime (e.g., [18-20]). In this regime, only stable and bi-stable states have been reported
in axisymmetric configuration when the top jet has lower density, i.e. stable stratification [8,9,23-25,31-
37].
Finally, it should be noted that in the presence of a flame, the heat release is intense resulting in high
temperatures and thus notable buoyancy-driven convection at large radii, which has been shown to
suppress the tendency of the flow-field to bifurcate [37]. This could be the reason that in the combustion
literature, the stability of counterflowing jets has attracted less attention despite its wide use. To the
authors’ knowledge, the only bi-stability of reacting flows in counterflow configuration was reported by
Ciani et al. [38] who observed a transition from a planar diffusion flame to an edge flame during extinction
where the behavior was attributed to flow field exchange of stability from one steady state to another.
Based on these considerations and the possible existence of bi-stable states in counterflow flame
ignition experiments in which a flame is not present, the main goal of the present chapter was to investigate
experimentally and computationally multi-dimensional effects and assess the validity of the 1-D model.
Specifically, the present investigation aimed to quantify the uncertainty that is introduced by bifurcation,
thermal boundary layers, and actual experimental geometry and boundary conditions. It should be noted
that such effects on the flow field have not been addressed adequately in previous studies, which have
been largely carried out in the absence of temperature gradients that are unavoidable in counterflow flame
ignition experiments.

18
2.2. Experimental approach
2.2.1. Method and apparatus
Experiments were carried out under atmospheric pressure in a counterflow burner configuration shown
schematically in Fig. 2.2. The apparatus consists of a straight upper quartz burner, from which a high
temperature jet is directed downward against an ambient temperature H 2/N2 jet that exits from the lower
straight burner; for both burners D = 14.5 mm. Air was used as the hot oxidizer to study the ignition of
non-premixed H2 flames. Coflowing N2 streams with an annulus gap of 3 mm surround both jets. All
gaseous flow rates were metered using sonic nozzles by controlling the pressure upstream of the orifice.

Fig. 2.2. Schematic of the counterflow ignition burner assembly.

The design of the upper burner that heats the air differs from that used in previous studies [31-34,39-
41]. More specifically, only external resistance heating wire (Kanthal A1) was used, and in order to obtain
temperatures at the center of upper burner exit, Tair, up to 1400 K, the length of the upper nozzle had to be
chosen significantly larger than conventional quartz burners that use internal heating elements. The burner
was an 85 cm long tube of which two segments were heated, each being 37 cm long. The power supplied
to each heated segment was controlled separately using variable AC transformers. Two bare B-type
thermocouples (Omega P30 R-008) were mounted on the outside surface of the tube at the end of each
D
L
δ
Cooling
Jacket
Air External
Thermocouple
Coflow(N
2
)
Coflow(N
2
)
Coflow(N
2
)
Coflow(N
2
)
Screen
Insulation
Wrap
Heating wire
Air
Heater Surface Thermocouple
Fuel Stream(Fuel/N
2
)

19
heated segment to control the temperature of the segment manually. The entire upper burner assembly
was wrapped with an alumina silica blanket to minimize heat loss.
Temperature measurements were carried out using an R-type thermocouple (Omega P13 R-002) with
a bead diameter of 120 ± 25 m determined by electron microscopy. Details of the thermocouple
assembly and measurements are elaborated in Appendix C.
The lower burner was made of stainless steel mounted onto a multistage translational stand capable of
translating about multiple axes in order to obtain desired alignment and separation distance between the
burners. L = D = 14.5 cm was chosen for all experiments (  = 1). A cooling coil was wrapped around
the lower burner to prevent any temperature rise due to radiation emitted from the heated quartz burner
and assure that the H2/N2 jet is maintained uniformly at an unburned temperature Tu = 298 K. Two fine-
wire screens (200 mesh) were placed in the lower burner as flow straighteners; one 8D and the other 10D
upstream of the burner exit. As a result, the lower burner exit velocity profile is of top-hat shape. The
measured burner exit profiles are shown in Appendix B.
Another major difference between the current design and the previous generation of burners is the fact
that no flow straightener was used in the upper nozzle. The burner is long enough for the flow field to
become fully developed both thermally and hydrodynamically at its exit. Note that temperature and axial
velocity profiles at the heated burner exit were reported to be parabolic when heated even when flow
straighteners are used [33,40,41]. This is largely due to the higher thermal diffusivity and kinematic
viscosity at high temperatures that result in thicker boundary layers. Given the non-uniformities at the
burner exits, balancing the momenta between the two jets is involved and the details of the process
followed in this study can be found in Appendix A.
Not using flow straighteners in the heated burner, several problems are eliminated. These include: (1)
agglomeration of seeding particles in the colder mesh caused by thermophoresis; (2) loss of flow
alignment due to thermal expansion and subsequent displacement of the mesh; (3) mechanical failure of
the quartz tube due to thermal stresses caused by difference between the thermal expansion coefficients
at the steel-quartz junction/interface; and (4) uncertainties caused by accounting for radiation correction
performed on the measured temperature due to glowing mesh at high temperatures.
The system was allowed to reach a steady temperature profile before each experiment. Upon
establishing the appropriate flow rates and temperatures at the two burner exits, the fuel flow rate was
gradually increased until ignition was observed (e.g., [31]). This approach was preferred over raising Tair
until ignition happens, as the uncertainty in the fuel side flow rate is notably lower compared to that of
Tair (e.g., [33,41]). Moreover, if the supplied power to the heating wire is changed in an attempt to raise

20
Tair, it was observed that the whole system has to be allowed to reach a new thermal equilibrium which
can take up to twenty minutes due to the high heat capacity of the heated burner.
The measured Tair at the state of ignition was defined as the ignition temperature (Tign) (e.g., [31]),
once corrected for radiative/convective heat transfer to/from the bead [33]. The uncertainty of the final
corrected temperature was determined to be less than ±35 K at the maximum achievable Tair incorporating
the methodology of Brady et al. [41]; the uncertainty was determined to be mostly due to the uncertainty
of thermocouple bead diameter measurements and its effect in radiation correction.

2.2.2. Flow and flame visualization
The flow was seeded with Al2O3 particles with a nominal diameter of 3 μm, and particle image
velocimetry (PIV) was used to measure the 2-D flow field on a plane that passes through the system
centerline. The PIV details can be found in previous studies by the authors (e.g., [30,43]). The axial flow
velocities along the stagnation streamline and the local strain rate, K, were derived based on the maximum
absolute value of the axial velocity gradient on the fuel-side [31]. The uncertainty of the PIV
measurements were determined by the methodology of Sarnacki et al. [20] and was found to be of the
order of ±0.04 m/s while the 95% confidence limit of K is ±30 s.
-1

High-speed recording of chemiluminescence from the flame along with shadowgraph imaging were
performed as well, in order to further evaluate the qualitative flow conditions during ignition and the
subsequent edge flame propagation (e.g., [44]). A direct Shadowgraph ensemble has been set up to
visualize the stagnation plane and the flame as shown in Fig. 2.3.

Fig. 2.3. Top view schematic of the visualization configuration.

Light was gathered from 13 mW epoxy encased white LED and focused on a pinhole to generate a
point source light. The diameter of the pinhole/aperture is adjustable and can be optimized to obtain the
best illumination. Light rays were further parallelized inside the test section using a series of spherical
lenses with a diameter of 2 inches and focused again at the camera lens after passing through the
established flow field. A high-speed charge-coupled device (CCD) camera (X-Sream, XS-3) was placed
very close to the focal point of the last lens. A frame rate of 500 Hz with the resolution of 1258  1280
Light
Source
CCD Camera
Flat Mirror
Pinhole

21
pixels and an exposure time of 2 ms was utilized to obtain the best balance between contrast and
resolution. This method was particularly efficient in determining 𝐿 and aligning the burners with a spatial
resolution better than 0.15 mm.
Flame luminosity was recorded to monitor the ignition location via the CCD camera with the same
frame rate and exposure of shadowgraph experiments. The fuel stream was doped with a trace amount of
methane to make the flame visible and with negligible effect on the reactivity of the system (e.g., [33,45]).
Experiments were repeated without methane, and ignition was confirmed at the same conditions using
shadowgraph. In order to verify that the ignition is taking place at the centerline, ignition had to be
recorded at two different angles. Therefore, a flat mirror was utilized to reflect the ignition event at an
angle with respect to the original position of the camera as shown in Fig. 2.3. Only after shots were taken
from these two views multiple times, the actual position of ignition could be found.
Tign was first determined for each experiment, following which the flow field was established at fuel
flow rates 3% below the ignition point at Tign where recording commenced. Only after triggering the
camera, the fuel flow rate was increased to the desired value for ignition. Five records were taken at the
ignition point to assure repeatability and quantify the uncertainty of the ignition position (± 3 mm).
It should be noted that PIV measurements, shadowgraph and flame visualization were performed
independently. The optical system for PIV measurements was mounted such that the laser light sheet was
perpendicular to the light arrays of shadowgraph.

22
2.3. Modeling approach
Transient numerical simulations of the 2-D axisymmetric problem involving coupled fluid dynamics,
detailed chemistry, and transport are carried out using the laminarSMOKE code [46], which is an open-
source computational package built on the OpenFOAM [47,48] platform and details are given in Appendix
B. The laminarSMOKE code has been specifically validated for counterflow flames [49] and also used
to simulate other laminar combustion phenomena in axisymmetric configuration such as low-pressure
flame sampling [50]. Numerical simulations were performed for atmospheric non-premixed H 2 flames at
pre-ignition state with a H2 mole fraction XF = 0.2 in the fuel jet. The H2/O2 sub-model of USC Mech-II
[51] was utilized to describe the kinetics. First, (quasi) 1-D simulations were performed to estimate Tign
for given XF and K using an opposed-jet code [1,52-54], and the maximum Tair for the 2-D simulations
was kept 20 K lower than the estimated Tign to make sure that ignition would not occur during the transient
simulations. This is a safe assumption since the conditions chosen correspond to the second explosion
limit of H2 where Tign is not sensitive to K [33].
In order to investigate the effects of the burner-flow-walls interactions, three axisymmetric geometries
were studied as depicted in Fig. 2.4. Geometry (a) corresponds to the experimental configuration used in
this study and involves realistic velocity and temperature BC’s, which are parabolic for the heated burner
and nearly top-hat for the lower burner as shown in Appendix B. Since parabolic jet exit profiles in
geometry (a) introduce undesirable 2-D complexities [4,5,7], it might be meaningless to compare the 1-D
and 2-D results. Hence, the simulations were repeated for geometry (b), which have BC’s with minimum
radial non-uniformities complying with the assumptions of 1-D formulation [1,10]. Solid surfaces
adjacent to coflow exits i.e., the insulation wrap for upper burner and the cooling jacket for lower burner
as shown in Fig. 2.2 are denoted as disk walls in Fig. 2.4a. Since it has been reported that recirculation
zones could be eliminated by removing disk walls and solid surfaces next to the jets exits [22], the
simulations were repeated for geometry (c) with the ideal BC’s of geometry (b) but with replacing disk
walls with open boundaries. Performing 2-D simulations for this geometry is particularly useful since the
mechanism leading to bifurcation is widely attributed to the recirculation zones [22,23,25,28].

23

Fig. 2.4. Schematics of the geometries used in 2-D simulations; Solid boundary lines represent walls; a)
geometry (a); b) geometry (c); c) geometry (b).

Table 2.1 summarizes the simulations performed and their correspondence and characteristics.
Additionally, BC’s, details of the various geometries, and detailed solution procedures are provided in
Appendix B.

24
Table 2.1. Summary of cases investigated computationally.
Geometry
2-D simulation
name
1-D simulation
Name
Boundary conditions for scalars, axial velocity
and its gradient
(a)
2-D, Upper Branch
1-D, Upper Branch Extracted from 2-D, Upper Branch at z = 0 & z = D
1-D (L = 0.4D),
Upper Branch
Extracted from 2-D, Upper Branch at z = 0.15D & z = 0.55D
(with domain size of L = 0.4D)
1-D (adjusted T),
Upper Branch
Extracted from 2-D, Upper Branch at z = 0 & z = D except
Tair = 917 K (value at z = 0.43D)
1-D (corrected),
Upper Branch
Extracted from 2-D, Upper Branch at z = 0 & z = D except
Tair = (Tair)eff (obtained from correction procedure)
2-D, Lower Branch 1-D, Lower Branch Extracted from 2-D, Lower Branch at z = 0 & z = D
(b)
2-D, geometry (b),
Upper Branch
1-D
Specified axial velocity and scalars, zero axial velocity
gradient
2-D, geometry (b),
Lower Branch
(c) 2-D, geometry (c)

To obtain both possible solution branches numerically in the presence of hysteresis an approach
similar to Ciani et al. [25] was incorporated; an initial non-reacting flow field was established in a coarse
grid at low enough Re at which only one stable state exist. At this state, the inlet velocities and
temperatures were fixed for the final desirable state but the kinematic viscosity in the momentum equation
was artificially increased to achieve a low 𝑅𝑒 solution. It was realized that the solution was already
asymmetric even at Re as low as 20, which is thought to be a result of profound influence of buoyancy at
low Re [8]. Then, Re was artificially increased in small steps by reducing the effective kinematic viscosity.
The existence of the other stable branch was subsequently tested as follows. If the original axisymmetric
solution was a lower (upper) branch then the flow field was perturbed by increasing (decreasing) the
average velocity of the fuel side to twice (half) of the original value and simulation was continued until
the flow establishes a new equilibrium. The new state can be categorized into one of the branches by
observing the separation bubble sizes at the upper and lower disk walls at exit of the jets; the bubble at the
lower disk wall for lower branch is always smaller than the corresponding one at upper wall (e.g., [25]).
The opposite is true for the upper branch. Subsequently, the fuel side velocity was adjusted to the value
at which the momenta of the two jets are balanced. If the transition was observed for the switched state,
extra care had to be taken to assure the solution would not jump to the initial branch. Re was then increased
further and the procedure was repeated until both stable solutions were obtained. Having the two

25
solutions, Re can be further increased to reach to the final state where refinement procedure is performed.
After refinement, the chemistry step solver was activated and the solution was tracked until steady state
was reached for all species in the domain.
2.4. Results and discussion
2.4.1. Bifurcation effects
All experiments were carried out for fixed Tair = 926 K and K = 360 s,
-1
and the ignition state was
characterized by the XF value that would lead to ignition. Ignition was found to occur consistently at two
different XF values namely XF = 0.18 and XF = 0.22, and the number of ignition occurrences at those two
XF values was found be approximately the same for more than 20 trials. The uncertainty of XF was
estimated using the methodology of Brady et al. [41] and it was determined to be less than 0.01. Therefore,
the observed differences of the XF values at ignition cannot be attributed to experimental uncertainty.
Shadowgraph images were taken at XF = 0.15, which corresponds to a pre-ignition state. For the same
BC’s, two stable states were found as shown in Fig. 2.5. The state of the flow-field was determined to be
susceptible to external disturbances. For example, inserting the thermocouple may cause the transition
between the two states. Also after extinguishing an ignited flame by turning off the fuel supply, the flow
could be in a state different from that before ignition. This is a further proof that the flow does indeed
switch between the aforementioned states. In the absence of the external disturbances, the two states are
obtainable by modifying the cold jet flow rate similar to the numerical procedure outlined in Section 2.3.

Fig. 2.5. Shadowgraph images of pre-ignition states showing two possible steady states for the same
boundary conditions; (a) Lower branch; (b) Upper branch.
The stagnation plane is expected to have a similar shape as the non-planar discontinuity shapes in the
light intensity shown in Fig. 2.5 and the stagnation point along the centerline is closer the lower burner
for the conditions of Fig. 2.5a. Therefore, the images shown in Figs. 2.5a and 2.5b correspond to the
lower and upper branches respectively; recall that fuel (air) is supplied by the lower (upper) burner.
For both the upper and lower branches, axial velocity profiles along the centerline were measured and
computed using 2-D simulations for geometry (a). The results are shown in Fig. 2.6 along with results
from 1-D simulations. The coordinate system was chosen to be at the center of lower burner exit and the

26
spatial distances were scaled by D, as outlined in Fig. 2.4. 2-D and original 1-D simulation results match
at the nozzles exits for both branches. However, the position of the stagnation plane and the adjacent
velocity field differ notably between the two formulations, such discrepancy might lead to misleading
conclusions if the experimental results are compared to the 1-D simulations.

Fig. 2.6. Comparison between measured and computed axial velocity profiles at the stagnation
streamline for geometry (a).
It has been reported that the 1-D approach can capture the velocity profiles along the centerline if BC’s
are measured at the burner exits and implemented in the simulations (e.g., [3,16,17,20,52]). Although,
this procedure was implemented (see Table 2.1) for the results shown in Fig. 2.6, the 1-D approach still
fails to reproduce the velocity field as determined by experiments and 2-D simulations. However, if BC’s
are extracted from the inner points of the domain (case denoted as “L = 0.4D” in Fig. 2.6), the velocity
profiles computed by 1-D and 2-D models match closely. Therefore, the domain of applicability of the 1-
D approach in predicting the velocity field is smaller than what has been previously reported [3,17-20].
Figure 2.7 depicts the measured and computed axial velocity radial profiles at the burners’ exits for
both branches. Given that the symmetry of the profiles with respect to the centerline was better than 94%
in all measurements, only half of the domain is shown; the lower branch velocities are plotted against
negative radius values to aid the comparison with the upper branch profiles. The lower branch profile
exhibits a larger dip at the lower burner exit because the stagnation point is closer to the lower burner; the
observed dip results from pressure feedback from the stagnation point (e.g., [7,20,55]).
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
U, m/s
z /D
PIV, Lower Branch
2-D, Lower Branch
1-D, Lower Branch
PIV, Upper Branch
2-D, Upper Branch
1-D, Upper Branch
1-D (L = 0.4D), Upper Branch

27

Fig. 2.7. Comparison between measured and computed radial profiles of axial velocity at the burners exits
for geometry (a); (▲) and (----) z = 0.93D, (●) and (─) z = 0.07D. Symbols: experimental data; Lines: 2-
D simulation results; Blue: Lower branch; Red: Upper branch.
Using high-speed shadowgraph to capture the details of the ignition process, it was found that ignition
at XF = 0.18 (0.22) occurs at the upper (lower) branch. The observed 0.04 difference in XF at which ignition
occurs is rather significant in terms of the reactivity of the ignition kernel. Figure 2.8 depicts the spatial
distributions of the hydrogen radical mole fractions ( 𝑋 ) along the centerline computed using 1-D
approach which shows that there is a 31% difference between the maximum 𝑋 values for the two
branches.

Fig. 2.8. Computed 𝑋 along the stagnation streamline of geometry (a) using 1-D approach; Red: Upper
Branch, XF = 0.18; Blue: Lower Branch, XF = 0.22.

It is also a common practice to perform the flame ignition experiment in the counterflow configuration
by keeping XF constant and increasing Tair until ignition is observed. To investigate if the observed
difference in ignition propensity of the two branches is sensitive to the method of achieving ignition, Tign
was measured for XF = 0.18 and XF = 0.22 with increasing Tair and it was observed that ignition occurs at
Tair = 926 K for both cases. Therefore, by changing Tair, the effects of bifurcation could not be captured.
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
r/D
U, m/s
U, m/s
r/D
0.E+00
2.E-12
4.E-12
6.E-12
8.E-12
1.E-11
0.0 0.2 0.4 0.6 0.8 1.0
X
H
z /D

28
The low sensitivity of Tign to XF is expected given that for the conditions considered in the present
investigation, ignition is not controlled by fuel diffusion [33]. Figure 2.9 shows the variation of Tign with
XF in which results of 1-D simulations are compared to the measurements. The 1-D results revealed that
there is only 4 K difference between Tign’s at XF = 0.18 and XF = 0.22.

Fig. 2.9. Variations of Tign as a function of XF for K = 360 s
-1
; (●) 1-D simulation results; (▲)
measurements for the upper branch; (■) measurements for the lower branch.

High-speed imaging of the flame luminosity revealed that ignition always initiates close to the
centerline for both branches as shown in Fig. 2.10a. Subsequently, the resulting edge flame propagates
radially through a path whose shape is similar to the stagnation plane before finally settling as a single-
state stable non-premixed flame shown in Fig. 2.10c. These observations suggest that the lower branch is
unstable in the presence of the flame, which is also evident from the transition of the lower branch flame
(Fig. 2.10b) to the upper stable flame (Fig. 2.10c). A possible explanation is that in the presence of a
flame, the lower branch flame becomes unstable due to buoyancy, which tends to push the stagnation
plane wings upward at larger radial distances away from the centerline [8,41].

Fig. 2.10. High-speed images of ignition for H2 non-premixed flame at different times after increasing 𝑋
to the ignition value for pre-ignition conditions corresponding to the lower branch; (a) ignition kernel
formation; (b) unstable lower branch flame; (c) stable upper branch flame.

It should be noted that by comparing the features of Figs. 2.5b and 2.10c, the flame appears to be more
planar compared to the shape of the light intensity discontinuity before ignition (Figs. 2.5b). This is due
910
920
930
940
950
960
970
980
0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28
T
ign
, K
X
F

29
to the fact that the flame changes the characteristics of the flow field due to the heat release and buoyancy.
Therefore, judging the ideality of the flow field before ignition by the flame topology upon ignition can
be misleading.
It was realized also that the reported bifurcation is sensitive to Tair and K. At low K and high Tair, only
the upper branch was attainable, while the lower branch was realized only for high K and low Tair values,
which could be caused by buoyancy [37]. Thus, given the wide range of Tair and K values that are
encountered in counterflow flame ignition experiments, it is not clear what flow conditions will be
prevailing at the pre-ignition states, and thus velocity measurements or shadowgraph visualization are
required.
2.4.2. Simulation results for geometry (a)
2-D simulations were performed for geometry (a) for XF = 0.2 and two stable states were found for the
same BC’s. Figure 2.11 depicts the temperature distribution along with the streamlines for the two stable
branches. For comparison, the results of the lower branch solution are shown on the left while the results
the results for upper branch are shown on the right, and they were separated by the stagnation streamline.
At the junction of the coflow exits and disk walls, the flow separates creating recirculation zones. The
upper separation bubble is smaller (bigger) than the lower one for the upper (lower) branch. The shape of
the light intensity discontinuity captured by shadowgraph and shown in Fig. 2.5 is reproduced qualitatively
by the temperature distribution in Fig. 2.11.

Fig. 2.11. Comparison between computed temperature maps superimposed by the streamlines for two
stable solutions of geometry (a); Left: lower branch; Right: upper branch.

The spatial distributions of 𝑋 and 𝑋 are shown in Fig. 2.12 for lower and upper branches. The
maximum 𝑋 occurs at the centerline indicating the higher ignition propensity at the centerline in
agreement with the results of Fig. 2.10a in which ignition was found experimentally to initiate at the
centerline. All other minor species and products also maximize at the centerline implying that the ignition
kernel is located at 𝑟 = 0 . However, the 𝑋 iso-contours curve at larger radii due to the parabolic nature

30
of the axial velocity profile exiting the upper burner and the effect of elliptic pressure field on the lower
burner exit. The most apparent difference between two branches is the shift of the ignition kernel due to
the displacement of the stagnation point.

Fig. 2.12. Computed 𝑿 𝑯 maps overlaid by 𝑿 𝐇 𝟐 iso-contours for the two stable solutions of geometry (a);
Left: lower branch; Right: upper branch.

Computed 𝑋 profiles along the stagnation streamline are shown in Fig. 2.13 to further investigate the
differences between the two solution branches. The maximum 𝑋 for the upper branch is 4% higher than
the corresponding value for the lower branch, which could explain why lower-branch ignition happens at
higher fuel concentrations experimentally. The 𝑋 profile computed using the 1-D opposed jet code are
shown also in Fig. 2.13 and the maximum value exceeds by nearly 2.33 times those predicted by the 2-D
simulations (for 1-D, Upper Branch case).
It should be noted also that central differencing schemes was used to discretize the 1-D domain as
order of discretization has shown to be important when comparing 1-D and 2-D results [11]. Moreover,
since comparisons of 1-D and 2-D results is only valid when rigorous quantification of discretization
uncertainty for 2-D results is provided [11], the error bar for maximum 𝑋 value (as explained in detail in
Appendix B) was quantified and shown in Fig. 2.13.

31

Fig. 2.13. Computed 𝑿 𝑯 along the stagnation streamline of geometry (a).

2.4.3. Effects of jet exit parabolic temperature profile
To further investigate the discrepancies between the results obtained by the 1-D and 2-D formulations,
the computed temperature profiles along the centerline are shown in Fig. 2.14.

Fig. 2.14. Computed temperature profiles along the stagnation streamline of geometry (a).

While the 1-D results predict a constant temperature from the exit of the top hot burner to the mixing
layer, the 2-D simulations reveal a linear decrease of temperature, which is due to radial heat loss due to
the parabolic temperature profile at the hot boundary (as shown in Fig. 2.15) which causes the heat to be
0.E+00
1.E-12
2.E-12
3.E-12
4.E-12
5.E-12
6.E-12
7.E-12
8.E-12
9.E-12
0.0 0.2 0.4 0.6 0.8 1.0
X
H
z /D
2-D, Lower Branch
2-D, Upper Branch
1-D, Upper Branch
1-D (adjusted T), Upper Branch
1-D (corrected), Upper Branch
1-D (L=0.4D), Upper Branch
900
905
910
915
920
925
930
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
T, K
z /D
2-D, Lower Branch
2-D, Upper Branch
1-D, Upper Branch
1-D (adjusted T), Upper Branch
1-D (corrected), Upper Branch
1-D (L=0.4D), Upper Branch

32
transferred from the hotter gases in the center of the jet to the outer periphery as they are being convected
axially.

Fig. 2.15. Comparison between uncorrected data and computed radial profiles of temperature of geometry
(a); profiles at z = 0.93D are obtained in the absence of the fuel-jet.

Since it is known that such losses can affect significantly the interpretation of counterflow results [9],
1-D simulations were repeated accounting for the temperature drop; the temperature at the axial position
at which the temperature drop deviates more than 1% from its linear trend (z = 0.43D), was chosen as a
new Tair for the “adjusted T” simulations. The results are shown in Figs. 2.13 and 2.14. The discrepancy
between the predicted maximum 𝑋 in the centerline was thus reduced from 232% to only 11%, indicating
that the most significant non-ideality of the experimental approach is the parabolicity of the temperature
profile at the hot boundary whose effect on the reported results has not quantified in past studies
[33,40,41]. Such parabolicity is required to ensure that ignition initiates at the centerline or else 1-D
assumptions are violated [39-41]. On the other hand, the use of the 1-D formulation is frequently justified
by asserting that although the measured velocity and temperature profiles have a parabolic shape in
general, they exhibit uniformity in a region around the centerline as small as few millimeters (e.g.,
[33,44]).
As already shown by previous studies [4-10], in the counterflow configuration, axial velocity profiles
must be of top-hat shape at the exit of the burners for the valid comparison of experimental data and results
of 1-D codes. The present results suggest also that that the burner exit temperature must be uniform as
well, and this is never the case in experiments.
The observed sensitivity to the temperature profile is expected considering the exponential dependence
of kinetics on temperature. Therefore, care should be taken when measuring temperatures at the hot
boundary and it was determined that the thermocouple readings are sensitive to the thermocouple wire
500
550
600
650
700
750
800
850
900
950
1000
0.0 0.1 0.2 0.3 0.4 0.5
T, K
r/D
2-D, inlet (z= 6D)
2-D, z= 0.93D
data, z= 0.93D
2-D, Upper Branch (z= D)
2-D, Upper Branch
corrected,
(z = z
kernel
)
z = z
kernel

33
size such that not utilizing thermocouples with a very small bead will introduce intrusive measurement
errors as big as 90 K as shown in Appendix C. It was found also that the intrusive effects can be minimized
by using micro-thermocouples and measuring in the single jet configuration with fuel jet turned off as
shown in Fig. 2.15.
The agreement between the uncorrected measured temperature and the 2-D numerical simulation
results is less compared to velocities specifically at larger radii. This is mostly due to the extra requirement
of the momentum balance at the burner exits for the numerical simulations, as discussed in Appendix B.
Moreover, the fact that the nozzle walls were assumed to be adiabatic is not realistic. However, at the
centerline, the temperatures after the correction do match and is equal to Tair = 926 K.
Given the parabolic temperature profiles at the burner exit [33,40,41] and the monotonic increase of
the curvature of the temperature profile at the centerline by increasing Tair [40], it is essential to introduce
corrections before comparing data to 1-D simulation results. Particularly, there is a need to estimate the
temperature at the ignition kernel in the experiment in the presence of radial heat loss, Tkernel, to
approximate the effective boundary temperature, (Tair)eff for the 1-D simulations similar to what was
performed in “adjusted T” simulations. Tkernel can be approximated by solving the 1-D heat conduction
problem in the radial direction. For given BC’s, the ignition kernel position as marked by the maximum
𝑋 , zkernel, can be determined through 1-D opposed jet simulations. By integrating the axial velocity profile
at the centerline obtained from 1-D simulations from z = D to z = zkernel, the time for the heat loss process,
tloss, can be computed. Subsequently, Tkernel can be obtained by letting the known temperature profile at
the hot boundary to evolve in time from its initial state to tloss. At last, since the boundary temperature is
always slightly higher than Tkernel, effective boundary temperature can be obtained such that
(Tair)eff = Tkernel × (
( )
) where (𝑇 )
is the computed kernel temperature from 1-D
opposed-jet simulation.
The success of this method was verified by applying on the upper branch solution. Evolution of the
radial temperature profile from the 2-D numerical results at the exit of the air jet (z = D) was tracked
utilizing the 1-D heat conduction simulations performed using laminarSMOKE code in a 1-D radial
geometry and “corrected” results are compared to the 2-D ones as shown in Fig. 2.15. The temperature
gradient was specified to be zero at r = 0 while at r = 0.5D, temperature was chosen to be equal to the
corresponding value from 2-D numerical results at the exit of the air jet. Results show only a difference
of 2 K at the centerline between the “corrected” and the 2-D results at z = zkernel = 0.39D. The discrepancy
increases at larger radii due to radial convection of gases with higher temperature from the centerline,
which causes the 2-D temperature profiles to reach higher values at larger radii. The 1-D opposed jet

34
simulations were repeated using (Tair)eff obtained from the “corrected” solution; results are also plotted in
Figs. 2.13. The “corrected” solution under-predicts the maximum 𝑋 only by 10%, therefore the
correction provides reasonable predictions considering the ~7% discretization uncertainty of the 2-D
numerical simulations.
Therefore, without having to perform expensive 2-D numerical simulations or extensive thermocouple
measurements that is subject to intrusive errors, 1-D conduction simulations could provide corrected air
boundary temperature to be used by 1-D opposed jet simulations whose results could be compared against
the data. It should be mentioned also that since the 1-D simulations denoted as “L = 0.4D” were performed
using Tair = 917 K, which is not a correct BC (as shown in Fig. 2.14), the predicted reactivity in the kernel
is still off by 58% (as shown in Fig. 2.13) although the velocity profiles are captured by the simulations
as demonstrated in Fig. 2.6. This is mostly due to the fact that the proper axial positions to apply the BC’s
for 1-D simulations is not known in priori because of the non-uniform temperature profile at the hot
boundary.
2.4.4. Simulation results for geometry (b)
2-D numerical simulations were performed also for geometry (b) and it was realized that bifurcation
also occurs despite having BC’s that conform to the assumptions of the 1-D formulation, which is mainly
due to the interaction of the jets with the adjacent disk walls. It should be noted that since the 1-D model
is derived for infinite-diameter burners, the introduction of finite diameters and adjacent walls could result
in recirculation bubbles and eventually bifurcation. Figure 2.16 depicts the 𝑋 distribution overlaid by
the streamlines for the two stable branches. Similar to Fig. 2.11, the recirculation bubbles at the junction
of coflow exits and disk walls have different sizes for the two branches indicating that the way the flow
orients itself is different for the two branches as it convects outward. It is also apparent from Fig. 2.16
that there is considerable amount of H2 both in upper and lower recirculation zones for both solutions
branches. This is mainly due to the large sizes of these recirculation bubbles, which cause considerable
amount of mixing [37,56-58].

35

Fig. 2.16. Comparison between computed 𝑿 𝐇 𝟐 maps overlaid by the streamlines for two stable solutions
of geometry (b); Left: lower branch; Right: upper branch.

Figure 2.17 depicts the 𝑋 variation along the centerline for both branches and for the 1-D
formulation. Results show that the bifurcation does not affect 𝑋 values and that the discrepancy between
the 1-D and 2-D results is about 12%. No noticeable difference was found between the axial velocity
profiles at the centerline (shown in Appendix D) for the two branches and the 1-D formulation reproduces
closely the flow velocities. This is expected as for geometry (b), axial velocities exhibit zero axial
gradients at the burner exits complying with the “plug flow” assumptions of the 1-D formulation (e.g.,
[1,3,10,20]).

36

Fig. 2.17. Computed 𝑿 𝑯 profiles along the centerline of geometry (b) and (c); inset: same profiles
plotted close to the stagnation point to aid the comparison.

To further investigate the 2-D effects that caused the discrepancy in reactivity shown in Fig. 2.13 and
2.17, the assumption of the constant radial pressure curvature eigenvalue J = (1/r)( P/ r) invoked in the
1-D formulation [1] was also tested. The typical results for obtained using 2-D simulations and those
obtained using the 1-D formulation are shown in Fig. 2.18.

Fig. 2.18. Comparison between the typical computed axial profiles of radial pressure curvature along the
stagnation streamline.
0.0E+00
5.0E-13
1.0E-12
1.5E-12
2.0E-12
2.5E-12
3.0E-12
3.5E-12
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
X
H
z /D
2-D, geometry (b), Lower Branch
2-D, geometry (b), Upper Branch
1-D
2-D, geometry (c)
2.9E-12
3.1E-12
3.3E-12
3.5E-12
0.58 0.60 0.62 0.64
-6
-5
-4
-3
-2
-1
0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
z /D
2-D, Upper Brach, geometry (a)
1-D, Upper Brance, geometry (a)
2-D, Lower Branch, geometry (b)
2-D, geometry (c)
1-D
𝟏 /𝒓 ( 𝝏 𝑷 /𝝏 𝒓 )×

37

J varies along the centerline for both geometries but more notably for geometry (a) for which 2-D
effects are more profound. Nevertheless, the average over the whole axial domain for upper branch of
geometry (a) is only 10% higher than the value predicted by the 1-D formulation the variations close to
the stagnation point are smaller in accordance with the findings of previous studies [5,16]. Similar order
of magnitude differences were observed for the upper branch of geometry (a). The discrepancy is much
smaller for geometry (b) such that the maximum deviation of J from the 1-D prediction is 5%. Johnson
et al. [11] showed that this maximum deviation for non-premixed H 2 flames is 3% with BC’s similar to
those used in geometry (b).
It was also found that the radial variations of the axial velocity and temperature profiles at an axial
position close to the stagnation point had considerably different shapes for the two branches. As a result,
the stagnation plane and kernel surface shapes are also different for the two branches. These results are
discussed in Appendix D.
While the 1-D and 2-D results of the velocities and scalars along the centerline are in close agreement,
the 𝑋 values predicted by 2-D simulations and shown in Fig. 2.19 indicate that the maximum reactivity
is realized well outside the region of influence of the main jets for both branches. Additionally, the
maximum 𝑋 value for the upper branch is 13% higher compared to the lower branch.

Fig. 2.19. Comparison between computed 𝑿 𝑯 maps for two stable solutions of geometry (b); Left: lower
branch; Right: upper branch.
Such 2-D effects are due to the presence of the recirculation zones next to jet exits. Although the size
of coflow annulus chosen for geometry (b) is larger than what is usually employed in experiments, it
cannot isolate the main jets from interacting with the reactants fed back to the high temperature region by
the recirculation bubbles as shown in Fig. 2.16. Specifically, since the bubble next to the upper burner
exit is bigger for the lower branch, high temperature gases lose more heat as they recirculate through the
bubble. Hence, the relatively hotter upper recirculation zone for the upper branch results in higher
reactivity in that region as confirmed by the 𝑋 maps illustrated in Fig. 2.20, which reveals a profound
difference in the distribution of HO2 (an important pre-ignition radical) between the two branches. While

38
the maximum 𝑋 occurs close to the wall at the upper recirculation zone for the upper branch, for the
lower branch it is located between the two recirculation zones with the peak value being 10 times smaller
than that of the upper branch. Hence, if geometry (b) is chosen in experiments, ignition will initiate at
large radii. In this case, the comparison of the experimental results with prediction using the 1-D
formulation may be meaningless.

Fig. 2.20. Comparison between computed 𝑿 𝐇𝐎 𝟐 maps for two stable solutions of geometry (b); Left: lower
branch; Right: upper branch.

The reason why such radical pool build up at large radii was not observed in the recirculation zones
of geometry (a) is that the temperature drops rapidly at larger radii in that geometry due to the parabolic
temperature profile of air jet in that geometry. The presence of the separation bubbles no matter how far
they are pushed away from the main jets by the coflow, affects the way the flow field develops and
modifies the species spatial distribution due to interactions of the main jest with the recirculation zones.
In fact, the distribution of minor species could be modified in geometry (b) if ambient temperature coflow
streams are utilized to eliminate the possibility of ignition initiation at larger radii, but on the other hand
radial temperature gradients will be introduced and which are inconsistent with the assumptions of
geometry (b).
2.4.5. Simulation results for geometry (c)
Only one stable state was found for geometry (c) as shown by temperature maps overlaid by
entrainment streamlines in Fig. 2.21. Moreover, only one large recirculation zone forms between the top
burner wall and upper outflow boundary. Due to larger size of this recirculation zone compared to the
ones observed in geometry (b), the temperature inside the bubble is approximately the ambient
temperature. High levels of buoyancy cause the main jet flow to rise upwards towards the upper open
boundary and suppress the bifurcation tendency [37].

39

Fig. 2.21. Computed temperature map superimposed by the entrainment streamlines for the only stable
solutions of geometry (c).

The variation of 𝑋 along the centerline of geometry (c) was also compared to the corresponding
profiles for geometry (b) and 1-D results in Fig. 2.17; 1-D results match the ones from geometry (c) within
the discretization uncertainty (7%). Same order of agreement was observed for the axial velocities along
the centerline as shown in Appendix D.
𝑋 isocontours overlaid by the 𝑋 maps for geometry (c) are shown in Fig. 2.22 revealing negligible
2-D effects. In contrast to geometry (b), 𝑋 maximizes at the centerline. The ignition kernel surface was
determined to be nearly planar with a slight curvature towards the upper burner. All other intermediate
radicals also peak in the region of jets interaction and iso-contours of the maxima are also flat along the
radius. Due to the absence of recirculation zones, the main reactivity region is isolated between the jets.
Moreover, the maximum variations of J along the centerline is only 2% for geometry (c) as shown in Fig.
2.18, which is in agreement with findings of Johnson et al. [11] for unconfined non-premixed flames.
Therefore, the results of 1-D codes match closely those for geometry (c).

40

Fig. 2.22. Computed 𝑿 𝑯 maps overlaid by 𝑿 𝐇 𝟐 iso-contours for geometry (c).

Therefore the design of a counterflow flame ignition experiment should be such that it matches the
conditions of geometry (c) as close as possible which is the burner design proposed by Niemann et al. [10]
with no solid surfaces close to jets exists. Solid surfaces close to the jet exits must be removed or trimmed
wherever possible to avoid recirculation [22], which can result in bi-stable states. It should also be noted
that the presence of vortices could introduce 3-D complications at higher Re. The axisymmetric vortex
rings can lose stability by side jet generation phenomena [59] and the flow can lose its axisymmetry by
azimuthal disturbances [60]. This is particularly important since the requirement of eliminating vortices
is not always fulfilled in the burner designs (e.g., [8,38]) and it has been reported to cause 2-D
complications in one case [38].
It is expected that for a given Tair and K and for momentum-balanced jets, there is a maximum
allowable thickness of disk walls i.e., cooling jacket thickness (denoted as δ in Fig. 2.2) above which
recirculation zones form close to the jet exists. These recirculation bubbles are susceptible to undesirable
instabilities. Therefore, it is worthwhile to perform a parametric study to investigate the effect of δ on the
bifurcation propensity. Moreover, since the existence of bi-stable states depend on both Tair and K it would
be useful if bifurcation maps could be generated to quantify the effects of each variable.

41
2.5. Concluding remarks
The one-dimensional assumptions invoked in counterflow flame ignition experiments were assessed
experimentally and computationally for non-premixed H 2 flames. The novelty of the present study is that
the combined effects of momentum and thermal boundary layers on the flow field established in
counterflow configurations have been considered.
For the first time, the existence of two stable states of the pre-ignition flow field was identified
experimentally for a given set of boundary conditions. Ignition was observed to occur at different fuel
mole fractions for the two stable states. Two-dimensional (axisymmetric) numerical simulations
reproduced the two stable solutions that are characterized by different distributions of velocities and
scalars and thus different ignition propensities. The solid surfaces next to the jets exits were identified to
cause the bifurcation due to flow separation.
It was determined also that one-dimensional simulations overestimate notably the reactivity around
the centerline compared to the results of the two-dimensional simulations. This difference is mainly
caused by the non-uniform radial temperature profiles that are unavoidable in ignition experiments. A
novel correction method based on simple numerical heat conduction calculations was developed to reduce
the discrepancies. It was shown that by utilizing an effective hot air boundary temperature as the boundary
condition for the one-dimensional approach, the two-dimensional results could be recovered.
To isolate two-dimensional effects, axisymmetric simulations were carried out also for ideal boundary
conditions that conform to those of the one-dimensional formulation with the only exception being the
use of finite burner diameter. While the two formulations result in similar velocities and scalar
distributions along the centerline as expected, it was found that radicals build up at larger radii due to the
jet-wall interaction. Therefore, ignition can initiate off-center. This undesirable effect can be eliminated
by removing solid surfaces from the jets vicinity.
It is recommended that the flow behavior observed only at pre-ignition states needs to be considered
in the design of the experimental apparatus, and that direct high-speed visualization of the entire flow-
field is required in order to evaluate potential errors when using standard one-dimensional codes to model
and interpret the data.

2.6. Appendix A: Momentum balance
One of the common practices in counterflow experiments is matching the momenta of the jets to adjust
the stagnation plane position to be roughly midway between the nozzle exits. Although one dimensional
codes are not limited to this requirement as long as the stagnation plane is not extremely close to one end
so that finite spatial gradients of concentration or temperature develop at the boundary, in experiments not

42
doing so, may result in radical losses at the burner rims or failure of propagation of the ignited flame due
to geometric constrains. More importantly, momentum imbalance may cause undesirable two-
dimensional non-uniformities.
However, due to presence of radial non-uniformities when the jets are heated, the momentum-
matching procedure is not as straightforward. First, velocity profiles at the exit have been reported to be
parabolic when heated even for nozzles with flow straighteners [33,40,41]. This is mostly due to high
thermal and viscous diffusion at high temperatures i.e., boundary layer thickness is larger at higher
temperatures. To further elaborate on the reason, it should be noted that by increasing 𝑇 , the density of
the air jet decreases; therefore, to maintain the momentum, the air velocity has to be increased. Moreover,
assuming that the velocity profile right upstream of the flow straightener is flat, the “entrance length”
relation in ducts can be used to estimate the length it takes for the flow to hydrodynamically fully develop
(e.g., [61]). Since the entrance length is directly proportional to 𝑅𝑒 and the momentum-balanced air
velocity is proportional to 𝑇 / and considering the dependency of the kinematic viscosity on temperature,
simple scaling analysis reveals that as 𝑇 increases, the entrance length decreases such that the length
required for the almost top hat velocity profile downstream of the flow straightener to develop into a
parabola decreases. In fact, this length was calculated to be six times smaller at 𝑇 = 1400 K than the
room temperature value. Moreover, flow straightener are not usually placed right at the exit of the nozzle
since immediately downstream of the screen the velocity profile is similar to small parallel jets strongly
modulated by the wakes of the wires rather than top hat [22] that is not desirable either. It is also known
that as 𝑇 increases, the flat portion at the center of the velocity profile shortens so the profile appears to
be more parabolic [40,41].
In this work, the momentum balancing procedure was performed first by measuring the temperature
at the center of the top nozzle and knowing that the cold jet is at room temperature i.e., 298 K, densities
of each stream were calculated accordingly by knowing the composition of the gases; using average exit
velocities i.e. assuming that velocities are equal to the ratio of their volumetric flow rates to the cross-
section areas of the nozzles, the momentum of the cold jet 𝜌 𝑉 was adjusted to match the value for the
hot stream 𝜌 𝑉 ; where 𝜌 , 𝑉 and 𝜌 , 𝑉 are average velocities and densities at the exit of the fuel and
air side respectively. Adjusting the velocity of the fuel stream was preferred over that of air since the
temperature of the air-jet will changes as a result of the heat release rate modification in the heated burner
and it was found that the whole system had to be allowed to reach to a new thermal equilibrium which
could take up to twenty minutes due to high heat capacity.
This method is common in the counterflow ignition experiments (e.g., [10,32]) but is approximate in
the sense that using the temperature at the center of the top nozzle exit to calculate the top jet momentum

43
is not accurate since the temperature profile is parabolic rather than uniform. More accurate method would
require to actually calculate the density profile by measuring the temperature profile at the exit and
integrate the momentum over the whole nozzle exit area taking into account both measured velocity and
calculated density profile shapes. Consequently, since 𝑉 is lower than the actual exit velocity at the center
due to the parabolic velocity profile, the momentum calculated using this velocity is underestimated
resulting in a smaller 𝑉 calculated based on momentum balance equation at the centerline. Furthermore,
since the velocity profile at the exit of the fuel side nozzle is “M shape” with a dip at the center due to
pressure feedback into the fuel side nozzle and the parabolic profile of the air-jet, the actual exit velocity
at the centerline is even smaller than what is required to fully balance the momentum at the center. The
consequence is that the stagnation point will be always closer to the fuel-jet exit.
Buoyancy driven velocities can affect the flow field through natural convection, which can in turn
change the position of the stagnation plane. Amantini et al. [8] showed that gravity has increasing effect
at larger radii and at low strain rates. It was also reported that the centerline location of the flame would
not be affected unless the effective strain rate [62] is lower than 60 s
-1
(e.g., [8,63]). Nevertheless, the
counterflow experiments are always performed at notably higher strain rates since stable planar flames
cannot be achieved otherwise.
Furthermore, N2 coflow is typically adjusted such that there is minimal degree of shear at the fluid
interface with the main flows (e.g., [33]) which usually translates to matching the average exit velocities
between the main and coflow (e.g., [64]) and since the heated coflow stream has higher temperature than
the cold one, there always exist momentum imbalance between the two coflow streams. This issue can
be resolved by eliminating the requirement for zero shear; this could however lead to shear instabilities
and it is usually avoided.
Moreover, due to the presence of mixing layer close to the stagnation plane, the stagnation point can
never be exactly positioned midway between the two burners, when the jets exhibit temperature and
composition differences, even when ignoring buoyancy and two-dimensional effects as it can be readily
demonstrated by one-dimensional codes.
Finally, even when dealing with unheated jets at room temperature, the stagnation plane will drift
away from the center after the bifurcation occurs at some critical Reynolds number, 𝑅𝑒 due to nozzle-
flow-walls proximity effect. Operational Reynolds numbers for counterflow experiments are always
higher than the criticality for isothermal jets meaning that the asymmetry always exists. For instance
Pawlowski et al. [23] showed that 𝑅𝑒 is a function of the aspect ratio, 𝛼 (ratio of separation distance to

44
burner diameter) and for 𝛼 < 2 is in the form of 𝑅𝑒 =
+ 30. Therefore, the symmetry breaks at
Reynolds numbers higher than 90 and 150 for aspect ratio of 1 and 0.5 respectively.
Based on the aforementioned considerations, momentum balance in experiments is rather a trial-and-
error practice to force the stagnation plane to settle midway between the nozzle exits (e.g., [10,65]). In
the present study, for consistency, the procedure outlined here was followed.

2.7. Appendix B: Details of modeling approach
Transient numerical simulations were performed for two-dimensional (2-D) axisymmetric problem.
The full set of conservation equations of mass, momentum, species, and energy along with reaction
chemistry as specified by a chemical mechanism were solved using the segregated approach by employing
the Strang splitting [66]. Pressure-velocity coupling was handled using PISO algorithm [67]. Since
buoyancy plays an important role in the counterflow flames particularly away from the centerline and at
low strain rates [8] it was specifically included in the momentum equation. Thermal radiation and Soret
diffusion were neglected, and the thermodynamic, mixture-averaged transport, and kinetic properties were
evaluated using the OpenSMOKE library [68]. Discretized equations were solved using finite volume
approach. Interpolation schemes used to evaluate variables at the control volume faces were linear
therefore schemes were second-order accurate in space. Implicit Euler time stepping was used to arrive
at the steady-state solution. While the Poisson equation for the pressure was solved using the geometric
algebraic multi-grid (GAMG) solver with Gauss-Seidel iteration, all of the remaining linear systems were
solved using the preconditioned biconjugate gradient (PBiCG) preconditioned by diagonal incomplete LU
(DILU). In order to resolve the small mass fractions of active radicals in pre-ignition state that can be
small as 10
(e.g., [69]) the tolerances of iterative sequences had to be decreased beyond what is usually
employed in flames [49]. For chemistry time integration which is handled by DLSODE Ordinary
differential equation solver in ODEPACK [70], absolute tolerance and relative tolerance were chosen to
be 10
and 10
respectively. For all other variables values of 10
and 0 were chosen
respectively. Such low level of tolerance increases the computational cost significantly such that the small
time steps as low as 10
had to be chosen for numerical stability. Such level of accuracy was found to
be necessary in order to obtain realizable values for active radical concentrations [71]. The code is highly
parallelized through the domain decomposition method on distributed-memory machines. The
communication between the processors occurs through the MPI protocol. Boundary conditions (BC) for
velocity, 𝑈 , temperature, 𝑇 , pressure, 𝑃 and species mole fractions, 𝑋 are summarized in Table 2.B.1.

45
Table 2.B.1. Boundary conditions employed for the 2-D numerical simulations.
Walls Inlets Outlet Slip
𝑈 zero Drichlet Drichlet zero Neumann for
normal component
zero Neumann for
tangential component,
zero Dirichlet for normal
component
𝑇 zero Neumann Drichlet zero Neumann zero Neumann
𝑃 zero Neumann zero Neumann Drichlet zero Neumann
𝑋 zero Neumann Drichlet zero Neumann zero Neumann

Since the current numerical formulation is solving the pressure Poisson equation (PPE) on collocated
grids, it requires the pressure to be specified at the outlet. This requirement necessitates the outlet to be
far away from the regions that have high gradient of variables. The radial position of the outlet was chosen
to be 30𝐿 from the symmetry axis. Atmospheric pressure was specified at the outlet. It was observed that
if the outlet is chosen to be at radial values less than 13𝐿 , the steady state solution could not be achieved
which was speculated to be the result of the presence of the recirculation zones up to that radial extend
due to high 𝑅𝑒 of the investigated conditions (𝑅𝑒 ≈ 1187 for the fuel stream). With the radial domain
extension of 30𝐿 the flow in between the discs fully develops thermally and hydrodynamically in the
absence of buoyancy as confirmed by velocity and temperature profiles at the outlet.
For wall-confined simulations, uniform structured mesh was used due to low performance of GAMG
method in solving the PPE with anisotropic grids [72]. Non-uniform structured mesh was chosen for
simulations without wall confinement due to large mesh size required to properly model the open
boundaries. The problem was first solved on a series of coarse grids and subsequently mapped on to a
grid with four times number of cells obtained by splitting each cell in two in each spatial direction. To
ensure mesh independent solution, the maximum 𝑋 of the minor species that has the highest concentration
in the flow-field was tracked through the mesh refinement steps. The solution was considered mesh
independent when the change of the mentioned variable is less than 2% upon refinement. The finest
spatial resolution is 50 m.
Numerical error quantification using Roaches’ method and Richardson extrapolation [73] was
performed to investigate the grid dependence of the solutions and specifying the error margin of the
results. The solutions of three consecutive mesh refinements were traced to estimate the order of accuracy
for maximum H radical mole fraction. Since the final order of accuracy is 1.54, the safety factor of 1.5

46
was chosen for grid convergence index calculations. Using these parameters, the discretization error of
maximum H radical mole fraction was determined to be 7% for all geometries considered.

2.7.1. Geometry (a) Specifications
This geometry is a representative of experimental setup as depicted in Fig. 2.4a. 𝛼 = 1 was chosen
to comply with experimental conditions. In order to capture the effect of pressure feedback into the tubes
(e.g., [20,55]) caused by high static pressure at the stagnation point in high strain rates, the domain is
extended by the amount of two diameters into the nozzles. Considering this extension is essential to
capture the behavior of the bifurcated state [25]. To ensure the adequacy of the chosen value for into the
tubes extension, the value of 4𝐷 was also tested in non-reactive flow and it was confirmed that the axial
velocity profiles at 𝑧 = 0 and 𝑧 = −2𝐷 are identical.
The experimental geometry had to be simplified for the numerical grid generation purposes. To be
specific, the exit of the main tubes and coflow streams were assumed to be at the same height. Burner
rims were assumed to be straight at the exit without any trims and both nozzle rims were assumed to have
a same size such that top and bottom nozzles are mirror images with respect to the center of the domain.
The cooling jacket for the lower burner and insulation blanket for the upper burner were extended all the
way to outlet.
The attempt to simulate the flow inside the tube introduces another complexity to the problem i.e., the
profile of the velocity inside the tube cannot be measured directly and used as an inlet condition to the
domain for numerical simulations. Bouvet et al. [15] used uniform profiles and extended the domain
inside the tube up to the position of the flow straightener. In this approach, since all the numerical cells
have the same value of velocity, numerical difficulties are raised at the cell right the next to the tube wall
at the inlet due to presence of large gradient. Moreover, right after the screen the velocity profile is similar
to parallel small jets strongly modulated by the wakes of the wires rather than top hat [22].
Alternatively, one can try to match the experimental exit profile of a single jet streaming into the
domain by trial-and-error on the inlet velocity profile inside the tube. First, radial profiles of the axial
velocity were measured at the exit of the burner in the absence of the opposing jet. Since the exit velocity
profile of the fuel jet is top-hat around the centerline with a finite boundary layer close to the burner rim,
by knowing the maximum velocity and the boundary layer thickness at the exit, the shape of the velocity
profile at the inlet of the domain can be estimated taking advantage of basic boundary layer theory (e.g.,
[61,74]). Then another set of simulations had to be performed using the estimated inlet velocity profile to
find the optimum. The functional form for the fuel stream inlet axial velocity was chosen according to Eq.
(2.1),

47

𝑈 (𝑟 ) =
⎩
⎪
⎨
⎪
⎧ 𝑈 , 0 ≤ 𝑟 ≤
𝐷 2
− 𝜃 𝑈 1 − 𝑟 −
𝐷 2
+ 𝜃 𝜃 ,
𝐷 2
− 𝜃 ≤ 𝑟 ≤
𝐷 2

(2.1)
where 𝑈 is the maximum velocity in the top-hat region and 𝜃 is a length scale representing the
boundary layer thickness. This profile represents the top-hat section by the uniform value and a finite
effect of the boundary layer by a parabolic drop of velocity. 𝑈 = 1.6 𝑚 /𝑠 and 𝜃 = 𝐷 4 ⁄ were found
to best match the measurements as shown in Fig. 2.B.1. Flatness of the profile close to the center helps to
conform to the 1-D assumptions.

Fig. 2.B.1. Comparison between measured and computed of radial profiles of axial velocity for the fuel-
jet streaming at the experimental conditions into the domain of geometry (a) in the absence of air-jet.
The same procedure was followed for the air-jet by considering both velocity and temperature profiles
at the exit of air nozzle. Since the flow is expected to be fully developed inside the nozzle, functional
form of the inlet velocity and temperature were assumed to be parabolic in the form of Eq. (2.2) and Eq.
(2.3) respectively.

𝑈 (𝑟 ) = 𝑈 (1 − 2𝑟 𝐷 )

(2.2)
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
0.0 0.1 0.2 0.3 0.4 0.5
U, m/s
r/D
inlet, z= -2D
z= 0.07D
PIV, z= 0.07D

48

𝑇 (𝑟 ) = (𝑇 − 𝑇 )(
2𝑟 𝐷 )
+ 𝑇
(2.3)
Optimization was performed on the values of the maximum velocity of the inlet profile, 𝑈 , the
maximum temperature (𝑇 ) and temperature at the wall at inlet (𝑇 ) in order to achieve profiles at the
exit of air nozzle that best match the experimental ones in the absence of the fuel-jet. This optimization
procedure for the air-jet had to be performed with a constraint on the system, namely the momentum
balance. Therefore, 𝑈 , 𝑇 , and 𝑇 were varied systematically such that the momentum of the
air-jet remains within 5% of the fuel-jet considering uniform temperature of 𝑇 = 298 K for the fuel-jet.
Therefore, a program had to be written to automate the process and simulations were subsequently
performed for the single jet configuration to obtain profiles at the exit of the air nozzle to be compared to
the measured profiles. 𝑈 = 3.05 𝑚 /𝑠 , 𝑇 = 675 K and 𝑇 = 960 K found to best replicate the
measurements as demonstrated in Figs. 2.B.2 and Fig. 2.15.

Fig. 2.B.2. Comparison between measured and computed radial profiles of axial velocity for the air-jet at
the experimental conditions streaming into the domain of geometry (a) in the absence of fuel-jet.
Simulation results capture the velocity exit profile of the air-jet with success over the most radial
positions expect very close to the nozzle rims where the predicted value is slightly higher than the PIV
data. This discrepancy is due to buoyancy-driven instabilities, which occur because of low flow velocities
in that region as confirmed by shadowgraph images.
It should be noted that with the above-mentioned optimized variables, the fuel-jet momentum is 4%
smaller than the air-jet in accordance with the experimental approach used to balance the jet momenta
-3.1
-2.6
-2.1
-1.6
-1.1
-0.6
-0.1
0.0 0.1 0.2 0.3 0.4 0.5
U, m/s
r/D
inlet, z= 6D
z= 0.93D
PIV, z= 0.93D

49
discussed in Appendix A. The inlet velocities of the coflow streams were assumed to have functional
form of the fully developed annulus flow [75] in the form of Eq. (2.4),

𝑈 (𝑟 ) = 𝑈 1 − 𝑟 𝑅
+
1 −
𝑅 𝑅 ln (
𝑅 𝑅 )
× ln (
𝑟 𝑅 )
(2.4)
where 𝑅 and 𝑅 are inner and outer radius of the annulus respectively as denoted in Fig. 2.4a. 𝑈 is the
velocity at the axis of the cylinder for a pipe Poiseuille flow with the same volumetric flow rate. For
𝑅 ≈ 𝑅 , this profile resembles a parabola with the slight shift of the maximum velocity position towards
the inner radius. In order to minimize shear between the main flow and coflow, the average velocities of
the two were set to be equal during the experiments and simulations. Dou et al. [76] showed that the
average velocity in the annulus, 𝑈 can be written in the form of Eq. (5).

𝑈 =
𝑈 2
1 + 𝑅 𝑅
−
1 −
𝑅 𝑅 ln (
𝑅 𝑅 )

(2.5)
Using 𝑈 = 64.03 𝑚 /𝑠 and 𝑈 = 97.49 𝑚 /𝑠 for the fuel- and air-jets respectively, the zero-shear
requirement was guaranteed. Since the coflows could not be seeded in the current experimental setup,
uniform temperature of 𝑇 = 675 K was assumed for the air-side coflow to make the momentum of the
coflow streams as close as possible while maintaining zero shear. The coflow surrounding the fuel-jet was
set to be at room temperature. With the mentioned conditions, the fuel coflow has a momentum, which is
2% smaller than the air coflow in accordance with main flows momenta.

2.7.2. Geometry (b) Specifications
This geometry is representing an ideal counterflow setup, which has the minimum non-uniformity
when compared to the assumptions of 1-D formulations. Inlet velocity profiles are top-hat in the form of
Eq. (2.1). The air- and fuel-jet have uniform temperatures of 920 K and 298 K respectively at the burners
exits; the inlet of the domain is placed at the burners exits so that the effects of the pressure-feedback are
eliminated. The attempt to place flow straighteners very close to the exit of the tube incorporated by
previous studies [22,32,77] was to replicate these conditions at the exit of the nozzles. Moreover, since
1-D codes only model the centerline, they implicitly assume that the jet has an infinite diameter that occurs
for very low aspect ratios and thus 𝛼 = 0.5 was selected. Supposedly, at such low level of separation
distance, the multi-dimensional effect has to be lower due to high levels of stretch. In order to avoid

50
feeding of unconsumed reactants back to the ignition kernel by the recirculation bubbles created by the
separation at the jet exit, N2 coflow had to be included next to the main inlets. To eliminate any possible
shear or wake between the coflow and main jets, the coflow boundary starts at the radial position where
the main inlet boundary ends and the velocities are equal as schematically demonstrated in Fig. 2.4c.
However, it was realized that this approach introduces numerical instabilities since there exist large
artificial density gradient at the boundary of main inlet and coflow due to difference in gas composition.
To resolve this issue, a finite length mixing layer was considered between the main jet and the coflow
such that the fuel concentration decreases almost linearly form its value at 𝑟 ≈ 𝐷 2 ⁄ where main fuel inlet
ends to zero at the end of the coflow boundary (𝑟 ≈ 1.06𝐷 ) as depicted in Fig. 2.B.3. It should be noted
that the size of the inlets was chosen to satisfy the momentum balance requirement as discussed below.
Figure 2.B.3. also depicts the procedure that was performed for the air-side. It should be noted that as the
H2 and O2 mole fractions decrease radially in the mixing layer, the value for N 2 increases correspondingly
as the main component of the coflow.

Fig. 2.B.3. Radial profiles of inlet mixture mole fractions for geometry (b) demonstrating the optimized
finite mixing layer thickness chosen to represent the coflow streams.

The N2 coflow annulus size (0.56𝐷 ) is wide enough to reduce the effects of reactant residence time
increase in the high temperature region between the jets due to back flow. At the same time, is not wide
enough for the annulus flow to become hydrodynamically unstable [76]. 𝑈 = 1.8 𝑚 /𝑠 was chosen
for the fuel-side in order to simulate conditions close to experimental conditions and those of geometry
(a) in terms of 𝑅𝑒 . By fixing this value the corresponding value for the air-side will be determined by the
momentum balance requirement at the centerline since temperatures and compositions are known. The
values of the boundary layer thicknesses had to also be determined by momentum balance requirement.
The procedure herein is not as straightforward as geometry (a) due to complex concentration profiles.
0.00
0.05
0.10
0.15
0.20
0.25
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
X
r /D
O2, z=1.06D
H2, z=0
0.75
0.80
0.85
0.90
0.95
1.00
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
X
r /D
N2, z=1.06D
N2, z=0

51
Another program was written to search for the optimized values in their allowable range and assure
momentum balance. Table 2.B.2 shows the summary of the results.

Table 2.B.2. Optimized values of maximum velocity and boundary layer thickness used in the inlet
velocity profiles of geometry (b) in the form of Eq. (1).
Fuel Side Air Side
𝑈 , 𝑚 /𝑠 1.8 2.81
𝜃 0.33𝐷 0.13𝐷

2.7.3. Geometry (c) Specifications
Geometry (c) is similar to geometry (b) with the only difference of replacing walls by open boundaries.
The velocity and scalar BC’s of jets inlets of geometry (b) are exactly adopted for jets inlets of geometry
(c). In order to implement the open boundaries properly, the BC’s utilized in the work by Oh et al. [63]
were adopted. Since the entire flow is convecting upwards towards the upper open boundary (Side C) due
to buoyancy, Side C is specified as outlet, the right boundary (Side B) is specified as slip, and the lower
boundary (Side A) is of inlet type [63] all of which are schematically shown in Fig. 2.4b. Since the inlet
entrainment N2 velocity, Uentrainment, of the lower open boundary is not known in priory in contrast to the
study performed by Oh et al. [63], an approximation had to be performed for finding the proper value for
this BC. Simulations were first performed with turning off the chemistry solver and setting the
acceleration of gravity to zero in the momentum equation along with specifying all of the open boundaries
as outlets. In doing so, it removes the buoyancy effects and lets the entire flow convecting radially
outwards from the right boundary and generating entrainment flows that enter the domain from the top
and bottom open boundaries (Side A and C). After achieving a steady state solution, the axial velocity at
the Side A boundary was extracted. The inlet entrainment velocity profile for the buoyant simulations
was specified as a uniform axial velocity whose value is twice of the average axial velocity at Side A for
non-buoyant simulations. This is justified since in the buoyant case the entrainment is generated only by
Side A. It should be noted that since Uentrainment is an order of magnitude smaller than the main jets
velocities, it is not expected to change the global stability of the flow field. This fact is further verified
by performing tests with increasing and decreasing Uentrainment by 50% to check the global stability
variations.
In order to damp the entropy and vorticity reflections from the right slip and top outflow boundaries
(e.g., [78,79]) and avoiding the natural convection instabilities (e.g., [80]) occurring in the buoyant rise of
the flow as it covects out of the Side C, four measures had to be undertaken. First, grading in the mesh

52
had to be considered outside the domain of jets interaction such that cells close to the open boundaries are
roughly 10 times bigger than the cells in the uniform region inside the jets interaction domain. Second,
tests had to be performed to find the minimum value of the radial domain size, W such that upon increasing
W by 10%, streamlines of the main jets are not affected. Third, another set of tests was required to find
the minimum distance of the surrounding inlet from the fuel jet inlet, Hbottom, such that the flow along the
bottom tube wall become fully developed at 𝑧 = 0 and similarly upon increasing Hbottom by 10%
streamlines of the main jets are not affected. Finally, a sponge-layer [79] had to be considered in the mesh
for 𝑧 ≥ 𝑧 by utilizing upwind-differencing scheme [81] for cell-face interpolations of the non-
linear term in the momentum equation. The last set of tests had to performed by varying the distance of
the top outlet boundary from the air jet inlet, Htop, and 𝑧 such that buoyant instabilities of the main
jets streamlines are avoided while their orientation remain intact in order to find the maximum value for
𝑧 and minimum for Htop. All the aforementioned tests were performed by turning off the chemistry
solver. Table 2.B.3 summarizes the optimized values.

Table 2.B.3. Optimized geometrical parameters and BC values used for geometry (c).
Uentrainment, m/s 𝑧 W Hbottom Htop
0.02 4.06D 12.5D 2D 7D

2.8. Appendix C: Effect of intrusive temperature measurements
Temperature is a key observable in flame ignition experiments and is measured using thermocouples
(e.g., [31,32]). The effect of the thermocouple probe was assessed using shadowgraph. Temperature
measurements are performed typically by inserting the thermocouple bead into the test section resulting
from two impinging jets (e.g., [31,32]). Since utilizing thermocouple is an intrusive procedure, and
considering the reports of inaccuracies caused by intrusive measurements in flames (e.g., [50]), to make
sure that the presence of the thermocouple wire is not affecting the flow and the measurement itself,
measurements were performed in a single jet configuration by switching off the fuel-jet. This approach
is justified since the correlations for convection correction are derived only for the uniform cross flow
about a cylinder or sphere [82], which is not the case when two impinging jets are involved due to velocity
gradients. Additionally, the thermocouple wire is assumed not to disturb the flow if the diameter of the
wire is small [64], and the wake influence of the wire is considered negligible for Re < 40 based on wire
diameter [83]. Although this requirement is satisfied for all experimental conditions of the present study,

53
it was observed that the wake is noticeable if measurements were performed in the presence of two
impinging jets, as shown in Fig. 2.C.1.

(a)

(b)
Fig. 2.C.1. Shadowgraph images of counterflow field showing the effect of intrusive temperature
measurement by thermocouple; the thermocouple rake was translated gradually from the right side into
the field such that the shadow of ceramic insulator is in the view; (a) micro-thermocouple; (b) ordinary
thermocouple where the rake is purposely tilted to make the thermocouple bead and wire visible.

When the thermocouple is inserted close to the air nozzle exit, the shape and position of the stagnation
plane are modified (Fig. 2.C.1a). It was observed that the flow decelerates downstream of the
thermocouple wire causing a cusp in the stagnation plane despite the fact that Re = 8. A constant tension
thermocouple rake [40] was utilized to avoid possible drooping of the thermocouple wire in hot
environment. The schematic of the thermocouple assembly is shown in Fig. 2.C.2. The bead was translated
in increments of 0.5 mm to resolve the temperature field. The cylindrical geometry was used instead of a
spherical one since for fine wire thermocouples the heat transfer from the thermocouple wires dominates
the junction energy balance [82]. Constant emissivity of 0.2 was used for the thermocouple bead [32].
Since the temperatures measured in this study are relatively low compared to the ones encountered in the
flames, the amount of radiation correction is not sensitive to the value chosen for thermocouple bead
emissivity. It was determined that 25% change in the emissivity will only change the corrected value by
2 K.

Fig. 2.C.2 Schematic of the thermocouple.

It should be noted that the length of the wire between the two ceramic insulators was chosen to be
twice that of the burner diameter to minimize the effect of ceramic insulators on the flow field. The effect
of the thermocouple wire diameter was also investigated by using a thermocouple with a wire five times
Bead
Ceramic insulator

54
thicker (250 m) than the original one. As shown in Fig. 2.C.1b, the flow first attaches to the wire and
creates a vortex street and vortex roll over in the shear layer causing an incomplete heat transfer and
incorrect temperature readings as shown in Fig. 2.C.3, in which uncorrected temperature readings for both
thermocouples in the presence of both jets are shown.

Fig. 2.C.3. Temperature readings at the exit of air jet (z = 0.93D) using two kinds of thermocouples
demonstrating the effect of intrusive measurements in the counterflow jet domain.

Error bars were computed based on the temperature controller resolution and fluctuations in readings.
More fluctuations observed at larger radii for micro-thermocouple while these values for the ordinary
thermocouple readings were within 2 K. The ordinary thermocouple readings were considerably lower
than the micro-thermocouple such that at the centerline it showed 91 K lower temperature. When the flow
attaches to the wire of ordinary thermocouple, as the bead gets closer to the centerline, the temperature
profile also becomes flat. This can result in the misleading conclusion that at least around the centerline,
the assumptions of 1-D formulation are met (e.g., [33,44]).
The attempt to use a thermocouple with an even smaller wire (25 m) failed since the wire is extremely
fragile at these scales and it readily breaks under any tension by the spring assembly of the rake. To the
author’s knowledge, all the counterflow ignition studies have been performed using thermocouples with
wire diameter larger than 25 m [31,32,40,33,41]. One might argue that since the correction to the
temperature based on convection/radiation balance of the bead is directly proportional to the bead diameter
[33], the readings might match after the correction. However, the correction is estimated to be less than
10 K if one assumes that the Nusselt number correlation employed in [33] still holds even for the unsteady
500
550
600
650
700
750
800
850
900
950
1000
0.0 0.1 0.2 0.3 0.4 0.5
T, K
r/D
micro-thermocouple
ordinary thermocouple

55
flow attached to the wire. This correction is an order of magnitude lower than the differences seen in Fig.
2.C.3. It should be noted also that the possibility of compromised temperature readings due to “aging” as
a result of the thermocouple overuse (e.g., [33,41]) was removed as the data presented were taken with
brand new thermocouples.
On the other hand, by measuring the temperature field in the single jet configuration the
aforementioned effects were eliminated, results of which were shown in Fig. 2.B.3. Thus, the typically
used argument regarding thermocouple flow disturbances that is based on uniform flow assumption over
a cylinder (e.g., [83]), does not apply.

2.9. Appendix D: Further results and discussion
Figure 2.D.1 depicts the axial velocity profiles at the centerline computed for both branches of
geometry (b) and for the 1-D formulation; the 1-D results are overlapping the upper branch solution from
the 2-D simulations. The differences between the lower branch solution and the one for 1-D are minor
compared to those for geometry (a) shown in Fig. 2.6. Figure 2.D.1 also includes the results of 2-D
simulation for geometry (c). No notable difference is observed between the upper branch axial velocity
profile of geometry (b) and the one for geometry (c).

56

Fig. 2.D.1. Computed axial velocity profiles along the stagnation streamline for geometry (b) and (c);
inset: same profiles plotted close to the stagnation point to aid the comparison

To investigate 2-D effects in geometry (b), the radial variations of the axial velocity and temperature
profiles were computed for both branches at an axial position close to the stagnation point shown in Figs.
2.D.2 and 2.D.3.
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
U, m/s
z /D
2-D, geometry (b), Upper Branch
2-D, geometry (b), Lower Branch
2-D, geometry (c)
1-D
-0.25
-0.05
0.15
0.50 0.55

57

Fig. 2.D.2. Comparison between results of 2-D numerical simulations for radial profiles of axial velocity
close to the stagnation plane (z = 0.57D) for two stable solutions of geometry (b).

While the axial velocity profile for the lower branch exhibits a notable dip at the centerline, the upper
branch profile is almost flat up to r/D = 0.8. The results of Fig. 2.D.3 reveal a non-uniform radial
temperature profile for the lower branch, caused by the radial convective transport of thermal energy. On
the other hand, the radial temperature profile for the upper branch is relatively flat mirroring the behavior
of the axial velocity profile. The differences between the two branches can affect the ignition behavior.

-0.9
-0.8
-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
U, m/s
r /D
Lower Branch
Upper Branch

58

Fig. 2.D.3. Comparison between results of 2-D numerical simulations for radial profiles temperature close
to the stagnation plane (z = 0.57D) for two stable solutions of geometry (b).

Figure 2.D.4 depicts the computed loci of the stagnation plane as well as the ignition kernel surface
identified by the maximum 𝑋 iso-surface. The stagnation plane for the lower branch is curved down
towards the lower burner, which is expected due to the dip observed for the axial velocity profiles in Fig.
2.D.2. On the other hand, the stagnation plane for the upper branch is almost flat with a slight downward
inclination. Kernels for both branches formed on the air-side due to presence of high temperature mixture
typical to non-premixed flame ignition studies (e.g., [33,41]). Both kernel surfaces are almost flat up to
a radius at which the coflow inert gas is injected.
850
855
860
865
870
875
880
885
890
895
900
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1
T,K
r /D
Lower Branch
Upper Branch

59

Fig. 2.D.4. Magnified view of the computed loci of the stagnation planes and kernel surfaces for two
stable solutions of geometry (b); halfway distance between the nozzle exits is marked for comparison
purposes; (--) kernel surface (─) stagnation plane; Red: upper branch, Blue: lower branch.

2.10. Acknowledgements
This material is based upon work supported as part of the CEFRC, an Energy Frontier Research Center
funded by the U.S. Department of Energy, Office of Science, and Office of Basic Energy Sciences under
Award Number DE-SC0001198. The numerical simulations were carried out using the Extreme Science
and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation
grant number TG-CTS140012. The authors would like to gratefully acknowledge the assistance of Dr.
Vyaas Gururajan for his feedback and technical support with the numerical simulations, and of Dr.
Francesco Carbone for the velocity measurements.

0.45
0.50
0.55
0.60
0.65
-1.1 -0.9 -0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 0.9 1.1
z /D
r/D
Midway position

60
2.11. References
[1] R.J. Kee, J.A. Miller, G.H. Evans, G. Dixon-Lewis, A computational model of the structure and
extinction of trained, opposed flow, premixed methane-air flames, Symp. (Int.) Combust. 22 (1988)
1479–1494.
[2] A. Ern, M.D. Smooke, Vorticity-velocity formulation for three-dimensional steady compressible
flow, J. Comput. Phys. 105 (1993) 58-71.
[3] H.K. Chelliah, C.K. Law, T. Ueda, M.D. Smooke, F.A. Williams, An experimental and theoretical
investigation of the dilution, pressure and flow-field effects on the extinction condition of methane-
air-nitrogen diffusion flames, Symp. (Int.) Combust. 23 (1990) 503-511.
[4] Y.M. Kim, H-J. Kim, Multidimensional effects on structure and extinction process of counterflow
nonpremixed hydrogen-air flames, Combust. Sci. Tech. 137(1998) 51-80.
[5] C.E. Frouzakis, J. Lee, A.G. Tomboulides, K. Boulouchos, Two-dimensional direct numerical
simulation of opposed-jet hydrogen-air diffusion flame, Symp. (Int.) Combust. 27 (1998) 571–577.
[6] W.-C. Park, A. Hamins, Investigation of velocity boundary conditions in counterflow flames, J.
Mech. Sci. Technol. 16 (2) (2002) 262 –269.
[7] E. Korusoy, J.H. Whitelaw, Inviscid, laminar and turbulent opposed flows, Int. J. Numer. Meth. Fl.
46 (2004) 1069-1098.
[8] G. Amantini, J.H. Frank, M.D. Smooke, A. Gomez, Computational and experimental study of steady
axisymmetric non-premixed methane counterflow flames, Combust. Theory and Model. 11 (2007)
47-72.
[9] V. Mittal, H. Pitsch, F.N. Egolfopoulos, Assessment of counterflow to measure laminar burning
velocities using direct numerical simulation, Combust. Theory and Model. 16 (2012) 419-433.
[10] U. Niemann, K. Seshadri, F.A. Williams, Accuracies of laminar counterflow flame experiments,
Combust. Flame 162 (2015) 1540-1549.
[11] R.F. Johnson, A.C. VanDine, G.L. Esposito, H.K. Chelliah, On the Axisymmetric Counterflow
Flame Simulations: Is There an Optimal Nozzle Diameter and Separation Distance to Apply Quasi
One-Dimensional Theory?, Combust. Sci. Tech. 187 (2015) 37-59.
[12] J.M. Bergthorson, D.G. Goodwin, P.E. Dimotakis, Particle streak velocimetry and CH laser-induced
fluorescence diagnostics in strained, premixed, methane–air flames, Proc. Combust. Inst. 30 (2005)
1637-1644.
[13] P. Versailles, J.M. Bergthorson, Optimized Laminar Axisymmetrical Nozzle Design Using a
Numerically Validated Thwaites Method, J. Fluids Eng. 134 (2012).

61
[14] P. Versailles, G.M. Watson, A.C. Lipardi, J.M. Bergthorson, Quantitative CH measurements in
atmospheric-pressure, premixed flames of C1–C4 alkanes, Combust. Flame 165 (2016) 109-124.
[15] N. Bouvet, D. Davidenko, C. Chauveau, L. Pillier, Y. Yoon, On the simulation of laminar strained
flames in stagnation flows: 1D and 2D approaches versus experiments, Combust. Flame 161 (2014)
438-452.
[16] J.M. Bergthorson, K. Sone, T.W. Mattner, P.E. Dimotakis, D.G. Goodwin, D.I. Meiron, Impinging
laminar jets at moderate Reynolds numbers and separation distances. Phys. Rev. E 72 (2005)
066307.
[17] J.M. Bergthorson, S.D. Salusbury, P.E. Dimotakis, Experiments and modelling of premixed laminar
stagnation flame hydrodynamics, J. Fluid Mech. 681 (2011) 340-369.
[18] L.W. Kostiuk, K.N. Bray, R.K. Cheng, Experimental study of premixed turbulent combustion in
opposed streams. Part I—Nonreacting flow field, Combust. Flame 92 (1993) 377-395.
[19] G. Pellett, K. Isaac, W. Humphreys, L. Garttrell, W. Roberts, C. Dancey, G. Northam, Velocity and
thermal structure, and strain-induced extinction of 14 to 100% hydrogen-air counterflow diffusion
flames, Combust. Flame 112 (1998) 575–592.
[20] B.G. Sarnacki, G. Esposito, R.H. Krauss, H.K. Chelliah, Extinction limits and associated
uncertainties of nonpremixed counterflow flames of methane, ethylene, propylene and n-butane in
air, Combust. Flame 159 (2012) 1026-1043.
[21] J.L. Fox, G.W. Morgan, On the stability of some flows of an ideal fluid with free surfaces, Brown
Univ. Providence RI DIV of Applied Mathematics (1952), paper TR-2.
[22] J.C. Rolon, D. Veynante, J.P. Martin, Counter jet stagnation flows, Exp. Fluids 11 (1991) 313-324.
[23] R.P. Pawlowski, A.G. Salinger, J.N. Shadid, T.J. Mountziaris, Bifurcation and stability analysis of
laminar isothermal counterflowing jets, J. Fluid Mech. 551 (2006) 117-139.
[24] R.P. Pawlowski, Numerical studies of complex reacting Flows, Ph.D. thesis, 2000.
[25] A. Ciani, W. Kreutner, C.E. Frouzakis, K. Lust, G. Coppola, K. Boulouchos, An experimental and
numerical study of the structure and stability of laminar opposed-jet flows, Comput. Fluids 39 (2010)
114-124.
[26] S.T. Strogatz, Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry, and
engineering, Westview Press, 2000, pp. 44-93.
[27] R.M. Feran, T. Mullin, K.A. Cliffe, Nonlinear flow phenomena in a symmetric sudden expansion,
J. Fluid Mech. 221 (1990) 596-608.
[28] D.A. Hammond, L.G. Redekopp, Local and global instability properties of separation bubbles, Eur.
J. Mech.,B/Fluids 17 (1998) 145-164.

62
[29] P.G. Drazin, W.H. Reid, Hydrodynamic stability, Cambridge University Press, 2004, pp. 126-144.
[30] O. Park, P.S. Veloo, H. Burbano, F.N. Egolfopoulos, Studies of premixed and non-premixed
hydrogen flames, Combust. Flame 162 (2015) 1078-1094.
[31] N. Liu, C. Ji, F.N. Egolfopoulos, Ignition of non-premixed C3–C12 n-alkane flames, Combust. Flame
159 (2012) 465-475.
[32] R. Seiser, K. Seshadri, E. Piskernik, A. Linán, Ignition in the viscous layer between counterflowing
streams: Asymptotic theory with comparison to experiments, Combust. Flame 122 (2000) 339-349.
[33] C.G. Fotache, Ignition in convestive-diffusive systems, Ph.D. thesis, 1997.
[34] C.G. Fotache, T.G. Kreutz, D.L. Zhu, C.K. Law, An Experimental Study of Ignition in Nonpremixed
Counterflowing Hydrogen versus Heated Air, Combust. Sci. Technol. 109 (1995) 373–393.
[35] V. Gupta, S.A. Safvi, T.J. Mountziaris, Gas-phase decomposition kinetics in a wall-less environment
using a counterflow jet reactor: design and feasibility studies, Ind. Eng. Chem. Res. 35 (1996) 3248-
3255.
[36] W. Zhang, Z. Chai, Z. Guo, B. Shi, Lattice Boltzmann Study of Flow and Temperature Structures
of Non-Isothermal Laminar Impinging Streams, Commun. Comput. Phys. 13 (2013) 835-850.
[37] N. Hassan, S.A. Khan, Two-dimensional interactions of non-isothermal counter-flowing streams in
an adiabatic channel with aiding and opposing buoyancy, Int. J. Heat Mass Transfer 54 (2011) 1150-
1167.
[38] A. Ciani, W. Kreutner, W. Hubschmid, C.E. Frouzakis, K. Boulouchos, Experimental investigation
of the morphology and stability of diffusion and edge flames in an opposed jet burner, Combust.
Flame 150 (2007) 188-200.
[39] S. Humer, K. Seshadri, R. Seiser, Combustion of jet fuels and its surrogates in laminar nonuniform
flows, 5
th
US Combustion Meeting (2007), paper E09.
[40] A. Ansari, F.N. Egolfopoulos, Further considerations of flame ignition in counterflow configuration,
Western States Section of the Combustion Institute Spring Meeting (2014), paper 087LF-0051.
[41] K.B. Brady, X. Hui, C.J. Sung, K.E. Niemeyer, Counterflow ignition of n-butanol at atmospheric
and elevated pressures, Combust. Flame 162 (2015) 3596-3611.
[42] O. Park, P.S. Veloo, N. Liu, F.N. Egolfopoulos, Combustion characteristics of alternative gaseous
fuels. Proc. Combust. Inst. 33 (2011) 887-894.
[43] Y. Dong, C.M. Vagelopoulos, G.R. Spedding, F.N. Egolfopoulos, Measurement of laminar flame
speeds through digital particle image velocimetry: mixtures of methane and ethane with hydrogen,
oxygen, nitrogen, and helium, Proc. Combust. Inst. 29 (2002) 1419–1426.

63
[44] C.S. Yoo, J.H. Chen, J.H. Frank, A numerical study of transient ignition and flame characteristics
of diluted hydrogen versus heated air in counterflow, Combust. Flame 156 (2009) 140-151.
[45] X.L. Zheng, J.D. Blouch, D.L. Zhu, T.G. Kreutz, C.K. Law, Ignition of premixed hydrogen/air by
heated counterflow, Proc. Combust. Inst. 29 (2002) 1637-1643.
[46] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, A computational tool for the detailed kinetic
modeling of laminar flames: Application to C 2H4/CH4 coflow flames, Combust. Flame 160 (2013)
870-886.
[47] http://www.openfoam.org.
[48] H.G. Weller, G. Tabor, H. Jasak, and C. Fureby, A tensorial approach to computational continuum
mechanics using object-oriented techniques, Computers in Physics 12 (1998) 620–631.
[49] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, Numerical modeling of laminar flames with detailed
kinetics based on the operator-splitting method, Energ. Fuel. 27 (2013) 7730-7753.
[50] V. Gururajan, F.N. Egolfopoulos, Direct numerical simulations of probe effects in low-pressure
flame sampling, Proc. Combust. Inst. 35 (2015) 821-829.
[51] H.W. Wang, X. You, A.V. Joshi, S.G. Davis, A. Laskin, F.N. Egolfopoulos, C.K. Law, USC-Mech
Version II. High-temperature combustion reaction model of H 2/CO/C1-C4 Compounds.
http://ignis.usc.edu/USC-Mech-II.htm
[52] A.E. Lutz, R.J. Kee, J.F. Grcar, F.M. Rupley, OPPDIF: A FORTRAN program for computing
opposed-flow diffusion flames, Snadia Report, SAND96-8243, Sandia National Laboratories, 1997.
[53] M. Nishioka, C.K. Law, T. Takeno, A flame-controlling continuation method for generating S-curve
responses with detailed chemistry, Combust. Flame 104 (1996) 328–342.
[54] F.N. Egolfopoulos, P.E. Dimotakis, Non-premixed hydrocarbon ignition at high strain rates, Symp.
(Int.) Combust. 27 (1998) 641–648.
[55] C.M. Vagelopoulos, F.N. Egolfopoulos, Direct experimental determination of laminar flame speeds,
Symp. (Int.) Combust. 27 (1998) 513-519.
[56] S.M. Hosseinalipour, A.S. Mujumdar, Flow and thermal characteristics of steady two dimensional
confined laminar opposing jets: Part I. Equal jets, Int. Commun. Heat Mass 24 (1997) 27–38.
[57] S. Devahastin, A.S. Mujumdar, A numerical study of flow and mixing characteristics of laminar
confined impinging streams, Chem. Eng. J. 85 (2002) 215-223.
[58] S.J. Wang, S. Devahastin, A.S. Mujumdar, Effect of temperature difference on flow and mixing
characteristics of laminar confined opposing jets, Appl. Therm. Eng. 26 (2006) 519–529.
[59] P.A. Monkewitz, D.W. Bechert, B. Barsikow, B. Lehmann, Self-excited oscillations and mixing in
a heated round jet. J. Fluid Mech. 213 (1990) 611-639.

64
[60] G. Coppola, A. Gomez, Experimental study of highly turbulent isothermal opposed-jet flows, Phys.
Fluids 22 (2010) 105101.
[61] R.K. Shah, A.L. London, Laminar Flow Forced Convection in Ducts: a source book for compact
heat exchanger analytical data. Vol. 1. Academic press, 1980, pp. 78-152.
[62] K. Seshadri, F.A. Williams, Laminar flow between parallel plates with injection of a reactant at high
Reynolds number, Int. J. Heat Mass Transf. 21 (1978) 251–253.
[63] C.B. Oh, A. Hamins, M. Bundy, J. Park, The two-dimensional structure of low strain rate
counterflow nonpremixed-methane flames in normal and microgravity, Combust. Theor Model 12
(2008) 283-302.
[64] C. Godrèche, P. Manneville, Hydrodynamics and Nonlinear Instabilities, Cambridge, UK, 2005, pp.
164-193.
[65] X.L. Zheng, J.D. Blouch, D.L. Zhu, T.G. Kreutz, C.K. Law, Ignition of premixed hydrogen/air by
heated counterflow, Proc. Combust. Inst. 29 (2002) 1637-1643.
[66] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5
(1968) 506-517.
[67] R.I. Issa, Solution of the implicitly discretized fluid flow equations by operator-splitting, J. Comput.
Phys. 62 (1986) 40-65.
[68] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, OpenSMOKE: numerical modeling of reacting
systems with detailed kinetic mechanisms, XXXIV Meeting of the Italian Section of the Combustion
Institute, Rome, Italy, 2011.
[69] F.N. Egolfopoulos, P.E. Dimotakis, Non-premixed hydrocarbon ignition at high strain rates, Symp.
(Int.) Combust. 27 (1998) 641–648.
[70] C. Hindmarsh, ODEPACK: A systematized collection of ODE solvers, North-Holland, Amsterdam,
1983, pp. 55 −64.
[71] A.E. Lutz, R.J. Kee, J.A. Miller, SENKIN: A FORTRAN program for predicting homogeneous gas
phase chemical kinetics with sensitivity analysis. SAND-87-8248, Sandia National Laboratories,
1988.
[72] W.L. Briggs, S.F. McCormick, A multigrid tutorial, SIAM, Denver, Colorado, 2000, pp. 119-133.
[73] W.L. Oberkampf, C.J. Roy, Verification and validation in scientific computing, Cambridge
University Press, 2010, pp. 309-329
[74] H. Schlichting, Boundary-layer theory, New York, 1968
[75] T.C. Papanastasiou, G.C. Georgiou, A.N. Alexandrou, Viscous fluid flow, CRC Press, 1999.

65
[76] H.S. Dou, B.C. Khoo, H.M. Tsai, Determining the Critical Condition for Flow Transition in a Full-
Developed Annulus Flow, arXiv preprint physic/0504193 (2005).
[77] H. Tsuji, Counterflow diffusion flames, Prog. Energy Combust. Sci 8 (1982) 93-119.
[78] C. Walchshofer, H. Steiner, G. Brenn, Robust outflow boundary conditions for strongly buoyant
turbulent jet flames. Flow, turbulence and combustion, 86 (2011) 713-734.
[79] A. Mani, Analysis and optimization of numerical sponge layers as a nonreflective boundary
treatment. J. Comput. Phys. 231 (2012) 704-716.
[80] B.M. Cetegen, Y. Dong, M.C. Soteriou, Experiments on stability and oscillatory behavior of planar
buoyant plumes. Physics of Fluids 10 (1998) 1658-1665.
[81] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid
Flows, Ph.D. Thesis, 1996.
[82] C.R. Shaddix, Correcting thermocouple measurements for radiation loss: a critical review,
Proceedings of the 33rd National Heat Transfer Conference (1999), Paper HTD99-282.
[83] R.D. Blevins, Flow induced vibrations classifications and lessons from practical experiences, New
York, 2008, pp. 1-44.

66
Chapter 3
Effects of confinement, geometry, inlet velocity profile, and Reynolds
number on the asymmetry of opposed-jet flows
3.1. Introduction
Opposed-flow, impinging, axisymmetric jets form a counterflow, as schematically shown in Fig. 3.1a.
The two important length scales are the jet nozzle diameter, D, and the nozzle separation distance, L,
which can be combined into an aspect ratio,   L/D. Between the nozzles exists a stagnation-type flow
field (e.g., [1]) whereby the stagnation plane is defined for all points with zero axial velocity, uz.
Furthermore, the stagnation point is defined as the single location with both zero uz and zero radial
velocity, ur. This configuration is used in various applications [2], for example to facilitate the
measurement of global flame properties such as ignition temperatures, laminar flame speeds, and
extinction strain rates (e.g., [3-5]), as well as the characterization of the flame structure [6]. Often, these
counterflow flame measurements are modeled using a quasi one-dimensional (1-D) formulation where the
flow field in the impingement region between the jets is assumed to be steady. The implicit assumption
of the 1-D formulation is that D  ∞ and thus  → 0 for any L, so that the interaction of the jets with the
surrounding can be ignored. If uz and all scalar quantities exhibit zero radial gradients along the jets axis
of symmetry, then the radial and axial momentum may be decoupled and only the domain along the axis
of symmetry needs to be modeled. The distribution of uz, depicted in Fig. 3.1b and the scalars along the
axis of symmetry can be obtained from 1-D calculations for a given set of boundary conditions (BCs) at
the jet exits [7]. The applicability of the 1-D formulation in flame experiments has been assessed in several
studies (e.g., [8-16]), which addressed the effects of BC, two-dimensionality, and nozzle design.

67

Fig. 3.1. (a) Schematic of a typical axisymmetric counterflow; (b) typical axial velocity distribution
along the symmetry axis; K is the strain rate at the stagnation point.

Since the assumption of a steady flow field underpins the 1-D formulation, it must be validated before
any comparison can be made to experimental results. The hydrodynamic stability of the counterflow
configuration has been the subject of previous studies (e.g., [17-23]), which revealed flow field regimes
with notably more complex dynamics than is permitted under the 1-D formulation. These studies can be
classified into two categories, confined configurations that generate recirculation zones and unconfined
configurations that are free of recirculation zones.
For counterflow configurations in confined chambers, such as flame experiments in a pressure
chamber (e.g., [5,24,25]), the presence of solid surfaces near the jet exits causes the flow to separate and
generate recirculation zones subject to instabilities (e.g., [26,27]), possibly leading to bifurcation of the
flow field (e.g., [16,19,20]) due to the Coanda effect [28]. In a recent study [16], it was shown that such
solid surfaces must be removed from the vicinity of the jets in order for the flow field to satisfy the 1-D
formulation assumptions. However, such surfaces are incorporated in both experiments (e.g., [16,18,29])
and numerical simulations (e.g., [10,19,30]). Rolon et al. [18] showed experimentally that for
counterflowing ambient temperature air jets with matched momenta, the stagnation plane might not be
located at the middle plane between the nozzles. Recirculation zones were observed that were eliminated
by removing the solid surfaces from the vicinity of jets exits. However, the stagnation plane could not be
stabilized at the middle plane in either case. It was further shown that the offset of the stagnation plane
location increases with the L. Pawlowski et al. [19] performed numerical bifurcation analysis on confined
isothermal counterflowing jets and found that the flow field bifurcates at a critical Reynolds number, Recrit
-1.0
-0.5
0.0
0.5
1.0
0 0.5 1
u
z
z
r
D
L
Bottom nozzle
Stagnation
Plane
Top nozzle
z
Symmetry
Axis
K = ( u
z
/ z)
(a) (b)
Stagnation
Point

68
based on D, and the average nozzle exit velocity, U. The instability was identified as “pitchfork
bifurcation” where exchange of stability occurs between the symmetric and asymmetric solutions [19].
Figure 3.2 depicts a schematic stability diagram for confined counterflow jets. For balanced jet
momenta producing a flat stagnation plane shape, the symmetric solution occurs when the stagnation point
rests at the mid-plane while bifurcated solutions occur for all other stable stagnation point positions. States
for which the stagnation point is closer to the lower nozzle is denoted as the “Lower Branch” while that
with the stagnation point closer to the upper nozzle is called the “Upper Branch” [19,20]. For the same
BCs, two stable states of the flow field were observed after pitchfork bifurcation, which were mirror
images with respect to the mid-plane. The magnitude of the stagnation point offset has been reported to
increase with increasing Re and  [19,20]. The flow structure has been reported also to be particularly
sensitive to the surrounding geometry [28] leading to coupling of the bifurcated state to small geometrical
asymmetries. In this case, the flow field displayed preference to one of the bifurcated states and the less
dominant state was observed only in a small range of Re [28]. Such a significant sensitivity to small
imperfections is indication of an “imperfect bifurcation” (e.g., [31]). Figure 3.2 also depicts schematics of
a typical imperfect bifurcation map, which shows that the two solutions are no longer mirror images with
respect to the mid-plane.

69

Fig. 3.2. Schematic of a typical stability diagram for confined counterflowing jets; black: axial location
of the stagnation point position in a system with a pitchfork bifurcation [19,20]; blue: axial location of the
stagnation point in a system with an imperfect bifurcation [28]; solid lines represent stable solutions; the
two stable solutions are distinct after the bifurcation; (----) unstable solution; (─ ∙ ─) nozzle axial
locations; red: upper nozzle; green: lower nozzle.
Counterflow configurations can be further divided into two categories based on the value of . For
 ≥ 2, the flow field exhibits characteristics of a free jet due to the large L and is called “free-floating”
(e.g., [12,14,32,33]). The 1-D formulation is not applicable in this regime, so flame studies are often
performed for  < 2 so that the 1-D formulation assumptions are approached experimentally.
For unconfined counterflowing jets free of influence by nearby surfaces, there are no recirculation
zones but there is entrainment [16,18,34]. Stability analysis for this type of flows has been performed
experimentally for a wide range of  but there are discrepancies among the results. Three distinct flow
regimes depending on Re and  have been identified based on the behavior of the stagnation point position
(e.g., [18,21-23,32,34-37]). For jets with equal momenta, the three regimes are: (1) the steady, symmetric
regime at the mid-plane; (2) the steady, offset regime away from the mid-plane; and (3) the unsteady
Upper Branch
Lower Branch
(Re
crit
)
pitchfork
Axial location, z
Upper Branch
Lower Branch
(Re
crit
)
imperfect
Re
Lower Nozzle
Upper Nozzle

70
regime characterized by oscillation about the mid-plane. Rolon et al. [18] reported a steady, offset regime
for  = 1 and Re ≈ 1026 with the stagnation point shifted away from the mid-plane by 6% of D. In another
study [23], the unsteady regime was reported for  = 1 and Re ≥ 786. Scribano et al. [38] reported a steady,
symmetric regime for 0.5 ≤  ≤ 2 and 200 ≤ Re ≤ 4800. The reported numerical results are inconsistent
with each other and disagree with the experiments. For example, outside the free-floating regime and in
the absence of a flame, a stagnation point offset has not been reported [1,13,38]. Kelvin-Helmholtz
instabilities have been observed at the jet edges for Re = 2000 and  = 1 but the stagnation point and the
flow field in the impingement region were unaffected [1]; such instabilities have been reported also
experimentally [21,22]. No instability has been reported numerically, either for ambient air jets with
 = 0.5 and 390 < Re < 2600 [13] nor in the presence of a flame for 0.5 ≤  ≤ 1.23 and 190 < Re < 2900
[9,14].
In some cases, the flow field was found to be asymmetric even when symmetric unconfined jets were
established [18,39], which indicates notable sensitivity of the flow field symmetry to the surrounding
geometry. Yet another nuance of the counterflow dynamics is the high sensitivity of the stagnation point
offset to imbalanced jet momenta [21,22] which decreases with Re and increases with . This offset has
been reported to depend on the jet exit velocity profile shape [1,22].
The discrepancies between experimental and numerical results can be justified based on the non-linear
dynamics of the counterflow configuration, which result in sensitivity to surrounding confinement, nozzle
design, and flow BCs. Since absolute perfection in system construction and alignment cannot be achieved
(e.g., [18,38,39]), real experiments likely suffer from design-specific dependencies.
Previous studies have not addressed the effects of different confinement levels on the stability of the
axisymmetric counterflow jets despite its pivotal role nor have they provided the sensitivity of the results
to imperfections in nozzle assembly, surroundings, and momentum balance procedure outside the free-
floating regime. Numerical stability analysis has been only performed in the confined case considering
only one or two variations of counterflow surrounding geometries [19,20]. To the author’s knowledge,
numerical stability analysis of axisymmetric ambient-temperature air opposed jets without the influence
of far field walls has never been performed. Inherent instabilities of axisymmetric counterflow jets [17]
are also rarely discussed in the context of flame experiments (e.g., [13,38]). This chapter focuses on the
stability of axisymmetric, laminar, non-reacting counterflow jets with respect to the underlying
assumptions of the 1-D modeling formulation. Investigating the non-reacting jets is an essential first step
in counterflow flame experiments (e.g., [12,13,32,39]) given that the presence of the flame could introduce
additional complexities that need to be addressed subsequently. The objectives of the present chapter are

71
to use experimental measurements and computational results to: (1) identify instabilities and asymmetric
flow fields; (2) identify the critical parameters influencing such behavior; (3) understand the underlying
physical mechanisms; (4) quantify the consequences on the validity of underlying assumptions of the 1-
D formulation; and (5) suggest guidelines in experimental design to maximize conformity with 1-D
modeling assumptions.

72
3.2. Experimental approach
Experiments were performed at ambient pressure and temperature in a counterflow assembly shown
schematically in Fig. 3.3 and which has been used in previous studies (e.g., [5,15]). It consists of two
identical and carefully aligned, axisymmetric, straight-tube pipes fixed within a cylindrical pressure
chamber. Windows mounted on the sidewalls of the chamber provided optical access.
L = D = 2.80 ± 0.01 cm was chosen for all experiments so that  = 1. Air was supplied through calibrated
sonic nozzles and directed into the feeding tubes. Two 80-mesh screens were placed inside the nozzles as
flow conditioners, one 5D and the other 7D upstream of the nozzle exit. The only solid surfaces near the
jet nozzles were removable caps, as shown in Figs. 3.3b and 3.3c, which provided flexibility to assess the
sensitivity of measurements to the jets’ surroundings [18]. With the caps removed, the only solid surfaces
to exert influence on the flow structure were the outer surface of the jet nozzles and very mild impact is
expected from the chamber walls that are located 10D away from the mid-plane. It should be emphasized
that no annular coflow was used in this study.
Particle Image Velocimetry (PIV) was used to characterize the two-dimensional (2-D) flow field.
Micrometer-scale silicone oil droplets were used as flow tracers and were generated by a glass pneumatic
nebulizer (Meinhard HEN-170-A0.1). A pair of ultra-compact folded resonator Nd:YAG lasers emitted
light at 532 nm. An optical system was set to shape and direct the laser beam to the measurement region.
Scattered light was recorded by a 12-bit charge-coupled device camera with 1392×1024 pixels of
resolution (Pixelfly 270XD) at the maximum frequency of 5.3 Hz per image pair. LaVision's DaVis 7.2
software was used for timing, synchronization, image acquisition, and processing. Nominal velocity
values are reported as the mean of 100 image pairs and velocity uncertainty as the two  standard deviation
from the mean value. More details can be found in previous studies (e.g., [5,40]).
The flow field in the impingement region was characterized in the range of 250 < Re < 1800, which is
typical in flame studies (e.g., [3,5,6,8,12,15]). Experiments were first performed with equal jet momenta
as verified by calculating the total volumetric flow rate through each nozzle obtained from calibrated flow
rates through the sonic nozzles. Since the uncertainty in the PIV measurements and the reconstruction of
the volumetric flow rate from PIV was higher than the one obtained by sonic nozzle calibration, the later
was used to control the momenta accurately. However, it should be noted that the volumetric flow rates
obtained by both methods matched within the uncertainties. At last, the average velocity of the bottom jet
is varied within 10% of the balanced value to investigate the importance of absolute momentum balance
on flow structure. Such imbalance can be encountered in flame experiments due to the uncertainties
involved [16]. The stagnation point and stagnation plane positions were used as global markers of the flow

73
field behavior (e.g., [18-23]). Results are presented in terms of the average exit velocity ratios U1/U2 where
U1 and U2 are the average velocity at the exit of the bottom and top nozzles respectively. U1 and U2 were
obtained by dividing the volumetric flow rate by the nozzle exit area.
As the stagnation point offset outside the free-floating regime can be an order of magnitude smaller
than D [16-19,32,35], a rigorous uncertainty quantification was required to resolve any offset given
measurements uncertainties. The momentum balance uncertainty can be expressed by the uncertainty of
U1/U2. U1 and U2 uncertainties were first calculated from four sources: (1) the uncertainty in the volumetric
flow rate measurements when calibrating the sonic nozzles; (2) the scatter in the calibration data; (3) the
uncertainty of the pressure reading from pressure gauges; and (4) the uncertainty in D measurement. The
relative uncertainties of these quantities were determined to be 0.05%, 1.00%, 0.30%, and 0.30%
respectively in the worst-case scenario. Since two feed lines were used for each nozzle, one for the main
flow and one for the seeding flow as demonstrated in Fig. 3.3a, the first three sources of uncertainty had
to be accounted for twice for each nozzle. Therefore, the total uncertainty of U1 and U2 was obtained to
be 1.54%, which was calculated based on the assumption that the four sources of uncertainty were
independent. Hence, the present uncertainty is the most conservative estimate. The uncertainty of U1/U2
was calculated to be 2.17% since U1 and U2 were controlled independently.
The uncertainty of the stagnation point location was derived based on the uncertainties of the axial
position and the axial velocity. The main contributors to the overall uncertainty of the stagnation point
were the large standard deviation of the velocity field near the stagnation point and the systematic errors
in determining the location of lower nozzle exit, which was chosen as the origin of axial coordinate. The
first source of uncertainty for the axial position originated from the fact that the lower nozzle exit appeared
to be not completely flat in the PIV images. Although the camera was aligned with respect to the nozzle
exits, the left corner of the bottom nozzle rim appeared to be higher than the right corner. Moreover, due
to the laser reflection from the nozzles, the nozzle rims appeared as intensity-saturated regions in the PIV
images. The axial location of the upper most saturated pixel at the left corner of the bottom nozzle rim
was chosen as the axial location of the origin. Such choice of a single pixel to represent an area of saturated
pixels (approximately 10 pixels) introduced uncertainties. Furthermore, the chosen pixel was found to
change location from one image pair to another. However, the fluctuations were quantifiable such that an
average axial location could be defined for the origin. These systematic errors are usually avoided in flame
experiments by choosing an arbitrary point in the measurement field as an origin (e.g., [3,5,12,15,38,40]).
In this study, having a fixed and well-defined reference point was essential in determining any asymmetry
and such errors could not be eliminated. The lower nozzle flatness uncertainty and the uncertainty due to
the reflection were determined to be 0.03 cm and 0.04 cm respectively which led to the total uncertainty

74
of 0.05 cm in the axial location. In order to determine the stagnation point location, the symmetry axis
was first determined using four-point interpolation based on the ur = 0 criterion. The axial location of the
stagnation point was then determined using another four-point interpolation based on the uz = 0 criterion
on the symmetry axis. The uncertainty of the stagnation point location was obtained by the weighted least-
square regression method [41] using the known uncertainty of the measured velocities in the vicinity of
the stagnation point. The final uncertainty was the combination of the uncertainties in determining the
lower nozzle exit axial location and the obtained value in the previous step assuming they were
independent quantities.

Fig. 3.3. Schematic of the counterflow system used in the experiment; (a) overview; (b) nozzle details
with the presence of the cap; (c) nozzle details after removing the cap.

75
3.3. Modeling approach
3.3.1. 1-D simulations
1-D simulations were performed using an opposed-jet code [7,42]. The distribution of uz along the
symmetry axis was obtained by specifying the jet exit axial velocities at the symmetry axis as BCs. The
current version of the code is also capable of applying “mixed BCs” [12]. In this formulation, the measured
axial velocity gradient in the axial direction, uz/ z, and the measured uz at nozzle exits and at the
symmetry axis are used directly as BCs. Thus, extending the applicability of the code for nozzles designs
that do not conform to “plug flow” assumption i.e., uz/ z = 0 at the nozzle exits (e.g., [8,12-14]). This
code has been used in many previous counterflow flame studies (e.g., [3-6,8,10,12,14-16,24]).

3.3.2. 2-D simulations
3.3.2.1. Transient flow simulations
Transient, 2-D, axisymmetric, incompressible flow simulations were performed for ambient
temperature air jets. The solver is based on the icoFoam application of the OpenFOAM suite [43] and has
been modified and used previously [44,45]. The mass and momentum conservation equations were solved
using a finite volume formulation and Pressure Implicit with Splitting of Operators (PISO) method [46].
A second-order implicit time differencing scheme was used for stability, and spatial derivatives were
second-order accurate. The time step was chosen such that the maximum Courant number [47] in the mesh
remained at or below 0.4. Within each time step, four outer iterations were performed to update the
velocity field. Within each outer iteration, three inner iterations were performed to correct the pressure
field. This process ensures that the solution satisfied both the Navier-Stokes and the continuity equations
and kept the introduced error within acceptable limits

3.3.2.2. Bifurcation analysis
For counterflowing jets with equal momenta, the characteristic flow speed is U and geometrical length
scale is D. These can be used to non-dimensionalize the Navier-Stokes equations and the Navier-Stokes
operator, 𝑁 , can then be defined as:
𝑁𝒖 = −𝒖 . ∇𝒖 − ∇𝑝 +
∇
𝒖 (3.1)
where 𝒖 (𝒙 , 𝑡 ) is the non-dimensional velocity vector field, and 𝑝 (𝒙 , 𝑡 ) is the pressure field determined by
the continuity equation, ∇. 𝒖 = 0. Thus, the Navier-Stokes equation is given by:
𝜕𝒖 𝜕𝑡 ⁄ = 𝑁𝒖 (3.2)

76
For a fixed geometry and set of BCs, Re is the only parameter dictating the flow physics. A velocity field
𝒖 (𝒙 ), along with its associated pressure field 𝑝 (𝒙 ), are an equilibrium or steady-state solution of 𝑁 if
and only if 𝑁𝒖 = 0. For a perturbation velocity field, 𝒖 ′(𝒙 , 𝑡 ), and perturbation pressure field, 𝑝 ′(𝒙 , 𝑡 ),
the linearization of 𝑁 about 𝒖 yields the linearized Navier-Stokes operator, ℒ, given by:
ℒ𝒖 = −𝒖 . ∇𝒖 − 𝒖 . ∇𝒖 − ∇𝑝 +
∇
𝒖 (3.3)
and 𝑝 ′ is determined by the continuity equation, ∇. 𝒖 ′ = 0. Furthermore, ℒ possesses a homogeneous BC
since the perturbed flow field given by 𝒖 + 𝒖 ′ and 𝑝 + 𝑝 ′ must obey the same BCs as the equilibrium
flow given by 𝒖 and 𝑝 . With these conditions, the linearized Navier-Stokes equation can be expressed
by:
𝜕𝒖 ′ 𝜕𝑡 ⁄ = ℒ𝒖 ′ (3.4)
Eigendecomposition of the linearized Navier-Stokes operator can then be performed to find the
eigenvalues and eigenmodes, which determine the long-time behavior of the linearized dynamics for a
given geometry and Re. The authors have developed the aforementioned approach for bifurcation analysis
details of which can be found in previous studies [44,45]. In short, the procedure begins by specifying
geometry and Re and then finding the steady-state solution, 𝒖 that satisfies 𝑁𝒖 = 0 . This is
accomplished by employing a doubly iterative procedure. In the outer iteration, Newton iteration with the
Armijo rule [48,49] was used to solve for a root of 𝑁 over the vector field of spatially discretized velocity
fields, 𝒖 . The iteration begins with an initial guess 𝒖 for the true root 𝒖 . For stable base flows, this guess
was determined from a transient simulation run for long times. The Generalized Minimal Residual
(GMRES) method [50,51] was used to solve the Jacobian-vector equation
𝒖 |
𝒖 𝒉 = −𝒖 for the Newton
step 𝒉 in the inner j
th
Newton iteration. A “time-stepping" method [51] was utilized to approximate
𝒖 |
𝒖 𝒉 , where the finite difference
𝒖 |
𝒖 𝒉 =
𝒖 𝒉 𝒖 + 𝑂 (𝜖 ) was employed with 𝜖 = 10
in
the GMRES iteration, and the 𝑂 (𝜖 ) term was further truncated. The solution was considered temporally
converged when the norm of 𝑁𝒖 scaled by the volume of the computational domain was less than
5 × 10
.
Once a satisfactory approximation of the steady-state solution 𝒖 was obtained, 100 Arnoldi iterations
[52] with discrete-time variants of ℒ were performed to estimate 100 leading eigenvalues and eigenmodes
of ℒ. By tracking the most unstable eigenvalue as Re increases, the bifurcation point and type were
identified. Recrit was obtained when the real part of the leading eigenvalue was zero. These OpenFOAM-
based operators interfaced with a separate in-house library that performed the Newton, Armijo, and

77
GMRES iterations. The code was highly parallelized through the domain decomposition method on
distributed-memory machines.

3.3.2.3. Geometry, initial and boundary conditions and cases investigated
2-D simulations and bifurcation analysis were performed in the counterflow configuration for four
different geometries as shown schematically in Fig. 3.4. The  = {0.5,1} values were used to represent a
typical range for flame experiments (e.g., [3,5,12,13,15,16,24]). Uniform and parabolic inlet velocity
profiles were chosen as two limiting cases encountered in experiments. Ideally, the nozzle exit velocity
profile must be uniform to conform to the assumptions of the 1-D formulation. However, in some designs
(e.g., [3-5,38]) a boundary layer (BL) can form, which tends to develop a parabolic profile as the flow
approaches the nozzle exit. Three levels of confinement were investigated by changing the vertical spacing
between disk walls and the nozzle exits, H. For “unconfined” simulation cases, the disk walls were
replaced by open boundaries. H values of 0 and 4 were used for “fully confined” and “semi-confined”
cases, respectively.

78

Fig. 3.4. Schematics of the geometries investigated numerically; solid boundary lines represent walls
Geometry (a) represents a fully confined opposed jet, H = 0, with a uniform inlet velocity profile.
Three levels of “confinement imperfection”, w, were considered; only the bottom disk wall had an offset
of w with respect to the bottom nozzle exit. Geometry (b) features less confinement than geometry (a) and
offers the flexibility to assess the effects of asymmetries in the nozzle rims that can be encountered
experimentally. Two values of the bottom nozzle rim thickness, δbottom, were considered in geometry (b),
that is δbottom = 0 and 5% of D while the top nozzle rim thickness, δtop, was assumed to be zero. Geometry
δ
bottom
(a) (b)
(c)
z
r
O utlet
0.5
w
𝒓 𝑶𝒖𝒕𝒍𝒆𝒕 1 O utlet
O utlet / Disk Wall
𝒓 𝑶𝒖𝒕𝒍𝒆𝒕 z
r
H
H
0.5
O utlet / Disk Wall
Disk Wall
Disk Wall
1
O utlet
O utlet
𝒓 𝑶𝒖𝒕𝒍𝒆𝒕 z
r
H
H
0.5
O utlet
Boundary Layer
z
r
0.5
0.5
O utlet / Disk Wall
𝒓 𝑶𝒖𝒕𝒍𝒆𝒕 H
O utlet / Disk Wall
O utlet
(d)
H
Symmetry
Axis
L
bottom
δ
top
L
top
δ
bottom
δ
top
δ
bottom
δ
top
L
bottom
L
top

79
(c) represents the actual experimental configuration. It is similar to geometry (b) with the difference that
the computational domain is extended into the feeding nozzles to allow investigation of the effects of BL
development as Re varies (e.g., [38]) and pressure field feedback into the nozzle (e.g., [12,15,16,18,53]).
In order to model the experimental geometry accurately, the nozzle rim thicknesses were first measured
using a micrometer (δtop ≈ δbottom = 0.14 ± 0.01 cm) in order to be directly implemented in the numerical
geometry. However, nozzle misalignments could generate geometrical asymmetries. Although the nozzles
were aligned with an accuracy of 0.05 cm in the PIV plane, they were aligned in the plane perpendicular
to the PIV plane without a visual aid. Therefore, misalignment could still persist in the later plane. The
misalignment and δ measurement uncertainties were modeled by specifying a value for δbottom that is larger
than δtop by 2.5% of D. Other geometrical asymmetries that could be encountered in the experiments were
also considered. Specifically, since the screens were placed manually by forcing them down the nozzle
from the jet exit, there existed an uncertainty in the axial location of the screens such that it would render
the nozzles asymmetric in terms of screen axial locations. The possibility of such error in manufacturing
was modeled by considering different computational domain extensions inside the tubes, Ltop and Lbottom
in geometry (c).
Finally, by replacing the inlet uniform velocity profiles of geometry (c) with parabolic ones, geometry
(d) is obtained. Table 3.1 summarizes all the cases investigated computationally and their main features.
A naming convention was utilized to distinguish between the various cases. The letters have the following
indication in the order that they appear in the name of various cases: (1) the confinement level; (2) the
pressure feedback effect; (3) the velocity profile shape; (4) ; (5) the symmetry of the geometry; and (6)
w. For example, in the CWU0.5S1.25 case, the letters have the following meaning: C ≡ Confined
geometry, W ≡ Without pressure feedback into the nozzles, U ≡ Uniform velocity profile, 0.5 ≡  = 0.5,
S ≡ Symmetric geometry, and 1.25 ≡ w = 1.25%. Similarly, in the UFP1A0 case: U ≡ Unconfined
geometry, F ≡ Feedback into the nozzles are included, P ≡ Parabolic velocity profile, 1 ≡  = 1,
A ≡ Asymmetric geometry, and 0 ≡ w = 0.

80
Table 3.1. Summary of cases investigated computationally.
Case Geometry α w H δbottom δtop Lbottom Ltop 𝒓 𝐎𝐮𝐭𝐥𝐞𝐭
Confinement
level
Bifurcation
type
Recrit
CWU0.5S0
(a)
0.5 0 0 N/A N/A 0 0 20 confined pitchfork 153
CWU1S0
1
0 0 N/A N/A 0 0 20 confined pitchfork 90
CWU1A1.25 1.25% 0 N/A N/A 0 0 20 confined imperfect 106
CWU1A12.5 12.5% 0 N/A N/A 0 0 20 confined imperfect 217
SWU1S0
(b) 1
0 4 0 0 0 0 30 semi-confined pitchfork 90
UWU1S0 0 ∞ 0 0 0 0 40 unconfined None N/A
UWU1A0 0 ∞ 5% 0 0 0 40 unconfined None N/A
UFU1S0
(c) 1
0 ∞ 0 0 5 5 40 unconfined None N/A
UFU1A0 0 ∞ 7.5% 5% 4.5 5 40 unconfined None N/A
CFP0.5S0
(d)
0.5 0 0 N/A 0 2 2 30 confined pitchfork 157
CFP1S0
1
0 0 N/A 0 2 2 18 confined pitchfork 54
SFP1S0 0 4 0 0 2 2 30 semi-confined pitchfork 54
UFP1S0 0 ∞ 0 0 2 2 40 unconfined None N/A
UFP1A0 0 ∞ 5% 0 2 2 40 unconfined None N/A

The velocities were specified at the inlets of the computational domain and set to zero at the walls.
The normal derivative of the velocity to the boundary was set to zero at open outlet boundaries. The
pressure was specified to a constant value of unity at open outlet boundaries. The “consistent BC” [54]
was used on all other boundaries where the pressure BC is derived from the normal gradient of the
momentum equation and the known velocities at the boundary without any restricting assumption.
Initially, velocities and pressure were set to zero throughout the domain.
For confined and semi-confined cases, the radial position of the open outlet boundary from the
symmetry axis, 𝑟 , was chosen such that the flow in between the disc walls fully develops ( ur/ r ≈ 0).
Otherwise, the results were found to depend on 𝑟 . Since the flow development depends on Re,  and
H, different 𝑟 had to be implemented for different cases as outlined in Table 3.1. For example, for
CFP0.5S0 at Re = 2000, specifying 𝑟 < 20 would not lead to a steady state solution which was
speculated to be the result of the presence of the recirculation zones up to r = 10. For each case, when a
steady-state solution was not achieved, the simulation was repeated for larger 𝑟 to assure that the
results were not an artifact of the open boundaries location on the results. Furthermore, the simulation was
repeated with larger 𝑟 for each case for an obtained steady state at the maximum Re investigated to
ensure less than 1% change of “quantities of interest” (the stagnation plane location and the maximum
Frobenius norm of the ∇𝒖 tensor in the computational domain). The difference between the quantities of
interest for UFU1S0 at Re = 500 were found to be 0.6% for a domain with 𝑟 = 45 and H = 25 and a
smaller domain with 𝑟 = 40 and H = 20. Hence, 𝑟 = 40 and H = 20 was chosen for all

81
unconfined cases. The computational cost was found to increase with decreasing confinement levels and
increasing Re. Thus, in order to simulate higher Re cases, the confinement level was kept as high as
possible.
For confined and semi-confined cases, a uniform structured mesh was used to ensure sufficient
resolution in the shear layer formed between the disk walls, which extends to large 𝑟 due to presence of
recirculation zones. Such resolution was required as shear instabilities usually trigger in the shear layer
(e.g., [26,27]) far from the jets’ impingement region. For unconfined cases, the grid was uniform up to
r = 10 and z ≈ ±5 beyond which a gradual decrease in mesh density resulted in cell lengths near the outlet
approximately 20 times that of the finest cell in the domain.
The maximum Frobenius norm of the ∇𝒖 tensor in the computational domain was tracked through the
mesh refinement steps and the solution was considered spatially converged when the change was less than
5% after dividing the cell lengths by two. Non-dimensional spatial resolution of 1/160 was found sufficient
to resolve the quantities of interest. For CWU1S0, Recrit = 90 was found computationally in the study by
Pawlowski et al. [19] and was repeated here as a validation target for the bifurcation analysis method.
Such comparison is particularly useful since a higher order spectral method was used for discretization in
Pawlowski’s study. Here, Recrit = 90 was obtained for momentum balanced condition in agreement with
Pawlowski et al. [19]. Therefore, the present method was considered adequate to capture the physics of
interest.

82
3.4. Results and discussion
3.4.1. Experimental results
For impinging jets with balanced momenta, the flow field was found to be asymmetric with respect to
the middle plane between the nozzles at low Re i.e., in the steady stagnation point offset regime [23]. Only
one steady asymmetric state with the stagnation plane closer to the lower nozzle was observed in contrast
to previous experimental studies [21-23]. This is demonstrated for Re = 290 in Fig. 3.5, which shows a
measured velocity vector field overlaying the uz contour along with the locations of the stagnation point,
stagnation plane, and middle plane. The coordinate system originates at the center of lower nozzle exit, as
outlined in Fig. 3.1, and all quantities are in the non-dimensional format.

Fig. 3.5. Measured velocity vector fields in the momentum balanced condition; (a) Re = 290; (b)
Re = 1550; overlaid by contour maps of axial velocities and locations of (----) the middle plane between
the nozzles; (─) the stagnation plane; (─ ∙ ─) the stagnation plane uncertainty band; (∙ ∙ ∙) the stagnation
point uncertainty ellipse.

The asymmetry was observed to be higher at low Re and the stagnation plane was closer to the middle
plane as Re increased. The stagnation point offset was defined as the difference between the axial location
of the stagnation point and the middle plane. This offset was found to be -13% ± 2% and -2% ± 10% of
D for Re = 290 and Re = 1550 respectively. Relative velocity fluctuations near the stagnation point were
determined also to be higher than the ones at nozzle exits, consistent with the findings of Lindstedt et al.

83
[39], thus contributing to the overall uncertainty. The stagnation plane was found to be flat within the
uncertainty of the PIV measurements.
Figure 3.6 depicts the uz profiles at the exit of the nozzles for the conditions of Fig. 3.5 both for single
jets and opposed jet configurations; the single jet configuration was obtained by turning off the opposing
stream. The velocity profiles differed in single and opposed jet cases and the differences were more
significant at Re = 290. As first described by Rolon et al. [18], the opposed jet flows generate a stagnation
pressure field (e.g., [12,15,16,53]) which modifies the flow field structure relative to the single jet case
even at the jet exits. The profiles for the opposed jet configuration and at the exit of the bottom nozzle
exhibit larger dip at the symmetry axis due to the closer proximity of the stagnation point to the lower
nozzle. Figure 3.6 further confirms that at a given Re, while the flows coming out of each nozzle have
equal average velocities in the single jet case, the opposed jet configuration changes the flow field such
that the whole system re-establishes itself into a new equilibrium.
Moreover, by comparing the velocity profiles of both jets in the single jet configuration, a slight
asymmetry was realized. The BL thickness of the top nozzle was larger than the corresponding value for
the bottom nozzle as apparent from single jet configuration profiles shown in Fig.6 at Re = 290. This was
speculated to be caused by the asymmetry in the flow conditioning screen positions in the nozzles as
discussed in Section 3.3.2.3. To be specific, the screen right downstream of the bottom nozzle exit was
probably located closer to the nozzle exit compared to the corresponding screen in the top nozzle.
Therefore, the uniform flow generated right downstream of the screen developed more in the top nozzle
before exiting the nozzle compared to the flow in the bottom nozzle [38]. As a result, the top nozzle exit
profile was more parabolic than the bottom profile, which caused the stagnation point to settle closer to
the bottom nozzle. The asymmetry in the velocity profiles was found to diminish by increasing Re such
that it was not noticeable at Re = 1550 as shown in Fig. 3.6.

84

Fig. 3.6. Measured and computed axial velocity profiles at the exit of the nozzles in the momentum
balanced condition; (─) UFU1A0, Re = 250; (●) PIV, Re = 290; (▲) PIV, Re = 1550; blue: z = 0.96; red:
z = 0.02; open symbols: single jet configuration; closed symbols: opposed jet configuration; representative
error bars are also shown to aid the comparison.

It was observed that by removing the nozzle caps as outlined in Figs. 3.3b and 3.3c and by repeating
the experiment, the stagnation plane settled at a different axial location. Hence, the flow field was sensitive
to its surrounding even though no recirculation zones or attachments to the surrounding surfaces were
observed. Figure 3.7 shows the stagnation point offset as a function of U1/U2, for various Re. The slopes
of various cases indicated the sensitivity of the stagnation point position to the momentum imbalance. The
offset was 6% higher in the momentum balance condition (U1/U2 = 1) when the caps were present. Since
the asymmetry and uncertainty was higher in the presence of caps, as indicated by the offset and scatter
in the data points of Fig. 3.7, the results presented in Figs. 3.5 and 3.6 were obtained without caps. The
higher uncertainty and asymmetry in the presence of the caps was mainly due to the lower accuracy of the
alignment procedure.
u
z
r r
u
z
u
z
u
z
r r
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
-1.8
-1.6
-1.4
-1.2
-1.0
-0.8
-0.6
-0.4
-0.2
0.0
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5
-1.8
-1.6
-1.4
-1.2
-1.0
-0.8
-0.6
-0.4
-0.2
0.0

85
Since U1 and U2 values have been used as controlling parameters in counterflow flame experiments
(e.g., [8,13,25]), the 1-D formulation predictions were also included in Fig. 3.7. The 1-D results were
obtained by assuming plug flow conditions and utilizing the U1 and U2 values as BCs. While the 1-D
simulation results showed dependency only on U1/U2 and not on the exact values of U1 and U2, the
experimental results revealed otherwise. As Re increases, the change in the stagnation point offset was
observed to decrease for small deviations from the balanced momentum state as evident by the decrease
in the trendline slope as Re increases. This is consistent with the findings of Li et al. [21,22] for turbulent
conditions. Furthermore, the 1-D predictions were closer to the experimental trends for high Re. It should
be noted that when mixed BCs were used in the 1-D simulations, the predicted stagnation point position
(not shown in Fig. 3.7) matched the experimental value as already demonstrated in previous studies [8,12].
Therefore, the values of U1 and U2 cannot fully describe the behavior of the flow field in the current design
and using mixed BCs is essential to accurately predict the stagnation point position in unconfined nozzle
designs in which the jet exit velocity profiles are affected by the stagnation pressure.

86

Fig. 3.7. Stagnation point offset from the middle plane between the nozzles as a function of the average
nozzle exit velocity ratio of the jets; symbols: experimental data; solid black line: simulation result;
trendlines and representative error bars are also plotted to aid the comparison; the momentum balanced
condition is marked by the vertical dashed line and symmetry position by the horizontal dashed line.

-40
-35
-30
-25
-20
-15
-10
-5
0
5
10
15
20
0.8 0.9 1.0 1.1 1.2
Offset [%]
U
1
/U
2
Re = 290 (with caps)
Re = 290
Re = 600
Re = 1550
1-D (plug flow)

87
3.4.2. Numerical results
3.4.2.1. Momentum balanced conditions
Results of the numerical simulations for momentum balanced cases are summarized in Table 3.1.
Bifurcation was only observed for confined and semi-confined cases. Such behavior is speculated to be
the result of recirculation zones, which form only in the presence of disk walls [16,18]. Figure 3.8 depicts
ur contour maps overlaid by streamlines, and the results illustrate the flow field differences between a
confined and an unconfined case at the same Re corresponding to subcritical condition. For CWU1S0
which is a confined case, two recirculation zones were formed due to flow separation from the top and
bottom disk walls as shown in Fig 3.8a. Since the generated radial flow outside the impingement region
when confined by the disk walls is exposed to more shear due to the back flow in the recirculation zones,
it has higher tendency for local shear instabilities (e.g., [26,55]).

Fig. 3.8. Radial velocity maps overlaid by the streamlines at Re = 89 showing the effect of confinement
on the flow field; (a) CWU1S0 (confined); (b) UWU1S0 (unconfined).

Pitchfork bifurcation was observed for confined and semi-confined cases in symmetric geometries
(δ, w = 0). For Re > Recrit, two stable states were realizable which were mirror images with respect to mid-
plane, consistent with previous studies [19,20]. The same Recrit was obtained for confined and semi-
confined cases, as shown in Table 3.1, suggesting that the presence of the back flow is the cause of
instability and not the magnitude of the shear. This is consistent with findings of Pawlowski et al. [19] in
which Recrit was found to be the same for H = 0 and H = 2; when far field disk walls are present, no matter

88
how far they are from the main impingement region, they create recirculation zones, which are subject to
instability.
Figure 3.9 depicts the stagnation point offset, both computed and measured, as a function of Re for the
cases in which an asymmetric flow field is present. In confined and semi-confined cases, the stagnation
point was found to deviate from the symmetry position for Re > Recrit. The offset first increases with Re
as found in previous studies [19,20]. However, the increase in asymmetry reversed after a certain Re
whereupon the offset decreases with Re. Similar non-monotonic behavior of the bifurcation maps with Re
were also reported for channel flows [28]. The decrease in asymmetry for large Re was due to the size
increase of the recirculation bubbles. As the bubbles grow with Re [20], the asymmetry is pushed further
radially outward and the influence on the symmetry axis is reduced. Based on the results shown in Fig.
3.9, the computed offset for all cases are within the experimental uncertainty at high Re except for CFP1S0
for which the offset is five times bigger than the experimental uncertainty at Re = Recrit + 300.
Furthermore, the simulations for UFU1A0 captured the measured offset in the experiment. The velocity
profiles at the nozzle exits for UFU1A0 were also shown in Fig. 3.6, which demonstrates the agreement
between the measured and computed results. Hence, UFU1A0 properly approximates the experimental
geometry.

89

Fig. 3.9. Computed and measured stagnation point offset of the stagnation point from the middle plane
between the nozzles for momentum balanced condition as a function of Re for various cases; (----)
symmetry position; insets: enlarged view of the same plots shown to aid comparison; dashed area: regions
that cannot be resolved experimentally due to measurement uncertainty; the uncertainty found to increase
with Re due to higher measurement standard deviation [39].

Considering that the 1-D formulation cannot capture the bifurcation, by invoking the plug flow
assumption it predicts zero stagnation point offset in the momentum balanced condition. It is expected
that by implementing mixed BCs obtained from 2-D simulation results the predictions of 1-D formulation
improve. Nevertheless, simulation results revealed that predictions by the 1-D formulation can be wrong
by as much as 23% in the worst-case scenario with maximum flow asymmetry, i.e., CFP1S0 and Re = 348;
for this case, even using mixed BCs increased the discrepancy to 27%. Such discrepancy is still significant
considering that the experimental uncertainty is of the order of 2% for these conditions. This was not the
-25
-20
-15
-10
-5
0
5
10
15
20
25
0 100 200 300 400 500 600 700
-0.06
-0.03
0
0.03
0.06
50 150 250 350
-25
-20
-15
-10
-5
0
5
10
15
20
25
0 500 1000 1500 2000
-0.6
-0.3
0
0.3
0.6
0 200 400
-25
-20
-15
-10
-5
0
5
10
15
20
25
0 500 1000 1500 2000
-3
-2
-1
0
1
2
50 150 250
(a) (b) (c)
Offset, Uncertainty [%]
Re Re Re
CWU0.5S0
CWU1S0
CWU1A1.25
CWU1A12.5
Experimental Uncertainty
SFP1S0
CFP0.5S0
CFP1S0
UFP1A0
Experimental Uncertainty
SWU1S0
UWU1A0
UFU1A0
Experimental Uncertainty
PIV Data

90
case in the experiments in which the flow field is unconfined. Hence, the validity of the 1-D formulation
results decreases by increasing confinement levels and comparing the experimental results to the 1-D
formulation predictions can be misleading when operating in the Re range of hundreds (e.g., [4,56]) for
confined geometries. In semi-confined and unconfined geometries, discrepancies are lower particularly if
the comparison is performed for higher Re as also discussed in Section 3.4.1. This Re effect can be
attributed to the nature of the Navier-Stokes equations that change from elliptic at low Re to hyperbolic at
high Re [57]. In other words, the solution of the Navier-Stokes equations is sensitive to BCs at low Re as
also demonstrated by Scribano et al. [38] such that small change in surrounding geometry and jet momenta
will significantly change the flow field in the impingement region. Moreover, since the 1-D formulation
assumptions are best satisfied in high Re [13], the compatibility of the 1-D formulation to experiments is
expected to be best at high Re.
The strain rate at the stagnation point, K ≡ uz/ z, (schematically shown in Fig. 3.1b) was also
compared to the 1-D formulation predictions as it is one of the most important parameters in characterizing
aerodynamic effects in counterflow flame experiments (e.g., [1,12,38]). The results are presented in Fig.
3.10 as the magnitude of the relative difference between 1-D and 2-D results and only for cases for which
the discrepancy is higher than 1%. Note that the K values for cases with geometries (a) and (b) can be
accurately predicted by the 1-D formulation and they are not shown in Fig. 3.10. For CFP0.5S0, the
discrepancy for Re > 1200 is below the measurement uncertainty. For all other cases shown in Fig. 3.10,
discrepancies up to 34% were identified although mixed BCs were implemented.

91

Fig. 3.10. Relative difference between strain rate at the stagnation point computed by 1-D and 2-D
approaches; dashed area: region that cannot be resolved experimentally due to measurement uncertainty.

The results of Figs. 3.9 and 3.10 suggest that 1-D formulation with mixed BCs can still provide
reasonable predictions in comparison with 2-D results if the nozzle design assures uniform exit velocities
at the nozzles exits, not only in the single jet configuration but also in the presence of the opposing jet.
Therefore, the stagnation pressure feedback into the nozzles should be eliminated for better agreements.
The agreement is still valid even when bifurcation occurs. However, bifurcation can further complicate
the data interpretation as shown in Chapter 2.
4.2.1.1 Confinement effects
As the confinement level decreases by increasing H, the effect of bifurcation on flow asymmetry and
the stagnation point offset decreases. This can be understood by comparing the bifurcation maps for
CWU1S0 (Fig. 3.9a) and SWU1S0 (Fig. 3.9b) for uniform inlet velocity profiles. The same trend was
observed for parabolic inlet velocity profiles (CFP0.5S0 and CFP1S0 in Fig. 3.9c). For unconfined cases
in symmetric geometries, the stagnation point rests at the middle plane for all Re considered in this study.
As the confinement level decreases the effect of recirculation bubbles on the impingement region
decreases. This is because decreasing confinement leads to increase in ur/ r and reduced velocity
magnitudes within the recirculation bubbles, thus weakening their effect on the impingement region.
0
5
10
15
20
25
30
35
0 500 1000 1500 2000
|K
2-D
- K
1-D
| / K
2-D
, Uncertainty [%]
Re
Experimental Uncertainty
UFU1S0
UFU1A0
CFP0.5S0
CFP1S0
SFP1S0
UFP1S0
UFP1A0

92
Comparing the results for CFP1S0 and SFP1S0 in Fig. 3.10 reveals that decreasing confinement also helps
improve the 1-D predictions of K.

3.4.2.1.2. Aspect ratio effects
As shown in Table 3.1, by increasing , Recrit decreases for both uniform velocity profiles, e.g.,
CWU0.5S0 and CWU1S0, and parabolic inlet velocity profiles, e.g., CFP0.5S0 and CFP1S0. This is
consistent with the findings of Pawlowski et al. [19] who showed that for uniform velocity profiles and
1 ≤  ≤ 2, Recrit = (60/ ) + 30. The validity of this correlation for  < 1 was investigated via CWU0.5S0
and it was found that it under predicts Recrit only by 2%. The flow asymmetry and the stagnation point
offset were found to increase with . as shown in Fig. 3.9 (comparing CWU0.5S0 and CWU1S0 or
CFP0.5S0 and CFP1S0), also consistent with the findings of Pawlowski et al. [19]. Similarly, the 1-D
predictions of K were found to improve by decreasing  as apparent by comparing curves for CFP0.5S0
and CFP1S0 in Fig. 3.10.

3.4.2.1.3. Inlet velocity profile effects
Ciani et al. [20] reported that replacing the uniform inlet velocity profiles with parabolic ones for the
same geometry reduces Recrit by almost half since the maximum velocity for parabolic profiles are twice
the average velocity. Such reduction is apparent by comparing Recrit for CWU0.5S0 and CFP0.5S0 or
comparing CWU1S0 and CFP1S0 in Table 3.1. It should be noted that Recrit found for CFP1S0 in [20] is
4% lower than corresponding value shown in Table 3.1.
The stagnation point offset was found to be significantly higher for cases with parabolic inlet velocity
profile, e.g., comparing CWU0.5S0 and CFP0.5S0 or CWU1S0 and CFP1S0 in Fig. 3.9. This was mainly
due to the fact that these parabolic inlet profiles exist only when no flow conditioning screens are present
at the nozzle exits, and when the flow field inside the nozzle tubes is accounted for in the modeling. In
this case, after the symmetric flow field becomes unstable at the bifurcation point, it goes through a
transient state before it settles at the final asymmetric bifurcated state. During this transition, the flow field
is free to adjust itself to the most stable state by modifying the nozzle exit velocity profiles. On the other
hand, by fixing the velocity at the nozzle exits by placing screens at the nozzle exits, the BCs are essentially
fixed at that location. Thus, the stagnation point for the bifurcated state cannot be very close to the nozzle
exits. Otherwise, the stagnation pressure would modify the nozzle exit velocities due to continuity which
violates the BCs for the current geometry and given flow conditions.
Nozzle designs suggested by Niemann et al. [13] and investigated in other studies [15,16] include flow
conditioner screens at the nozzle exit to eliminate the pressure feedback effect. Moreover, Fig. 3.10 shows

93
that the discrepancies for UFU1S0 are smaller compared to CFP1S0. Hence, the more parabolic the
velocity profiles at nozzle exits, the higher the discrepancies in K as obtained from the 1-D and 2-D
formulations.
In order to further investigate the degree of asymmetry and two-dimensionality after the bifurcation,
the location of the stagnation plane represented by the offset at various radii is shown in Fig. 3.11 for
different BCs. For confined cases, the stagnation plane is flat and located at the middle plane for
Re < Recrit. After the bifurcation, the stagnation plane becomes concave towards the lower nozzle for the
lower branch and the curvature of the stagnation plane initially increases with Re. Further increasing Re
reduces the curvature and forces the stagnation plane back to its original position before the bifurcation.
The curvature of the stagnation plane and the fact that curvature is opposite for the two solution branches
suggest that the velocities and other scalar distributions near the stagnation point are not 1-D and that the
solution branches have different properties, as shown in Chapter 2.
Furthermore, the shape of the stagnation plane depends on the BCs. While the stagnation plane has a
parabolic shape for uniform inlet velocity profiles shown in Fig. 3.11a, the curvature was only significant
for 𝑟 > 0.3 for parabolic inlet velocity profiles shown in Fig. 3.11b. These differences suggest also that
different inlet velocity profiles result in different flow fields in the impingement region as also reported
in previous studies [9,10,53].

94

Fig. 3.11. Stagnation plane locus in terms of offset from the middle plane between the nozzles at various
Re for different boundary conditions in the momentum balanced condition; (a) lower branch of CWU1S0;
(b) lower branch of CFP1S0; (c) UFU1A0.

3.4.2.1.4. Geometrical asymmetry and imperfection effects
As shown in Fig. 3.9a and Table 3.1, having a confinement imperfection (w) as small as 1.25% of D,
not only changes the bifurcation type from pitchfork to imperfect but also increases Recrit by 18%.
Therefore, even small imperfections in the surrounding geometry of the jets, i.e., w, result in disappearance
of the symmetric flow field irrespective of Re. Moreover, the flow field has more tendency towards the
branch on which the imperfection is imposed such that the other solution branch is realizable only for
large Re. This can be seen by comparing the bifurcation maps of CWU1A1.25 and CWU1A12.5 to
CWU1S0 in Fig. 3.9a. The stagnation point offset and flow asymmetry were also found to be higher in
the presence of confinement imperfection. Increasing w results in increased Recrit, stagnation point offset,
and flow asymmetry.
Asymmetry in the nozzle rim assembly (δ) was also determined to influence the symmetry of the flow
field. While no bifurcation was found for unconfined cases, the symmetry of the flow field is compromised
and the stagnation point shifts towards the nozzle with the larger δ as shown in Fig. 3.9 for UWU1A0,
-25
-20
-15
-10
-5
0
0.0 0.1 0.2 0.3 0.4 0.5
-25
-20
-15
-10
-5
0
0.0 0.1 0.2 0.3 0.4 0.5
-25
-20
-15
-10
-5
0
0.0 0.1 0.2 0.3 0.4 0.5
Offset [%]
Re = 53
Re = 57
Re = 157
Re = 348
Re = 1101
Re = 125
Re = 250
Re = 375
Re = 500
Re = 89
Re = 91
Re = 155
Re = 275
Re = 1425
r r r
(a) (b) (c)

95
UFU1A0, and UFP1A0. The sensitivity of the flow field to δ can partly explain the experimental
observations in Figs. 3.5 and 3.7. Particularly, the results of Fig. 3.9b for UWU1A0 reveal that the
stagnation point offset decreases with Re similar to the experimental results shown in Fig. 3.5. However,
the measured shifts were significantly higher than those caused by the asymmetry in δ alone. The
simulations could replicate the experimental results only when the asymmetry in the flow conditioning
screens was also included in UFU1A0. Therefore, the flow field was found to be more sensitive to the
symmetry in the screens locations than the symmetry in δ. The stagnation plane shape for UFU1A0 is also
shown in Fig. 3.11c. The results indicate that increasing Re decreases the stagnation plane curvature and
forces it towards the mid-plane position in agreement with the result of Fig. 3.5.
Moreover, for non-zero δ, the stagnation point is located at different positions for different BCs at the
same Re. By comparing the results for UWU1A0, UFU1A0, and UFP1A0 in Fig. 3.9, it is apparent that
having a uniform inlet velocity profile without allowing the feedback inside the nozzle, i.e., UWU1A0,
had the smallest offset.
It should be noted that no stagnation point offset was identified for the unconfined cases with
symmetrical geometries, i.e., UWU1S0, UFU1S0, and UFP1S0. The fact that flow fields that are otherwise
insensitive to Re and BCs exhibit notable sensitivity upon the introduction of a small asymmetry in the
nozzle assembly, indicates that care is required in manufacturing and designing the nozzles, especially
when experiments are to be performed at low Re (e.g., [4,56]).

3.4.2.2. Sensitivity to momentum imbalance
Figure 3.12 depicts the stagnation point offset versus U1/U2 for selected unconfined cases in order to
investigate the dependence on Re and the validity of the 1-D formulation in the absence of bifurcation.
UWU1S0 exhibits the smallest sensitivity to momentum imbalance and was best predicted by the 1-D
formulation with plug flow BCs. Therefore, in the case of uniform nozzle exit velocity profile and no
stagnation pressure feedback into the nozzles, only the measured U1 and U2 are needed in order to assure
compatibility with the 1-D formulation. This was not the case in the experimental results shown in Figs.
3.7 and 3.12b. Letting the pressure field influence the internal flow inside the nozzles (UFU1S0, Fig.
3.12b) resulted in more deviation from the 1-D predictions, i.e. higher sensitivity of the stagnation point
offset to momentum imbalance. However, this sensitivity was found to decrease with increasing Re.
Similar behavior was observed for the cases with asymmetric geometry as shown for UFU1A0 in Fig.
3.12b with the difference that asymmetries in the nozzle assembly cause the 1-D predictions to deviate
further from the 2-D simulation results. On the other hand, by implementing mixed BCs, the measured
offset is recovered by the 1-D formulation as shown in Fig. 3.12b. Hence, for nozzle designs in which

96
stagnation pressure effects are not negligible, mixed BCs are required to predict the stagnation point
position by 1-D simulations. However, even mixed BCs cannot assure success in predicting K as the results
revealed that deviations up to 34% can still exist as demonstrated in Fig. 3.10.
Figure 3.12b also reveals that the experimental results match their numerical counterparts i.e.,
UFU1A0. The highest sensitivity of the stagnation point position to momentum imbalance was found for
parabolic velocity profiles, as shown in Fig. 3.12c for UFP1S0, such that 10% miss-match caused
stagnation point shifts up to 28%. Additionally, increasing Re causes more asymmetry and stagnation
point offset for this case. In the worst-case scenario considered, UFP1S0 at Re = 500 and U1/U2 = 1.1,
even when utilizing mixed BCs the stagnation point location and K predicted by the 1-D formulation are
9% and 38% different respectively compared to the 2-D simulation results. Therefore, parabolic inlet
velocity profiles introduce the most deviation from one-dimensionality as also reported in previous studies
[9,10,16,53]. The results of Figs. 3.9 and 3.12 for geometry (d) indicate that the stagnation point offset is
more sensitive to U1/U2 than the asymmetry in δ.

97

Fig. 3.12. Stagnation point offset from the middle plane between the nozzles for various unconfined cases
as a function of the average exit velocity ratio at various Re; (a) UWU1S0; (b) closed symbols: UFU1S0,
open symbols with dotted lines: UFU1A0; (c) UFP1S0; momentum balanced condition is marked by a
vertical dashed line and the middle plane by a horizontal dashed line.

The evidence presented suggests that the stagnation point offset observed in the unconfined isothermal
opposed-jet experiments outside the free-floating regime [18,21-23,32,35] is due to the high sensitivity to
nozzle assembly asymmetries and momentum miss-match rather than the bifurcation suggested by
Pawlowski et al. [19].

U
1
/U
2
(a) (b)
U
1
/U
2
U
1
/U
2
(c)
Offset [%]
1-D (plug flow)
Re = 50
Re = 250
Re = 500
1-D (mixed BCs from Re = 50) PIV data , Re = 600
,
,
,
1-D (mixed BCs from Re = 600)
1-D (plug flow)
Re = 125
Re = 250
Re = 500
1-D (plug flow)
Re = 50
Re = 250
Re = 500
-30
-25
-20
-15
-10
-5
0
5
10
15
20
25
30
0.90 0.95 1.00 1.05 1.10
-30
-25
-20
-15
-10
-5
0
5
10
15
20
25
30
0.90 0.95 1.00 1.05 1.10
-30
-25
-20
-15
-10
-5
0
5
10
15
20
25
30
0.90 0.95 1.00 1.05 1.10

98
3.5. Concluding remarks
A counterflow configuration consisting of two laminar opposed-flow air jets at ambient temperature
and pressure was investigated experimentally and computationally in search of instabilities and
asymmetries. These non-idealities can cause complications especially when performing flame
experiments and comparing data with results from one-dimensional simulations. The effects of different
confinement levels, geometry, nozzle separation distance to nozzle diameter ratio, inlet velocity profiles,
and momentum imbalance were quantified over a wide range of Reynolds numbers. The results were
presented in terms of the deviation between experimental and/or two-dimensional simulation results for
stagnation point position, stagnation plane shape, and strain rates from those obtained from one-
dimensional calculations.
Bifurcation analysis revealed that the flow field between opposed jets having equal momenta
bifurcates after a critical Reynolds number, but only when the radial flow generated outside the jets
impingement region is confined between two disk walls. The bifurcation results in two co-existing
asymmetric stable states which display non-monotonic behavior with Reynolds number; the deviation of
one-dimensional predictions from their two-dimensional counterparts initially increases by increasing the
Reynolds number but decreases by further increasing the Reynolds number. It was also determined that
the deviations decrease by decreasing the confinement level and nozzle separation distance to nozzle
diameter ratio. The degree of asymmetry after the bifurcation was found to be different for uniform and
parabolic inlet velocity profiles and it increased with increasing confinement level and nozzle separation
distance to nozzle diameter ratio. Introducing a small asymmetry in the distance between the confining
walls and the nozzle exits changed the bifurcation type from pitchfork to imperfect and increased the
critical Reynolds number and flow asymmetry. To minimize these non-idealities in laminar counterflow
flame experiments, surrounding surfaces and confining walls should be removed from the vicinity of the
jets exits.
Experimentally, only unconfined flows with nozzle separation distance to nozzle diameter ratio of one
was explored. It was observed that even in the absence of surfaces and confining walls next to the jet exits,
and for momentum balanced condition, the flow field is asymmetric in low Reynolds numbers of the order
of 300 such that the stagnation plane is not located at the middle plane between the nozzles. This is caused
by the asymmetric nozzle exit velocity profiles, which are due to unavoidable geometrical asymmetries
and manufacturing imperfections in the nozzle assembly. The flow asymmetry was found to decrease by
increasing the Reynolds number and eventually cannot be detected. The observed behavior was replicated
by introducing an asymmetric location of the flow conditioning screens inside the nozzles in the numerical
simulation where the stagnation plane was found to shift towards the nozzle that had a screen closer to the

99
nozzle exit, i.e. the nozzle with a more uniform exit velocity profile. It was determined that the deviation
of one-dimensional predictions of the stagnation point location from their two-dimensional counterparts
decreases monotonically by increasing Reynolds number in unconfined geometries. This deviation for the
strain rate however was found to have a non-monotonic behavior with Reynolds number. The deviations
in stagnation point location and strain rate were found to be lowest for uniform nozzle exit velocity profiles
that are immune to the feedback of the stagnation pressure field. Moreover, it was confirmed that using
mixed BCs in one-dimensional simulations decreases the discrepancy in the stagnation point location that
is within the experimental uncertainty for unconfined geometries. The discrepancy in strain rate at the
stagnation point however was not resolved for some conditions by implementing the mixed BCs. The
persistence of the strain rate discrepancy underlines the importance of the nozzle design when
compatibility with the assumptions of one-dimensional formulation is desirable.
Among all the cases considered, parabolic inlet velocity profiles resulted in the largest flow asymmetry
and two-dimensionality whether due to bifurcation or asymmetries in the boundary conditions. While
bifurcation appeared to be the main reason for flow field asymmetry in confined cases, it was determined
that imperfections, momentum mismatch, and geometrical asymmetries mainly influence the unconfined
cases. Furthermore, results presented here suggest that except for parabolic jet exit velocity profiles, as
Reynolds number increases, there is a decrease in flow field asymmetry, sensitivity to geometrical
imperfections, and to errors in momentum balance.
In order to minimize the sensitivity of the flow field to experimental imperfections two approaches
were identified: (1) small nozzle separation distance to nozzle diameter ratios, and (2) nozzle assemblies
with flow conditioning screens that producing uniform inlet velocity profiles which minimize the feedback
of the stagnation pressure field.
It was found both experimentally and numerically that the flow field is more sensitive to geometrical
imperfections and asymmetries at low Reynolds numbers than at high Reynolds numbers. Therefore,
experimental observables from laminar flame experiments can be compared to the predictions of one-
dimensional formulation with more confidence at Reynolds number values of approximately 1000.

3.6. Acknowledgements
The Air Force Office of Scientific Research (Grant No. FA9550-15-1-0409) supported this work under
the technical supervision of Dr. Chiping Li. The numerical simulations were carried out using the Extreme
Science and Engineering Discovery Environment (XSEDE), which is supported by National Science
Foundation grant number TG-CTS140012. The authors would like to gratefully acknowledge the
assistance of Dr. Vyaas Gururajan for his feedback and technical support with the numerical simulations.

100
3.7. References
[1] J. Carpioa, A. Liñán, A.L. Sánchez, F.A. Williams, Aerodynamics of Axisymmetric Counterflowing
Jets, Combust. Flame 177 (2017) 137-143.
[2] A. Tamir, Impinging-stream reactors: Fundamentals and applications. Elsevier, Beer-Sheva 2014.
[3] O. Park, P.S. Veloo, H. Burbano, F.N. Egolfopoulos, Studies of premixed and non-premixed
hydrogen flames, Combust. Flame 162 (2015) 1078-1094.
[4] S.H. Won, B. Jiang, P. Diévart, C.H. Sohn, Y. Ju, Self-sustaining n-heptane cool diffusion flames
activated by ozone, Proc. Combust. Inst. 35 (2015) 881-888.
[5] R.R Burrell, Studies of Methane counterflow flames at low pressures, Ph.D. thesis, Mechanical
Engineering, University of Southern California, Los Angeles, California 2017.
[6] F. Carbone, A. Gomez, The structure of toluene-doped counterflow gaseous diffusion
flames, Combust. Flame 159 (2012 )3040-3055.
[7] R.J. Kee, J.A. Miller, G.H. Evans, G. Dixon-Lewis, A computational model of the structure and
extinction of strained, opposed flow, premixed methane-air flames, Symp. (Int.) Combust. 22
(1988)1479–1494.
[8] H.K. Chelliah, C.K. Law, T. Ueda, M.D. Smooke, F.A. Williams, An experimental and theoretical
investigation of the dilution, pressure and flow-field effects on the extinction condition of methane-
air-nitrogen diffusion flames, Symp. (Int.) Combust. 23 (1990) 503-511.
[9] Y.M. Kim, H-J. Kim, Multidimensional effects on structure and extinction process of counterflow
nonpremixed hydrogen-air flames, Combust. Sci. Tech. 137(1998) 51-80.
[10] C.E. Frouzakis, J. Lee, A.G. Tomboulides, K. Boulouchos, Two-dimensional direct numerical
simulation of opposed-jet hydrogen-air diffusion flame, Symp. (Int.) Combust. 27 (1998) 571–577.
[11] V. Mittal, H. Pitsch, F.N. Egolfopoulos, Assessment of counterflow to measure laminar burning
velocities using direct numerical simulation, Combust. Theory and Model. 16 (2012) 419-433.
[12] B.G. Sarnacki, G. Esposito, R.H. Krauss, H.K. Chelliah, Extinction limits and associated
uncertainties of nonpremixed counterflow flames of methane, ethylene, propylene and n-butane in
air, Combust. Flame 159 (2012) 1026-1043.
[13] U. Niemann, K. Seshadri, F.A. Williams, Accuracies of laminar counterflow flame experiments,
Combust. Flame 162 (2015) 1540-1549.
[14] R.F. Johnson, A.C. VanDine, G.L. Esposito, H.K. Chelliah, On the Axisymmetric Counterflow
Flame Simulations: Is There an Optimal Nozzle Diameter and Separation Distance to Apply Quasi
One-Dimensional Theory?, Combust. Sci. Tech. 187 (2015) 37-59.

101
[15] R.R. Burrell, R. Zhao, D.J. Lee, H. Burbano, F.N. Egolfopoulos, Two-dimensional effects in
counterflow methane flames, Proc. Combust. Inst. 36 (2016) 1387-1394.
[16] A. Ansari, F.N. Egolfopoulos, Flame ignition in the counterflow configuration: Reassessing the
experimental assumptions, Combust. Flame 174 (2016) 37-49.
[17] J.L. Fox, G.W. Morgan, On the stability of some flows of an ideal fluid with free surfaces, Brown
Univ. Providence RI DIV of Applied Mathematics (1952), paper TR-2.
[18] J.C. Rolon, D. Veynante, J.P. Martin, Counter jet stagnation flows, Exp. Fluids 11 (1991) 313-324.
[19] R.P. Pawlowski, A.G. Salinger, J.N. Shadid, T.J. Mountziaris, Bifurcation and stability analysis of
laminar isothermal counterflowing jets, J. Fluid Mech. 551 (2006) 117-139.
[20] A. Ciani, W. Kreutner, C.E. Frouzakis, K. Lust, G. Coppola, K. Boulouchos, An experimental and
numerical study of the structure and stability of laminar opposed-jet flows, Comput. Fluids 39 (2010)
114-124.
[21] W.F. Li, Z. Sun, H.F. Liu, F.C. Wang, Z. Yu, Experimental and numerical study on stagnation point
offset of turbulent opposed jets, Chem. Eng. J. 138 (2008) 283-294.
[22] W.F. Li, T.L. Yao, F.C. Wang, Study on factors influencing stagnation point offset of turbulent
opposed jets, AIChE J. 56 (2010) 2513-2522.
[23] W.F. Li, T.L. Yao, H.F. Liu, F.C. Wang, Experimental investigation of flow regimes of
axisymmetric and planar opposed jets, AIChE J. 57 (2011) 1434-1445.
[24] C.G. Fotache, T.G. Kreutz, D.L. Zhu, C.K. Law, An Experimental Study of Ignition in Nonpremixed
Counterflowing Hydrogen versus Heated Air, Combust. Sci. Technol. 109 (1995) 373–393.
[25] R. Seiser, K. Seshadri, E. Piskernik, A. Linán, Ignition in the viscous layer between counterflowing
streams: Asymptotic theory with comparison to experiments, Combust. Flame 122 (2000) 339-349.
[26] D.A. Hammond, L.G. Redekopp, Local and global instability properties of separation bubbles, Eur.
J. Mech.,B/Fluids 17 (1998) 145-164.
[27] S. Cherubini, J. Robinet, P. De Palma, The effects of non-normality and nonlinearity of the Navier–
Stokes operator on the dynamics of a large laminar separation bubble, Phys. Fluids 22 (2010)
014102.
[28] R.M. Feran, T. Mullin, K.A. Cliffe, Nonlinear flow phenomena in a symmetric sudden expansion,
J. Fluid Mech. 221 (1990) 596-608.
[29] G. Amantini, J.H. Frank, A. Gomez, Experiments on standing and traveling edge flames around
flame holes, Proc. Combust. Inst. 30 (2005) 313–21.

102
[30] G. Amantini, J.H. Frank, M.D. Smooke, A. Gomez, Computational and experimental study of steady
axisymmetric non-premixed methane counterflow flames, Combust. Theory and Model. 11 (2007)
47-72.
[31] S.T. Strogatz, Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry, and
engineering, Westview Press, 2000, pp. 44-93.
[32] L.W. Kostiuk, K.N. Bray, R.K. Cheng, Experimental study of premixed turbulent combustion in
opposed streams. Part I—Nonreacting flow field, Combust. Flame 92 (1993) 377-395.
[33] G. Pellett, K. Isaac, W. Humphreys, L. Garttrell, W. Roberts, C. Dancey, G. Northam, Velocity and
thermal structure, and strain-induced extinction of 14 to 100% hydrogen-air counterflow diffusion
flames, Combust. Flame 112 (1998) 575–592.
[34] G. Stan, D.A. Johnson, Experimental and numerical analysis of turbulent opposed impinging jets,
AIAA J. 39 (2001) 1901-1908.
[35] L.W. Kostiuk, K.N. Bray, R.K. Cheng, Experimental study of premixed turbulent combustion in
opposed streams. Part II—reacting flow field and extinction. Combust. Flame 92 (1993) 396-409.
[36] N. Ogawa, H. Maki, Studies on opposed turbulent jets: influences of a body on the axis of opposed
turbulent jets, Bulletin of JSME 29 (1986) 2872-2877.
[37] N. Ogawa, H. Maki, K. Hijikata, Studies on opposed turbulent jets: impact position and turbulent
component in jet center, JSME Int. J. 35 (1992) 205-211.
[38] G. Scribano, F. Bisetti, Reynolds number and geometry effects in laminar axisymmetric isothermal
counterflows, Phys. Fluids 28 123605 (2016).
[39] R.P. Lindstedt, D.S. Luff, J. H. Whitelaw, Velocity and strain-rate characteristics of opposed
isothermal flows, Flow Turbul. Combust. 74 (2005) 169-194.
[40] Y. Dong, C.M. Vagelopoulos, G.R. Spedding, F.N. Egolfopoulos, Measurement of laminar flame
speeds through digital particle image velocimetry: mixtures of methane and ethane with hydrogen,
oxygen, nitrogen, and helium, Proc. Combust. Inst. 29 (2002) 1419–1426.
[41] I. Hughes, T. Hase, Measurements and their uncertainties: a practical guide to modern error analysis,
Oxford University Press, New York 2010.
[42] A.E. Lutz, R.J. Kee, J.F. Grcar, F.M. Rupley, OPPDIF: A FORTRAN program for computing
opposed-flow diffusion flames, Snadia Report, SAND96-8243, Sandia National Laboratories, 1997.
[43] H. H. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum
mechanics using object-oriented techniques, Comput. Phys. 12 (1998) 620-631.
[44] K.K. Chen, C.W. Rowley, H.A. Stone, Vortex dynamics in a pipe T-junction: Recirculation and
sensitivity, Phys. Fluids 27 (2015) 034107.

103
[45] K.K. Chen, C.W. Rowley, H.A. Stone, Vortex breakdown, linear global instability and sensitivity
of pipe bifurcation flows, J. Fluid Mech. 815 (2017) 257-294.
[46] R.I. Issa, Solution of the implicitly discretized fluid flow equations by operator-splitting, J. Comput.
Phys. 62 (1986) 40-65.
[47] R. Courant, K. Friedrichs, H. Lewy, On the partial difference equations of mathematical physics,
IBM J. 11 (1967) 215-234.
[48] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pac. J.
Math. 16 (1966) 1–3.
[49] C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations. Society for Industrial and
Applied Mathematics, Philadelphia 1995.
[50] Y. Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving
nonsymmetric linear systems, SIAM J. Sci. Stat. Computing 7 (1986) 856–869.
[51] L. Tuckerman, D. Barkley, Bifurcation analysis for time steppers, In: E. Doedel, L. S. Tuckerman,
Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems, Springer-
Verlag, New York, 2000, pp. 453–466.
[52] L.N. Trefethen, D. Bau III, Numerical Linear Algebra, Siam, Philadelphia, 1997.
[53] E. Korusoy, J.H. Whitelaw, Inviscid, laminar and turbulent opposed flows, Int. J. Numer. Methods
Fluids 46 (2004) 1069-1098.
[54] P.M. Gresho, R.L. Sani, On pressure boundary conditions for the incompressible Navier-Stokes
equations, Int. J. Numer. Meth. Fluids 7 (1987) 1111-1145.
[55] C. Godrèche, P. Manneville, B. Castaing, Hydrodynamics and nonlinear instabilities, Cambridge
University Press, Cambridge, 2005.
[56] C.B. Oh, A. Hamins, M. Bundy, J. Park, The two-dimensional structure of low strain rate
counterflow nonpremixed-methane flames in normal and microgravity, Combust. Theory and
Model. 12 (2008) 283-302.
[57] C.A.J. Fletcher, Computational techniques for fluid dynamics 1: Fundamental and general
techniques, Berlin and New York, Springer-Verlag, 1988.

104
Chapter 4
Parameters influencing the burning rate of laminar flames propagating
into a reacting mixture

4.1. Introduction
The laminar flame speed ( ) is a fundamental property of a combustible mixture that is a measure
of its exothermicity and reactivity in a given diffusive medium. This quantity is widely used as a measure
of the burning rate in practical devices (e.g., [1,2]). The physics of flame propagation and the parameters
controlling has been thoroughly investigated in the typical range of 1 atm < P < 10 atm and
300 K < Tu < 500 K (e.g., [2]); P and Tu denote the thermodynamic pressure and the unburned mixture
temperature respectively. At those thermodynamic conditions, the mixture auto-ignition delay time (𝜏 )
is considerably larger than the characteristic residence time of the unburned mixture before the flame
arrival (𝜏 ) [3,4]. As a result, the unburned mixture is non-reactive for all practical purposes. In order to
quantify the extent of reactivity in the unburned mixture, a progress variable (C) is defined as
C = 1 - YF/(YF)o; YF is the instantaneous fuel mass fraction and (YF)o the initial supply value. Thus, C
varies in the range {0,1} with C = 0 representing no fuel consumption and C = 1 representing full fuel
consumption by auto-ignition.
In practical devices such as internal combustion (IC) engines, elevated Tu and P values can be
encountered resulting potentially in C > 0 (e.g., [5,6]). Under such conditions, 𝜏 >> 𝜏 may not be
applicable, and the physics of flame propagation and the parameters influencing the burning rate of the
reacting front may differ notably compared to non-reactive conditions (C = 0).
The topic of flame propagation into a reacting mixture has recently gained more attention due its
relevance to knock in IC engines (e.g., [5,6]) and the control of homogenous charge compression ignition
(HCCI) engines (e.g., [7,8]). These studies can be broadly classified into two categories.
The first category involves studies of flames that are statistically steady with respect to an inertial
reference frame. Experimentally, this has been achieved by first mixing heated fuel and oxidizer and then
providing sufficient f so that C > 0 before the “reactants” enter the flame that is stabilized at a fixed
location (e.g., [9-11]). These steady flames have recently been investigated numerically to reveal the
flame structure for a wide range of C values [4].
The second category includes studies of unsteady flame propagation into a mixture with C > 0 both
computationally (e.g., [7,8,12-20]) and experimentally (e.g., [21]). More specifically, the experiments

S
u
o

S
u
o

105
involved flames propagating in a constant volume chamber. During the later stages of flame propagation,
the unburned mixture is compressed to the point that C starts increasing leading potentially to auto-ignition
and knock (e.g., [21]).
Previous studies have revealed that the burning rate increases as C increases [3,4,8,12,22]. It is well
known that the flame is controlled by the balance of heat conduction and chemical heat release at C = 0
[1]. Hence the flame is called a “conduction-reaction” front. However, as C increases, the effect of
transport phenomena decreases to a point that the convection-reaction balance is controlling the flame at
C ≈ 1 [3,4,8,23]. Different regimes and burning rates have been observed experimentally depending on
the fuel type and the C values (e.g., [9-11])
Steady computations have demonstrated significant sensitivity of the burning rate to low temperature
chemistry (LTC) [17]. However, reaction rate perturbations can affect both the mass burning rate of the
reacting front (𝑚 ̇ ) and the ignition dynamics of the unburned mixture upstream of the flame.
Understanding the parameters influencing 𝑚 ̇ when C > 0 is essential in understanding the combustion
processes in next-generation engines that are expected to operate at higher pressures. Previous studies
under unsteady conditions mainly focused either on flame-pressure wave interactions (e.g., [13-18,24]) or
cool flames (e.g., [7,19,20]). Studying cool flames propagation is essential in validating LTC [17,19,20].
It should be noted though that cool flames could be established only under a narrow window of
thermodynamic conditions and time scales that may or may not be realized in engines [19].
On the other hand, it is worthwhile to investigate the physics of the transient “hot” flame in the
conduction-reaction controlled regime when C > 0. Specifically, it is not obvious what is causing 𝑚 ̇ to
increase right after the onset of the unburned mixture reactivity and before transitioning to the convection-
reaction controlled regime. Although the Tu increase seems to be the major cause of the 𝑚 ̇ increase
[3,8,22], it is not clear whether it is due to the subsequent enhancements of reactivity and/or transport
rates. These effects on 𝑚 ̇ are generally coupled [1] and it is essential to quantify their relative importance.
Moreover, composition changes in the unburned mixture are also coupled with the Tu increase but the
attendant effects on 𝑚 ̇ have yet to be quantified.
This chapter aims to further explore computationally the parameters responsible for the increase of 𝑚 ̇
when an unsteady flame propagates into a reacting mixture before transitioning to the convection-reaction
controlled regime and perform sensitivity analysis to identify the importance of each parameter.

106
4.2. Modeling approach
The transient one-dimensional reacting flow code (TORC) was used to integrate implicitly the
spatially discretized low Mach number conservation equations for mass, species and energy [25-27];
details can be found in Appendix A. A semi-infinite Cartesian domain was chosen to simulate an unsteady
flame propagation into a reacting mixture at constant pressure. The left side of the domain, x = 0, was
bounded by an adiabatic impermeable wall where a flame was initiated. An outflow boundary condition
was implemented at x = xmax. The appropriate xmax value was found by trial and error such that the
unburned mixture at x = xmax auto-ignites before the arrival of the flame. The diffusion velocities were set
to zero at the boundaries. Initially (t = 0), the composition in the domain was fuel/oxidizer/inert at (Tu)o
and constant P with C = 0. A hot pocket of equilibrium products was used to initiate a flame propagating
from left to right. Velocities were initially set to zero throughout the domain. An adaptive mesh
refinement algorithm was employed to provide sufficient grid resolution within the flame front (see
Appendix A).
The flame displacement speed with respect to the unburned mixture defined as Sd = (dxf/dt) 𝜌 /𝜌 , was
extracted at each time step; xf is the location of the peak heat release, while 𝜌 and 𝜌 represent the
densities of the unburned and burned mixture calculated at x = xmax and x = 0 respectively. As a result,
𝑚 ̇ = 𝜌 Sd. Note that for a reacting mixture, 𝜌 varies with time as both the composition and Tu change
with time.
Table 4.1 depicts the conditions investigated that are relevant to hydrogen IC engines [5] and HCCI
engines [3,7,8,14]. Stoichiometric H2/O2/N2 and n-C7H16/air mixtures were chosen as representative
mixtures. For H2/O2/N2 mixtures, H2:O2:N2 = 2:1:9 molar ratios were chosen to sufficiently dilute the
mixture and decrease the heat release rate in the flame therefore avoiding compressibility effects [13,23],
and the H2 kinetic model of Davis et al. [28] was utilized. For n-C7H16/air mixtures, the skeletal model
developed by Liu et al. [29] was incorporated that appropriately predicts the low and high temperature
chemistry. Both models have been tested successfully for the conditions investigated herein [14,18,27].
The simulations for H2/O2/N2 mixtures were performed for two values of (Tu)o to investigate the effect of
𝜏 on the results; 𝜏 was calculated using SENKIN based on the maximum rate of temperature increase
[30]. For the n-C7H16/air mixtures, the simulations were performed by varying (Tu)o in the vicinity of
negative temperature coefficients region such that the cases investigated have equal 𝜏 .

107
Table 4.1. Summary of mixtures and thermodynamic conditions investigated
Case Mixture
Pressure,
P (atm)
Initial temperature,
(Tu)o
(K)
Auto-ignition
delay time, 𝜏
(s)
Initial flame
displacement speed,
(Sd)o (cm/s)
H1 H2/O2/N2 55 900 6.08 × 10
-2
85
H2 H2/O2/N2 55 925 3.17 × 10
-2
101
C0 n-C7H16/air 1 298 ∞ 55
C1 n-C7H16/air 40 808 8.00 × 10
-4
110
C2 n-C7H16/air 40 900 8.00 × 10
-4
153
C3 n-C7H16/air 40 987 8.00 × 10
-4
207

4.3. Results and discussion
Results indicate that shortly after t = 0, the initial hot pocket of products generates a hot flame that
starts propagating into the unburned mixture. No cool flame was observed due to the high temperatures
used to initiate the flame [19]. It was determined that the flames display an initial quasi-steady propagation
where 𝑚 ̇ hardly changes until about t/𝜏 ≈ 0.6 and 0.2 for the H2 and n-C7H16 cases respectively.
Therefore, an initial constant (𝑚 ̇ )o was defined and used to scale the results as shown in Fig. 4.1, which
depicts also the Tu/(Tu)o as a function of the normalized time (t/𝜏 ). Table 1 also contains attendant (Sd)o
values for comparison purposes.

Figure 4.1. Evolution of the scaled burning rate and the scaled unburned mixture temperature as a
function of the normalized time.

After the initial quasi-steady propagation, 𝑚 ̇ starts to increase for the reacting cases. The state of the
unburned mixture was also found to vary as a function of time due to pre-flame reactivity as apparent
1.00
1.05
1.10
1.15
1.20
0 0.2 0.4 0.6 0.8 1
1.0
1.1
1.2
1.3
1.4
1.5
1.6
1.7
0 0.2 0.4 0.6 0.8 1
T
u
/(T
u
)
o
, Case H2
T
u
/(T
u
)
o
, Case H1
𝒎 ̇ /(𝒎 ̇ )
o
, Case H2
𝒎 ̇ /(𝒎 ̇ )
o
, Case H1
t/ 
ign
𝒎 ̇ /(𝒎 ̇ )
o
, Case C3
𝒎 ̇ /(𝒎 ̇ )
o
, Case C2
𝒎 ̇ /(𝒎 ̇ )
o
, Case C1
T
u
/(T
u
)
o
, Case C3
T
u
/(T
u
)
o
, Case C2
T
u
/(T
u
)
o
, Case C1

108
from the Tu increase in Fig. 4.1. Figures 4.2 and 4.3 illustrate the time evolution of C and the properly
scaled mass fractions of major species in the unburned mixture.

Figure 4.2. Evolution of the amount of fuel consumed, C, and the scaled water mass fraction in the
unburned mixture as a function of the normalized time for H 2/O2/N2 mixtures; equilibrium mass fraction
was used to scale the results for water.

Figure 4.3. Evolution of n-C7H16 (black) and formaldehyde (green) mass fractions in the unburned
mixture as a function of the normalized time for n-C7H16/air mixtures; (─) Case C1; (----) Case C2; (∙ ∙ ∙)
Case C3; initial and maximum values were used to scale n-C7H16 and formaldehyde mass fractions
respectively.

Note that the results were analyzed only in the conduction-reaction controlled regime. The threshold
for transition to the convection-reaction controlled regime was determined by tracking the instantaneous
volumetric heat conduction and heat release rate in the flame. As an illustration, the scaled volumetric
heat conduction rate (Qcond) and the scaled volumetric heat release rate (Qreact) for Case C2 are shown in
Fig. 4.4 at different times as functions of the temperature in the domain; both Qcond and Qreact were scaled
by the maximum value of Qreact. It is observed that maximum Qcond decreases with time indicating that
the flame starts transitioning to regime that is not controlled by conduction-reaction. The flame becomes
convection-reaction controlled at the time corresponding to the maximum heat release rate of the unburned
gas (tqmax) in agreement with the results of Martz et al. [8]; note that the flame is sustained by the
conduction-reaction balance up to tqmax that is approximately equal to 0.97𝜏 for the conditions of Fig.
4.4.
0
0.04
0.08
0.12
0.16
0.2
0
0.02
0.04
0.06
0.08
0.1
0 0.2 0.4 0.6 0.8 1
t/ 
ign
C, Case H2
C, Case H1
𝒀 𝐇 𝟐 𝐎 /(𝒀 𝐇 𝟐 𝐎 )
equil
, Case H2
𝒀 𝐇 𝟐 𝐎 /(𝒀 𝐇 𝟐 𝐎 )
equil
, Case H1
0.0
0.2
0.4
0.6
0.8
1.0
0 0.2 0.4 0.6 0.8 1
t/ 
ign
𝒀 𝐧 𝐂 𝟕 𝐇 𝟏𝟔 (𝒀 𝐧 𝐂 𝟕 𝐇 𝟏𝟔 )
o
𝒀 𝐂 𝐇 𝟐 𝐎 (𝒀 𝐂𝐇 𝟐 𝐎 )
Max

109
It should be noted also that C ≈ 0.1 at t = tqmax for H2 flames as shown in Fig. 4.2. Hence, the flame
transitions to the convection-reaction controlled regime before significant amount of fuel has been
consumed in the unburned mixture. For n-C7H16 flames on the other hand, C ≈ 0.9 at the threshold of
transition.

Figure 4.4. Scaled volumetric heat conduction rate (blue) and heat release rate (red) as a function of
temperature in the domain at various times for Case C2; (─) t/ ign = 0.13; (----) t/ ign = 0.60; (∙ ∙ ∙)
t/ ign = 0.97.

For H2 flames, 𝑚 ̇ exponential increases after the initial quasi-steady propagation. The value of (Tu)o
does not affect this trend as the results collapse after rescaling by the initial values as shown in Fig. 4.1;
note that (𝜏 )
≈ 2×(𝜏 )
. On the other hand, the results for n-C7H16 flames do not
collapse, suggesting the importance of LTC. The 𝑚 ̇ /𝑚 ̇ o curve follows the shape of Tu/(Tu)o that has a the
two-stage ignition behavior. In contrast to H2 flames, the variations of 𝑚 ̇ cannot be explained entirely by
Tu variations for n-C7H16 flames. For instance, for Case C3, the Tu increase is less than 1% until
t/𝜏 = 0.6 while 𝑚 ̇ increases by 12%. This initial increase is due to the variation in the unburned mixture
composition by the LTC pre-flame reactivity without significant attendant heat release. As demonstrated
in Fig. 4.3, 40% of the initial n-C7H16 is consumed at t/𝜏 = 0.6 for Case C3 such that a mixture of n-
C7H16 and smaller hydrocarbons (mostly formaldehyde and ethylene) enters the flame.
Note that the hump shown for Case C2 in Fig. 4.1 is also caused by the hump in the evolution of the
major species in the unburned mixture as shown in Fig. 4.3. As the reaction progresses in the unburned
mixture, the composition varies such that the rates of formaldehyde and ethylene production, along with
the n-C7H16 destruction rate for 0.3 < t/ 𝜏 < 0.5 are higher compared to 0.5 < t/ 𝜏 < 0.8. Since
formaldehyde and ethylene have higher than n-C7H16 [31], the rate of 𝑚 ̇ increase for 0.3 < t/𝜏 < 0.5
is higher compared to 0.5 < t/𝜏 < 0.8.
Hence, the 𝑚 ̇ evolution generally depends on the composition of the mixture entering the flame and
its temperature in a coupled manner, as noted also in previous studies [9-11]. However, when the results
-1.0
-0.5
0.0
0.5
1.0
900 1200 1500 1800 2100 2400 2700
Temperature (K)
Q
react
Q
cond

S
u
o

110
are plotted against C as shown in Fig. 4.5, for a given C, 𝑚 ̇ /𝑚 ̇ o is the largest for the case with the highest
Tu/(Tu)o.

Figure 4.5. Scaled flame burning rate and unburned mixture temperature as a function of the amount of
fuel consumed in the unburned mixture for n-C7H16/air mixtures.

In order to delineate the contributions of different parameters on 𝑚 ̇ , a sensitivity analysis was
performed by first identifying the important parameters. At a given time, 𝑚 ̇ can be expressed as:
𝑚 ̇ = 𝑚 ̇ (Tad, Yk(t), Tu(t), (t), Dk(t)) (4.1)
where Tad represents the adiabatic flame temperature, Yk the mass fraction of the species k in the unburned
state,  the thermal conductivity of the mixture, and Dk the diffusion coefficient of species k. As the flame
propagates into a reacting mixture, the pre-flame reactivity causes Tu and Yk to change with time. These
changes subsequently modify  and Dk and the combination of these effects lead to 𝑚 ̇ variations.
Although, for the constant pressure adiabatic cases considered here, Tad remains constant, the effects of
Tad on 𝑚 ̇ was quantified and compared to the other parameters.
In order to perform a meaningful sensitivity analysis, a decoupling approach was developed allowing
for the perturbation of one parameter at a time in Eq. (4.1). Computational experiments were then
performed to quantify the effects on 𝑚 ̇ . First, the parameters in Eq. (4.1) were decoupled for the baseline
cases shown in Table 4.1. Then, another TORC simulation was performed utilizing a slightly perturbed
value of each parameter and a new 𝑚 ̇ value was obtained. The results ware presented in terms of
logarithmic sensitivity coefficients ( 𝐿𝑆𝐶 ). For example, the effect of Tad on 𝑚 ̇ was quantified by
(𝐿𝑆𝐶 )
≈ (∆𝑚 ̇ /∆Tad) (Tad/𝑚 ̇ ), etc. On the other hand, the Yk effects were quantified through C by
computing (𝐿𝑆𝐶 )
, given the direct correlation between Yk and C.
It should be noted that the Tu and Yk perturbations relate primarily to “chemical effects” since they
affect the reaction rates. While Tu variations can only change the reaction rates for temperatures lower
than Tad, perturbing Tad only modifies the reaction rates in the high temperature range i.e., “high
1
1.1
1.2
1.3
1.4
1.5
0 0.2 0.4 0.6 0.8 1
𝒎 ̇ /(𝒎 ̇ )
o
, Case C3
𝒎 ̇ /(𝒎 ̇ )
o
, Case C2
𝒎 ̇ /(𝒎 ̇ )
o
, Case C1
T
u
/(T
u
)
o
, Case C3
T
u
/(T
u
)
o
, Case C2
T
u
/(T
u
)
o
, Case C1
Amount of fuel consumed, C

111
temperature chemical effect”.  and Dk perturbations assess “thermal effects” as by suppressing the other
effects utilizing the developed decoupling approach, the only effect of the temperature increase will be
the  and Dk modifications. Perturbing the pre-exponential factor of each reaction, Ai, standard sensitivity
analysis was also performed, and (𝐿𝑆𝐶 )
was compared to (𝐿𝑆𝐶 )
and (𝐿𝑆𝐶 )
. The effect of Dk is
insignificant in comparison to the other parameters (see Appendix B). Further details regarding the
perturbation procedure and 𝐿𝑆𝐶 computation can be found in Appendix B.
Figure 4.6 depicts some examples of the comparison between the baseline 𝑚 ̇ and those obtained by
perturbations. It was found that perturbing Tu, Yk, and Ai not only modifies 𝑚 ̇ but also the unburned
mixture evolution in time. For example, perturbing Ai of H + O2  O + OH (R1) at t = 0 results in
different C values compared to the baseline for t/𝜏 > 0.2 as apparent by the inset plot of Fig. 4.6 that
depicts the evolution of C.
Therefore, comparisons of 𝑚 ̇ values outside the quasi-steady regime are not justified since more than
one parameter is modified, and they are appropriate only when C ≈ 0 at t/𝜏 ≈ 0. However, quantifying
the perturbation effects for t/𝜏 > 0 is crucial in understanding the underlying physics. Hence, another
approach was developed to perform local perturbations at t/𝜏 > 0. First, the baseline simulation was
stopped at the time at which the perturbation was about to be introduced (tpert). Second, the parameters in
Eq. (4.1) were decoupled. Then, the perturbation was performed utilizing the unburned mixture state and
composition at t = tpert. Finally, a new simulation utilizing the perturbed parameters was conducted by
initiating a new flame using the hot pocket of equilibrium products at t = tpert. The results of the new
simulation were compared to the baseline by a simple time adjustment tnew = t + tpert. Figure 4.6 also
depicts the results obtained for non-zero tpert.
Stopping the baseline simulation and conducting a new one from t = 0 was essential in order to
properly decouple the parameters in Eq. (4.1) as elaborated in Appendix B. For example, since Tad
perturbations were achieved by modifying the N 2 diluent amount in the mixture, the perturbations could
not be simply imposed at t = tpert during the baseline simulation. Only by utilizing this procedure, local
isolated perturbations of the flame were achieved and the undesirable ignition dynamics modifications
discussed earlier were eliminated. Figure 4.6 confirms that the C values remain unchanged at
t = tpert = 0.59𝜏 by utilizing the local perturbation approach.

112

Figure 4.6. Comparison between the time evolution of the flame burning rate of the baseline and the
ones obtained by perturbing pre-exponential factors of selected reactions for Case C2; (R1):
H + O2  O + OH; (R2): C7H15O2  C7H14OOH; arrows indicate the time at which the perturbation
was imposed; inset: time evolution of the amount of fuel consumed in the unburned mixture.

Moreover, as apparent by the sudden initial jump in 𝑚 ̇ value for the perturbed cases in Fig. 4.6, the
effects of initial transient flame initiation are minimal irrespective of tpert. This was further confirmed by
comparing the baseline 𝑚 ̇ to its values obtained by restarting the baseline simulation from non-zero tpert
values without imposing any perturbations (see Appendix C). Hence, the comparison against the baseline
𝑚 ̇ could be performed with confidence at t ≈ tpert.
Comparison of the curves obtained by perturbing Ai of (R1) at tpert = 0 to the ones obtained by imposing
the perturbation at tpert = 0.59𝜏 in Fig. 4.6, reveals that the time history of 𝑚 ̇ and C depend also on tpert.
This observation is further justifying the local perturbation approach. Therefore, 𝐿𝑆𝐶 values were
computed only at t = tpert.
Figure 4.7 summarizes the sensitivity analysis results obtained at the initial quasi-steady propagation
(C = 0). The results for non-reactive Case C0 are shown also and compared to the reacting cases. At
room temperature and pressure, 𝑚 ̇ is not sensitive to Tu and the largest sensitivity was found to be, as
expected, to Tad.
For reacting cases, finite sensitivities to Tu and C are observed. However, (𝐿𝑆𝐶 )
remains the
largest among all 𝐿𝑆𝐶 values while (𝐿𝑆𝐶 )
is the lowest. For H2 flames, (𝐿𝑆𝐶 )
is negative due to the
considerable H2O production in the unburned mixture as demonstrated in Fig. 4.2. The dilution effect of
H2O has been also observed in H2/O2 flames [22]. However, positive (𝐿𝑆𝐶 )
values are observed for n-
C7H16 flames due to production of smaller hydrocarbons in the unburned mixture as shown in Fig. 4.3;
1.8
1.9
2.0
2.1
2.2
2.3
2.4
2.5
0.0 0.2 0.4 0.6 0.8 1.0
𝒎 ̇ (g/(cm
2
×s))
t/ 
ign
baseline
(R1), t
pert
/ 
ign
= 0
(R1), t
pert
/ 
ign
= 0.59
(R2), t
pert
/ 
ign
= 0
0.50
0.52
0.54
0.58 0.60 0.62
(R2), t
pert
/ 
ign
= 0.59
t/ 
ign
C

113
keeping all other parameters the same, the generated smaller hydrocarbons (mostly formaldehyde and
ethylene) have higher 𝑚 ̇ compared to n-C7H16 due to different chemistry involved [1,31]. Moreover, since
the generated H2O and CO2 are negligible in the early stages of unburned mixture reactivity, the dilution
effects are not present for n-C7H16/air flames at C ≈ 0.

Figure 4.7. Logarithmic sensitivity coefficients of flame burning rate to various parameters computed
for mixtures with no reactivity in the unburned mixture.

The most interesting outcome of the local sensitivity approach is the fact that 𝑚 ̇ is insensitive to LTC
for the conditions considered herein. The important LTC reactions were first identified by performing
sensitivity analysis of 𝜏 utilizing SENKIN (see Appendix D). As demonstrated in Fig. 4.6, the Ai
perturbation for C7H15O2  C7H14OOH (R2) leads to negligible 𝑚 ̇ variation, i.e. (𝐿𝑆𝐶 )
= 0.06; R2 is
one of the most important LTC pathways affecting 𝜏 [17]. 𝑚 ̇ was found to be mostly sensitive to the
high temperature radical branching and termination reactions and those controlling the heat release as
shown in Figs. 4.7 and 4.8. This is consistent with the established understanding of flame chemistry under
non-reactive conditions [31]. It should be noted that if the Ai perturbations are to be performed in the
steady computations [17], adhering to 𝜏 >> 𝜏 criterion is essential to obtain meaningful 𝐿𝑆𝐶 values.
Otherwise, Ai perturbation, would not only modify the flame chemistry but also the reactants entering the
flame as discussed earlier. Large sensitivity to LTC observed in the study of Pan et al. [17] is an artifact
of such ignition dynamics modifications. However, since (𝐿𝑆𝐶 )
was computed locally in the present
study, only the flame chemistry was affected. Therefore, the approach proposed here is more appropriate
to identify the reactions that affect 𝑚 ̇ the most.
6 7 -1 0 1 2 3 4
Case H1
𝑳𝑺𝑪 Case H2
Case C0
Case C1
Case C2
Case C3
T
ad
C
T
u
H + O
2
+ M → HO
2
+ M
HO
2
+ O H → O
2
+ H
2
O
H + O
2
→ O + OH
CO + O H → CO
2
+ H
O H + H
2
→ H + H
2
O
λ

114

Figure 4.8. Logarithmic sensitivity coefficient of flame burning rate to various parameters computed at
different times for Case C3.

All 𝐿𝑆𝐶 values vary as the flame and the unburned mixture evolve in time, and typical results are
shown in Fig. 4.8 for Case C3. With the exception of (𝐿𝑆𝐶 )
, the variations of 𝐿𝑆𝐶 are insignificant even
up to the threshold of transitioning to the convection-reaction controlled regime (tqmax ≈ 0.95𝜏 for the
conditions of Fig. 4.8), which indicates that the parameters controlling the flame do not change drastically
in the conduction-reaction controlled regime.
Finally, (𝐿𝑆𝐶 )
is positive for t/ ign = 0, 0.6. However, (𝐿𝑆𝐶 )
decreases over time and becomes
negative at t/𝜏 = 0.9 due to the H2O production in the unburned mixture. Similar trends are observed
for 𝐿𝑆𝐶 values at non-zero t/𝜏 for other cases investigated (see Appendix E). It should be also noted that
the decoupling approach incorporated enables simplifying the non-linear dependence of 𝑚 ̇ to the
controlling parameters. In another word, the change in 𝑚 ̇ values can be estimated by first calculating the
contribution of each parameters in Eq. (4.1) individually and then adding the contributions (see Appendix
F).

3.5 4 -0.5 0 0.5 1
𝑳𝑺𝑪 T
ad
C
T
u
H + O
2
+ M → HO
2
+ M
HO
2
+ OH → O
2
+ H
2
O
H + O
2
→ O + OH
CO + O H → CO
2
+ H
OH + H
2
→ H + H
2
O
λ
t/ 
ign
= 0.6 t/ 
ign
= 0.9 t/ 
ign
= 0

115
4.4. Concluding remarks
The propagation of transient one-dimensional laminar flames into a reacting mixture was investigated
numerically under constant pressure conditions. The study focused on the parameters influencing the
burning rate in the conduction-reaction controlled regime. Simulations were performed for H 2/O2/N2 and
n-C7H16/air mixtures under engine-relevant conditions.
It was demonstrated that in H2/O2/N2 mixtures, as the pre-flame reaction progresses the burning rate
increases. However, the evolution of the reactivity both in the unburned gas and the flame does not depend
on the initial temperature given that the results collapse upon proper scaling.
On the other hand, for n-C7H16/air mixtures that exhibit low temperature chemistry, the evolution of
the burning rate depends on the initial mixture temperature even for mixtures with the same auto-ignition
delay time. More specifically, the evolution of the burning rate was found to depend on the instantaneous
unburned mixture temperature and composition in a coupled way.
The parameters influencing the burning rate were identified and an approach was developed to isolate
them and to determine their relative importance. The novelty of the proposed approach is the ability to
perform local perturbations such that flame chemistry can be decoupled from the ignition dynamics. It
was shown that the burning rate does not directly depend on the fuel specific reactions even in the presence
of low temperature chemistry. The low temperature chemistry influences the composition of reactants
entering the flame. The burning rate was found to be sensitive only to reactions related to high temperature
chain branching/termination and heat release.

4.5. Appendix A: Details of transient numerical simulations
The mixture-averaged formulation [32] was used to calculate transport properties in TORC. The Soret
effect [33] was ignored as its effect on the burning rate was found to be minimal for the cases investigated.
As depicted in Fig. 4.A.1, including the Soret diffusion minimally affects the time evolution of the flame.
Although, the absolute value of (𝑚 ̇ )o decreased by 5% by including the Soret diffusion in this case, the
transient behavior and the time evolution of the flame which were the main focus of this study did not
change. Moreover, the 𝐿𝑆𝐶 values were found to be identical for both cases.

116

Figure 4.A.1. The effect of Soret diffusion on the evolution of scaled mass burning rate as a function of
the normalized time for Case H2.
Grid independent solutions were obtained by using a minimum spatial grid size ( x) of 0.25 and
0.125 m for H2 and n-C7H16 flames respectively. Figure 4.A.2 demonstrates the grid independence for
two representative cases.

Figure 4.A.2. The effect of minimum spatial grid size on the evolution of the flame displacement speed
as a function of the normalized time.

In order to verify the numerical tolerances incorporated, the evolution of the unburned mixture ahead
of the propagating flame was compared to the results obtained by the zero-dimensional SENKIN code
𝒎 ̇ /(𝒎 ̇ )
o
1
1.1
1.2
1.3
1.4
1.5
0 0.2 0.4 0.6 0.8 1
t/ 
ign
With Soret
Without Soret
100
140
180
220
260
300
340
380
0 0.2 0.4 0.6 0.8 1
t/ 
ign
80
100
120
140
160
180
0 0.2 0.4 0.6 0.8 1
t/ 
ign
S
d
, cm/s
x = 0.0625 m
x = 0.1250 m
x = 0. 25 m
x = 0. 50 m
Case H2 Case C1

117
[30]. Figure 4.A.3 demonstrates the H radical mass fraction in the unburned mixture as a function of the
unburned mixture temperature for Case C1. The SENKIN calculations were performed using the constant
pressure and enthalpy constraints. The excellent agreement between the results using both codes
demonstrates the ability of TORC to capture ignition accurately.

Figure 4.A.3. H radical mass fraction in the unburned mixture as a function of the unburned mixture
temperature for Case C1; results obtained by two codes are shown.

4.6. Appendix B: Details of the developed sensitivity analysis approach
The sensitivity analysis was performed by utilizing an approach similar to a recent study [31], by
perturbing one parameter at a time. In order to perturb Tad, the amount of N2 diluent in the mixture was
modified. Such modification results in change in  and Dk, which was not desirable. The  and Dk
dependence was eliminated by adopting a constant  and Dk formulation in TORC. In this formulation,
the initial mixture conductivity at (Tu)o, o, and an effective diffusion coefficient, Deff, equal to the fuel
diffusion coefficient at (Tu)o were chosen as constants for performing the simulations. The 𝑚 ̇ values
obtained by this formulation were found to deviate by as much as 37% from the original variable transport
formulation. However, since the main focus of the present study was to study the time evolution of the
flame and the parameters influencing 𝑚 ̇ , the absolute differences in 𝑚 ̇ would not change the general
conclusions; the relative differences were found to be negligible. In fact, the 𝑚 ̇ /(𝑚 ̇ )o value obtained by
the constant transport formulation was found to deviate from the variable transport formulation by less
than 1% and 5% for H2 and n-C7H16 flames respectively as demonstrated in Fig. 4.B.1. The larger
Y
H
 10
8
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
800 850 900 950 1000 1050 1100
T
u
, K
SENKIN
TORC

118
discrepancy for n-C7H16 flames was due to larger variation of  for these flames. The variation of the
unburned mixture conductivity ( u) is as also shown in Fig. 4.B.1.

Figure 4.B.1. Evolution of the scaled burning rate and the scaled unburned mixture conductivity as a
function of the normalized time for two representative cases.

By utilizing the constant transport formulation, the new parameter space will be:
𝑚 ̇ = 𝑚 ̇ (Tad, Yk(t), Tu(t), o, Deff) (4.2)
in which only Yk and Tu are varying in time. With the elimination of spatial and temporal dependence of
transport properties, Tad perturbation could be performed by only modifying the amount of N 2 diluent; the
perturbation does not change other parameters in Eq. (4.2). Subsequently, computational experiments
could be performed to assess the effect of Tad on 𝑚 ̇ . First, Tad was computed for the baseline mixtures
using the EQUIL code [34]. Then, transient numerical simulations were performed for the baseline
mixtures utilizing the constant transport formulation. Tad of the mixtures were then slightly perturbed by
N2 molar ratio modifications. At last, transient numerical simulations with the same o and Deff were
performed for mixtures with a perturbed Tad and compared against the baseline. The constant transport
formulation was specifically useful to conduct the parametric study since the non-linear dependence of
diffusion coefficients on temperature [1,33] could be eliminated. It should be noted that the perturbations
performed are all local as elaborated in Section 4.3. Figure 4.B.2 illustrates the computational approach
for computing (𝐿𝑆𝐶 )
.
1
1.05
1.1
1.15
1.2
0 0.2 0.4 0.6 0.8 1
t/ 
ign
1
1.1
1.2
1.3
1.4
1.5
1.6
0 0.2 0.4 0.6 0.8 1
t/ 
ign
𝒎 ̇ /(𝒎 ̇ )
o
, Constant transport

u
/ 
o
, Variable transport
𝒎 ̇ /(𝒎 ̇ )
o
, Variable transport
𝒎 ̇ /(𝒎 ̇ )
o
, Constant transport
𝒎 ̇ /(𝒎 ̇ )
o
, Variable transport
Case H2 Case C2

119

Figure 4.B.2. An illustrative diagram of the computational steps performed for obtaining (𝐿𝑆𝐶 )
; first,
the unburned mixture at t = tpert was extracted from the TORC baseline solution; this mixture has a
temperature of Tu(tpert), fuel mass fraction of YF(tpert) which corresponds to C(tpert), n% volumetric N2, and
an enthalpy which corresponds to Tad; the perturbation comprised of two processes: 1) the mixture was
diluted with more N2 to reduce Tad by 1%, 2) a new TORC simulation was performed using the transport
properties of the baseline mixture; 𝑚 ̇ obtained from both TORC calculations was used for calculating
(𝐿𝑆𝐶 )
.

𝐿𝑆𝐶 values for other parameters were also obtained using a similar approach. Since perturbing Tu also
modifies Tad, N2 diluent amount modifications was utilized again to restore Tad to its baseline value. Figure
4.B.3 illustrates the computational approach for obtaining (𝐿𝑆𝐶 )
.

Figure 4.B.3. An illustrative diagram of the computational steps performed for obtaining (𝐿𝑆𝐶 )
; first,
the unburned mixture at t = tpert was extracted from the TORC baseline solution; this mixture has a
temperature of Tu(tpert), fuel mass fraction of YF(tpert) which corresponds to C(tpert), n% volumetric N2, and
an enthalpy which corresponds to Tad; the perturbation comprised of three processes: 1) the mixture was
heated instantaneously in order to increase Tu by 1%; the heating increases Tad, 2) the mixture was diluted
with more N2 to restore Tad to its baseline value, 3) a new TORC simulation was performed using the
transport properties of the baseline mixture; 𝑚 ̇ obtained from both TORC calculations was used for
calculating (𝐿𝑆𝐶 )
.

T
u
(t
pert
)
C(t
pert
)
T
ad
- T
ad
=0.99 T
ad
(n+ n)%N
2
T
u
(t
pert
)
C(t
pert
)
T
ad
n% N
2
constant transport TORC (1)
2. constant transport TORC (2)
1. diluting
𝒎 ̇ base
𝒎 ̇ pert
(𝐿𝑆𝐶)
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
. T
u
(t
pert
)+ T
u
=1.01T
u
(t
pert
)
C(t
pert
)
T
ad
+ T
ad
n% N
2
T
u
(t
pert
)
C(t
pert
)
T
ad
n% N
2
constant transport TORC (1)
3. constant transport TORC (2)
1. heating
𝒎 ̇ base
𝒎 ̇ pert
(𝐿𝑆𝐶)
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
( )
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
. T
u
(t
pert
)+ T
u
=1.01T
u
(t
pert
)
C(t
pert
)
T
ad
(n+ n)%N
2
2. diluting

120
In order to perturb Yk, the mixture was allowed to react adiabatically for a residence time of t. The
mixture composition and temperature was followed as function of time in an adiabatic reactor at constant
pressure using SENKIN. Tu(tpert + t) and Yk(tpert + t) were then extracted and the mixture was
instantaneously cooled to Tu(tpert) in order to match the baseline temperature. The cooling however resulted
in Tad reduction. Hence, the N2 diluent amount modification was performed again to neutralize Tad
variation effects. At last, a new transient simulation was conducted using the updated composition and
Tu(tpert) as an initial condition. This approach assures perturbation of only Yk in Eq. (4.2). Figure 4.B.4
illustrates the computational approach for obtaining (𝐿𝑆𝐶 )
.

Figure 4.B.4. An illustrative diagram of the computational steps performed for obtaining (𝐿𝑆𝐶 )
; first,
the unburned mixture at t = tpert was extracted from the TORC baseline solution; this mixture has a
temperature of Tu(tpert), fuel mass fraction of YF(tpert) which corresponds to C(tpert), n% volumetric N2, and
an enthalpy which corresponds to Tad; the perturbation comprised of four processes: 1) the mixture was
allowed to react adiabatically at a constant pressure for a residence time of t to obtain a mixture that has
10% less fuel (10% higher C); this mixture has a different temperature and composition (as a result
different percentage N2) compared to the baseline, 2) the mixture was cooled instantaneously to restore
Tu to its baseline value; the cooling decreases Tad, 3) the volumetric percentage N2 in the mixture was
decreased to restore Tad to its baseline value, 4) a new TORC simulation was performed using the transport
properties of the baseline mixture; 𝑚 ̇ obtained from both TORC calculations was used for calculating
(𝐿𝑆𝐶 )
.

T
u
(t
pert
)
C(t
pert
) + C = 1.1C(t
pert
)
T
ad
- T
ad
(n+ n)%N
2
T
u
(t
pert
)
C(t
pert
)
T
ad
n% N
2
T
u
(t
pert
+ t)
C(t
pert
) + C = 1.1C(t
pert
)
T
ad
(n+ n)%N
2
constant transport TORC (1)
3. constant transport TORC (2)
1. Adiabatic reactor at
constant pressure (SENKIN)
𝒎 ̇ base
𝒎 ̇ pert
(𝐿𝑆𝐶)
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
( )
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
. 2. cooling
T
u
(t
pert
)
C(t
pert
) + C = 1.1C(t
pert
)
T
ad
(n+ n -  n)%N
2
3. diluting

121
o perturbation did not require any other modifications as the perturbation do not change other
parameters in (2). Figure 4.B.5 illustrates the computational approach for obtaining (𝐿𝑆𝐶 )
∘
.

Figure 4.B.5. An illustrative diagram of the computational steps performed for obtaining (𝐿𝑆𝐶 )
∘
; first,
the unburned mixture at the time at which the perturbation was about to be performed, tpert, was extracted
from the TORC baseline solution; the perturbation comprised of two processes: 1) o was decreased, 2) a
new TORC simulation was performed using 0.9 o and Deff as constants; 𝑚 ̇ obtained from both TORC
calculations was used for calculating (𝐿𝑆𝐶 )
∘
.

Due to non-linear dependence of Dk on temperature and non-linear dependence of 𝑚 ̇ on Dk [1],
utilizing the constant transport formulation was not justified to quantify the effects of Dk on 𝑚 ̇ . This is
mainly due to the fact that in this formulation, a constant effective diffusivity was used to represent
diffusion coefficient of all species. However, from the existing knowledge of laminar premixed flames in
non-reactive conditions, it is known that Dk can vary in wide range of values depending on the species
[25,31]. Moreover, each Dk has different effect on 𝑚 ̇ and often the effects are opposite such that for
species l  k, conditions exist such that (𝐿𝑆𝐶 )
> 0 while (𝐿𝑆𝐶 )
< 0 [25,31]. Hence, using Deff to
represent the effects of various Dk can lead to wrong conclusions.
Due to failure of constant transport formulation to properly address the dependence of 𝑚 ̇ on Dk, a
more rigorous approach was developed to calculate (𝐿𝑆𝐶 )
. In the new approach, Dk was calculated
utilizing the mixture-averaged formulation originally incorporated in TORC while constant  was
assumed i.e., constant  formulation. o was utilized as the appropriate constant throughout the simulations
similar to the constant transport approach. Constant  formulation is specifically immune to the possible
 variations that might occur as a result of Dk perturbations; Dk perturbations can essentially modify the
temperature and Yk distribution in the flame and therefore modify . Therefore, constant  formulation
assures variation of one parameter at a time.
T
u
(t
pert
)
C(t
pert
)

o
constant transport TORC (1)
2. constant transport TORC (2)
1. 
o
perturbation
T
u
(t
pert
)
C(t
pert
)

o
-  
o
=0.9 
o
𝒎 ̇ base
𝒎 ̇ pert
(𝐿𝑆𝐶)
∘
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
.

122
The Dk perturbations in the constant  formulation was performed by generating a new version of
constant  TORC for each species by artificially multiplying Dk(x,t) obtained by the mixture-averaged
calculation by a constant multiplication factor. Therefore, k number of versions were generated to be used
to quantify the effect of each Dk on 𝑚 ̇ . First, a constant  TORC simulation was conducted for the baseline
parameters. Then, for each specific k, the perturbed Dk version of constant  TORC was utilized to obtain
a new 𝑚 ̇ that was subsequently compared against the baseline value. Similar approaches have been
previously incorporated to examine the effects of transports in the flame propagation in reacting mixtures
[8,13]. Figure 4.B.6 illustrates the computational approach for obtaining (𝐿𝑆𝐶 )
.

Figure 4.B.6. An illustrative diagram of the computational steps performed for obtaining (𝐿𝑆𝐶 )
; first,
the unburned mixture at the time at which the perturbation was about to be performed, tpert, was extracted
from the baseline solution obtained by the constant  TORC formulation; the perturbation comprised of
two processes: 1) Dk(x,t) was decreased by a constant multiplication factor of 0.9, 2) a new simulation
was performed using the generated version of constant  TORC for species k that uses 0.9Dk(x,t) and o
as constants; 𝑚 ̇ obtained from both TORC calculations was used for calculating (𝐿𝑆𝐶 )
.

Figure 4.B.7 summarizes the sensitivity analysis results for Dk obtained at the initial quasi-steady
propagation (C = 0). Apart from (𝐿𝑆𝐶 )
for H2 flames, the other (𝐿𝑆𝐶 )
values were negligible
compared to the other 𝐿𝑆𝐶 values shown in Fig. 4.7. Moreover, (𝐿𝑆𝐶 )
was found to decreases over
time (as C increase) for all reacting cases. This is mainly due to flame transition to the regime that is not
controlled by conduction-reaction.
T
u
(t
pert
)
C(t
pert
)
D
k
(x, t
pert
)
constant  TORC (1)
2. constant  TORC (2)
1. D
k
perturbation
T
u
(t
pert
)
C(t
pert
)
D
k
(x, t
pert
) -  D
k
(x, t
pert
) = 0.9D
k
(x, t
pert
)
𝒎 ̇ base
𝒎 ̇ pert
(𝐿𝑆𝐶)
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
=
𝒎 ̇ pert
𝒎 ̇ base
𝒎 ̇ base
×
.

123

Figure 4.B.7. Logarithmic sensitivity coefficients of mass-burning rate to diffusion coefficients for
mixtures with no reactivity in the unburned mixture.

It should be noted that due to the presence of numerical uncertainty on the order of 2% in 𝑚 ̇
determination, the perturbations had to be chosen such that the resulting change in 𝑚 ̇ was larger than 4%.
Perturbation values of 1% for Tad and Tu and 10% for C, o, Dk and Ai were found to be sufficient to obtain
realizable 𝑚 ̇ variations.

4.7. Appendix C: Transient effects of stopping and restarting the simulation
The transient simulations were initiated by igniting the mixture to generate a flame. First, the
equilibrium composition and Tad were calculated using EQUIL for the thermodynamics conditions of the
unburned mixture. These quantities were further utilized to specify the initial hot pocket of gas as an
ignition source. Initially, the temperature of the entire domain was set to (Tu)o except at the hot pocket
where the temperature was set to Tad. In order to avoid discontinuity in the initial temperature profile, an
initial hyperbolic tangent profile was incorporated. The thickness of the profile was chosen to be close to
the flame thickness at the corresponding initial conditions. The initial flame thickness was estimated using
the PREMIX code [35]. However, since obtaining a PREMIX solution for reactive mixtures requires
particular care and interpretation [3], PREMIX simulations were only performed for unreactive mixtures
at reduced (Tu)o corresponding to C ≈ 0. Then o was estimated by linearly extrapolating the flame
thickness values obtained from the previous step. Similarly, the initial Yk profile contained two flat parts
joined by a hyperbolic tangent profile. The initial hot pocket filled the computational domain from x = 0
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4
Case H1
𝑳𝑺𝑪 Case H2
Case C0
Case C1
Case C2
Case C3
𝑫 𝐇 𝟐 𝑫 𝐎 𝟐 𝑫 𝐇 𝑫 𝐎𝐇

124
to x = xpocket = 10 o. Tests were conducted to evaluate the sensitivity of the results to xpocket. No significant
differences were observed by using xpocket = 5 o.
The perturbation procedure required stopping the baseline simulation and starting a new simulation
by initiating a new flame using a hot pocket of equilibrium products. Therefore, the transient effects of
the flame initiation process had to be assessed to assure that they do not influence the results. This was
achieved by stopping the baseline simulations and restarting a new one without imposing any
perturbations. In this test, the mixture was re-ignited at the same conditions that the simulation was
stopped. Figure 4.C.1 depicts the 𝑚 ̇ time evolution of the restarted cases along with their baseline 𝑚 ̇ for
couple of representative cases. The initial transient process is apparent by the sudden drop of 𝑚 ̇ value in
Fig. 4.C.1. After the sudden drop, 𝑚 ̇ matches its baseline value. Moreover, the total time of the transient
flame initiation process was less than 1% of the total simulation time for the restarted cases. Hence, the
transient effects of restarting the simulations were negligible.

Figure 4.C.1. Comparison between the time evolution of the scaled mass-burning rate of the baseline
and the cases obtained by stopping and restarting the baseline simulations for representative cases;
arrows indicate the time at which the perturbation was imposed.

In order to evaluate the effect of flame initiation method on the results, the simulation was repeated
for Case H2 with a different method of flame initiation. Instead of equilibrium products, unreactive
mixture composition at Tad was used as a hot pocket to initiate the flame. Results are shown in Fig. 4.C.2
which demonstrates insignificant change in flame time history.
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
0 0.2 0.4 0.6 0.8 1
𝒎 ̇ /(𝒎 ̇ )
o
t/ 
ign
Case C1, baseline
Case C1, t
pert
/ 
ign
= 0.80
Case H2, baseline
Case H2, t
pert
/ 
ign
= 0.91

125

Figure 4.C.2. The effect of flame initiation method on the evolution of the scaled mass-burning rate as
a function of the normalized time for Case H2.

4.8. Appendix D: Sensitivity analysis results for 𝝉 𝐢𝐠𝐧
Figure 4.D.1 depicts the sensitivity analysis results for 𝜏 computed using SENKIN which shows
that 𝜏 is sensitive to reactions that are different from the ones that 𝑚 ̇ is sensitive to (shown in Figs.
4.E.1 and 4.E.2).

Figure 4.D.1. Logarithmic sensitivity coefficient of the auto-ignition delay time to the reactions pre-
exponential factors.

1
1.05
1.1
1.15
1.2
0 0.2 0.4 0.6 0.8 1
t/ 
ign
𝒎 ̇ /(𝒎 ̇ )
o
Equilibrium products
Unreactive mixture
-1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
(𝑳𝑺𝑪 )
𝑨 𝒊 ≈ ( ∆𝝉 𝒊𝒈𝒏 / ∆𝑨 𝒊 ) (𝑨 𝒊 /𝝉 𝒊𝒈𝒏 )
Case H1
HO
2
+ HO
2
→ H
2
O
2
+ O
2
H + O
2
+ M → HO
2
+ M
Case H2
H + O
2
→ O + OH
O H + O H + M → H
2
O
2
+ M
H + H
2
O
2
→ H
2
+ HO
2
Case C1
Case C2
Case C3
O H + O H + M  H
2
O
2
+ M
C
7
H
15
O
2
 C
7
H
14
OOH
O
2
C
7
H
15
O  C
7
H
13
O
2
+ OH
n-C
7
H
15
+ O
2
 C
7
H
15
O
2
C
7
H
15
O
2
 n-C
7
H
15
+ O
2
i-C
7
H
15
+ O
2
 C
7
H
15
O
2
C
7
H
15
O
2
 i-C
7
H
15
+ O
2
(𝑳𝑺𝑪 )
𝑨 𝒊 ≈ ( ∆𝝉 𝒊𝒈𝒏 / ∆𝑨 𝒊 ) (𝑨 𝒊 /𝝉 𝒊𝒈𝒏 )

126
4.9. Appendix E: Sensitivity analysis results for 𝒎 ̇
Figures 4.E.1 and 4.E.2 summarize the sensitivity analysis results for 𝑚 ̇ for typical reacting cases
considered.

Figure 4.E.1. Logarithmic sensitivity coefficient of mass-burning rate to various parameters computed
at different times for Case H2.

Figure 4.E.2. Logarithmic sensitivity coefficient of mass burning rate to various parameters computed
at different times for Case C1.

-2 -1 0 1 2 3 4 5 6 7
𝑳𝑺𝑪 T
ad
C
T
u
H + O
2
+ M → HO
2
+ M
HO
2
+ O H → O
2
+ H
2
O
H + O
2
→ O + OH
O H + H
2
→ H + H
2
O
λ
t/ 
ign
= 0.91
t/ 
ign
= 0.96
t/ 
ign
= 0
-0.3 0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3 3.3 3.6 3.9 4.2
𝑳𝑺𝑪 T
ad
C
T
u
H + O
2
+ M → HO
2
+ M
HO
2
+ O H → O
2
+ H
2
O
H + O
2
→ O + OH
OH + H
2
→ H + H
2
O
λ
t/ 
ign
= 0.85
t/ 
ign
= 0.93
t/ 
ign
= 0
CO + O H → CO
2
+ H

127
4.10. Appendix F: Predicting the 𝒎 ̇ behavior using LSC values
The decoupling approach incorporated enables simplifying the non-linear dependence of 𝑚 ̇ to the
controlling parameters. In another word, the change in 𝑚 ̇ values (∆𝑚 ̇ ) can be estimated by first calculating
the contribution of each parameters in Eq. (4.2) individually and then adding the contributions. The
process is essentially a Taylor expansion of 𝑚 ̇ since 𝐿𝑆𝐶 can be considered as a first order partial
derivative. Hence, ∆𝑚 ̇ at each instant is approximately the superposition of effect of individual controlling
parameters which can be expressed as:
(∆ ̇ )
(𝑚 ̇ )
°
≈ [(𝐿𝑆𝐶 )
× ∆𝑇 /𝑇 ] + [(𝐿𝑆𝐶 )
× ∆𝑇 /𝑇 ] + [(𝐿𝑆𝐶 )
∘
× ∆𝜆 /𝜆 ] + [(𝐿𝑆𝐶 )
× ∆𝐶 /𝐶 ] (4.3)
In order to demonstrate the concept, (∆𝑚 ̇ )predicted/(𝑚 ̇ )o was calculated and compared to the actual
∆𝑚 ̇ /(𝑚 ̇ )o value for Case C1 at different time intervals. For each time interval, the change in the controlling
parameters (∆Tad, ∆Tu, ∆ , and ∆C) were calculated using the mixture states at the available solutions for
the end points of the interval. Then a representative state (start point or end point of the time interval) was
used to extract the 𝐿𝑆𝐶 values and controlling parameters to be plugged into Eq. (4.3) and obtain
(∆𝑚 ̇ )predicted/(𝑚 ̇ )o. The results are shown in Table 4.F.1 which demonstrates that there are only small
differences between the predicted and actual values of ∆𝑚 ̇ . Hence, the decoupling approach successfully
separates the physics involved and appropriately explains the behavior of 𝑚 ̇ .

Table 4.F.1. Comparison between the actual change in mass-burning rate to the predicted values for Case
C1.
Time interval
Time at which 𝐿𝑆𝐶 values and controlling
parameters are extracted to be used in prediction
equation (4.3)
(∆𝑚 ̇ )predicted/(𝑚 ̇ )o ∆𝑚 ̇ /(𝑚 ̇ )o
0 < t/𝜏 < 0.6 t/𝜏 = 0 0.05 0.03
0.6 < t/𝜏 < 0.85 t/𝜏 = 0.85 0.28 0.32
0.85 < t/𝜏 < 0.93 t/𝜏 = 0.93 0.11 0.08

4.11. Acknowledgements
The Air Force Office of Scientific Research (Grant No. FA9550-15-1-0409) supported this work under
the technical supervision of Dr. Chiping Li.

128
4.12. References
[1] C.K. Law, Combustion physics, Cambridge university press, 2010.
[2] F.N. Egolfopoulos, N. Hansen, Y. Ju, K. Kohse-Höinghaus, C.K. Law, Advances and challenges in
laminar flame experiments and implications for combustion chemistry, F. Qi, Prog. Energy
Combust. 43 (2014) 36-67.
[3] J.B. Martz, H. Kwak, H.G. Im, G.A. Lavoie, D.N. Assanis, Combustion regime of a reacting front
propagating into an auto-igniting mixture, Proc. Combust. Inst. 33 (2011) 3001-3006.
[4] A. Krisman, E.R. Hawkes, J.H, Chen, The structure and propagation of laminar flames under
autoignitive conditions, Combust. Flame 188 (2018) 399-411.
[5] C.M. White, R.R. Steeper, A.E. Lutz, The hydrogen-fueled internal combustion engine: a technical
review, Int. J. of Hydrogen Energy 31 (2006) 1292-1305.
[6] Z. Wang, H. Liu, H., R.D. Reitz, Knocking combustion in spark-ignition engines, Prog. Energy
Combust. 61 (2017) 78-112.
[7] Y. Ju, W. Sun, M.P. Burke, X. Gou, Z. Chen, Ju, Y., Sun, W., Burke, M. P., Gou, X., & Chen, Z.
(2011). Multi-timescale modeling of ignition and flame regimes of n-heptane-air mixtures near spark
assisted homogeneous charge compression ignition conditions, Proc. Combust. Inst. 33 (2011) 1245-
1251.
[8] J.B. Martz, G.A. Lavoie, H.G. Im, R.J. Middleton, A. Babajimopoulos, D.N. Assanis, The
propagation of a laminar reaction front during end-gas auto-ignition, Combust. Flame 159 (2012)
2077-2086.
[9] H. Oshibe, H. Nakamura, T. Tezuka, S. Hasegawa, K. Maruta, Stabilized three-stage oxidation of
DME/air mixture in a micro flow reactor with a controlled temperature profile, Combust. Flame 157
(2010) 1572-1580.
[10] S.H. Won, B. Windom, B. Jiang, Y. Ju, The role of low temperature fuel chemistry on turbulent
flame propagation, Combust. Flame 161 (2014) 475-483.
[11] C.B. Reuter, O. Yehia, Y. Ju, S.H. Won, Effect of Low-Temperature Reactivity on the Turbulent
Combustion of n-Octane/iso-Octane Mixtures in a Reactor-Assisted Turbulent Slot Burner, in: 55th
AIAA Aerospace Sciences Meeting, Grapevine, Texas, USA, 2017, p. 1781.
[12] W.J. Pitz, C.K. Westbrook, Interactions between a laminar flame and end gas autoignition, No.
UCRL-92295 CONF-850830-1, Lawrence Livermore National Lab, Livermore, CA, USA, 1985.
[13] H. Yu, Z. Chen, End-gas autoignition and detonation development in a closed chamber, Combust.
Flame 162 (2015) 4102-4111.

129
[14] P. Dai, Z. Chen, S. Chen, Y. Ju, Numerical experiments on reaction front propagation in n-
heptane/air mixture with temperature gradient, Proc. Combust. Inst. 35 (2015) 3045-3052.
[15] H. Terashima, M. Koshi, Mechanisms of strong pressure wave generation in end-gas autoignition
during knocking combustion, Combust. Flame 162 (2015) 1944-1956.
[16] J. Pan, G. Shu, P. Zhao, H. Wei, Z. Chen, Interactions of flame propagation, auto-ignition and
pressure wave during knocking combustion, Combust. Flame 164 (2016) 319-328.
[17] J. Pan, H. Wei, G. Shu, Z. Chen, P. Zhao, The role of low temperature chemistry in combustion
mode development under elevated pressures, Combust. Flame 174 (2016) 179-193.
[18] C. Qi, P. Dai, H. Yu, Z. Chen, Different modes of reaction front propagation in n-heptane/air mixture
with concentration non-uniformity, Proc. Combust. Inst. 36 (2017) 3633-3641.
[19] W. Zhang, M. Faqih, X. Gou, Z. Chen, Numerical study on the transient evolution of a premixed
cool flame, Combust. Flame 187 (2018) 129-136.
[20] P. Zhao, W. Liang, S. Deng, C.K. Law, Initiation and propagation of laminar premixed cool flames,
Fuel 166 (2016) 477-487.
[21] H. Hu, J. Keck, Autoignition of adiabatically compressed combustible gas mixtures, SAE technical
paper 872110 (1987).
[22] D. Kaufmann, P. Roth, Numerical simulation of one-dimensional laminar flames propagating into
reacting premixed gases, Combust. Flame 80 (1990) 385-394.
[23] J.H. Chen, E.R. Hawkes, R. Sankaran, S.D. Mason, H.G. Im, Direct numerical simulation of ignition
front propagation in a constant volume with temperature inhomogeneities: I. Fundamental analysis
and diagnostics, Combust. Flame 145 (2006) 128-144.
[24] X.J. Gu, D.R. Emerson, D. Bradley, Modes of reaction front propagation from hot spots, Combust.
Flame 133 (2003) 63-74.
[25] J. Jayachandran, R. Zhao, Determination of laminar flame speeds using stagnation and spherically
expanding flames: molecular transport and radiation effects, F.N. Egolfopoulos, Combust. Flame
(2014) 2305-2316.
[26] J. Jayachandran, F.N. Egolfopoulos, Thermal and Ludwig–Soret diffusion effects on near-boundary
ignition behavior of reacting mixtures, Proc. Combust. Inst. 36 (2017) 1505-1511.
[27] C. Xiouris, T. Ye, J. Jayachandran, F.N. Egolfopoulos, Laminar flame speeds under engine-relevant
conditions: Uncertainty quantification and minimization in spherically expanding flame
experiments, Combust. Flame 163 (2016) 270-283.
[28] S.G. Davis, A.V. Joshi, H. Wang, F.N. Egolfopoulos, An optimized kinetic model of H2/CO
combustion, Proc. Combust. Inst. 30 (2005) 1283-1292.

130
[29] S. Liu, J.C. Hewson, J.H. Chen, H. Pitsch, Effects of strain rate on high-pressure nonpremixed n-
heptane autoignition in counterflow, Combust. Flame 137 (2004) 320-339.
[30] A.E. Lutz, R.J. Kee, J.A. Miller, SENKIN: A FORTRAN program for predicting homogeneous gas
phase chemical kinetics with sensitivity analysis, SANDIA National Laboratories Report, SAND87-
8248, Livermore, CA, USA, 1990.
[31] J. Smolke, F. Carbone, F.N. Egolfopoulos, H. Wang, Effect of n-dodecane decomposition on its
fundamental flame properties, Combust. Flame 190 (2018) 65-73.
[32] J. Warnatz, Ber. Bunsenges. Calculation of the structure of laminar flat flames I: Flame velocity of
freely propagating ozone decomposition flames, Phys. Chem. 82 (1978) 193.
[33] D.E. Rosner, Transport Processes in Chemically Reacting Flow Systems, Dover Publications, 2000.
[34] A.E. Lutz, F.M. Rupley, R.J. Kee, W.C. Reynolds, E. Meeks, EQUIL: A CHEMKIN
implementation of STANJAN for computing chemical equilibria, Sandia Technical Report No.,
Sandia National Laboratories, Livermore, CA, USA, 1998.
[35] R.J. Kee, J.F. Grcar, M.D. Smooke, J.A. Miller, Premix: a FORTRAN program for modeling steady
laminar one-dimensional premixed flames, Sandia Report, SAND85- 8240, Sandia National
Laboratories, 1985.

131
Chapter 5
Concluding Remarks and Recommendations
5.1. Counterflow configuration
The counterflow configuration was investigated experimentally and computationally in search of
conditions and nozzle designs that best conform to the assumptions of quasi one-dimensional formulation.
It was found that solid surfaces must be removed from the vicinity of the jets exit in order to avoid flow
bifurcation that leads to two-dimensionality and violating the 1-D assumptions. Furthermore, it was shown
that the nozzle designs suggested by Niemann et al. [1] with the flow conditioning screens at the nozzle
exits best conform to the 1-D assumptions. Moreover, the experiments must be performed at Reynolds
number values of approximately 1000 to eliminate the sensitivity of the flow field to boundary conditions.
The flow field of isothermal counterflows in the absence of confining solid surfaces were found to be
hydrodynamically stable. However, both transient axisymmetric simulations and experiments revealed
that the flow field is extremely sensitive to perturbations for Reynolds numbers higher than 1500. The
sensitivity was observed by introducing a small change in the average velocity of one of the jets. Such
perturbation forced the system to go through a transient unsteady state before settling in a new equilibrium.
It was observed that the time to reach an equilibrium was extremely large (approximately 200 non-
dimensional time). The large transition time has been previously observed in other shear flows (e.g.,
[2][3]). Since the eigenmodes of Navier-Stokes equations are non-normal at high Reynolds numbers, any
perturbation will initially grow due to interactions of non-normality and non-linearity [2][3]. The
perturbations will later damp and the system will go back to the initial equilibrium state.
Hence, it is postulated that the axisymmetric counterflows are convectively unstable at high Reynolds
numbers (e.g., [4]) which means that any perturbation will result in an initial instability. The instability
however, will convect out of the jets impingement region for large times. This fact raises questions about
the validity of any transient measurement in the counterflow configuration. Particularly since ignition
(extinction) experiment is performed by gradually increasing (decreasing) the fuel content, it is inherently
a transient experiment. Therefore, ignition and extinction experiments in the counterflow configuration
can be unsteady processes in nature at high enough Reynolds numbers.
It is useful to conduct full transient axisymmetric ignition and extinction simulations at high Reynolds
numbers in order to investigate the behavior of the system during the transition between burning and non-
burning states. Such investigation will reveal the dynamics of the transient events during the extinction

132
and ignition. The study can also assess if the experiment can be modeled by 1-D codes which assumes a
steady flow field.
5.2. End gas ignition by spherically expanding flames in a constant volume chamber
The system of spherically expanding flame propagating in a constant volume chamber was
investigated in this dissertation in terms of the physical phenomena occurring when the end gas reaches
high enough pressures and temperatures such that it starts to react. The main focus of the study was the
flame itself along with the parameters controlling its burning rate. It was found that the flame is not
sensitive to the important reactions for ignition. The ignition dynamics only modifies the reactants that
are being consumed by the flame.
Identifying the parameters and physics controlling the flame propagation into a reacting mixture is
only the first step to understand the rich physics involved in the end gas ignition in the constant volume.
Hopefully, the results presented in this dissertation will help understanding the various physics involved
in end gas ignition experiments in a constant volume chamber.
The feasibility of achieving end gas ignition has been already demonstrated by initial tests. Moreover,
the Direct Numerical Simulation (DNS) capability [5][6] enables capturing all the relevant physics in the
experiment and comparisons can be performed with confidence with the models’ predictions. End gas
ignition experiments in a constant volume chamber is relatively cheaper and less complex compared to
shock tube and rapid compression machine experiments.
However, one question has to answered to validate the experimental approach: Is the observed pressure
oscillation a marker of end gas ignition only? Or it is due to flame instability? The propagation of a stable
flame in a constant volume generates a rate change of pressure that varies linearly as the pressure in the
chamber increases [7]. Hence, the question can be rephrased as: Is the change in slope for rate change of
pressure in the chamber due to reactivity of the end gas or flame instability?
It is essential to answer the mentioned questions since it has been shown that pressure oscillation can
be generated by flame instabilities alone [8]. Preliminary experimental results indicate that there is no
flame instability in the late stage of end gas compression if the hydrodynamic and thermo-diffusional
stability criterion is followed (e.g., [7]). Furthermore, initial tests revealed that the slope for rate change
of pressure can be modified by only doping the mixture with small amount of ignition promoters. Such
slope modification is in fact due to reactivity of the end gas as the trace amount of impurities do not change
the flame properties or its stability [9]. Therefore, circumstantial evidence exists to confirm that the
observed pressure oscillations are generated by the end gas ignition.

133
It will be useful if two-dimensional simulations of a spherical flame are performed in a constant
volume chamber to confirm that un-doped mixtures create flames that propagate without any instability
generated acceleration. Only then it can be claimed that the change in slope of rate change of pressure is
indeed due to reactivity.

5.3. References
[1] U. Niemann, K. Seshadri, F.A. Williams, Accuracies of laminar counterflow flame experiments,
Combust. Flame 162 (2015) 1540-1549.
[2] K. Balasubramanian, R.I. Sujith, Thermoacoustic instability in a Rijke tube: Non-normality and
nonlinearity, Physics of Fluids 20 (2008) 044103.
[3] S. Cherubini, J. Robinet, P. De Palma, The effects of non-normality and nonlinearity of the Navier–
Stokes operator on the dynamics of a large laminar separation bubble, Physics of Fluids 22 (2010)
014102.
[4] Godrèche, C., Manneville, P., Castaing, B.: Hydrodynamics and nonlinear instabilities. Cambridge
University Press, Cambridge (2005)
[5] J. Jayachandran, R. Zhao, F.N. Egolfopoulos, Determination of laminar flame speeds using
stagnation and spherically expanding flames: molecular transport and radiation effects, Combust.
Flame 161 (2014) 2305-2316.
[6] J. Jayachandran, F.N. Egolfopoulos, Thermal and Ludwig–Soret diffusion effects on near-boundary
ignition behavior of reacting mixtures, Proc. Combust. Inst. 36 (2017) 1505-1511.
[7] C. Xiouris, T. Ye, J. Jayachandran, F.N. Egolfopoulos, Laminar flame speeds under engine-relevant
conditions: Uncertainty quantification and minimization in spherically expanding flame
experiments, Combust. Flame 163 (2016) 270-283.
[8] A. S. Al-Shahrany, D. Bradley, M. Lawes, K. Liu, R. Woolley, Darrieus–Landau and thermo-
acoustic instabilities in closed vessel explosions, Combust. Sci. Tech. 178 (2006) 1771-1802.
[9] O. Mathieu, J. W. Hargis, E. L. Petersen, J. Bugler, H. J. Curran, F. Güthe, The Effect of Impurities
on Ignition Delay Times And Laminar Flame Speeds Of Syngas Mixtures At Gas Turbine
Conditions, ASME Turbo Expo Turbine Technical Conference and Exposition (2014).

Abstract (if available)

Abstract Worldwide dependence on energy supply and its ever increasing demand is a known and over-discussed topic. Burning fossil fuels to produce energy is particularly of interest due to their high energy density. Accurate chemistry models are required to understand and predict the combustion phenomena occurring when burning fuels. The chemistry models are developed and validated using measurements performed in “legacy” experiments. Such experiments must have controllable and well-defined observables such that the physics can be modeled with confidence. Otherwise, the measurements taken in the experiments cannot be used to validate chemistry models. In this dissertation, the counterflow configuration, as one of the legacy combustion experiments, was evaluated in terms of its controllability. The counterflow configuration is composed of impingement of two gaseous jets issuing into a stagnant atmosphere from coaxial cylindrical ducts. The accuracy of a widely used model to express the counterflow configuration (quais one-dimensional formulation) was assessed experimentally and computationally. Conditions were identified at which the underlying assumptions of the model fail. More specifically, the conditions at which the flow field becomes unstable and two-dimensional were identified. The consequences on the validity of the model were quantified and guidelines in experimental design were suggested to maximize conformity with the assumptions of the model. In another effort, this dissertation addresses the physics involved in combustion under engine relevant conditions. The flame propagation into a reacting mixture was studied computationally to understand the underlying physical mechanisms. The purpose of the investigation was to identify the parameters influencing the burning rate of the flame in extreme conditions that the unburned mixture is auto-igniting. A new approach was developed allowing for the decoupling the flame chemistry from the ignition dynamics as well as for the decoupling of parameters influencing the burning rate, so that meaningful sensitivity analysis could be performed.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Modeling investigations of fundamental combustion phenomena in low-dimensional configurations

PDF

Determination of laminar flame speeds under engine relevant conditions

PDF

Studies of combustion characteristics of heavy hydrocarbons in simple and complex flows

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Studies of methane counterflow flames at low pressures

PDF

Flame ignition studies of conventional and alternative jet fuels and surrogate components

PDF

Experimental and kinetic modeling studies of flames of H₂, CO, and C₁-C₄ hydrocarbons

PDF

Experimental studies of high pressure combustion using spherically expanding flames

PDF

End-gas autoignition investigations using confined spherically expanding flames

PDF

Flame characteristics in quasi-2D channels: stability, rates and scaling

PDF

Pressure effects on C₁-C₂ hydrocarbon laminar flames

PDF

Investigations of fuel effects on turbulent premixed jet flames

PDF

Studies of siloxane decomposition in biomethane combustion

PDF

Investigations of fuel and hydrodynamic effects in highly turbulent premixed jet flames

PDF

A theoretical study of normal alkane combustion

PDF

Studies on the flame dynamics and kinetics of alcohols and liquid hydrocarbon fuels

PDF

A theoretical investigation in multi-dimensional and non-thermal plasma effects on combustion processes and pollutant remediation

PDF

Mesoscale SOFC-based power generator system: modeling and experiments

PDF

Re-assessing local structures of turbulent flames via vortex-flame interaction

PDF

CFD design of jet-stirred chambers for turbulent flame and chemical kinetics experiments

Asset Metadata

Creator Ansari, Abtin (author)

Core Title Accuracy and feasibility of combustion studies under engine relevant conditions

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publication Date 07/20/2018

Defense Date 06/05/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag chemistry models,combustion,counterflow configuration,ignition,laminar flame,model validation,numerical simulation,OAI-PMH Harvest,stability analysis

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Egolfopoulos, Fokion (committee chair), Bermejo-Moreno, Ivan (committee member), Lee, Jiin-Jen (committee member), Ronney, Paul (committee member)

Creator Email abtin.ansari88@gmail.com,abtinans@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-21132

Unique identifier UC11668684

Identifier etd-AnsariAbti-6439.pdf (filename),usctheses-c89-21132 (legacy record id)

Legacy Identifier etd-AnsariAbti-6439.pdf

Dmrecord 21132

Document Type Dissertation

Format application/pdf (imt)

Rights Ansari, Abtin

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA