ABSTRACT 
 
 
 
 
Title of Document: EFFECT OF STRUCTURAL HEAT 
CONDUCTION ON  THE PERFORMANCE 
OF MICRO-COMBUSTORS AND MICRO-
THRUSTERS 
  
 Timothy Thierry Leach, Doctor of Philosophy, 
2005 
Directed By: Professor Christopher P. Cadou, Department of 
Aerospace Engineering 
 
 
This thesis investigates the effect of gas-structure interaction on the design and 
performance of miniaturized combustors with characteristic dimensions less than a 
few millimeters.  These are termed ?micro-combustors? and are intended for use in 
devices ranging from micro-scale rocket motors for micro, nano, and pico-satellite 
propulsion, to micro-scale engines for micro-Unmanned Air Vehicle (UAV) 
propulsion and compact power generation.  Analytical models for the propagation of 
a premixed laminar flame in a micro-channel are developed.  The models? predictions 
are compared to the results of more detailed numerical simulations that incorporate 
multi-step chemistry, distributed heat transfer between the reacting gas and the 
combustor structure, heat transfer between the combustor and the environment, and 
heat transfer within the combustor structure.  The results of the modeling and 
simulation efforts are found to be in good qualitative agreement and demonstrate that 
  
the behavior of premixed laminar flames in micro-channels is governed by heat 
transfer within the combustor structure and heat loss to the environment. 
The key findings of this work are as follows:  First, heat transfer through the 
micro-combustor?s structure tends to increase the flame speed and flame thickness.  
The increase in flame thickness with decreasing passage height suggests that micro-
scale combustors will need to be longer than their conventional-scale counterparts.  
However, the increase in flame speed more than compensates for this effect and the 
net effect is that miniaturizing a combustor can increase its power density 
substantially.  Second, miniaturizing chemical rocket thrusters can substantially 
increase thrust/weight ratio but comes at the price of reduced specific impulse (i.e. 
overall efficiency).  Third, heat transfer through the combustor?s structure increases 
steady-state and transient flame stability.  This means that micro-scale combustors 
will be more stable than their conventional-scale counterparts.  Fourth and finally, the 
extended temperature profile associated with the broadened flame causes a different 
set of elementary reactions to dominate the operation of the overall reaction 
mechanism at the micro-scale.  This suggests that new chemical mechanisms may 
need to be developed in order to accurately simulate combustion at small-scales.  It 
also calls into question the efficacy of single-step mechanisms presently used by other 
researchers.  
 
 
 
  
 
 
 
 
 
 
 
 
 
 
 
EFFECT OF STRUCTURAL HEAT CONDUCTION ON  THE PERFORMANCE 
OF MICRO-COMBUSTORS AND MICRO-THRUSTERS 
 
 
 
By 
 
 
Timothy Thierry Leach 
 
 
 
 
 
Dissertation submitted to the Faculty of the Graduate School of the  
University of Maryland, College Park, in partial fulfillment 
of the requirements for the degree of 
Doctor of Philosophy 
2005 
 
 
 
 
 
 
 
 
 
Advisory Committee: 
Professor Christopher P. Cadou, Chair/Advisor 
Professor James Baeder 
Professor Georg Dolzmann 
Professor Andre Marshall 
Professor Kenneth Yu 
  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
? Copyright by 
Timothy Thierry Leach 
2005 
 
 
 
 
 
 ii 
 
Dedication 
 
 
 
To Jennifer and Anne-Sophie. 
 
 iii 
 
Acknowledgements 
 
First of all, I would like to acknowledge the Air Force Office of Scientific 
Research (AFOSR) and Dr. Mitat Birkan who supported this work under AFOSR 
F496200110435. 
I would like to express my sincere gratitude to my advisor Chris Cadou for all he 
did for me during my graduate education at the University of Maryland.  Thank you 
for supporting me, teaching me how to do research, giving me the opportunity to take 
initiative, believing in me, and for all the interesting discussions. 
I would also like to thank the members of my dissertation committee for their time 
and insightful discussions: Professors Baeder, Dolzmann, Marshall, and Yu. 
I would as well like to acknowledge the great help of Dr. Greg Jackson from the 
department of Mechanical Engineering at the University of Maryland, who provided 
me with the initial version of the numerical model used in this work and assisted in its 
modification. 
Also, I would like to thank all of my friends at the University of Maryland: Anand, 
Daniel, John, Kiran, Scott, Shyam. 
I would also like to thank Anand for all the interesting discussions and for his 
experimental data. 
I would also like to express my deepest gratitude to Dr. Edna Hokenson for her 
incredible generosity which helped Jennifer and me so much.  You will always be in 
our thoughts.  RIP. 
 
 iv 
 
 I would like to express my greatest appreciation to my parents Pascale and 
Edward and parents in law Judy and Ray for all their support.  Thank you Ray for 
supporting me, believing in me, for all the always interesting discussions (and for the 
Rum and cigars).   
Last, but not least, I would like to thank my wife Jennifer for all her support, love, 
and patience during these past four years and a half.  You made me the man I am and 
I would not have achieved what I did without you.  And, thank you also for giving me 
the most beautiful baby in the world! 
 
 v 
 
Table of Contents 
 
 
Dedication .................................................................................................................ii 
Acknowledgements ..................................................................................................iii 
Table of Contents ......................................................................................................v 
List of Tables............................................................................................................ix 
List of Figures ...........................................................................................................x 
List of Figures ...........................................................................................................x 
List of Abbreviations..............................................................................................xix 
List of Abbreviations..............................................................................................xix 
Chapter 1: Introduction........................................................................................1 
1.1 Motivation .......................................................................................................1 
1.1.1. Nano and pico-satellite propulsion............................................................1 
1.1.2. Replacement for batteries .........................................................................2 
1.1.3. Micro-engines for UAVs..........................................................................3 
1.2 Challenges .......................................................................................................4 
1.3. Previous Work ................................................................................................7 
1.3.1. Modeling..................................................................................................7 
1.3.2. Experimental work.................................................................................12 
1.4. Objectives and Approach ..............................................................................14 
Chapter 2: The Flame Broadening Phenomenon ................................................17 
2.1. Introduction ..................................................................................................17 
2.2. Analytical Model ..........................................................................................18 
 
 vi 
 
2.2.1. Adiabatic theory.....................................................................................19 
2.2.2. Non-adiabatic theory..............................................................................24 
2.3. Numerical Model ..........................................................................................27 
2.4. Results and Discussion..................................................................................32 
2.4.1. Adiabatic results.....................................................................................32 
2.4.2. Non-adiabatic results..............................................................................43 
2.4.3. Preliminary implications on power density.............................................45 
2.4.4. Pressure loss...........................................................................................49 
2.5. Conclusion....................................................................................................52 
Chapter 3:  Power Density and Efficiency of Micro-Combustors ........................54 
3.1. Introduction ..................................................................................................54 
3.2. Approach ......................................................................................................55 
3.3. Results and Discussion..................................................................................56 
3.3.1. Hydrogen Fuel .......................................................................................56 
3.3.2. Methane Fuel .........................................................................................69 
3.4. Conclusion....................................................................................................81 
Chapter 4:  Micro-Thruster Design and Performance ..........................................83 
4.1. Introduction ..................................................................................................83 
4.2. Model for the micro-nozzle...........................................................................84 
4.3. Results and discussion...................................................................................89 
4.4. Conclusions ..................................................................................................94 
Chapter 5:  Flame Stability .................................................................................96 
5.1. Introduction ..................................................................................................96 
 
 vii 
 
5.2. Analytical Model ..........................................................................................97 
5.3. Numerical Model ........................................................................................107 
5.4. Steady Stability: Results and Discussion .....................................................108 
5.4.1. Analytical model results .......................................................................108 
5.4.2. Numerical model results.......................................................................113 
5.4.3. Comparison to results of other researchers............................................119 
5.4.4. Experimental validation........................................................................122 
5.5. Transient Stability: Results and Discussion .................................................127 
5.5.1 Introduction...........................................................................................127 
5.5.2. Sample calculations..............................................................................128 
5.5.3. Effect of passage height........................................................................131 
5.5.3. Effect of oscillatory period ...................................................................132 
5.6. Conclusions ................................................................................................133 
Chapter 6:   Effect of Flame-Structure Thermal Interaction on Chemistry..........135 
6.1. Introduction ................................................................................................135 
6.2. Alternate definitions of reaction zone thickmess..........................................136 
6.3. Effect of thermal coupling on the kinetics of a multi-step reaction mechanism
..........................................................................................................................142 
6.4. Conclusions ................................................................................................147 
Chapter 7:   Conclusions and Perspectives .........................................................149 
7.1. Contribution and key findings .....................................................................149 
7.2. Future work ................................................................................................150 
Appendix 1............................................................................................................153 
 
 viii 
 
Appendix II ...........................................................................................................158 
II.1. Hydrogen Combustion Mechanism.............................................................158 
II.2. Methane Combustion Mechanism...............................................................159 
Bibliography..........................................................................................................163 
 
 ix 
 
List of Tables 
 
Table 1:  Flow regimes associated with different Knudsen number ranges from.......36 
Table 2: Summary of optimum power density configurations for a silicon micro-
combustor.  The fuel is H2 burning in air and the equivalence ratio is 0.5. .......68 
Table 3:  Major reactions occurring in the flame. ...................................................143 
Table 4: Hydrogen combustion mechanism ...........................................................158 
Table 5: The Leeds reaction mechanism for methane combustion..........................162 
 
 x 
 
List of Figures 
 
 
Figure 1-1: Estimated range of a micro-UAV as a function of energy storage 
efficiency and energy conversion efficiency.......................................................3 
Figure 1-2: Diagram of a simple channel/micro-combustor........................................4 
Figure 1-3: Evolution of the surface to area ratio as a function of passage height in a 
parallel plate reactor...........................................................................................5 
Figure 1-4: Influence of the passage height on the Reynolds number for a 1 m/s flow 
of air between two parallel plates at 300K..........................................................7 
Figure 2-1: Flame stabilized in a passage where H~?
r
..............................................19 
Figure 2-2:  Thermal exchange in a flame in close contact with a structure.  
gas
q&  and 
cond
q&  represent the heat conducted from the reaction zone to the pre-heat zone 
through the gas and through the structure respectively, 
out
q& and 
in
q&  represent 
convective heat transfer between the gas and the structure, and 
loss
q&  represents 
the heat lost to the environment........................................................................20 
Figure 2-3: Thermal resistor network for visualizing heat exchange between gas and 
structure and heat loss to the environment.  The two parallepipeds represent 
streamtubes for heat transfer in the gas and the structure.  
gen
q&  represents the heat 
generated by the combustion process, 
loss
q&  represents the heat lost to the 
environment.....................................................................................................22 
Figure 2-4: Schematic diagram of a flame stabilized in a micro-channel showing heat 
exchange between the reacting fluid and the structure, axial heat exchange 
 
 xi 
 
through the structure, and heat loss to the environment.  The computational grid 
is superimposed on the lower half of the domain..............................................30 
Figure 2-5: Fractional increase in reaction zone length d
r
 as a function of the 
normalized flow passage height H/d
r,fr 
for a range of Nusselt numbers.  The lines 
without symbols correspond to the analytical model and the lines with symbols 
correspond to results from the numerical simulations.  The predictions of the 
simulation and model are no longer valid in the shaded region because the 
Knudsen number is greater than 0.1. ................................................................33 
Figure 2-6: Fractional increase in reaction zone length ?
r
 as a function of the 
normalized flow passage height H/?
r,fr
  for a range of normalized structure 
thermal conductivities ? . Lines without symbols correspond to the analytical 
model and the lines with symbols correspond to numerical model results and 
k
r
=0.21 W/mK. ................................................................................................40 
Figure 2-7: Fractional increase in reaction zone length ?
r
 as a function of non-
dimensional thermal conductance (?? ) of the structure for a range of normalized 
passage heights.  ?
r,fr
 is the reaction zone length associated with a freely 
propagating flame, ? =1, Nu=3, and k
r
=0.21 W/mK.  Lines without symbols 
correspond to the analytical model and lines with symbols correspond to 
numerical model results. ..................................................................................42 
Figure 2-8: Fractional increase in reaction zone length ?
r
 as a function of the 
normalized flow passage height H/?
r,fr
  for a range of heat transfer coefficients in 
W/m
2
K. Lines without symbols correspond to the analytical model and the lines 
with symbols correspond to numerical model results........................................44 
 
 xii 
 
Figure 2-9: Comparison between the evolution of the normalized flame thickness, 
?
r
/?
r,fr
, and the normalized laminar flame speed, S
L
.  Results are from the 1D 
numerical model. .............................................................................................47 
Figure 2-10: Normalized power density as a function of normalized passage height 
for two values of the normalized passage length.  Results are from the 1D 
numerical model. .............................................................................................49 
Figure 2-11: Normalized pressure drop as a function of normalized passage height.  
The evolution of the non-dimensional flame thickness is superimposed on the 
plot for comparison..........................................................................................51 
Figure 3-1: Non-dimensional flashback and blow-off limits for H
2
-air combustion in a 
silicon micro-channel as a function of non-dimensional channel (combustor) 
height and for two values of the channel length L.   The equivalence ratio is 0.5 
and the pressure is 1 atm. .................................................................................57 
Figure 3-2: Axial profiles of temperature and major species concentration for a H
2
-
Air flame stabilized in a silicon micro-channel with 
frr
H
,
? = 0.002 and 
frr
L
,
? = 
9.  The equivalence ratio is 0.5 and the pressure is 1 atm..................................58 
Figure 3-3:  Adiabatic H
2
-air combustion in a silicon micro-channel; ?=0.5, P=1atm.  
Solid lines: Contours of non dimensional power density ( )
refDD
ww
,
&&  as a 
function of non-dimensional channel height and length.  Dashed lines: Contours 
of overall efficiency. ........................................................................................62 
Figure 3-4: Adiabatic (solid lines) and non-adiabatic (dashed lines) H
2
-air combustion 
in a silicon micro-channel; ?=0.5, P=1atm, h
env
=1 W/m
2
K, L/?
r,fr
=20.  The 
circles show non-dimensional laminar flame speed as a function of non-
 
 xiii 
 
dimensional channel height while the squares show non-dimensional reaction 
zone thickness..................................................................................................63 
Figure 3-5: Non-adiabatic H
2
-air combustion in a silicon micro-channel; ?=0.5, 
P=1atm, h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional power 
density ( )
refDD
ww
,
&&  as a function of non-dimensional channel height and length.  
Dashed lines: Contours of overall efficiency. ...................................................64 
Figure 3-6: Non-adiabatic H
2
-air combustion in a silicon micro-channel; ?=0.5, 
P=5atm, h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional power 
density ( )
refDD
ww
,
&&  as a function of non-dimensional channel height and length.  
Dashed Lines: Contours of overall efficiency. ..................................................65 
Figure 3-7: Evolution of power density with increasing pressure for a meso-scale 
combustor and a micro-combustor. ..................................................................67 
Figure 3-8. Evolution of normalized flame speed and non-dimensional temperature 
overshoot as a function of normalized passage height.......................................71 
Figure 3-9: Evolution of normalized flame speed as a function of normalized passage 
height for different values of the normalized passage length.............................72 
Figure 3-10. Adiabatic CH
4
-air combustion in a silicon micro-channel; ?=0.5, 
P=1atm.  Contours of non dimensional power density ( )
refDD
ww
,
&&  as a function 
of non-dimensional channel height and length..................................................73 
Figure 3-11: Non-adiabatic CH
4
-air combustion in a silicon micro-channel; ?=0.5,  
h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional power density 
 
 xiv 
 
( )
refDD
ww
,
&&  as a function of non-dimensional channel height and length.  
Dashed Lines: Contours of overall efficiency. ..................................................75 
Figure 3-12: Comparison of normalized temperature overshoot for methane and 
hydrogen fuels. ................................................................................................76 
Figure 3-13: Comparison of the evolution of flame thickness for hydrogen and 
methane fuel as passage height is reduced........................................................77 
Figure 3-14:  Thermal exchange in a flame in close contact with a structure.  
gas
q&  and 
cond
q&  represent the heat conducted from the reaction zone to the pre-heat zone 
through the gas and through the structure respectively, 
out
q& and 
in
q&  represent 
convective heat transfer between the gas and the structure, and 
loss
q&  represents 
the heat lost to the environment........................................................................78 
Figure 3-15: Evolution of chemical efficiency as the passage length is reduced at 
fixed passage height.  The fuel is methane at equivalence ratio of 0.5...............80 
Figure 4-1: Illustration of the strategy for estimating the performance of a micro-
thruster.............................................................................................................84 
Figure 4-2: Schematic diagram of a model micro-nozzle..........................................85 
Figure 4-3:  Contours of nozzle inlet Mach number as a function of combustor 
configuration....................................................................................................86 
Figure 4-4: Contours of r, the required throat to inlet height ratio to achieve a 
subsonic-supersonic isentropic flow as a function of combustor configuration. 87 
Figure 4-5: Normalized thrust-to-weight ratio as a function of normalized passage 
height and length..............................................................................................90 
 
 xv 
 
Figure 4-6: Solid lines: Specific impulse as a function of normalized passage height 
and length; Dotted lines: Normalized thrust-to-weight ratio as a function of 
normalized passage height and length. .............................................................92 
Figure 4-7: Solid lines: Relative pressure drop (in %) as a function of normalized 
passage height and length; Dotted lines: Normalized thrust-to-weight ratio as a 
function of normalized passage height and length. ...........................................94 
Figure 5-1: Schematic diagram of a parallel plate reactor in which a flame is 
stabilized showing structure and gas temperatures, the three zones described in 
the text and the heat flux from the structure to the preheat zone........................98 
Figure 5-2: Thermal resistance network associated with the flame stabilization 
problem. ..........................................................................................................99 
Figure 5-3: Burning velocity as a function of flame position for different values of the 
thermal conductivity of the structure in W/mK. The families of curves 
correspond to solutions of the thermal network where h
Te
=50 W/m
2
K and U=1.8 
m/s.................................................................................................................104 
Figure 5-4: Burning velocity as a function of flame position for different values of the 
thermal conductivity of the structure. The families of curves correspond to 
solutions of the thermal network where h
Te
=50 W/m
2
K and U=0.8 m/s. .........106 
Figure 5-5:  Flame position as a function of inlet velocity for different values of the 
thermal conductivity of the structure computed using the analytical model.  The 
fuel considered is hydrogen at ? =0.5. ............................................................108 
Figure 5-6: Effect of the thermal conductivity of the structure on the stability limits of 
the flame for H=1mm and L=1cm computed using the analytical model.........110 
 
 xvi 
 
Figure 5-7:  Effect of combustor length on the stability limits of the flame computed 
using the analytical model..............................................................................111 
Figure 5-8:  Effect of combustor height on the stability limits of the flame computed 
using the analytical model..............................................................................112 
Figure 5-9:  Flame position as a function of inlet velocity for different values of the 
thermal conductivity of the structure computed using the numerical model. ...113 
Figure 5-10:  Gas and structure temperature profiles in the fast-burning and slow-
burning regimes identified above. ..................................................................114 
Figure 5-11: Profile of normalized structural temperature gradient for the fast-burning 
and slow-burning regimes. .............................................................................115 
Figure 5-12: Comparison of the results of the analytical and numerical approaches:  
effect of thermal conductivity for H=1cm and L=1cm. ...................................117 
Figure 5-13:  Effect of combustor length on the stability limits of the flame computed 
using the numerical model. ............................................................................118 
Figure 5-14:  Effect of combustor height on the stability limits of the flame computed 
using the numerical model. ............................................................................119 
Figure 5-15:  Flame position as a function of inlet velocity for different values of the 
thermal conductivity of the structure from
15
...................................................120 
Figure 5-16:  Effect of the thermal conductivity of the structure on the stability limits 
of the flame from
15
.........................................................................................121 
Figure 5-17:  Schematic diagram of the silicon-walled micro-burner......................123 
Figure 5-18:  Infrared images of the plates for different values of the inlet velocity of 
the fuel-air mixture.  The reactants are methane and air at ? =0.5. ..................124 
 
 xvii 
 
Figure 5-19:  Dependence of flame location on inlet velocity predicted using 
analytical, numerical and experimental methods. ...........................................126 
Figure 5-20: Schematic diagram illustrating the response of a flame stabilized in a 
micro-channel to perturbations in inlet velocity and equivalence ratio............127 
Figure 5-21:  Transient behavior of H
2
-Air flames (F=0.5) subjected to perturbations 
in inlet velocity.  The upper figure shows the velocity perturbation and the lower 
figure shows the response of the flame for two different passage heights........129 
Figure 5-22:  Transient behavior of H
2
-Air flames subjected to perturbations in 
equivalence ratio.  The upper figure shows the velocity perturbation and the 
lower figure shows the response of the flame for two different passage heights.
......................................................................................................................130 
Figure 5-23: Evolution of normalized maximum flame position with normalized 
perturbation amplitude for different passage heights. .....................................132 
Figure 5-24: Evolution of normalized maximum flame position with normalized 
perturbation amplitude for different oscillation periods. .................................133 
Figure 6-1: H and O radical mass fraction profiles as a function of non-
dimensionalized streamwise coordinate for L=1cm, ? =0.5 and P=1atm.........137 
Figure 6-2: Definition of the flame thickness based on the evolution of a radical 
species mass fraction......................................................................................138 
Figure 6-3: Definition of the flame thickness based on the H
2
O mass fraction........139 
Figure 6-4: Evolution of non-dimensional flame thickness as a function of non-
dimensional passage height for several definitions of the flame thickness. The 
 
 xviii 
 
legend indicates the parameter upon which the flame thickness is computation is 
based..............................................................................................................140 
Figure 6-5: Evolution of chemical efficiency as a function of non-dimensionalized 
passage height for L=1cm and L=0.5cm.  The fuel is H
2
 burning in air at an 
equivalence ratio of 0.5 and a pressure P=1 atm.............................................142 
Figure 6-6: Evolution of rates of progress of the five most important reactions as a 
function of downstream distance for L=5mm and H=1cm. .............................144 
Figure 6-7: Evolution of rates of progress of the five most important reactions as a 
function of downstream distance for L=5mm and H=100?m. .........................144 
Figure 6-8: Comparison of the relative importance of the five most important 
reactions for two different values of the passage height..................................146 
Figure A1-1: Freely propagating laminar premixed flame......................................153 
Figure A1-2: Temperature profile through the freely propagating flame.................154 
Figure A1-3: Thermal circuit for a freely-propagating flame..................................156 
 
 xix 
 
List of Abbreviations 
 
A Cross-sectional area 
a Surface area per unit volume 
a Non-dimensional oscillation amplitude 
C
P 
Heat capacity at constant pressure 
D Diffusion coefficient 
E
a
 Activation energy 
H Channel height 
h Enthalpy 
h
T
 Heat transfer coefficient 
I
sp
 Specific impulse 
Kn Knudsen number 
k Thermal conductivity 
L Channel length 
M Mach number 
Nu Nusselt number 
P Pressure 
R Gas constant 
R Thermal resistance 
RR Reaction rate 
Re Reynolds number 
r Ratio of throat to inlet nozzle heights 
 
 xx 
 
S
L
 Laminar flame speed 
S/V Surface-to-volume ratio 
T
 
Temperature 
T Thrust 
T Oscillation period 
T/W Thrust-to-weight ratio 
t Thickness of the plate 
U, u Velocity 
V Volume 
W Molecular weight 
W Width of the plates 
W
&
 Power generated by the combustion process 
D
w&  Power density 
x Streamwise coordinate 
Y Mass fraction 
 
Greek letters 
?  See Equation (2-4) 
Dp Pressure loss 
?
ph
 Pre-heat zone thickness 
?
r
 Reaction zone thickness 
e Small arbitrary velocity 
? Efficiency 
 
 xxi 
 
?  Equivalence ratio 
g Ratio of specific heats 
?  Ratio of the thermal conductivity of the structure to the thermal conductivity 
of the gas 
m Viscosity 
?
&  Chemical production rate 
 
Subscripts 
ad Adiabatic 
bl Boundary layer 
ch Channel 
e, env Environment 
eff Effective 
f Flame 
i Ignition 
k Quantity associated with species k 
fr Freely propagating value 
m Mean 
r, g Gas 
s Structure 
th Thermal 
0,i Inlet 
o Overall
 
 1 
 
Chapter 1: Introduction 
 
1.1 Motivation 
1.1.1. Nano and pico-satellite propulsion 
Over the past several years, the Aerospace community has become increasingly 
interested in building smaller satellites weighing between 0.1 kg and 10 kg. These 
types of satellites are termed pico (0.1-1kg) and nano-satellites (1-10kg).  Several 
factors are driving the satellite miniaturization trend.  First, the miniaturization of 
electronics and computer systems makes it possible to reduce the weight and the 
volume of the avionics previously required to operate a satellite. Consequently, it is 
now possible to achieve the mission of a normal size satellite using a spacecraft 
whose size and payload weight are much smaller, thereby to lowering the mission 
cost. Second, constellations of many smaller satellites operating cooperatively offer 
increased reliability and enable new missions. Third, small satellites could serve as 
mobile platforms for satellite inspection, refueling, and repair. 
However, the technology associated with these small satellites is in its infancy and 
much more research needs to be conducted in order to be able to fully take advantage 
of the possibilities offered by nano and pico-satellites. One major challenge is the 
development of small, light-weight propulsion systems that are also have adequate 
levels of efficiency.  Existing micro-propulsion systems are very inefficient 
(Mueller
1
). 
 
 2 
 
There are several reasons why chemical propellants are attractive energy storage 
medium for small satellite propulsion systems. First, the mass-specific energy density 
of liquid fuels is at least on the order of 50 MJ/kg, which is relatively high compared 
to, for example, the energy density of batteries which is on the order of 1 MJ/kg. The 
volume-specific energy density of such fuels is also large, which is a great advantage 
for aerospace applications where volume as well as mass is a decisive parameter. 
Second, combustion is a very well established and reliable means for energy release. 
Since it is well understood at regular scale, one can expect that it can be well-
understood at microscale and this should make it very reliable. Third, it has been 
shown
2
 that combustion in small volumes is possible in spite of the greatly increased 
heat loss to the environment. 
   
1.1.2. Replacement for batteries 
Owing to the large energy density of liquid hydrocarbon (HC) fuels compared to 
electro-chemical materials (~100MJ/kg vs. ~1MJ/kg), micro-combustors fueled by 
liquid HCs could replace batteries, provided that the energy in the hydrocarbon fuels 
is converted with at least 1% efficiency.  These power systems could be used in a 
wide range of devices that are presently powered by batteries.  If conversion 
efficiencies greater that 10% are achieved, the operating lifetime of these systems 
could be increased by a factor greater than 10.  The use of a liquid fuel also offers the 
possibility of easy refueling. 
 
 
 3 
 
1.1.3. Micro-engines for UAVs 
The development of micro-UAVs (Unmanned Air Vehicles) depends on the 
availability of efficient and reliable powerplants.  Figure 1-1 illustrates how the range 
of a micro-UAV depends on the method used to store and release energy in the 
vehicle.   
1
10
100
10 100 1000 10000 100000 1000000
Storage Efficiency Q
r
 (W-Hr/Kg)
P
o
w
e
r
p
l
a
nt
 E
f
f
i
ci
ency 
?
pw
r
 (%
)
JP
1
0
Batteries
Model Aircraft Engine 
(150g; 'glow fuel')
DARPA DMFC
(1000 W-Hr/Kg)
Objective
Endurance = 1min
(Range = 1.8 km)
10 min
(18 km)
1000 min
(1800 km)
100 min
(180 km)
10000 min
(18000 km)
?
f
 = 0.45
?
p
 = 0.7
L/D = 8
v = 30 m/ s
20
 
Figure 1-1: Estimated range of a micro-UAV as a function of energy storage 
efficiency and energy conversion efficiency. 
 
The range of the vehicle is computed using the Brequet range equation
3
 
( ()
f
r
po
x
g
Q
D
L
R +
?
?
?
?
?
?
= 1ln?? ) for a fixed-wing UAV with lift-to-drag ratio L/D of 8, 
propulsion efficiency ?
p
 of 70%, ground velocity of 30 m/s, and fuel mass fraction x
f
 
 
 4 
 
of 45%.  The figure suggests that in order to achieve the range objective, an internal 
combustion engine burning a heavy fuel (JP10) with efficiency greater than 20% 
appears to be the only option. 
 
1.2 Challenges 
One major barrier to developing practical H.C. fueled micro-rockets and micro-
power systems is heat loss to the environment. Micro-combustors are very susceptible 
to this because their surface to volume ratio decreases with decreasing combustor 
size. As an example, consider the flow in the simple channel pictured in figure 1-2. 
 
Figure 1-2: Diagram of a simple channel/micro-combustor. 
 
The interior surface area can be expressed as: 
WLS 2=                                                  (1-1) 
 while the interior volume is: 
Direction of flow
W ?  ? 
L 
H 
t
 
 5 
 
HWLV =                                                 (1-2) 
Consequently the surface to area ratio is: 
HV
S
S
v
2
=?                                                (1-3) 
This surface-to-volume ratio increases as the size of the combustor is reduced. 
Figure 1-2 shows the evolution of the surface to volume ratio in such a reactor as a 
function of the passage height H. It shows that as the size of combustors decreases, 
the surface area increases with respect to the enclosed volume. Since the heat loss is 
proportional to surface area while the heat evolved during the chemical reaction is 
proportional to the enclosed volume, figure 1-2 suggests that as the size of 
combustors decreases, heat loss to the environment decreases more slowly than the 
amount of heat released by the combustion process. 
10
-4
10
-3
10
-2
10
-1
10
0
10
0
10
1
10
2
10
3
10
4
10
5
H (m)
S /
 V
 (
m
-1
)
 
Figure 1-3: Evolution of the surface to area ratio as a function of passage 
height in a parallel plate reactor. 
Conventional scale 
Microscale
 
 6 
 
 
When heat loss becomes larger than heat release during combustion, chemical 
reaction cannot be sustained and the flame quenches. 
Another major barrier to building efficient micro-combustors arises from recent 
experimental work
4
 that indicates that flames in micro-combustors may be broader 
than in conventional devices. This, in turn, suggests that the distance needed to 
convert the fuel to products may have to be larger than in conventional combustors 
and that the design and construction of efficient micro-combustors may be more 
complicated than originally anticipated.  The topic of flame broadening will be 
investigated in greater detail in Chapter 2. 
Another challenge arises from the need to achieve high power density.  This 
requires small length scales and large flow velocities which lead to small residence 
times.  However, this also lowers the Damkohler number (i.e. the ratio of residence 
time to chemical time) making it potentially difficult to sustain a combustion wave. 
Finally, the small length scales also make it difficult to use turbulent processes to 
increase the reaction rate and hence the power density.  This is illustrated in figure 1-
4 that shows the influence of the passage height on the Reynolds number for a flow of 
air at 300K in the simple channel of figure 1-1 at a velocity of 1m/s.  For passage 
heights less than 2 cm, the flow is laminar and consequently turbulent enhancement 
of mixing and reaction rate cannot be achieved.  Note that all the combustor 
dimensions considered in this work fall under this 2 cm threshold such that turbulence 
does not come into play. 
 
 
 7 
 
10
-4
10
-3
10
-2
10
-1
10
0
10
1
10
2
10
3
10
4
10
5
10
6
H (m)
Re
Turbulent transition (Re~3000)
2 cm 
LAMINAR TURBULENT 
 
Figure 1-4: Influence of the passage height on the Reynolds number for a 1 
m/s flow of air between two parallel plates at 300K.  
 
1.3. Previous Work 
1.3.1. Modeling 
The problem of the effect of the presence of a structure in close proximity to a 
flame has been addressed by many researchers. One of the first was Zeldovich
5
 in 
1941. He studied the propagation of flames in a series of circular tubes and showed 
that there was a minimum value for the radius of the tube below which the flame was 
unable to propagate. He introduced the concept of thermal quenching to explain this 
 
 8 
 
observation of a minimum tube size for flame propagation.  It was based on the idea 
that if the heat loss to the environment is larger than the heat generated by 
combustion, the combustion cannot continue and therefore the flame cannot 
propagate. 
In 1956, Spalding
6
 developed an analytical theory to explain Zeldovich?s 
observations.  He considered the problem of a flame propagating non-adiabatically, 
with heat loss to the environment accounted for using a single (fixed) loss term in the 
gas-phase energy equation.  One of the most important results of his analysis was the 
fact that two flame speeds are associated with the reactive mixture. The lower one is 
unstable and the higher one is always less than the freely propagating flame speed.  
As the magnitude of heat loss increases, the stable and unstable flame speeds 
approach each other and eventually coincide at the inflammability limit of the 
mixture.  While Spalding?s model was simplistic, the results and behaviors which he 
predicted are qualitatively consistent with modern experimental measurements in 
micro-combustors.  Nevertheless, an improved theory is required that accounts for 
heat exchange within the combustor structure and the scale and geometry 
dependencies of the heat loss term.   
While the thermal interaction between flame and structure can result in heat loss to 
the environment thereby lowering the flame speed of the mixture and degrading 
combustor perfomance, it can also be beneficial by providing an additional path for 
heat transfer from post-flame to pre-flame via the structure.  This concept was first 
recognized and investigated by Weinberg
7
 who introduced the term ?heat 
recirculating burner? to describe the latter situation.  Weinberg?s work was 
subsequently followed by Takeno
8
, who performed detailed theoretical and 
 
 9 
 
experimental studies of the effect of heat recirculation.  One of his important 
observations was that heat recirculation could result in gas temperatures that are 
higher than the adiabatic flame temperature.  These have been termed ?excess 
enthalpy? flames and the porous media burners that produce them are often called 
excess enthalpy burners. 
More recently, researchers have investigated the influence of thermal exchange in 
the structure on the thermodynamic performance of miniaturized power systems. In 
1998, Peterson
9
 identified the minimum sizes for MEMS-based regenerative heat 
engines from a purely thermodynamic viewpoint. Using models for conduction within 
the structure of a generic Stirling engine and the heat loss to the environment, he was 
able to determine the efficiency of the heat engine as a function of the cycle 
temperature ratio, regenerator effectiveness and a non-dimensional ?conduction 
parameter? that represented the ratio of heat lost to the environment via conduction 
through the structure to the heat released by combustion in the engine. He showed 
that as the conduction parameter increases, the size of the heat engine would need to 
increase in order to achieve a ?reasonable? level of efficiency. The net result was that 
for given thermal properties of the structure of a combustor, a heat engine needs to 
have a certain minimum size to achieve acceptable operating efficiency. However, he 
did not incorporate any type of model for the combustion process and he neglected 
thermal exchange in the gas. 
Aichlmayer
10
 performed a more detailed investigation that included more realistic 
thermal coupling between the gas and the structure.  He did so in the context of a 
miniature free-piston Homogeneous Charge Compression-Ignition (HCCI) engine 
where a premixed fuel-air mixture is compressed until ignition. HCCI engines differ 
 
 10 
 
from diesel engines in that the fuel and air are premixed before entering the cylinder. 
They differ from spark-ignited engines because the mixture auto-ignites under 
compression and therefore requires no spark. In an HCCI engine, the mixture burns 
rapidly and homogeneously throughout the combustor volume. As a result, one would 
expect wall quenching to be less important than in regular engines and that HCCI 
engines would be a very advantageous operating mode for micro-combustors. 
Aichlmayer?s calculation included full chemistry and also showed that there is a 
minimum size for HCCI engines that arises from the heat transfer to the structure. 
Reducing the aspect ratio of the engine improved efficiency by reducing the surface-
to-volume ratio and consequently heat loss to the environment. However his model 
did not incorporate the effect of thermal exchange (i.e. conduction) within the 
structure itself. 
In 2002
11
, Daou and Matalon investigated combustion in a 2-D rectangular channel 
with fully developed Poiseuille flow and heat loss to the walls occurring 
perpendicularly to the flow direction. They focused on the effect of the heat loss to 
the channel walls and the inlet velocity on the flame speed and the burning rate of the 
fuel. They showed that as the passage height becomes smaller, the burning rate 
decreases because of heat loss to the walls. They also identified two modes of flame 
extinction: total extinction in small channels and local extinction at the walls in larger 
channels. In 2003
12
, Daou and Matalon extended their work by studying the detailed 
structure of the flame. They showed that the wall quenching zone enlarged as heat 
transfer to the wall increased. However, both works relied on a simplified single step 
model for the chemical reaction and did not consider any axial (i.e. in the flow 
direction) conduction within the combustor walls. 
 
 11 
 
In 2004, Ronney
13
 developed a model for a ?Swiss roll? heat recirculating burner 
where heat from the post flame region is transferred upstream to pre-heat the 
incoming reactants using a counter-flow heat exchanger. The combustion process was 
modeled as a perfectly stirred reactor (PSR).  This work showed that heat conduction 
within the structure of the combustor in the flow direction is a dominating effect at 
microscale. 
Recent work by Ju
14
 in 2004 also highlighted the importance of heat conduction in 
the combustor structure.  He extended Spalding?s 1D theory to include the effect of 
heat recirculation within the structure of a micro-combustor in addiction to heat loss 
to the environment.  Like Spalding, Ju observed two distinct flame speeds where the 
lower one was unstable.  However, one major difference was the possibility for the 
higher flame speed to be in excess of the laminar flame speed due to heat 
recirculation.  
Also in 2004, Vlachos
15
 reported results of a 2D numerical simulation showing 
how axial heat conduction within the structure of a micro-channel can enhance flame 
stability.  The work was motivated by earlier experimental observation of ?Thermally 
Stabilized Combustion? by Churchill
16
.  Vlachos also showed that the enhancement in 
stability is strongly affected by the thermal conductivity of the structure as well as 
heat loss to the environment.  Vlachos showed that there is an optimum value of the 
structural thermal conductivity that maximizes flame stability with respect to heat 
loss to the environment. 
Relatively little modeling of flows within micro-thrusters has been reported in the 
open literature where the effect of heat transfer to the thruster structure are 
incorporated.  Levin
17
 performed a study of the thermal coupling between reacting gas 
 
 12 
 
and structure in such micro-thrusters at low pressure.  The 2D gas-phase flow was 
resolved using DSMC while the structure temperature was computed using a Finite 
Element Method.  Levin?s work was motivated by the idea that high temperature 
flows in microscale systems cause the structure to reach high temperature which 
could in turns degrade its mechanical strength.  Levin?s work focused on the transient 
heating of the structure and how it affects the maximum impulse of micro-thrusters.   
Other simulations of flows in micro-nozzles performed by Ketsdever
18
 focused on 
how friction and viscous losses affect performance.  This work aimed to determine 
whether a simple orifice configuration was more efficient than a contoured nozzle at 
micro-scale.  Ketsdever?s conclusion was that the nozzle to orifice thrust ratio 
depends on the Reynolds number and is less than one at low Re.  However, nozzle to 
orifice specific impulse ratio was always greater than one, indicating that contoured 
nozzles provide better propulsive efficiency. 
 
1.3.2. Experimental work 
There has been a relatively large amount of work involving the design and 
construction of micro-combustors and micro-power systems that has met varying 
levels of success. 
The MIT micro-engine project
19
 is the oldest of the combustion based micro-scale 
power generation projects.  The basic idea is to construct a shirt-button size gas 
turbine engine that could be attached to a lunchbox-size fuel tank (full of JP4) and be 
used to generate electric power for soldier-bourne navigation and communication 
gear.  After more than ten years of effort, operation of the different elements of the 
 
 13 
 
micro-engine have been demonstrated in isolation.  However very substantial 
challenges remain in order to integrate those elements into a working (and suitably 
efficient) system.  One important outcome of the research on the micro-gas turbine 
engine combustor was the demonstration that despite the challenges highlighted 
above, combustion in millimeter-scale silicon-based devices is possible (Mehra
20
).  
Additional work by Spadaccini
21
 characterized the performance of the combustor and 
showed that catalytic micro-combustion of hydrocarbons was also possible.  
However, the overall efficiency of the micro-engine combustor is very low (<60% on 
HC fuels) and this alone could prevent the engine?s thermodynamic cycle from 
closing. 
The Swiss-roll combustor discussed previously is also part of an experimental 
program to develop a compact power system based on catalytic combustion and 
thermo-electrical generation.  The Swiss-roll geometry minimizes heat loss to the 
environment while providing the opportunity for super-adiabatic flame temperatures 
(excess enthalpy) via heat recirculation through the walls of the combustor.  The 
catalyst allows for combustion to be sustained at very low equivalence ratio, thereby 
limiting structure temperature.  The use of a catalyst deposited on the combustor 
walls also capitalizes on the aforementioned high surface-to-volume typical of micro-
combustors.   Thermo-electric materials separate the hot and cold sides of the counter 
current heat exchanger and are used to generate electric power.  Preliminary results
22
 
suggest that the approach is viable. 
Another interesting approach is that taken by Masel?s group at UIUC.  The focus 
of this work
23
 is on surface treatment and how it can enable combustion in very small 
 
 14 
 
channels by preventing processes associated with flame quenching such as radical 
recombination at the wall.  The results obtained demonstrate the possibility of sub-
millimeter combustion and suggest the importance of surface properties and 
preparation on the design of efficient micro-combustors. 
The first MEMS bi-propellant micro-thruster was developed by London
24
 at MIT, 
in 2000.  The performance of the thruster was surprisingly good despite the very 
small scale of the thruster in the transverse direction (Thrust~1N, I
sp
~180s).  Pressure 
losses through the device were smaller than expected because of the non-continuum 
effects (i.e. velocity slip at the wall). 
Hitt?s group
25
 at the University of Vermont has demonstrated a MEMS-based 
H
2
O
2
 thruster.  Hydrogen peroxide was injected in a sub-millimeter scale combustion 
chamber filled with platinum-coated rods.  Reaction occurred in the combustion 
chamber and hot products were expanded through a nozzle.  The device produced 
500?N of thrust with a specific impulse of up to 180s. 
Finally, Ketsdever?s group
26
 performed an experimental study of the supersonic 
flow in micro-orifices and nozzles and validated the results obtained during the 
computational phase of its investigation. 
 
1.4. Objectives and Approach 
The review of the literature has shown that axial heat transfer through a micro-
combustor?s structure has a significant impact on its overall performance and viability 
as part of a micro-chemical thruster or a compact power system.  Therefore, the 
principal objective of this thesis is to develop a fundamental understanding of the role 
 
 15 
 
of axial heat transfer on the performance of micro-scale combustors for these 
applications.  This is accomplished through a combined program of analytical 
modeling and numerical simulation that focuses on the propagation of a laminar 
premixed flame in a rectangular channel where the flame exchanges heat with the 
structure, heat flows within the structure itself, and heat is lost to the environment.  
This ?model problem? has the advantage of being simple while at the same time 
incorporating most of the critical aspects of the problem.  As a result, it can represent 
combustion in a wide range of micro-power systems. 
The approach taken is as follows:  First, a simple analytical model for flame 
propagation in a micro-channel is developed in order to gain a basic understanding of 
the physics that govern the problem and lead to the unusual observations of flame 
broadening and increased burning velocity at small scales.  These results are 
compared to the predictions of more detailed one-dimensional numerical simulations 
that include conjugate heat transfer, heat loss to the environment, and multi-step 
chemistry.  Second, the results of these calculations are used to identify the effects of 
varying the combustor geometry and the degree of heat loss to the environment on the 
power density and efficiency of a micro-combustor.  Third, these results are used to 
predict how the performance of a micro-chemical thruster varies with its geometry 
and heat losses to the environment.  Fourth, flame stabilization due to axial heat 
transfer is investigated.  Fifth, and finally, the implications of axial heat transfer on 
the operation and selection of multi-step chemical kinetic mechanisms used to 
represent the combustion process in micro-channels are discussed. 
 
 
 16 
 
 
 
 
 17 
 
Chapter 2: The Flame Broadening Phenomenon 
 
 
2.1. Introduction 
The work reported in this chapter focuses the challenge identified in the previous 
chapter that is associated with flame broadening in micro-combustors. The objective 
is to identify the basic governing physics of the phenomenon and to use this 
understanding to make recommendations regarding the design of efficient micro-
combustors.  To this end, a mechanism is proposed by which flame broadening 
occurs via the exchange of heat from the post-flame region to the pre-heat region 
through the structure of the combustor.  
A simple analytical model is developed to describe the basic physics of the fluid-
structure coupling problem. The derivation of the model employs reasoning similar to 
that used by Mallard and Le Chatelier to predict the velocity of a freely propagating 
flame but extends it to include the effect of the thermal exchange within the structure 
of the combustor. Both adiabatic and non-adiabatic situations are considered. The 
model is used to understand the basic physics of the problem by predicting the 
variation of the flame thickness with different the design parameters that characterize 
the combustor.  These include the passage height, the thermal conductivity, the 
Nusselt number for heat transfer between the reacting flow and the structure, and the 
heat transfer coefficient between structure and the environment.  
 
 18 
 
The predictions of the analytical model are subsequently compared to those of a 
more detailed one-dimensional numerical simulation with full chemistry, species 
diffusion and heat transfer in the structure. 
Preliminary work in the area was reported by Leach
27
 but significant 
improvements in the analytical model and numerical simulation, as well as an 
additional discussion of power density and pressure loss are presented here.  Some of 
the results presented here are in press
28
 
 
2.2. Analytical Model 
Mallard and Le Chatelier?s thermal model for flame propagation
29
 predicts that the 
laminar flame speed S
L
 is proportional to the square root of the product of the thermal 
diffusivity (
p
r
C
k
?
) of the unburned reactants and the chemical reaction rate RR:  
RR
T
i
T
i
T
f
T
p
C
r
k
L
S
0
?
?
=
?
                                                (2-1) 
T
o
, T
i
 and T
f
 are respectively the temperatures of the unburned mixture, the 
temperature at ignition, and the temperature of the post-flame gases.  The theory 
applies to free-burning (i.e. unconfined) laminar flames and is predicated on the 
notion that the flame speed is set by the rate at which heat diffuses upstream into the 
unburned mixture.  Appendix I provides a theoretical overview of Mallard and Le 
Chatelier?s thermal model for flame propagation, and an explanation of how a 
thermo-electrical analogy can be used to derive equation (2-1). 
 
 
 19 
 
2.2.1. Adiabatic theory 
The thermal theory of Mallard and Le Chatelier applies to freely propagating 
flame that does not exchange heat with the combustor structure or the surroundings.  
These effects are expected to become important when the flame thickness approaches 
the characteristic size of the combustion chamber as illustrated in figure 2-1.   A more 
complex picture emerges when the theory is extended to include an additional heat 
exchange path through the structure from the post-flame region upstream to the pre-
flame region.  This is illustrated in figure 2-2 which shows that two additional heat 
transfer paths emerge in addition to conduction upstream through the gas ? namely   
through the structure from the post-flame to the pre-flame, and from the structure to 
the environment.  The path from the post-flame to the pre-flame has three elements: a 
path from the post-flame to the structure, a path through the structure, and a path from 
the structure to the pre-flame. 
 
 
S
L
 
?
r
H 
CV 
 
Figure 2-1: Flame stabilized in a passage where H~?
r
.  
 
 20 
 
 
Reaction
T 
 
loss
q
&
in
q
&
gas
q
&
x 
T
i
T
f
Pre-Heat 
T
0 
?
ph
?
r
out
q
&
cond
q&
Structure
 
Figure 2-2:  Thermal exchange in a flame in close contact with a structure.  
gas
q&  and 
cond
q&  represent the heat conducted from the reaction zone to the 
pre-heat zone through the gas and through the structure respectively, 
out
q& and 
in
q&  represent convective heat transfer between the gas and the 
structure, and 
loss
q&  represents the heat lost to the environment. 
 
The thermal resistor network of figure 2-3 can be used to visualize the heat 
exchange between the reaction and the preheat zones as well as the heat loss to the 
environment. In the figure, 
gen
q&  is the rate of heat generation by the reaction, 
loss
q&
 
is 
the rate at which heat is lost to the environment by convective heat transfer from the 
plates, T
e
 is the environment temperature, T
s
 is the structure temperature and T
f
 is the 
flame temperature.  R
r
 and R
s
 are the resistances to conductive heat transfer in the gas 
 
 21 
 
and the structure respectively and R
bl
 is the resistance to convective heat transfer from 
the gas to the structure. These thermal resistances can be expressed as follows: 
convr
bl
ss
r
s
rr
r
r
NuAk
H
R
Ak
R
Ak
R
2
,, ===
??
                              (2-2) 
where k
r
 and k
s
 are the thermal conductivity of the gas and of the structure 
respectively, ?
r
 is the reaction zone thickness, H is the passage height, Nu is the 
Nusselt number and A
r
, A
s 
and A
conv
 are the cross-sectional surface areas for 
conduction in the gas and the structure, and for convection between gas and structure 
respectively.  A
conv
 is assumed to be based on the flame thickness ?
r
.   
Finally, it should be noted that the Mallard and Le Chatelier theory for freely-
propagating laminar flames does not include a finite-dimensional pre-heat zone 
upstream of the reaction zone and so a model for heat transfer in the pre-heat zone is 
not included in the extension of the theory to micro-combustion.  This is not intended, 
however, to imply that modeling the pre-heat zone is not important.  It is simply to 
say that modeling this mechanism is beyond the scope of the simplified theory for 
micro-combustors presented here.  As we will see, this omission does not prevent the 
model from producing qualitatively correct results.  
 
 22 
 
 
?
r
R
bl
 R
bl
R
s
R
r
T
i
 T
o
 
r
 
T
s 
T
e
gen
q
&
T
f
Structure
Gas
 
Figure 2-3: Thermal resistor network for visualizing heat exchange between 
gas and structure and heat loss to the environment.  The two parallepipeds 
represent streamtubes for heat transfer in the gas and the structure.  
gen
q&  
represents the heat generated by the combustion process, 
loss
q&  represents the 
heat lost to the environment. 
 
In the adiabatic case (corresponding to a tube that is perfectly insulated from the 
environment), 
loss
q&  is set to 0.  Since the flame temperature is a constant equal to the 
adiabatic flame temperature, the rate at which heat is generated by the flame (
gen
q& ) 
does not need to be defined. Applying conservation of energy (the equivalent of 
Kirchoff?s law) to the thermal circuit illustrated in figure 2-3 and assuming that 
RRS
Lr
=?  as is done in the Mallard-Le Chatelier theory leads to the following 
expression for the reaction zone thickness: 
 
 23 
 
 
RRTT
TT
C
k
i
if
p
r
r
1
0
?
?
=
?
??                                           (2-3) 
Beta (? ) is a parameter that accounts for the variation of the overall thermal 
resistance of the gas-structure combination as a function of the passage height H,  the 
ratio ? of the cross sectional areas for axial conduction in the structure and the 
reactant stream-tube 
?
?
?
?
?
?
=
r
s
A
A
?  and the ratio ?  of the thermal conductivity of the 
structure to that of the reactive mixture 
?
?
?
?
?
?
=
r
s
k
k
? . 
Nu
H
Nu
H
r
r
1
21
12
11
2
2
2
2
?
??
?
??
?
+
?
?
?
?
?
?
?
?
++
=                                           (2-4) 
Combining equations 2-3 and 2-4 yields a fourth order polynomial equation in ?
r
 for 
which four solutions can be found.  Of these, only one value of ?
r
 is real positive and 
hence physically realistic.  For a rectangular channel, 
H
t2
=?  where t is the 
thickness of the combustor wall.  Combining equations 2-3 and 2-4 to solve for ?
r
 
results in a fourth order algebraic equation with four possible roots.  Only one of 
these is real and positive and hence physical. 
Equations 2-3 and 2-4 indicate that the reaction zone thickness depends on the 
normalized thermal conductivity of the structure, the area ratio for axial conduction ?, 
the passage height H, and the Nusselt number in addition to the usual parameters 
associated with the Mallard-Le Chatelier theory.  While equation 2-3 is an implicit 
expression for ?
r
 because ? is a function of ?
r
, we can check that it has the proper 
 
 24 
 
asymptotic behavior using equations 2-3 and 2-4. At zero Nusselt number or infinite 
passage height, equation 2-3 indicates that 1?? .  Inserting this result into equation 
2-4 recovers the traditional formulation in which the reaction zone thickness only 
depends on the thermal diffusivity of the gas.  Conversely, at large Nusselt numbers 
or small passage heights, equation 2-3 indicates that ??? +? 1 , which is always 
greater than one. This indicates that the reaction zone thickness increases, leading to 
flame broadening.  At large Nusselt numbers, the thermal resistance between the gas 
and the structure is small and the transport of heat (and hence the reaction zone 
length) is determined by the thermal resistance of the structure acting in parallel with 
the thermal resistance of the gas.  As a result, in either the large Nusselt number or 
the small passage height limits, the reaction zone thickness is determined by the 
thermal conductivity of the structure and the gas, and the relative areas for axial heat 
conduction.  
Thus, the simple model has shown that a conductive heat transfer conductive 
structure can increase heat transfer from the reaction zone to the upstream reactants.  
This can also lead to an increase in burning rate due to faster upstream conduction 
and pre-heating of the reactants.  This is consistent with the effects of pre-heating 
reactants in laminar flames
30
, and the pre-heating mechanism is very similar to that of 
?excess enthalpy? burners
31,32
 where increases in burning rate are also observed. 
 
2.2.2. Non-adiabatic theory 
The effect of heat loss to the environment, an important consideration in micro-
combustors because of their large surface/volume ratios, can be incorporated by 
 
 25 
 
adding a heat sink to the thermal network and by recognizing that the chemical 
reaction acts like a heat  (or analogous current) source.  In this case, 
gen
q&  is computed 
by adjusting the rate of heat generation in the adiabatic case 
ad
q&  by a factor that 
accounts for the effect of an arbitrary flame temperature T
f
 
?
?
?
?
?
?
?
?
+?=
adf
a
f
a
adgen
RT
E
RT
E
qq
,
exp&&                                    (2-5) 
In this expression, E
a
 is the activation energy, R is the gas constant and T
f,ad
 is the 
adiabatic flame temperature. The value for E
a
 corresponds to a one-step overall 
reaction mechanism for stoichiometric H
2
-air combustion and was taken from
33
. 
loss
q&
 
is expressed as a convective heat loss: 
()
eslossTeloss
TTAhq ?=&                                             (2-6)               
where h
T,e
 is the coefficient for heat transfer from the plates to the environment and 
A
loss
 is the surface area for convective heat loss. Since heat loss occurs mainly in the 
region where the gas temperature is high, the area for heat loss is assumed to 
correspond to a length equal to the flame thickness. The heat transferred from the 
reaction zone to the preheat zone is written as the difference between the heat 
generated by the reaction and the heat lost to the environment: 
lossgen
qqq &&& ?=                                                  (2-7) 
The value of q
&
 obtained using equation 2-7 is used to compute an ?effective? flame 
temperature which includes the effect of heat loss. 
effif
RqTT &+=                                              (2-8) 
 
 26 
 
In this expression, T
i
 is the ignition temperature and R
eff
 is the effective thermal 
resistance of the thermal network defined as: 
111
)2(
???
++=
blsreff
RRRR                                  (2-9) 
Balancing the heat fluxes through the boundary layer, structure, and to the 
environment leads to the following expression for the peak temperature of the 
structure T
s
: 
()
()
blsr
blsr
loss
bls
bls
ifis
RRR
RRR
q
RR
RR
TTTT
++
+
?
+
+
?+= &
2
       (2-10) 
As before, performing a thermal balance on the entire network gives an expression for 
the flame thickness: 
()RRTTC
q
ip
r
0
?
=
?
?
&
                                       (2-11) 
The reaction rate in equation 2-11 is the adiabatic reaction rate adjusted by a factor 
accounting for the ?effective? flame temperature T
f
  
?
?
?
?
?
?
?
?
+?=
adf
a
f
a
ad
RT
E
RT
E
RRRR
,
exp                          (2-12) 
The system of equations {2-5 - 2-12} can be solved using an iterative procedure and 
an initial guess that corresponds to the adiabatic case: 
?
?
?
=
=
ad
adff
qq
TT
&&
,
                                        (2-13) 
 
 27 
 
 
2.3. Numerical Model 
A transient 1-D model developed by Zhu and Jackson for studying the 
performance of catalytic reactors is adapted to study the behavior of a flame 
stabilized in a channel.  A complete description of the code may be found in
34
  and 
only an abbreviated description is presented here.  For the numerical model in this 
study, all surface chemistry is excluded from the calculations and as such there is no 
mass exchange between the surface and the gas phase.  Further, the pressure in the 
channel flow is assumed constant, which implies that the momentum equation is 
neglected entirely and velocity is calculated from mass conservation.  While this is 
appropriate for the purposes of comparison with the analytical model, the pressure 
drop in micro-channels can be a significant performance parameter because of its 
impact on the efficiency of thermodynamic cycles.  In this case, however, we are 
interested in minimum size combustors with lengths on the order of the flame 
thickness.  As a result, pressure losses are not expected to be significant.  This will be 
further discussed in the ?Results and Discussion? part of this chapter. 
The numerical model divides the channel into axial cells and solves for the 
following variables in each cell: ?, u, T, T
s
, and Y
k
 where Y
k
 is the vector of mass 
fractions for the species in the gas phase. The governing equations in each cell for 
these variables include the following: 
? Ideal gas law: 
 
RT
PW
m
=?                                             (2-14) 
 
 28 
 
? Conservation of mass: 
()
0=
?
?
x
u?
                                              (2-15) 
? Conservation of energy in the gas-phase: 
()
()( )
?
?
???
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
?
?
?
?
?
?
?
?
?
?
?=
?
?
k
kkksTchr
k
kk
P
hW TTha
x
T
k
xx
Yh
u
t
T
C ???
&
(2-16) 
? Conservation of species in the gas-phase: 
kk
kk
k
k
W
x
Y
u
x
Y
D
xt
Y
???? &+
?
?
??
?
?
?
?
?
?
?
?
?
=
?
?
                  (2-17) 
? Conservation of energy in the solid structure: 
() ( )
eseTessTs
s
s
s
sPs
TThaTTha
x
T
k
xt
T
C ???+
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
?
?
,,,
?       (2-18) 
In the above equations, 
m
W is the mean molecular weight of the gas flow, 
k
?
&
 is the 
production rate of the k
th
 species in the gas phase. a
ch
 represents the internal structural 
surface area per unit gas phase flow volume, and a
s
 and a
s,e
 represent the internal and 
external structural surface area per unit structure volume respectively.  Since we 
consider a parallel plate geometry these last parameters simplify to a
ch
=2/H and a
s
= 
a
s,e
 = 1/t. 
For the numerical model, a 9-species, 19-reaction mechanism
35
 for H
2
-air 
combustion is employed.  Table 4 in Appendix II shows the reactions considered in 
the mechanism.  Properties for the gas-phase species varied with temperature and 
were computed based on fitted polynomials (h
k
 and C
P
) and on kinetic theory 
 
 29 
 
approximations (k
r
 and D
k
).  The properties for the solid structure (?
s
, C
Ps
 and k
s
) did 
not vary with temperature in the numerical model. 
 
 
The boundary conditions for equations 2-15 ? 2-18 in the numerical model are as 
follows: 
? at the inlet (x = x
in
,): 
() ()
()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=?
?=?
=
=
x
T
kTTh
TThTTCu
YY
uu
s
ssinT
sinTinP
inkk
in
,
,
,
??
                 (2-19) 
? and at the oulet (x = x
out
,): 
?
?
?
?
?
?
?
?
?
=
?
?
=
?
?
=
?
?
0
0
0
x
T
x
T
x
Y
s
k
                                  (2-20) 
For the inlet boundary at x = x
in
, h
T,in
 represents the effective heat transfer coefficient 
to the upstream portion of the gas flow, and for the current study h
T,in
 = 0 such that
 
T  
= T
in
 for the first cell boundary.  For all runs in the current study, T
in
 = 300 K, and 
Y
k,in
 are set for a H
2
-air mixture with ? =  1.0.  Pressure was held constant at 1.0 atm 
for all simulations.  The value for u
in
 is determined during a model run as discussed 
below to ensure that the flame stabilization occurred away from the reactor inlet and 
thereby flame thicknesses were found independent of the inlet boundary conditions. 
 
 30 
 
 
Figure 2-4: Schematic diagram of a flame stabilized in a micro-channel 
showing heat exchange between the reacting fluid and the structure, axial 
heat exchange through the structure, and heat loss to the environment.  The 
computational grid is superimposed on the lower half of the domain. 
 
A discretized version of equations {2-14 ? 2-18} is solved simultaneously on an 
evenly spaced grid consisting of 201 elements with a domain length of 1.0 cm and the 
boundary conditions shown in equations 2-19 ? 2-20.  Adequate grid resolution has 
been verified by doubling the number of grid points and finding that flame 
thicknesses did not change from the 201-cell simulation.  The 1.0 cm computational 
domain was sufficient to contain the reaction zone for all runs presented here.  The 
code uses the method of lines to implement an adaptive time step.  The convective 
terms are discretized using a Godunov-type upwind finite volume method while a 
Nu 
T
e
d
r
 
t
U, P, T
o
 
L
h
env
=1 W/m
2
K 
H 
Computational 
grid 
?  =2t/H=1 
 
 31 
 
central difference method was used to discretize the diffusive terms.  The resulting set 
of discrete equations was integrated using the stiff PDE solver LIMEX
36
. 
The numerical solution depends strongly on user-specified parameters including 
H, Nu, k
s
 and h
T,e 
, and for a given combination of user-defined parameters, the 
necessary u
in
 (? S
L
) to stabilize the flame in the domain away from the inlet boundary 
is unknown.  Conventional flame models have used the approach of solving the 
steady-state version of equations {2-14 ? 2-18}, fixing a temperature within the 
computational domain, and solving for ?u as an eigenvalue for the system
37
.  
However, in the current study, it is desired to model conditions when the flame 
undergoes extinction, and thus rather than anchoring the flame with a fixed 
temperature in the mesh, an alternative approach based on the transient integration of 
equations {2-14 ? 2-18} is implemented.  An initial condition is set with T and T
s
 = 
1000 K throughout the domain except the first gas point of the domain where T =300 
K.  The high initial values of T and T
s
 in the rest of the domain ensure that the flow 
?lights?.  The inlet velocity u
in
 is set to a value below the expected S
L
.   A steady 
solution is found typically after 10 s of transient integration and the flame is anchored 
at the inlet of the reactor.  u
in
 is then increased by 0.1% and the solution is 
recalculated.  This process is repeated until the induction zone length (based on the 
decay of the O
2
 mass fraction in the flame zone) equals that calculated for the freely 
propagating flame.  The value for the laminar flame speed associated with a freely 
propagating flame was found in
38
 and equaled 2.1 m/s at the equivalence ratio 
considered here.  At the final u
in
, the integration continues for 100 s to ensure that the 
 
 32 
 
final solution is at a steady-state condition for both the solid and gas phases and that 
the flame position is fixed in relation to the walls. 
Upon determining the final solution at a value of u
in
 that provides a reaction zone 
away from the inlet, the flame thickness can be determined.  The definition of flame 
thickness proposed by Sung and Law
39
 is used where the flame thickness is the ratio 
of the total temperature difference to the maximum value of the temperature 
derivative: 
max
0
?
?
?
?
?
?
?
=
dx
dT
TT
f
r
?                                                  (1-21) 
For H
2
-air mixtures with ?  = 1.0 and T
0
 = 300 K, the flame thickness for a freely 
propagating flame computed using equation 1-21 is 0.42 mm and is consistent with 
values reported elsewhere
40
.  For the same flame, the total temperature change across 
the flame is T
f
 ?T
0
 = 2036 K.  This is 61 K lower than the difference associated with 
complete combustion by an adiabatic flame.  This reduced temperature difference is 
the result of dissociation. 
 
2.4. Results and Discussion 
2.4.1. Adiabatic results 
Figure 2-4 is a plot of the non-dimensional flame broadening (i.e. the ratio of the 
reaction zone of thickness ?
r
 to the reaction zone thickness associated with a freely 
propagating flame in the absence of the structure ?
r,fr
) as a function of the normalized 
passage height H/?
r,fr
 for four values of Nu under adiabatic structure conditions (i.e., 
 
 33 
 
h
T,e
 = 0).  The value of ?
r,fr
 found using the numerical simulation is 0.42 mm.  The 
Nusselt number is a constant for fully-developed flows that depends on the geometry 
of the flow passage.  Therefore, curves corresponding to different values of the 
Nusselt number are included to illustrate the effects of changing the thermal coupling 
between the gas and the structure by changing the flow geometry.   
H/?
r,fr
Kn
?
r
/
?
r,
f
r
10
-2
10
-1
10
0
10
1
10
2
10
-5
10
-4
10
-3
10
-2
10
-1
0
1
2
3
4
5
6
7
8
Nu=0
Nu=3
Nu=5
Nu=10
Nu=0
Nu=3
Nu=5
Nu=10
Model
?=48, ?=1
Simulation
 
Figure 2-5: Fractional increase in reaction zone length d
r
 as a function of the 
normalized flow passage height H/d
r,fr 
for a range of Nusselt numbers.  The 
lines without symbols correspond to the analytical model and the lines with 
symbols correspond to results from the numerical simulations.  The 
predictions of the simulation and model are no longer valid in the shaded 
region because the Knudsen number is greater than 0.1. 
 
 34 
 
 
The lines without symbols correspond to the simple analytical model while the 
dotted lines with symbols correspond to the numerical simulations. The analytical and 
numerical models have good qualitative agreement as both indicate that ?
r
 undergoes 
a limit cycle behavior with variations in H.  At large H, ?
r
/?
r,fr
 = 1 and ?
r
 is 
determined by k
r
 (~0.2 W/mK).  As H decreases, ?
r
 becomes more strongly 
influenced by the structure thermal conductivity k
s
 (10 W/mK).  In both models, the 
strength of the thermal coupling between the fluid and the structure (determined in 
part by Nu) determines at what H the transition from gas to structure-dominated heat 
transfer occurs.  The transition begins to occur as H decreases below 10?
r,fr
 (~ 5 mm) 
and  the  broadening  has  reached  its  maximum limit  when H < 0.01?
r,fr
 the 
unconfined flame thickness.  Many factors contribute to the discrepancy in the 
amplitude of the broadening predicted by the numerical and analytical models.  Some 
important differences include the use of multi-step versus single step chemistry, 
variable versus fixed gas-phase properties, and the use of discrete versus distributed 
heat transfer between the gas and the structure.  Considering these major differences, 
the correspondence between the model and simulation is quite remarkable and 
indicates that the model captures the dominant aspects of the governing physics.   
The results also show that inclusion of heat transfer through the structure leads to 
an important difference in how the burning velocity is related to the axial temperature 
gradient in freely propagating flames versus those stabilized in micro-channels. Both 
the original theory of Mallard and Le Chatelier and the theory for microchannels 
presented here are predicated on the notion that in order for a flame to propagate, heat 
 
 35 
 
from the post-flame must be transferred upstream to heat the incoming reactants to 
the ignition temperature.  Therefore, increasing the velocity of the reactants in both 
theories requires an increase in the heat flux from the post-flame to the pre-flame.  In 
the original theory of Mallard and Le Chatelier, the only mechanism for 
accomplishing this is to increase the strength of the temperature gradient through the 
flame by decreasing the flame thickness.  Thus, in freely-propagating flames, the 
burning velocity must be inversely related to the flame thickness.  In micro-channels, 
however, an additional heat transfer path through the structure is available.  This has 
two consequences:  First, the presence of the second path takes some of the load off 
axial conduction in the gas.  This means that less heat needs to be transferred through 
this path and the temperature gradient through the flame does not have to be as steep.  
This, in turn, means that the reaction zone can be longer.  Second, the heat transfer 
through the structure increases with increasing contact area with the gas.  Since the 
contact area increases with increasing reaction zone thickness, another mechanism for 
increasing heat transfer upstream is to increase the reaction zone thickness.   
Therefore for a flame stabilized in a micro-channel, the burning velocity is directly 
related to the flame thickness.  
Finally, it is important to note that some of the data points in the left portion of 
figure 2-4 correspond to passage heights that are so small that the continuum 
hypothesis, and thus our numerical representation of the reacting fluid, may no longer 
apply.  These sections have been shaded gray in this figure and the ones that follow.  
The model and simulation results are only projected into these areas to illustrate the 
complete limit-cycle behavior.  The common metric for identifying length scales 
 
 36 
 
where non-continuum effects become important is the Knudsen number.  It is defined 
as the ratio of the mean free path ? of a gas molecule to the characteristic length scale 
of the device.  In this work, the characteristic length scale is the passage height so: 
H
Kn
?
=           (2-22) 
Table 1 summarizes the commonly accepted boundaries for the different flow 
regimes as a function of the Knudsen number
41
.   
 
Flow Regimes 
 
Kn < 0.01 Continuum 
0.01 < Kn < 0.1 Slip 
0.1 < Kn < 10 Transitional 
10 < Kn Free Molecular
 
Table 1:  Flow regimes associated with different Knudsen number ranges 
from
42
 
 
The secondary x-axis on figure 2-4 shows that the smallest passage heights fall 
into the slip regime.  While this is not necessarily a problem because the Navier-
Stokes equations are valid through the entire slip regime
43
, slip boundary conditions 
should be applied in the slip regime to account for possible discontinuities in velocity 
and temperature at the wall.  Since the momentum equation is neglected in this work, 
 
 37 
 
the only discontinuity that needs to be accounted for is the gas temperature at the 
wall.  A fit to results of direct simulation Monte Carlo calculations for fully 
developed flow between infinite parallel plates
44
 shows that the Nusselt number 
varies with Knudsen number as follows: 
 
Kn
eNu
86.3
86.7
?
=                        (2-23) 
Figure 2-5 compares predictions of flame broadening made assuming the Nusselt 
number is constant (7.86, the value associated with fully developed continuum flow 
between parallel plates) to those made using equation 2-23.  It shows that the near-
wall temperature discontinuity in the slip regime has a negligible impact on the 
broadening phenomenon.  Therefore, it does not appear that including this effect is 
particularly important for predicting the performance of micro-combustors.  
Furthermore, the smallest dimensions shown in figure 2-5 are probably not very 
realistic for practical devices.  Instead, the lower size limit will most-likely be set by 
heat loss from the structure to the environment.  This will be discussed later. 
 
 
 38 
 
H/?
r,fr
Kn
?
r
/
?
r,
fr
10
-2
10
-1
10
0
10
1
10
2
10
-5
10
-4
10
-3
10
-2
10
-1
0
1
2
3
4
5
6
7
8
Nu(Kn)
Nu=const
Nu(Kn)
Nu=const
Simulation
Model
?=48, ?=1
Figure 2-5: Comparison of the fractional increase in reaction zone length d
r
 as a 
function of the normalized flow passage height H/d
r,fr
 when temperature slip at 
the wall is and is not included.  The lines without symbols correspond to the 
analytical model and the lines with symbols correspond to the numerical 
simulations. The solid lines show the results when temperature slip is included 
while the dashed lines correspond to a fixed Nusselt number of 7.86 
corresponding to fully-developed flow between two parallel plates. 
 
 
 39 
 
Figure 2-6 compares analytical and numerical predictions of the non-dimensional 
flame broadening ?
r
/?
r,fr
 for flow between infinite parallel plates as a function of the 
normalized passage height H/?
r,fr
  for three values of k
s
.  The change in Nusselt 
number with Knudsen number is included in this and all subsequent figures.  Both the 
model and the simulation show that the broadening depends on k
s
.  Increasing k
s
 
increases the heat transfer rate through the structure, which thereby reduces the 
necessary heat transfer rate through the gas phase. This enables the temperature 
gradient through the gas to be shallower and allows the reaction zone thickness to be 
longer than that required for a freely propagating flame. While both models show that 
the passage height corresponding to the transition from gas to structure-dominated 
thermal exchange does not depend strongly on k
s
, the amplitude of the broadening 
phenomenon appears to be a stronger function of the thermal conductivity of the 
structure in the analytical model than in the numerical model.   
 
 40 
 
H/?
r,fr
Kn
?
r
/
?
r,
fr
10
-2
10
-1
10
0
10
1
10
2
10
-5
10
-4
10
-3
10
-2
10
-1
0
1
2
3
4
5
6
7
8
?=10
?=24
?=48
?=10
?=24
?=48
{
Model
Simulation
{
Nu(Kn), ?=1
 
Figure 2-6: Fractional increase in reaction zone length ?
r
 as a function of the 
normalized flow passage height H/?
r,fr
  for a range of normalized structure 
thermal conductivities ? . Lines without symbols correspond to the analytical 
model and the lines with symbols correspond to numerical model results and 
k
r
=0.21 W/mK. 
 
So far, the ratio of the cross-sectional areas for axial conduction of the structure to 
that of the gas has been help constant as the passage height is varied.  However, 
 
 41 
 
increasing the relative thickness of the structure compared to the gas will allow more 
heat to be conducted through the structure in the same way than increasing structure?s 
thermal conductivity would.  Therefore, the truly important parameter is the 
structure?s conductance which is defined as ?? .  Figure 2-7 shows the effect of the 
structure?s conductance on the flame broadening phenomenon for three values of 
H/?
r,fr
.  Both the analytical model and numerical simulation indicate that the increase 
in the reaction zone thickness approaches a limit with increasing conductance of the 
structure.  This limit depends on Nu as well as H.  For a highly conductive structure, 
the resistance to heat transfer at the gas/structure interface limits flame broadening. 
The magnitude of broadening predicted at the high structure conductance limit by the 
analytical model is approximately 50% larger than that predicted by the numerical 
simulation, but both show similar limit cycles. 
 
 42 
 
??
?
r
/
?
r,
fr
10
-2
10
-1
10
0
10
1
10
2
10
3
10
4
10
5
1
2
3
4
5
6
H/?
r,fr
=2.4
H/?
r,fr
=1.2
H/?
r,fr
=0.6
H/?
r,fr
=2.4
H/?
r,fr
=1.2
H/?
r,fr
=0.6
{
{
Nu(Kn), ?=1
Model
Simulation
 
Figure 2-7: Fractional increase in reaction zone length ?
r
 as a function of 
non-dimensional thermal conductance (?? ) of the structure for a range of 
normalized passage heights.  ?
r,fr
 is the reaction zone length associated with a 
freely propagating flame, ? =1, Nu=3, and k
r
=0.21 W/mK.  Lines without 
symbols correspond to the analytical model and lines with symbols 
correspond to numerical model results. 
 
 
 43 
 
So far, we have assumed that the combustor is adiabatic and that no heat is lost 
from the structure to the environment.  The results for the adiabatic structure cases in 
figures 2-5 ? 2-7 indicate that conduction of heat from the post-flame through the 
structure to the upstream preheat zone can broaden the reaction zone by a factor of 4 
or more.  The results also show that the thermal resistance between the gas flow and 
the structure limits heat transfer upstream and in turn limits the magnitude of the 
broadening phenomenon. 
 
2.4.2. Non-adiabatic results 
Heat loss from the structure to the environment is also important in micro-scale 
combustors because it dissipates energy that could be used to pre-heat the reactant 
mixture and can also lead to flame quenching.  Figure 2-8 shows the non-dimensional 
flame broadening as a function of the normalized passage height for three values of 
h
T,e
 for both the analytical model and the numerical simulation.  For h
T,e
 of O(10) or 
less, reaction zone broadening arises as before, but as H decreases the heat loss to the 
environment eventually leads to a reduction in the flame zone thickness followed by 
extinction.  Extinction occurs when the surface to volume ratio of the micro-channel 
becomes so large that heat loss to the environment prevents sufficient thermal 
feedback to the preheat zone to sustain continuous ignition of the flame.  As h
T,e
 
increases to large values, the broadening phenomenon is replaced by a rapid reduction 
of the flame thickness with decreasing H and flame extinction occurs at relatively 
large H.  The reduction in reaction zone thickness with decreasing H for structures 
with large thermal losses arises because steeper temperature gradients in the gas and 
 
 44 
 
structure must be sustained in order to transfer an adequate fraction of the reaction 
heat release back to the preheat zone.  These results suggest that the ability to prevent 
heat loss to the environment sets one limit on the minimum size of a micro-
combustor.  Finally, it is important to note that the heat transfer coefficient has been 
varied over an extremely wide range in figure 2-8 in order to illustrate what happens 
when heat loss to the environment is relatively weak and relatively strong.  It would 
probably be quite difficult to achieve a heat transfer coefficient of 10 W/m
2
K or 
smaller in a practical device. 
H/?
r,fr
Kn
?
r
/
?
r,
f
r
10
-2
10
-1
10
0
10
1
10
2
10
-5
10
-4
10
-3
10
-2
10
-1
0
1
2
3
4
5
6
7
8
h
T,e
=0
h
T,e
=10
h
T,e
=1000
h
T,e
=0
h
T,e
=10
h
T,e
=1000
{
{
Simulation
Model
Nu(Kn), ?=48, ?=1
 
Figure 2-8: Fractional increase in reaction zone length ?
r
 as a function of the 
normalized flow passage height H/?
r,fr
  for a range of heat transfer 
coefficients in W/m
2
K. Lines without symbols correspond to the analytical 
model and the lines with symbols correspond to numerical model results. 
 
 45 
 
 
While the qualitative agreement between the analytical and numerical models 
lends credence to the principles underlying the modified Mallard-Le Chatelier theory, 
the quantitative discrepancy between the analytical model and numerical simulation is 
hardly surprising.  As indicated earlier, the discrepancy can arise from numerous 
factors including the implementation of variable properties and multi-step chemistry 
in the numerical model.  In addition, the single-step chemistry of the analytical model 
assumes that the reaction goes to completion while the multi-step chemistry of the 
analytical model accounts for dissociation and the accompanying reduction in flame 
temperature.  Finally, the numerical model distributes the heat feedback from the 
structure throughout the reaction and preheat regions.  This reduces to some extent 
the impact of the structural heat feedback in comparison to the analytical model 
which sends all of the heat transferred to the structure in the post-flame back to the 
preheat zone.  Accordingly, the numerical model is less sensitive to changes in 
passage height, structure thermal conductivity, and Nusselt number than the analytical 
model. 
 
2.4.3. Preliminary implications on power density 
While the broadening phenomenon is of interest theoretically, it also has very 
practical implications on the development of high power density combustors.  For the 
purpose of this work, the power density is defined as follows: 
V
W
w
D
&
& =                                                     (2-22) 
 
 46 
 
where W
&
 is the power generated by the combustion process and V is the volume of 
the combustor in which the chemical energy is released.  We have already shown that 
thermal exchange with the structure increases the reaction zone thickness and this 
observation along with equation 2-22 indicates that the overall power density of the 
combustor might decrease as the size of the combustor is reduced.  This, of course, 
would run counter to the objective of miniaturizing combustors in the first place, 
which is to increase power density.  At the same time, however, we have already 
shown that the increased heat feedback from the structure to the preheat zone results 
in a substantial increase in the allowable burning velocity so the net effect is not clear.   
The relative increase in the non-dimensional burning velocity and the non-
dimensional flame thickness determined using the numerical model is illustrated in 
figure 2-9.  It is important to note that the simple model predicts that the power 
density should remain constant as the passage height is varied because changes in 
flame speed are directly related to changes in flame thickness (since 
RRS
Lr
=?
 is 
assumed).  The results from the numerical model presented in figure 2-9, however, 
show that the evolution of those parameters is slightly different in both the amplitude 
of the increase and the passage height H corresponding to the transition point from a 
thin to a thick flame. This discrepancy can again be attributed to the mismatch in 
complexity between the simple analytical model and the numerical simulation.   
 
 47 
 
H/?
r,fr
Kn
10
-2
10
-1
10
0
10
1
10
-4
10
-3
10
-2
10
-1
1
2
3
4
5
?
r
/ ?
r,fr
S
L
/S
L,fr
Nu(Kn), ?=48, ?=1
 
Figure 2-9: Comparison between the evolution of the normalized flame 
thickness, ?
r
/?
r,fr
, and the normalized laminar flame speed, S
L
.  Results are 
from the 1D numerical model. 
 
The results from the numerical model presented in figure 2-9 also suggest that 
maximizing power density may require a careful optimization over geometric (as well 
as thermophysical) properties of the structure.  While the results of such an 
optimization are reported in the next chapter of this work, some insight can be gained 
 
 48 
 
by computing the variation of power density with passage height for two different 
combustor lengths.   The power density is computed from the numerical results using 
the following equation:  
()
L
dTTCU
w
f
T
T
p
D
?
=
0
?
&                                      (2-23) 
Note that the passage height and width appear in both the numerator and denominator 
and thus cancel each other out.  As a result, the power density depends only on the 
enthalpy change per unit area across the flame and the passage length L.  
Figure 2-10 shows the evolution of non-dimensional power density 
D
w&  as a 
function of normalized passage height H/?
r,fr
 for two different values of the 
normalized passage length L/?
r,fr
. The structure is assumed to be adiabatic, the values 
of Nu, ?, and ? are fixed, and the power density is made non-dimensional by 
normalizing by the power density associated with an arbitrarily defined ?macro-scale? 
configuration corresponding to L=1cm and H=1cm.  ?he plot shows that reducing the 
passage height and reducing the passage length both increase the power density.  The 
increase in power density with decreasing passage height is a direct consequence of 
the increased burning rate caused by enhanced heat transfer to the pre-flame at small 
passage heights while the increase in power density with decreasing combustor length 
is a consequence of miniaturization only.  
 
 49 
 
H/?
r,fr
Kn
W
D
/W
D,
m
a
c
r
o
10
-1
10
0
10
1
10
-4
10
-3
10
-2
0
5
10
15
20
25
Nu(Kn), ?=48, ?=1
L/?
r,fr
=10
L/?
r,fr
=2
 
Figure 2-10: Normalized power density as a function of normalized passage 
height for two values of the normalized passage length.  Results are from the 
1D numerical model. 
 
 
2.4.4. Pressure loss 
Pressure loss is another important factor affecting the design, performance, and 
ultimate practicality of micro-combustors.  A first order estimate of its magnitude in 
 
 50 
 
the parallel plate geometry investigated here is made assuming that the flow is fully-
developed
45
: 
2
8
H
UL
p
?
=?                                                (2-24) 
In this expression ?p is the pressure drop, ? is the dynamic viscosity (~50.10
-6
 kg/m.s 
for air at 1500K) and U is the velocity (taken to be the laminar flame speed S
L
 as an 
approximation).  The pressure drop corresponding to the analytical model is estimated 
by taking the length L to be equal to the minimum combustor length (or the flame 
thickness) while the pressure drop corresponding to the numerical simulation is 
estimated by taking L to be equal to 1 cm (the length of the computational domain).  
The solid lines in figure 2-11 show pressure drop normalized by the pressure at the 
channel entrance (1 atm) as a function of non-dimensional passage height.  The flame 
broadening results are also plotted (dashed lines) so that the passage height where the 
broadening phenomenon begins to occur can be compared to the passage height at 
which pressure losses start to become large.  It is clear from the figure that the 
broadening phenomenon occurs well before the pressure losses between the plates 
become significant (i.e. greater than 5% of the total pressure).  As a result, flame 
broadening appears to be a phenomenon that could be encountered in a ?practical? 
micro-combustor.  
 
 51 
 
H/?
r,fr
Kn
?
p/
p
?
r
/
?
r,
fr
10
-2
10
-1
10
0
10
1
10
2
10
-5
10
-4
10
-3
10
-2
10
-1
0
0.2
0.4
0.6
0.8
1
0
1
2
3
4
5
6
7
8
?p / p - Analytical
?p/p-Numerical
?
r
/ ?
r,fr
-Analytical
?
r
/ ?
r,fr
- Numerical
Nu(Kn), ?=48, ?=1
 
Figure 2-11: Normalized pressure drop as a function of normalized passage 
height.  The evolution of the non-dimensional flame thickness is 
superimposed on the plot for comparison. 
 
Finally, it should be noted that the increase in burning velocity with decreasing 
passage height is consistent with other work that uses somewhat different methods to 
model the chemical reaction and account for the effects of axial heat conduction in 
the structure.  In particular, Ronney
46
 used a perfectly stirred reactor (PSR) coupled 
 
 52 
 
with a counter-current heat exchanger to show that higher mass flow rates were 
required to stable combustion when heat recirculation was present or, equivalently, 
when the thermal conductivity of the tubes was not small.  This is qualitatively 
equivalent to the predictions made by the analytical model as well as those made by 
the numerical model where the distributed nature of the reaction zone is accounted 
for. 
 
2.5. Conclusion 
A simple model based on a thermal resistance network has been developed to 
extend the Mallard-Le Chatelier theory for a pre-mixed flame freely propagating into 
a reactant mixture, to a premixed flame propagating in and exchanging heat with a 
thermally conducting channel. In spite of its simplicity, comparisons to the results of 
more detailed numerical simulations that include axially distributed heat transfer and 
full chemistry indicate that the simple model captures many of the important features 
of a combustion wave stabilized in a micro-channel.  Both the analytical and 
numerical models show that heat exchange between the reacting fluid and the 
structure tends to broaden the reaction zone.  Both models predict that the passage 
height at which axial heat transfer becomes important depends on the Nusselt 
number, and for a given Nusselt number predict similar values for the transition 
passage height.  Both models also predict limiting behaviors at large and small 
channel heights, and both indicate that increasing the conductance of the structure 
(either by increasing the thermal conductivity of the structure material or the structure 
thickness) increases the amplitude of the broadening phenomenon.  The maximum 
 
 53 
 
limit of flame broadening at very small heights is predicted within a factor of two by 
the analytical model in comparison to the numerical simulation.  The simple model 
has been extended to include the effects of quenching via heat loss to the 
environment.  The model extension does not predict quenching diameters well, 
however they still lie within a factor of 10 of the numerical simulations? predictions 
and show the same qualitative trends with passage height.   
While the increase in flame thickness at micro-scale may initially appear to be a 
drawback for designing high power density combustors because it suggests that 
proportionally more combustor volume is required, the increase in burning rate 
associated with pre-heating the reactants more than compensates for this effect and 
the net result is an increase in power density greater than what one would expect as a 
result of miniaturization alone.  Finally, it should be noted that mass transport to the 
walls as well as surface reactions there are other important aspects of the micro-
combustion problem.  These are not considered here and are excellent topics for 
future work. 
Taken together, the results demonstrate the utility of the simple extension to the 
Mallard-Le Chatelier theory and the importance of fluid-structure coupling in micro-
combustion.  Experimental data is required to verify the predictions of the model and 
additional work in the modeling area is needed to explore optimum combustor 
configurations and the influence of the choice of the fuel on power density and the 
broadening phenomenon.  The results of these investigations will provide valuable 
information to designers of microscale combustion systems.  Work focused on 
identifying optimum combustor configuration is presented in the next chapter. 
 
 54 
 
Chapter 3:  Power Density and Efficiency of Micro-
Combustors 
 
3.1. Introduction 
In the previous chapter of this work, an analytical model was developed to explain 
the fact that reaction zones appear broader in small-scale combustors.  This model 
also predicted that the burning velocity should increase because of axial heat 
conduction through the structure.  The numerical simulations provided preliminary 
evidence that the net effect of miniaturization was an increase in power density. 
However, the results presented were computed for only a few configurations of the 
combustor and the general behavior of power density as both the height and length of 
the micro-combustor were decreased was unclear.  In addition, the effect of heat loss 
on power density was not included.  Nevertheless, these preliminary results as well as 
intuition suggested the possibility for an optimum combustor configuration that 
maximizes power density. 
The objective of the work reported in this chapter is to systematically investigate 
the dependence of power density on the combustor design which, in this case, is 
determined by its length and height.  The overall combustion efficiency is also 
computed as it is another important performance metric.  Finally, the influence of fuel 
type is considered by performing the calculations with both H
2
 and CH
4
 fuels and 
comparing the results.   
 
 55 
 
The results shown for hydrogen combustion were presented at the 30
th
 
International Symposium on Combustion held at UIC in March 2004.  They were 
published in Proceedings of the Combustion Institute 30 (2005) pp. 2437-2444. 
The results for methane combustion were presented at the 43rd AIAA Aerospace 
Sciences Meeting and Exhibit, Reno, Nevada, January 10?13, 2005 and partly 
published as ?Effect of Fuel Type on Optimum Micro-Combustor Design and 
Performance?, AIAA Paper 2005-939. 
 
3.2. Approach 
The numerical model presented in Chapter 2 of this work is used to model micro-
combustion between two parallel plates.  A schematic diagram of the problem 
considered can be found in the previous chapter of this work (Figure 2-4).  Surface 
reactions are not included and hence there is only thermal (no mass) exchange 
between the gas and the structure.  The Nusselt number for heat transfer between the 
gas and structure is assumed constant at 8, the value for fully developed flow between 
two parallel plates
47
.    
Since a ?realistic? combustor had to be considered in order for the results obtained 
to be of practical use, the plates were assumed to be made out of Silicon, a typical 
MEMS fabrication material.  The use of silicon as a combustor material imposes two 
additional constraints on the simulation.  The first is that the equivalence ratio be low 
enough so that the adiabatic flame temperature remains below the melting point of 
silicon (1680K).  The second is that the temperature dependence of silicon?s thermal 
 
 56 
 
conductivity be included.  This is important because silicon?s thermal conductivity 
changes by more than a factor of five over the expected range of temperatures and is 
computed using a power law fit to tabulated data
48
.   
Using the numerical simulation it is then possible to compute the maximum inlet 
velocity for specific values of the geometric dimensions of the combustor L and H 
and the inlet velocity U.  The maximum inlet velocity is chosen because we are 
interested in determining the maximum power density associated with a particular 
configuration. 
 
3.3. Results and Discussion 
3.3.1. Hydrogen Fuel 
In this first part of the results section, we consider hydrogen-air combustion at 
equivalence ratio of 0.5.  The gas-phase combustion chemistry is modeled using a 
nine-species, nineteen-reaction mechanism
49
.  Table 4 in Appendix II shows the 
reactions considered in the mechanism.   
Figure 3-1 is a plot of non-dimensional inlet velocity versus non-dimensional 
passage height for two non-dimensional channel lengths. First, note the increase in 
flame speed (burning velocity) as the channel height is reduced.  Second, it shows 
that there are two limiting values of the flame speed at a particular design point.  The 
minimum value corresponds to the flashback limit where a flame is stabilized near the 
channel entrance.  The maximum value corresponds to the blow-off where the flame 
is stabilized near the channel exit.  A continuum of flame speeds is possible in 
 
 57 
 
between these limits.   Third, it shows that at large passage heights where the thermal 
coupling with the structure is weak, the blow-off and flashback velocities converge to 
the flame speed associated with a freely propagating flame.  As the passage height is 
reduced and the strength of the coupling between the reacting fluid and the micro-
channel structure increases, the blow-off and flashback limits diverge and it becomes 
possible to stabilize a combustion wave over a range of inlet velocities.  Fourth, and 
finally, decreasing the passage length tends to decrease the magnitude of the flame 
speed increase associated with decreasing passage height, and reduces the difference 
between the blow-off and flashback limits.    
10
-3
10
-1
10
1
10
3
0
5
10
15
H / 
?
r,fr
U /
 
S
L,f
r
Min., L / 
?
r,fr
=10
Max., L / 
?
r,fr
=10
Min., L / 
?
r,fr
=1
Max., L / 
?
r,fr
=1
 
Figure 3-1: Non-dimensional flashback and blow-off limits for H
2
-air 
combustion in a silicon micro-channel as a function of non-dimensional 
channel (combustor) height and for two values of the channel length L.   The 
equivalence ratio is 0.5 and the pressure is 1 atm. 
 
 
 58 
 
The reasons for these observations can be inferred from figure 3-2 which shows 
axial profiles of non-dimensional temperature and major species mass fractions at 
inlet velocities corresponding to blow-off and flashback.  
0 0.2 0.4 0.6 0.8 1
0
0.5
1
x / L
T/T
ad
Y
O
2
Y
H
2
O
0 0.2 0.4 0.6 0.8 1
0
0.5
1
x / L
T/T
ad
Y
O
2
Y
H
2
O
Blow off Limit 
Flashback Limit 
U / S
L,fr
=14 
U / S
L,fr
=11 
 
Figure 3-2: Axial profiles of temperature and major species concentration for 
a H
2
-Air flame stabilized in a silicon micro-channel with 
frr
H
,
? = 0.002 and 
frr
L
,
? = 9.  The equivalence ratio is 0.5 and the pressure is 1 atm. 
 
 59 
 
 
 The profiles indicate that the flame is stabilized near the channel exit at the blow-
off limit while it is stabilized near the channel entrance at the flashback limit.  The 
flame speed is higher when the flame is stabilized at the channel exit because there is 
more surface area for heat exchange between the structure and the reactants, and 
hence more pre-heating of the mixture before reaction occurs.  Decreasing the overall 
length of the combustor reduces the surface area for heat exchange between the 
structure and the reactants, reduces the pre-heating of the mixture, and thus reduces 
the flame speed.  The blow-off and flashback limits come together as the channel 
length is reduced because decreasing the length of the micro-channel decreases the 
range of possible stable flame locations within the channel.  This, in turn, decreases 
the range of stable inlet velocities.    
It is interesting to note that the temperature profiles presented in figure 3-2 do not 
show the temperature overshoot characteristic of excess enthalpy burners.  This is 
because the profiles in figure 3-2 correspond to relatively small passage heights 
where the heat transfer coefficient is large and thermal equilibrium between the gas 
and the structure is established very quickly.  This is consistent with the results of 
Takeno
50
 which show a reduction in the thermal overshoot when the thermal coupling 
between the gas and the structure is strong.  Thermal overshoot is observed at larger 
passage heights where the thermal coupling is weaker and thermal equilibrium takes 
longer to establish.  The physical bases for the ?smoothing out? of temperature and 
species gradients in micro-scale systems have been discussed by Pello
51
.  
Thermal coupling between the reacting gas and the structure has important 
implications for the design of efficient, high power density combustors.  The overall 
 
 60 
 
efficiency of a micro-channel combustor can be written as the product of a chemical 
efficiency and a thermal efficiency: 
thcho
??? =                           (3-1) 
The chemical efficiency is defined as the ratio of the total heat evolved in the reaction 
to the total heat available if the reaction were to go to completion.  In this work, we 
assume that this is equivalent to the fraction of H
2
 that gets consumed before the 
mixture exits the micro-channel: 
 
i
out
H
H
ch
Y
Y
availableheatMaximum
rxnbyevolvedheatTotal
,2
,2
1
.
?==?     (3-2) 
The thermal efficiency is defined as the ratio of the total heat actually delivered to 
the flow exiting the channel divided by the total heat evolved by the reaction.  This is 
computed from the axial temperature profile through the micro-channel, the heat loss 
to the environment, and the mass flow rate: 
mqdTTC
dTTC
rxnbyevolvedheatTotal
flowtodeliveredheatTotal
loss
T
T
p
T
T
p
th
out
i
out
i
&&+
=
=
?
?
)(
)(
?
     (3-3) 
The power density is defined as follows: 
V
W
w
D
&
& =                                  (3-4) 
 
 61 
 
where W
&
 is the energy release rate due to combustion (computed as 
?
out
i
T
T
p
dTTCm )(& ) 
and V  is the combustor volume.   
The solid contours in figure 3-3 show maximum non-dimensional power density 
as a function of micro-channel length and height for H
2
-air combustion at an 
equivalence ratio of 0.5, atmospheric pressure, and no heat loss to the environment.  
Channel height and length are non-dimensionalized by the reaction zone thickness of 
a freely propagating H
2
-air flame.  Power density is non-dimensionalized with respect 
to a reference value associated with combustion in a channel where L and H equal 1 
cm.  This reference channel is sufficiently large that there is essentially no influence 
of the structure on the combustion process.  Note that the maximum power density 
associated with a particular combustor design (i.e. L and H ) occurs when the inlet 
velocity is at the blow-off limit.   
Figure 3-3 shows that there is an optimum micro-channel length where power 
density is maximized.  The maximum is a direct consequence of thermal exchange 
within the structure.  It arises from a competition between the decrease in channel 
volume and the decrease in flame speed (and hence energy release rate) as the length 
of the combustor is reduced.  The figure shows that thermal exchange within the 
channel wall makes it possible to achieve power densities in excess of 60 times the 
maximum associated with the reference channel (H=1cm, L=1cm).   Note that the 
power density also decreases as the channel height is increased.  This is because 
increasing the channel height reduces the thermal coupling between the fluid and the 
micro-channel wall.   
 
 62 
 
10
0
10
1
10
-2
10
-1
10
0
L / 
?
r,fr
H
 
/
 
?
r
,fr
60 
50 
40 
30 
20 
10 95 % 
90 % 
85 % 
80 % 
 
Figure 3-3:  Adiabatic H
2
-air combustion in a silicon micro-channel; ?=0.5, 
P=1atm.  Solid lines: Contours of non dimensional power density ( )
refDD
ww
,
&&  
as a function of non-dimensional channel height and length.  Dashed lines: 
Contours of overall efficiency.    
 
The dashed contours in figure 3-3 show the overall efficiency of the micro-
channel combustor as a function of channel height and length.  Since the combustor is 
adiabatic, the thermal efficiency is 100% and changes in overall efficiency are 
entirely the result of changes in chemical efficiency.  The figure shows that power 
density and chemical efficiency trade against each other:  Attempts to increase power 
density by decreasing the length of the channel reduce the degree of chemical 
conversion by decreasing the time available to complete chemical reaction.  The 
figure also shows that heat exchange with the structure drives the configurations 
associated with maximum efficiency and maximum power density farther apart.   
 
 63 
 
Thermal efficiency also changes with combustor configuration when heat loss to 
the environment is included.  Figure 3-4 shows that heat loss, simulated by 
introducing a heat transfer coefficient of 1 W/m
2
K between the micro-channel and the 
environment, eventually quenches reaction as channel height is reduced.  One would 
expect similar behavior as channel length is reduced because of reduced chemical 
efficiency.   
10
-3
10
-1
10
1
10
3
0
5
10
15
H / 
?
rfr
?
r
 / 
?
r,fr
 - Ad.
S
L
 / S
L,fr
 - Ad.
?
r
 / 
?
r,fr
 - Non ad.
S
L
 / S
L,fr
 - Non ad.
 
Figure 3-4: Adiabatic (solid lines) and non-adiabatic (dashed lines) H
2
-air 
combustion in a silicon micro-channel; ?=0.5, P=1atm, h
env
=1 W/m
2
K, 
L/?
r,fr
=20.  The circles show non-dimensional laminar flame speed as a 
function of non-dimensional channel height while the squares show non-
dimensional reaction zone thickness. 
 
 
 64 
 
These effects are observed in figure 3-5 as the contours of overall efficiency 
(dashed lines) shift upward and to the right.  The subsequent reduction in overall 
efficiency that arises from heat loss at small passage heights reduces the power 
density and creates the global maximum in power density shown in the solid 
contours.  There is now a single, optimum, combination of L and H that maximizes 
power density.  Not surprisingly, the maximum power density and the efficiency at 
this point are somewhat lower than those in the adiabatic case.  It is important to note 
that the non-adiabatic power density is only slightly less than the adiabatic value 
because the heat transfer coefficient between the combustor and the environment is 
very small.  In practice, one would expect heat losses to be much larger and the 
maximum realizable power density to be much lower. 
10
0
10
1
10
-2
10
0
L / 
?
r,fr
H /
 
?
r,
f
r
60 
50 
40 
30 
20 
10 
95 % 
90 %
85 % 
80 % 
20 
10 
75 % 
 
Figure 3-5: Non-adiabatic H
2
-air combustion in a silicon micro-channel; 
?=0.5, P=1atm, h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional 
power density ( )
refDD
ww
,
&&  as a function of non-dimensional channel height 
and length.  Dashed lines: Contours of overall efficiency.  
 
 65 
 
 
Figure 3-6 is a plot of power density and efficiency at a combustor pressure of 5 
atmospheres under non-adiabatic conditions.   
10
0
10
1
10
-3
10
-1
10
1
L / 
?
rfr
H /
 
?
r,
f
r
300 
250 
200 
150 
100 
50 
20 
10 
95 % 
90 % 
85 % 
80 % 
 
Figure 3-6: Non-adiabatic H
2
-air combustion in a silicon micro-channel; 
?=0.5, P=5atm, h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional 
power density ( )
refDD
ww
,
&&  as a function of non-dimensional channel height 
and length.  Dashed Lines: Contours of overall efficiency. 
 
It shows that increasing the pressure increases the maximum power density and 
reduces the values of L and H associated with the optimum. The reason for this is that 
the reaction zone thickness decreases with increasing pressure.  The power density 
increases because the flame speed is higher and because an acceptable level of 
chemical conversion can be achieved in a smaller volume.  Note that an increase in 
 
 66 
 
flame speed with pressure is consistent with mixtures having burning velocities 
greater than 50 cm/s
52
.  In these flames, temperatures are high enough that 
dissociation is important.  Increasing the pressure represses dissociation, steepens the 
temperature profile, and increases the burning velocity.  Since thermal feedback to the 
pre-flame through the micro-channel structure tends to increase the mixture 
temperature, one would expect flames stabilized in micro-channels to be more 
sensitive to pressure than freely propagating flames. 
Figure 3-6 also shows that increasing the pressure increases the overall efficiency 
at the maximum power density point.  The reason for this is that the wall temperature, 
and hence the heat loss from the micro-channel to the environment, remains 
approximately constant as the pressure is increased.  This increases the thermal 
efficiency by decreasing the ratio of the heat lost to the heat produced.   
Figure 3-7 illustrates the general effect of pressure on the the performance of 
combustors at conventional scale and micro-scale.  The figure shows that for our 
?reference? meso-scale combustor for which H=L=1cm, power density increases as 
P
1.2
.  This is to compare to the increase in maximum power density of a micro-scale 
combustor, which increases as P
2
. 
 
 67 
 
 
Figure 3-7: Evolution of power density with increasing pressure for a meso-
scale combustor and a micro-combustor. 
 
The results of the power density and efficiency calculations are summarized in 
Table 2 and show that the reaction zone thickness doesn?t change much as pressure is 
increased from 5 to 10 atm.  As a result, the optimum channel length and height 
associated with the maximum power density remains approximately constant.  It is 
important to point out that the optimum values of L and H listed in table 2 are 
estimated using the contour plots and are not the result of an optimization in L and H.  
Since the peaks in power density are relatively flat, the uncertainty in the optimum 
values of L and H is relatively large.  As a result, we are unable to make precise 
10
0
10
1
10
2
10
3
10
4
10
5
10
6
10
7
P (atm)
W
D
 (
M
W / m
3
)
W
D,ref
W
D,max
W
D
,r
e
f
 ~
 P
1
.2
 
W
D
,
m
a
x
 
~
 
P
2
 
 
 68 
 
estimates of the optimum combustor volume per unit channel depth and the values in 
the last row of table 2 are included only as general references. 
 
 1 atm  5 atm  10 atm 
Parameter Ad. 
Non-
Ad. 
 Ad. Non-Ad.  Ad. Non-Ad. 
ref,DD
ww &&  
64 60  326 323  475 471 
fr,r
H ?  
<0.002 0.01  <0.012 0.006  <0.1 0.016 
fr,r
L ?  
1 1  1.31 1.31  1.6 1.6 
ref,L
SU  
5.44 5.16  4.44 4.42  5.06 5.01 
Efficiency (%) 90 89  92 91  96 96 
ref,D
w& (MW/cm
3
) 
177 177  1230 1230  3180 3180 
fr,r
? (cm) 
0.055 0.055  0.008 0.008  0.006 0.006 
fr,L
S (cm/s) 
90 90  136 136  187 187 
V (cm
3
/cm) <6x10
-6
 30x10
-6
  <1x10
-6
 0.5x10
-6
  <5.8x10
-6
 0.9x10
-6
 
 
  
Table 2: Summary of optimum power density configurations for a silicon 
micro-combustor.  The fuel is H2 burning in air and the equivalence ratio is 
0.5. 
 
Finally, is important to note there are other factors that are important to the design 
of optimum micro-combustors that have not been discussed here.  The first is the 
conductance (the product of the thermal conductivity and the cross-sectional area) 
which limits upstream heat conduction through the structure. In this work, the ratio of 
 
 69 
 
the wall thickness to the passage height is fixed at ?.  If this value were made very 
small, no thermal feedback would occur and hence no optimum power density 
configuration would be observed.  Similarly, making this parameter very large would 
increase the power density at the optimum until a limit imposed by the heat transfer 
coefficient between the gas and the structure is reached.  The second factor is the 
pressure drop which is a very important consideration when integrating a combustor 
into a power cycle.  The pressure drop cannot be determined using the present one-
dimensional formulation that does not include the momentum equation.  However, we 
have used a simple Poiseuille flow approximation to provide an estimate.  Pressure 
drops near the optimum in figure 3-6 are about 10% or less.  This indicates that 
combustors designed to be at or near the optimum power density configuration could 
be integrated into practical power systems.   
 
3.3.2. Methane Fuel 
The results presented in the previous part of this work considered hydrogen-air 
combustion.  To determine how the results obtained change with fuel, methane-air 
combustion at an equivalence ratio of 0.5 is now considered.  The Leeds 
mechanism
53
, which includes 37 species and 175 reactions, is used to represent the 
gas phase chemistry.  The Leeds mechanism was chosen because it does not include 
nitrogen chemistry, which is not of interest in the present study and greatly increases 
the speed of the calculations.  Using a larger mechanism like GRI-Mech would have 
been computationally prohibitive.  Table 5 in Appendix II shows the reactions 
considered in the mechanism.   
 
 70 
 
Temperature overshoot is a well-known characteristic of heat recirculating 
burners
24
. Some of the enthalpy generated in the post-flame can be transferred to the 
pre-heat zone, allowing combustion to occur at a higher preheat temperature and thus 
increasing the flame temperature and burning velocity.  A very comprehensive 
analysis of heat recirculating burners can be found in Weinberg
24
.  
A non-dimensional temperature overshoot is defined as: 
 
ad
ad
T
TT ?
=
max
?   (3-5) 
where T
max
 is the maximum temperature achieved within the flame, and T
ad
 is the 
adiabatic flame temperature.  It is important to note that while temperature overshoot 
is observed in hydrogen-fueled heat recirculating burners, the effect is much stronger 
in methane-fueled burners.   
 Figure 3-8 shows the evolution of the non-dimensional temperature overshoot and 
non-dimensional burning velocity 
frL
L
S
S
,
as a function of non-dimensional passage 
height 
frr
H
,
?
 for a fixed value of the non-dimensional passage length ( 10
,
=
frr
L
?
) 
and under adiabatic conditions. The figure shows that for large values of H, the 
thermal coupling between the gas and the structure is weak, and consequently the 
amount of heat recirculation is small. This leads to a small temperature overshoot as 
well as a burning velocity close to that of a freely-propagating flame.  
 
 
 71 
 
10
-2
10
-1
10
0
10
1
0
5
10
15
20
25
30
35
H / ?
r,fr
S
L
 / S
L,fr
? x 100
 
Figure 3-8. Evolution of normalized flame speed and non-dimensional 
temperature overshoot as a function of normalized passage height. 
 
As the passage height is decreased, the thermal coupling becomes stronger and this 
creates a significant temperature overshoot.  This, in turn, elevates the burning 
velocity.  At some point, however, continuing to reduce the passage height starts to 
decrease the temperature overshoot again.  This is a direct consequence of the fact 
that heat recirculation broadens the reaction zone:  Broadening the reaction zone 
reduces the temperature gradient driving heat from the post-flame to the pre-flame, 
and decreases the overall heat transfer to the reactants.  With less heat being 
transferred upstream, the chemical kinetics run slower again.  However, since the 
flame speed is the reaction zone thickness divided by the overall reaction rate, the 
flame speed remains elevated.  The limiting value of the flame speed as H becomes 
very small is set by the amount of heat that can be conducted through the structure, or 
 
 72 
 
by its thermal conductance.  The net result, however, is that the burning velocity is 
maximized for a certain value of the passage height.  This maximum is a direct result 
of a tradeoff between the strength of the thermal coupling between the gas and the 
structure and the ability to support temperature gradients in the axial (streamwise) 
direction.  Making the passage length shorter magnifies the effect of flame 
broadening on the axial temperature gradients responsible for the peak in flame speed 
and eventually causes the peak to disappear.  This is illustrated in figure 3-9. 
10
-2
10
-1
10
0
10
1
0
5
10
15
20
25
30
35
H / 
?
r,fr
S
L
 / S
L
,fr
L / 
?
r,fr
 = 10
L / 
?
r,fr
 = 5
L / 
?
r,fr
 = 2
L / 
?
r,fr
 = 1
 
Figure 3-9: Evolution of normalized flame speed as a function of normalized 
passage height for different values of the normalized passage length. 
 
Figure 3-10 shows contours of non-dimensional power density as a function of 
normalized passage height and length in the adiabatic case.   
 
 73 
 
1
0
1
0
1
5
1
5
15
15
1
5
1
5
1
5
2
0
20
2
0
20
20
2
0
2
0
2
0
2
5
25
25
2
5
2
5
2
5
2
5
3
0
L / ?
r,fr
H /
 
?
r
,fr
1 2 3 4 5
10
-2
10
-1
10
0
 
Figure 3-10. Adiabatic CH
4
-air combustion in a silicon micro-channel; 
?=0.5, P=1atm.  Contours of non dimensional power density ( )
refDD
ww
,
&&  as a 
function of non-dimensional channel height and length.  
 
Since the flow is adiabatic, the overall efficiency is determined solely by the 
chemical efficiency.  In the range of configurations considered in figure 3-10, the 
chemical efficiency is very close to 1.  A possible explanation for this will be 
presented later in this section.   The figure shows that there are actually two maxima 
in power density (as opposed to only one for H
2
-air combustion). The first maximum 
is global and corresponds to a particular configuration of the combustor, more 
precisely H=0.5mm and L=2.5mm.  It arises from the two phenomena identified 
earlier: the maximum in burning rate at a particular height, and the trade-off between 
compactness and heat recirculation potential as passage length decreases.  The second 
maximum is analogous to the maximum observed in earlier work on adiabatic 
 
 74 
 
hydrogen-air micro-combustors
 
where there is a combustor length that maximizes 
power density for passage heights smaller than a particular value.  It results from the 
elevation of burning rate at small passage heights and the increase in compactness as 
the passage length is reduced. 
Heat loss to the environment imposes important additional constraints on the 
performance of micro-combustors.  Figure 3-11 shows contours of non-dimensional 
power density as a function of normalized passage height and length when the 
structure of the combustor loses heat to the environment with a convective heat 
transfer coefficient of 1 W/m
2
K.  Contours of thermal efficiency are also 
superimposed on the graph.   
Two distinct maxima are apparent.  The first is similar to the first one discussed in 
figure 3-10 as it corresponds to the combined effects of temperature overshoot and 
axial gradient suppression.  The second arises from the fact that heat loss and the 
resulting thermal quenching places a lower bound on the passage height where 
sustained combustion is possible.  This is also analogous to the results obtained in 
previous analyses of hydrogen-air micro-combustors
 
such as those presented in figure 
3-5 of this work.   
 
 75 
 
55
55
1
0
1
0
101010
15
15
1
5
1
5
15
15
1515
20
20
20
2
0
2
0
2
0
20
20
21
21
2
1
2
1
2
1
21
2
1
25
25
25
25
2
5
3
0
505050 606060
808080
85
85
85
9090
90
98
98
98
9
9
9
9
9
9
9
9
L / ?
r,fr
H /
 
?
r,
f
r
1 2 3 4 5
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
 
Figure 3-11: Non-adiabatic CH
4
-air combustion in a silicon micro-channel; 
?=0.5,  h
env
=1 W/m
2
K.  Solid lines: Contours of non dimensional power 
density ( )
refDD
ww
,
&&  as a function of non-dimensional channel height and 
length.  Dashed Lines: Contours of overall efficiency. 
 
The optimum power density configuration occurs at H=0.3mm and L=2mm (i.e. 
H/?
r,fr
=0.3, L/?
r,fr
=2 in non-dimensional units) and is associated with the temperature 
overshoot.  The second maximum is not as beneficial as the first one for two reasons.  
First, the value for power density obtained at the second maximum is lower than that 
obtained at the first one (22 vs. 30).  Second, since the passage height associated with 
the second maximum is lower than that for the first one, heat loss is more significant 
and consequently thermal efficiency is lower (90% vs. 98%).  
 
 76 
 
Figures 3-12 and 3-13 show why there are two optimum configurations for a 
methane combustor whereas there is only one for a hydrogen combustor.  First, the 
temperature overshoot is much higher for methane than for hydrogen.  Second, for 
hydrogen, the flame thickness increases monotonically with decreasing H whereas it 
decreases then increases for methane.   
10
-2
10
-1
10
0
10
1
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
H / 
?
r,fr
?
CH
4
H
2
 
Figure 3-12: Comparison of normalized temperature overshoot for methane 
and hydrogen fuels. 
 
 
 
 77 
 
10
-2
10
-1
10
0
10
1
0
5
10
15
20
25
30
H / ?
r,fr
?
r
 / 
?
r,f
r
H
2
CH
4
 
Figure 3-13: Comparison of the evolution of flame thickness for hydrogen 
and methane fuel as passage height is reduced. 
 
The reasons for this difference remain unclear at this point, but one may speculate 
that it due either to differences between methane?s and hydrogen?s chemistries or 
differences in the thermal properties of the reactive mixtures.   
 
 78 
 
Reaction
T 
 
loss
q
&
in
q
&
gas
q
&
x 
T
i
T
f
Pre-Heat 
T
0 
?
ph
?
r
out
q
&
cond
q&
Structure
 
Figure 3-14:  Thermal exchange in a flame in close contact with a structure.  
gas
q&  and 
cond
q&  represent the heat conducted from the reaction zone to the 
pre-heat zone through the gas and through the structure respectively, 
out
q& and 
in
q&  represent convective heat transfer between the gas and the 
structure, and 
loss
q&  represents the heat lost to the environment. 
 
Nevertheless, the fact that the temperature overshoot is larger for methane than 
hydrogen can be found to be a direct consequence of the non-monotone behavior of 
methane?s flame thickness with decreasing passage height.  Indeed, one of the 
consequences of the fact that flame thickness decreases can be found by computing 
an approximation of the total heat flux from reaction zone to pre-flame in the 
 
 79 
 
situation shown in figure 3-14.  Assuming perfect heat transfer between gas and 
structure, this can be written: 
r
s
s
g
gflameprereaction
x
T
k
x
T
kq
?
1
?
?
?
?
?
?
?=
??
&                           (3-6) 
The above expression shows that this heat flux is to first order inversely 
proportional to flame thickness.  Thus, if flame thickness decreases, the amount of 
heat transferred to the pre-heat zone increases, thereby increasing the amount of pre-
heating and contributing to a large temperature overshoot. 
The second maximum shared by CH
4
 and H
2
 appears at different configurations. 
For methane, it occurs at H=300 ?m and L=2 mm (i.e. H/?
r,fr
=0.3, L/?
r,fr
=2 in non-
dimensional units),  while for hydrogen the optimum configuration is H=5 ?m and 
L=500 ?m (i.e. H/?
r,fr
=0.01, L/?
r,fr
=1 in non-dimensional units). While the optimum 
non-dimensional passage length is similar for the two fuels considered, there is a 
large discrepancy between the optimum non-dimensional passage heights.  This is 
due to the fact that quenching of the hydrogen flame occurs at a much smaller 
passage height than for the methane flame. 
Another notable difference is that for CH
4
, the chemical efficiency appears to be 
approximately constant and equal to 1 (i.e. 100%) over the range of configurations 
considered.  This contrasts with the H
2
 results that showed there were some regions of 
the design space where chemical efficiencies were less than 1.  A possible 
explanation lies in the fact that while the flame thickness is of the same order of 
magnitude for both fuels, the burning velocity is much higher for hydrogen than for 
methane (~1 m/s vs. 7 cm/s).  As a result, for a given configuration, the residence 
time in a micro-channel where a methane flame is stabilized is much larger than in a 
 
 80 
 
channel where a hydrogen flame is stabilized.  If the length of the combustor is 
reduced beyond the range covered in Fig. 3-11, the chemical efficiency eventually 
does decrease as illustrated in Fig. 3-15 but it does so at much smaller passage lengths 
than one would expect for hydrogen.  
 
10
-2
10
-1
10
0
10
1
92
93
94
95
96
97
98
99
100
L / ?
r,fr
?
c
hem
H / 
?
f,fr
adiabatic
 
Figure 3-15: Evolution of chemical efficiency as the passage length is reduced 
at fixed passage height.  The fuel is methane at equivalence ratio of 0.5. 
 
Taken together, these results indicate that there may be some advantages to using 
methane as opposed to hydrogen in micro-combustors.  First, since an optimum 
power density configuration is realized at a relatively large combustor size, optimized 
methane micro-combustors ought to be easier to construct and have smaller pressure 
losses. Second, the overall efficiency associated with this configuration is relatively 
H/?
r,fr
=0.5 
Adiabatic 
 
 81 
 
high so that one is not forced to trade chemical and thermal efficiency for power 
density as in hydrogen-powered micro-combustors.  Third, the shapes of the 
efficiency contours indicate that it should be possible to construct smaller methane 
combustors that still have high chemical efficiency.  The disadvantage, however, is 
that the maximum improvement in power density that can be achieved via 
miniaturization of a methane combustor is about half that of a hydrogen combustor:  
Peak normalized power density is on the order of 30 for methane as opposed to 60 for 
hydrogen. 
 
3.4. Conclusion 
A one-dimensional numerical model for reacting flow in a micro-channel has 
been developed to identify optimum configurations for micro-combustors.  The 
model shows that as the size of a combustor is reduced, axial heat transfer through the 
structure from the post-flame to the pre-flame plays an increasingly important role in 
determining its performance.  While axial heat transfer tends to broaden the reaction 
zone, it also increases the burning rate by pre-heating the incoming reactants.  Since 
the increase in reaction rate is larger than the increase in reaction zone thickness, the 
net result is an increase in power density.  Axial heat transfer through the structure 
also provides a mechanism for flame stabilization and the maximum power density 
for a particular combustor design (values of L and H) occurs at the blow-off limit.  
Heat transfer from the structure to the environment places a lower limit on the 
combustor volume that tends to decrease with increasing pressure.  The combustor 
 
 82 
 
design that maximizes power density generally does not maximize efficiency.  
However, operation at elevated pressures tends to improve overall efficiency by 
increasing the burning rate and decreasing the reaction zone thickness. Similar results 
were obtained when methane fuel was considered.  However, the presence of a 
second optimum configuration in power density (at larger dimensions) suggests the 
importance of chemistry in the optimization process.  In summary, axial heat transfer 
through the structure, in addition to heat loss to the environment, has important 
effects on the stability and performance of micro-combustors and should be included 
when estimating micro-combustor performance. 
 
 
 83 
 
Chapter 4:  Micro-Thruster Design and Performance 
 
4.1. Introduction 
The results presented in the previous chapters have focused on an idealized 
problem - a parallel plate micro-combustor.  They suggested that the configuration of 
a micro-combustor could be optimized for maximum power density, but at the price 
of lower combustion efficiency.  The ultimate objective of the research project, 
however, is to evaluate the effect of structural heat conduction on the design and 
performance of micro-thrusters.  Therefore, the objective of this chapter is to apply 
what we have learned from the study of the model micro-combustor problem to the 
design of a micro-thruster.   
To do so the following approach is used:  The idea is to couple an arbitrary 
geometry micro-nozzle to the idealized micro-combustor to simulate a complete 
micro-thruster.  This is accomplished by using the results obtained in the previous 
chapter of this work to determine the flow conditions exiting the micro-thruster?s 
combustion chamber, and the passing those results on to a simplified numerical 
model for the compressible flow within the micro-nozzle.  The strategy is illustrated 
in figure 4-1. 
 
 
 84 
 
 
Figure 4-1: Illustration of the strategy for estimating the performance of a 
micro-thruster. 
 
4.2. Model for the micro-nozzle 
For the purpose of this analysis, the arbitrary geometry illustrated in figure 4-2 is 
used to represent a planar converging-diverging nozzle.  Since the geometry is pnar, 
the flow is assumed to be quasi-1D and the device is assumed to have unit width.  The 
nozzle has an entrance height of H
0
, a throat height of rH
0
 and an exit height of 2H
0
.  
The thickness of the walls of the nozzle is H
0
/2 to match the similar assumption made 
for the micro-combustor considered in Chapter 3.  The length of the nozzle is L
N
=2H
0
 
and the throat is located at L
N
/4.  H(x) is the x-varying nozzle height.  Since the 
problem considered is quasi-1D, H(x) defines the nozzle area A(x). 
H
Outlet flow conditions of combustion chamber 
 ?    Inlet conditions of micro-nozzle 
Micro-thruster outlet
Computed using micro-
nozzle numerical code 
Combustion chamber 
Reactants 
(entering at S
L
)
L
Micro-nozzle 
 
 85 
 
 
 
 
 
Figure 4-2: Schematic diagram of a model micro-nozzle.   
 
The results reported in the previous chapter of this work are used to determine the 
conditions at the nozzle entrance.  The exit temperature and velocity from the 
combustor define the Mach number of the flow entering the nozzle: 
i
i
i
RT
u
M
?
=                                                   (4-1) 
Figure 4-3 shows the contours of Mach number as a function of L and H. 
rH
0
t=H
0
/2 
M=M
i
 
T=T
i
 
? = ?
i
 
2.H
0
x 
2 H
0 
H
0
 / 4 
0 
H
0
 
 
 86 
 
0
.
0
1
0
.
0
1
0
.
0
1
0.01
0.01 0.01
0.020.020.02
0
.
0
2
0
.
0
2
0
.
0
2
0
.
0
2
0.05
0.05
0
.
0
5
0
.0
5
0
.
0
5
0
.
0
7
L / ?
r,fr
H /
 
?
r
,fr
10
0
10
-2
10
-1
10
0
 
Figure 4-3:  Contours of nozzle inlet Mach number as a function of 
combustor configuration. 
 
The required value of r (throat to inlet height ratio) necessary to achieve a 
subsonic-supersonic isentropic flow is determined using the data presented in figure 
4-3 and the well-known area-Mach number relationship
54
: 
 
5.0
2
i
2
i
M
2
1
1
1
2
M
1
r
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+
+
=
?
?
                             (4-2) 
Figure 4-4 shows contours of r as a function of combustor configuration. 
 
 87 
 
0
.
0
2
0
.
0
2
0
.
0
2
0
.
0
2
0.02 0.02
0.0
0.050.050.05
0
.
0
5
0
.
0
5
0
.0
5
0.1
0.1
0
.1
0
.1
L / ?
r,fr
H /
 
?
r
,fr
10
0
10
-2
10
-1
10
0
 
Figure 4-4: Contours of r, the required throat to inlet height ratio to achieve 
a subsonic-supersonic isentropic flow as a function of combustor 
configuration. 
 
Once the value of r required to achieve proper subsonic to supersonic transition is 
determined, the geometry of the micro-nozzle is set and the area-Mach number 
relationship is used to determine the evolution of M with downstream distance.  
Knowledge of M and the assumption that the flow in the nozzle is isentropic allows 
the flow properties T, r, p and U to be determined using the following relationships: 
1
2
2
1
1
?
?
?
?
?
?
? ?
+= MTT
i
?
                                             (4-3) 
 
 88 
 
1
1
2
2
1
1
?
?
?
?
?
?
?
?
?
+=
?
?
?? M
i
                                         (4-4) 
1
2
2
1
1
?
?
?
?
?
?
?
?
?
+=
?
?
?
Mpp
i
                                        (4-5) 
RTMu ?=                                                        (4-6) 
The thrust produced per unit width is: 
()
ie
uumT ?= &                                                   (4-7) 
where  
ee
Hum 2?=&                                                    (4-8) 
However, a more appropriate performance metric is thrust-to-weight ratio.  To 
evaluate this, the mass of the thruster is estimated using geometric considerations.  
The thruster mass per unit width is: 
)32(
2
HHm
Si
+= ?                                               (4-9) 
where r
Si
 is the mass density of silicon (2330 kg/m
3
).  The thrust-to-weight ratio 
becomes: 
mg
T
WT =/                                                     (4-10) 
Another important performance metric is the specific impulse.  It is defined as:  
g
u
I
e
sp
=
                                                       (4-11) 
Another important performance characteristic of micro-scale systems is pressure 
losses.  At small scales, wall velocity gradients are large, and these lead to large 
pressure loss.   In order to estimate these losses, we assume that a fully-developped 
 
 89 
 
Poiseuille flow is established in the nozzle.  Thus, the relative pressure loss can be 
evaluated as: 
 
?
= dx
)x(H
)x(u)T,x(8
p
1
p
2
0
?
?                           (4-12) 
where x is the downstream distance, u is the velocity in the nozzle, ? is the 
temperature dependent viscosity of the gas (which is assumed to be N
2
), H(x) is the 
distance between the nozzle walls, and p
0
 is the combustion chamber pressure (1atm).  
 
4.3. Results and discussion 
Figure 4-5 is a contour plot of thrust-to-weight ratio as a function of combustor 
passage height H and passage length L.  These dimensions have been normalized by 
?
r,fr
, the flame thickness of a freely propagating flame.  The thrust-to-weight ratio 
T/W has been non-dimensionalized by that of the space shuttle main engine for 
comparison.  This value was evaluated from data found in
55
.  The reactants are H
2
 and 
air at an equivalence ratio of 0.5.  While air is not a realistic rocket motor oxidizer, 
this mixture was selected so that the results from the previous part could be re-used. 
 
 90 
 
11
1
1
111
1
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
4
4
4
4
4
5 5
5
5
6
6
6
7
7
8
8
9
10
L / ?
r,fr
H /
 
?
r
,fr
10
0
10
-2
10
-1
10
0
 
Figure 4-5: Normalized thrust-to-weight ratio as a function of normalized 
passage height and length. 
 
 
Just like power density, figure 4-5 suggests that T/W is maximized for a specific 
configuration of the combustor.  The optimum occurs at the same configuration 
corresponding to optimum power density and illustrates the strong influence 
combustor configuration can have on a micro-thruster?s performance.  The figure also 
suggests that the peak non-dimensional T/W is significantly greater than 1.  This 
indicates that it may be possible to achieve significantly higher T/W using micro-
thrusters than can be achieved in macro-scale devices.  Since the actual fuel used in 
the reference engine (i.e. the space shuttle main engine) is high-pressure liquid 
 
 91 
 
hydrogen and oxygen which is much more energetic than the atmospheric pressure 
H
2
-air  mixture considered in this study we might expect the maximum value of T/W 
to be much higher than 10 for realistic rocket propellants.  Consequently the benefits 
of miniaturizing chemical thrusters could be even greater than what is described here. 
While thrust-to-weight ratio is an important performance metric, specific impulse 
is also important as it reflects the efficiency of the combustion process.  The solid 
lines in figure 4-6 show contours of specific impulse as a function of normalized 
passage height H and length L.  The dotted lines are contours of thrust-to weight ratio 
which are superimposed on the plot for reference.  The figure shows that at the 
optimum configuration with respect to T/W, the specific impulse is on the order of 
140s.  At large H and L where there is essentially no interaction between the reacting 
gas and the structure of the combustor, the specific impulse is larger than 180s.  This 
suggests that, just like for the results presented in the previous chapter, it is possible 
to design micro-thrusters with large thrust-to-weight ratios, but that this comes at the 
price of specific impulse.  However, as mentioned earlier, specific impulse would 
have been larger had we considered a more realistic oxidized for rockets.  
 
 92 
 
11
1
1
111
1
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
4
4
4
4
4
5 5
5
5
6
6
6
7
7
8
8
9
10
L / ?
r,fr
H /
 
?
r
,fr
150150
150
60
1
6
0
1
6
0
160
1
7
0
1
7
0
1
7
0
1
7
0
1
8
0
1
8
0
10
0
10
-2
10
-1
10
0
 
Figure 4-6: Solid lines: Specific impulse as a function of normalized passage 
height and length; Dotted lines: Normalized thrust-to-weight ratio as a 
function of normalized passage height and length. 
 
One key limiter of miniaturization in these high-speed devices is pressure drop.  
Figure 4-7 is a contour plot of normalized pressure drop, defined by equation 4-12, as 
a function normalized passage height and length.  Again, contours of normalized 
trust-to-weight ratio have been superimposed as dotted lines for the purpose of 
comparison.  As the passage height decreases, the velocity gradient at the wall 
increases and thus so does the pressure drop.  The variation with passage length is 
more complex.  The figure shows that as the length of the combustion chamber 
 
 93 
 
decreases, the conditions at the outlet of the combustor change: The combustion 
temperature decreases because the degree of conversion of the fuel decreases, and the 
flow velocity decreases because there is less pre-heating of the mixture via heat 
recirculation through the combustor structure.  This leads to a decrease in Mach 
number at the outlet of the combustor (or the inlet of the nozzle) with L.  Thus, the 
value of r required to reach sonic conditions at the throat increases, leading to a 
smaller throat height and larger pressure losses.  Therefore, the least amount of 
pressure drop occurs at the largest passage height and length.  As the configuration 
moves away from this optimum, pressure drop becomes significant. Beyond a 
normalized pressure drop of 100%, the device becomes totally unviable but this 
would really occur much earlier.  The figure suggests that the range of possible 
configurations associated with relative pressure drops that are less than 100% is 
actually quite narrow.  This limits the maximum achievable normalized thrust to 
weight ratio to 5.  However the combustion chamber pressure considered in this study 
is 1atm which is very low for a rocket.  Considering a higher pressure would reduce 
the significance of the pressure losses and expand the range of viable configurations. 
 
 94 
 
11
1
1
111
1
2
2
2
2
2
2
2
2
3
3
3
3
3
3
3
3
4
4
4
4
4
5 5
5
5
6
6
6
7
7
8
8
9
10
L / ?
r,fr
H /
 
?
r
,fr
2
0
2
0
5
0
5
0
50 50
0
0
100
100 100
10
0
10
-2
10
-1
10
0
 
Figure 4-7: Solid lines: Relative pressure drop (in %) as a function of 
normalized passage height and length; Dotted lines: Normalized thrust-to-
weight ratio as a function of normalized passage height and length. 
 
4.4. Conclusions 
The performance of hydrogen-air fueled micro-thrusters was investigated by 
appending a planar, ideal nozzle to the parallel plate micro-combustor geometry 
investigated in the previous chapter. This simulation took into account heat 
recirculation within the structure of the combustor.  The results suggest that the 
combustion chamber configuration strongly affects the performance of the micro-
 
 95 
 
thruster.  In particular, there exists a combustion chamber configuration which 
maximizes thrust-to weight ratio.   However, the optimum configuration has low 
specific impulse and high pressure loss, making it unlikely that the ?real? device 
could actually achieve the predicted T/W and I
sp
.  However, the results suggest that 
thrust-to-weight ratios much greater than those of macro-scale thrusters could be 
achievable, provided that good thermal management of the combustion chamber is 
achieved and that the devices can be operated at high pressures.  This demonstrates 
the potential benefit of miniaturizing chemical thrusters 
However, more work is required in order to evaluate the performance of more 
realistic thrusters.  In particular, higher combustion chamber pressures and more 
realistic propellants ought to be considered.  These more realistic thrusters should 
perform better than the ones described here.  In addition, the effect of the thermal 
coupling between gas and structure of the micro-nozzle itself should be taken into 
account to obtain more realistic performance estimates.   
Nevertheless, this chapter does provide some preliminary insight regarding the 
design and performance of micro-thrusters and the potential benefits of miniaturizing 
chemical thrusters.  However, maximizing T/W would only be beneficial if the 
acceleration of the spacecraft considered is the driving design parameter.  If, on the 
other hand, the lifetime/efficiency of the spacecraft is of greater importance, then 
conventional scale thrusters are more advantageous as they have greater specific 
impulse.  In general, the optimum thruster configuration ultimately depends on the 
particular application. 
 
 96 
 
Chapter 5:  Flame Stability 
 
5.1. Introduction 
In the second chapter of this work, numerical and analytical models were developed 
for flames propagating in micro-channels.  They showed that flames stabilized in 
micro-channels behaved differently than their conventional-scale counterparts: 
Reaction zones seemed to be broader and burning velocities seemed to be higher than 
what one would expect for freely propagating flames.  It was shown that these 
behaviors are the result of heat recirculation from the post-flame to the pre-flame 
through a micro-combustor?s structure that becomes more important as the size of the 
combustor is reduced.  Some of the consequences of this behavior were explored in 
the third and fourth chapters of the thesis, which concentrated on predicting the 
performance of micro-combustors and micro-chemical thrusters.  The flame 
broadening phenomenon and the increase in burning velocity of the mixture were 
shown to have a great impact on power density of micro-combustors and the thrust-
to-weight ratio of micro-thrusters.  Brief mention was also made of the fact that 
thermal coupling also seemed to widen the stability limits of the flame:  Figure 3-2 
showed that there seems to be a range of velocities over which stable combustion can 
be sustained.  This observation is consistent with early analytical work which showed 
that the interaction between a flame and its surroundings can enhance its stability
56
.  
Subsequent experimental work confirmed this result, and the process was given the 
name ?Thermally Stabilized Combustion?
57
.  This contrasts however with other 
 
 97 
 
analytical theories for freely propagating flames that predict a single value of the 
burning velocity that is an eigenvalue of the gas-phase problem.   
In this chapter, we attempt to improve our understanding of the basic physical 
mechanisms responsible for the enhancement of flame stability in micro-combustors.  
This is accomplished by developing a simple analytical model for flame propagation 
in a narrow channel and comparing its predictions to those obtained using more 
detailed numerical simulations and experiments.  Both steady-state and transient 
stability are investigated.  Note that combustion instability associated with the 
coupling of pressure and heat release oscillations is not under consideration here. 
 
5.2. Analytical Model 
Figure 5-1 is an idealized representation of a flame stabilized in a micro-channel 
formed by two closely-spaced parallel plates that extend to infinity in the directions 
normal to the plane of the drawing.  The distance between the plates is H and their 
thickness t is H/2. This ensures that the respective cross-sectional areas for heat 
conduction in the gas and the structure remain constant as the distance between the 
plates is varied. The length of the micro-channel is L and the position of the flame 
from the flow entrance is x
flame
. Premixed hydrogen fuel and air enter the combustor 
from the left at an equivalence ratio ? =0.5, pressure P=1 atm and velocity U.  Since 
the plates are semi-infinite, the width of the channel does not have any effect on the 
problem and it is taken to be unity.  
 
 98 
 
 
Figure 5-1: Schematic diagram of a parallel plate reactor in which a flame is 
stabilized showing structure and gas temperatures, the three zones described 
in the text and the heat flux from the structure to the preheat zone. 
 
A one-dimensional approach is used to model exchange processes in the flow 
direction.  The channel is divided into three parts: the preheat zone, in which the gas 
temperature is raised from the inlet temperature T
0
 to a preheat temperature T
ph
, the 
reaction zone in which gas temperature increases sharply from the preheat 
temperature to the flame temperature T
ad
, and the post-flame zone, where the gas 
temperature decreases because heat is lost to the passage walls and possibly to the 
environment. The final exit temperature is T
f
 which can be written using energy 
conservation as 
T
ph 
T
0 
0 L x
flame 
x 
 
Q
s 
T
f 
Gas temperature 
Structure temperature 
H 
t=H/2 
Pre-heat region Post-flame region 
T
ad 
Reaction zone 
 
 99 
 
M
Q
T
HUC
Q
TT
loss
ad
p
loss
adf
?=?=
?
                           (5-1) 
where T
ad
 is the adiabatic flame temperature (1680K in this study) and 
loss
Q  is the 
heat loss to the environment per unit width of the plate. 
It is possible to represent thermal exchange processes between gases and solid 
structures using thermal resistance networks
58
.  Figure 5-2 shows the network used to 
represent thermal exchange between the gas, the structure, and the environment in 
this problem.  Note that it is only necessary to represent the upper half of the micro-
channel.   
 
 
Figure 5-2: Thermal resistance network associated with the flame 
stabilization problem. 
 
T
0 
Q
s 
T
0 
T
s1 
T
s2 
T
f 
L-x
flame 
x
flame 
T
0 
R
b2 
R
b1 
R
s 
0 L 
x
flame 
H/2 
R
e1 
R
e2 
Post-flame region 
Pre-heat region 
R
r 
L
 
 100 
 
In the network, R
s
 is the resistance to conductive heat transfer in the combustor 
structure, R
b1
 and R
b2
 are the resistances to convective heat transfer in the boundary 
layers in the pre-heat and post-flame regions respectively and, R
e1
 and R
e2
 are the 
resistances to convective heat transfer from the structure to the environment in the 
pre-heat and post-flame regions respectively. These resistances can be written as: 
()
()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
==
==
?
==
==
==
==
flameeTpreeT
e
flameeTpreeT
e
flamepostgpostpostT
b
flamepregprepreT
b
s
r
gpreg
r
r
sss
s
xLhAh
R
xhAh
R
xLkNu
H
Ah
R
xkNu
H
Ah
R
HkAk
R
tk
L
Ak
L
R
,,
2
,,
1
,,
2
,,
1
,
11
11
21
21
??
                  (5-2) 
where k
s
, k
g,pre
, and k
g,post 
are respectively the thermal conductivities of the structure, 
the gas in the pre-heat region, and the gas in the post-flame region.  A
s
, A
g,pre
, and 
A
g,post 
are the corresponding cross-sectional areas for conduction, Nu is the Nusselt 
number for heat transfer between the gas and the structure, and ?
r
 is the flame 
thickness.  
While it has been shown earlier that ?
r
 can increase due to structural heat 
conduction, in this work we assume that ?
r
 is fixed to 0.5mm ? the value for a freely 
propagating flame.  This is justified for two reasons.  First, we are most concerned 
with the position of the flame which is taken to be the middle of the reaction zone in 
this work.  Second, the thermal resistance of the gas column is already so much larger 
 
 101 
 
than that of the structure and boundary layers that it passes a small portion of the total 
heat flux in the system.  As a result, any changes in the thermal resistance of the gas 
column that arise from changes in the reaction zone thickness should have a 
negligible impact on the solution of the thermal network problem. Third, flame 
broadening is not expected to occur at the lengths scales considered here. 
The structure temperatures T
s1
 and T
s2
 are computed by performing heat flux 
balances at the nodes of the thermal resistance network.  This leads to: 
C
BTAT
R
RR
1
RR
RR
T
RR
R
R
1
R
1
RT
T
1f10
2
s
LR
S2b
LR
f
2es
R
1b1e
L0
1s
+
=
?
+
?
?
?
?
?
?
++
=           (5-3) 
C
BTAT
R
RR
1
R
R
T
RR
R
RR
R
R
1
RT
T
2f20
2
s
LR
2b
R
f
s1b
L
s1e
L
2e
R0
2s
+
=
?
+
?
?
?
?
?
?
++
=          (5-4) 
where 
     
s2b2e
R
R
1
R
1
R
1
1
R
++
=                                 (5-5)   
s1b1e
L
R
1
R
1
R
1
1
R
++
=                                      (5-6) 
2es
R
1b1e
1
RR
R
R
1
R
1
A ++=                                  (5-7) 
 
s1b
L
s1e
L
2e
2
RR
R
RR
R
R
1
A ++=                              (5-8) 
S2b
LR
1
RR
RR
B =                                       (5-9) 
2b
R
2
R
R
B =                                         (5-10) 
 
 102 
 
2
1
s
LR
R
RR
C ?=                                      (5-11) 
The heat loss to the environment is written as: 
()()( )[ ]
flame02sflame01se,Tloss
xLTTxTTh2Q ??+?=
     (5-12)       
where h
T,e 
is the heat transfer coefficient between the structure and the environment. 
We now have four equations (1, 3, 4 and 12) with 4 unknowns (T
f
, T
s1
, T
s2
, and Q
loss
).  
A closed form solution can be found for T
f
: 
M
B
1
M
A
TT
T
0ad
f
+
?
=                                               (5-13) 
where 
2e1e2e
2
1e
1
R
1
R
1
CR
A
CR
A
A ??+=
                                     (5-14) 
2e
2
1e
1
CR
B
CR
B
B +=                                                (5-15) 
Note that T
f
 depends on x
flame
 through the different resistances R
e1
, R
e2
, R
b1
 and R
b2
.   
Now that T
f
 is known, T
s1
, T
s2
 and Q
loss
 can subsequently be evaluated using equations 
5-3, 5-4 and 5-12.  An additional constraint is imposed by requiring that the heat 
generated by combustion must be higher than the heat loss to the environment.  This 
is written: 
loss0fpL
Q )T-H(TCS >?                                    (5-16) 
The net heat flux from the post-flame to the pre-heat zone is determined using: 
r
0f
1b
01s
s
R
TT
R
TT
2Q
?
+
?
=                                  (5-17) 
 
 103 
 
Assuming that heat transfer from the structure to the pre-heat zone occurs uniformly, 
the heat flux per unit area can be written as: 
flame
s
s
x
Q
q =                                               (5-18) 
If the heat flux is uniform, it can be shown that the temperature rises linearly in the 
preheat region
59
 so 
flame0ph
x
dx
dT
TT +=                                    (5-19) 
where  
HCU
q
dx
dT
p
s
?
=                                         (5-20) 
Heat recirculation increases the unburned (or preheat) gas temperature by transferring 
a certain amount of enthalpy from the post-flame to the pre-heat zone.  This, in turn, 
increases the burning velocity of the mixture.  At the same time, however, heat loss to 
the environment decreases the flame temperature and thus the burning rate. To 
account for these dependencies, we use the following simplified expression for the 
flame speed that includes the effects of pre-heating
60
: 
?
?
?
?
?
??
?
?
f
A
2/n
fph
375.0
L
RT2
E
expTTTS
r
                   (5-21) 
where 
 
2
TT
T
fph
+
=                                           (5-22) 
The global order of the reaction, n
r
, is 3 for the hydrogen-air reaction considered here 
and the activation energy E
A
 is 12.5kcal/mole
61
. 
 
 104 
 
In order to investigate the stability limits of the combustor, we look for values of 
the flame position x
flame
 that lead to burning velocities that are equal to the inlet 
velocity.   Figure 5-3 shows a sample plot of non-dimensional burning velocity for 
hydrogen fuel at ? =0.5 computed using equation (5-21) as a function of non-
dimensional flame position for H=0.5mm, L=1cm, U=1.8m/s, h
T,e
=50 W/m
2
K and for 
three different values of the thermal conductivity of the structure (0 W/mK, 5 W/mK, 
and 50 W/mK). In the figure, the flame position is non-dimensionalized by the 
combustor length and the flame speed is non-dimensionalized by the speed of a freely 
propagating flame ? i.e. one that propagates without exchanging heat with the 
structure (S
L,fr
=0.9 m/s). 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
x
flame
 / L
S
L
 / 
S
L,
f
r
0
5
50
Stable
Unstable
x
1
x
2
S
L
=U=1.8m/s
 
Figure 5-3: Burning velocity as a function of flame position for different 
values of the thermal conductivity of the structure in W/mK. The families of 
curves correspond to solutions of the thermal network where h
Te
=50 W/m
2
K 
and U=1.8 m/s. 
 
 105 
 
 
Figure 5-3 shows that if the thermal conductivity of the structure is high enough, the 
horizontal dashed line corresponding to the inlet velocity of 1.8 m/s intersects the plot 
of S
L
 vs. x
flame
 in two places indicating that two solutions for the flame position are 
possible.  Of these, only the upstream one is stable: This can be seen by realizing that 
if x
flame
 increases at constant U, the burning velocity increases.  Since the burning 
velocity now exceeds the inlet velocity U, the flame tends to move upstream thereby 
returning to its original position.  Similar lines of reasoning show that the upstream 
solution is also stable if x
flame
 decreases at constant U, and that the downstream 
solution is unstable.  It should be noted that this stabilization mechanism is similar to 
that found in Bunsen burner flames
62
.  Figure 5-3 also shows that there is a minimum 
threshold value for the structure thermal conductivity below which no solution is 
possible.  This illustrates the stabilizing role played by structural heat conduction.  
When the thermal conductivity of the structure is zero and in the absence of heat loss, 
the non-dimensional flame speed is 1 and the only possible inlet velocity is the freely 
propagating flame speed.  This is consistent with earlier theories for flame 
propagation that do not account for heat transfer to the structure.  
While the behavior described above is observed for inlet velocities larger than the 
freely propagating flame speed S
L,fr
, a rather different behavior is observed at small 
velocities.  This is because as the inlet velocity is reduced the importance of heat loss 
becomes more significant.  Figure 5-4 is similar to Figure 5-3 in that it shows a 
sample plot of non-dimensional burning velocity as a function of non-dimensional 
flame position for H=0.5mm, L=1cm, h
T,e
=50 W/m
2
K and for three different values of 
 
 106 
 
the thermal conductivity of the structure (0 W/mK, 5 W/mK, and 50 W/mK).  
However the inlet velocity considered is now 0.8m/s, compared to 1.8m/s in Figure 5-
3.   
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
x
flame
 / L
S
L
 / S
L
,fr
0
5
50
S
L
=U=0.8m/s
No flame due to heat loss
 
Figure 5-4: Burning velocity as a function of flame position for different 
values of the thermal conductivity of the structure. The families of curves 
correspond to solutions of the thermal network where h
Te
=50 W/m
2
K and 
U=0.8 m/s. 
 
Figure 5-4 shows that there is a domain of x
flame
 for which no solutions of the set of 
equations described above is possible because of the heat loss constraint expressed in 
equation (5-16).  For small values of x
flame
, there is a large part of the structure 
downstream the flame and whose temperature is relatively high.  Thus heat loss is 
significant.  The figure also shows however that a stable solution does exist if the 
 
 107 
 
thermal conductivity of the structure is small (0 or 5 W/mK).  This contrasts with the 
sample case described above where the inlet velocity is high.   
Finally, the degree of heat recirculation is controlled in this work by varying the 
thermal conductivity of the structure.  The reason that thermal conductivity is chosen 
as the design parameter is that it is convenient both for understanding how the choice 
of combustor structural material affects flame stability and for comparing our results 
to those of others.  However, it should be reiterated that this is also somewhat 
misleading as the degree of heat recirculation also depends on the relative thickness 
of the structural material to the gas column.  As a result, the conductance ? which is 
the product of the thermal conductivity and the cross-sectional area for conduction in 
the structure ? is a more physically appropriate parameter to vary in order to 
understand the effect of heat conduction in the structure on the stability of the flame. 
 
5.3. Numerical Model 
The numerical model introduced in Chapter 2 of the thesis can be used to check the 
predictions of the analytical model regarding steady-state stability.  In addition the 
transient capability of the code allows the study of unsteady stability, i.e. the response 
of the flame to external time-dependant perturbations.  For the purposes of this 
analysis, flame position in the numerical model is defined as the axial position when 
the H radical mass fraction peaks.  This is justified by the fact that the presence of a 
large H radical concentration typically indicates that combustion is taking place 
within the mixture. 
 
 
 108 
 
5.4. Steady Stability: Results and Discussion 
5.4.1. Analytical model results 
Figure 5-5 shows the evolution of the normalized flame position as a function of 
the normalized inlet velocity computed using the analytical model presented in 
section 5.2..  In the figure, the flame position is non-dimensionalized by the 
combustor length L and the burning velocity is non-dimensionalized by the freely 
propagating flame speed S
L,fr
=0.9 m/s. 
0.5 1 1.5 2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
U / S
L,fr
F
l
am
e pos
i
t
i
on /
 L
k
s
=10W/mK
k
s
=50W/mK
H=1mm, L=1cm, h
T,e
=50W/m
2
K
Low Vel. BO
High Vel. BO
 
Figure 5-5:  Flame position as a function of inlet velocity for different values 
of the thermal conductivity of the structure computed using the analytical 
model.  The fuel considered is hydrogen at ? =0.5. 
 
 
 109 
 
The results in the figure correspond to H=1mm, L=1cm and h
T,e
=50W/m
2
K and to 
two different values of the thermal conductivity of the structure (10 W/mK and 50 
W/mK).  The figure shows that as the inlet velocity increases, the flame 
monotonously moves farther downstream the combustor, as intuitively expected.  The 
U-shaped curves show that there are two stability limits corresponding to high and 
low-velocity blow-off.  The high velocity blow-off limit occurs when heat 
recirculation through the structure is unable to provide adequate preheating of the 
mixture to yield a burning velocity equal to the inlet velocity.    
The low velocity blow-off limit is a consequence of the fact that the total heat 
release in the combustor is directly proportional to the inlet velocity whereas the heat 
loss to the environment scales with the structure temperature which remains 
approximately constant as the inlet velocity is reduced.  As a result, reducing the inlet 
velocity U reduces the amount of heat available for pre-heating and this, in turn, 
reduces the burning velocity.  At some point, U exceeds the burning velocity for all 
possible flame positions and the flame blows off.  
Figure 5-5 also illustrates the effect of changing k
s
, the thermal conductivity of the 
structure.  Increasing k
s
 shifts the stability limits to the right or to higher velocities.  
This makes sense because a larger structure thermal conductivity allows for greater 
heat recirculation from the post-flame to the pre-flame and thus facilitates the pre-
heat of incoming reactants.  This observation is consistent with the effect of k
s
 
described in Chapter 2 of this work.  However, a higher structural thermal 
conductivity also provides less thermal insulation by allowing less structure 
temperature to reach a high temperature, thus increasing heat loss to the environment 
and causing the flame?s low velocity blow-off limit to also shift to higher velocities. 
 
 110 
 
Figure 5-6 illustrates the general effect of structural thermal conductivity on the 
stability limits of the flame.   
10
-2
10
-1
10
0
10
1
10
2
10
-0.3
10
-0.2
10
-0.1
10
0
10
0.1
10
0.2
10
0.3
k
s
 (W/mK)
U /
 
S
L
,fr
Low Vel. BO
High Vel. BO
H=1mm, L=1cm, h
T,e
=50W/m
2
K
 
Figure 5-6: Effect of the thermal conductivity of the structure on the stability 
limits of the flame for H=1mm and L=1cm computed using the analytical 
model. 
 
The figure shows that for the reasons identified above increasing the thermal 
conductivity of the structure shifts both the high and low-velocity blow-off limits to 
higher velocities.  The change in burning velocity shows a limit-cycle behavior as k
s
 
becomes very large or very small.  The explanation of this behavior is similar to that 
of Figure 2-7:  At large thermal conductivities of the structure, thermal feedback from 
post-flame to pre-flame becomes limited by convective heat transfer from gas to the 
 
 111 
 
structure, i.e. by resistances R
b1
 and R
b2
 in figure 5-2.  Conversely, for small k
s
 there 
is no more thermal feedback through the structure and the behavior of the flame is 
governed entirely by heat loss to the environment.   
Figure 5-7 illustrates the effect of reducing the length of the combustor on the low 
and high-velocity stability limits.  
10
-1
10
0
10
1
10
2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
k
s
 (W/mK)
U /
 
S
L
,fr
Low-Vel BO
High-Vel BO
Low-Vel BO
High-Vel BO
}
H=1mm, L=1cm
H=1mm, L=5mm
}
 
Figure 5-7:  Effect of combustor length on the stability limits of the flame 
computed using the analytical model. 
 
 It shows that decreasing the length of the combustor shifts both blow-off limits to 
lower velocities when the thermal conductivity of the structure is relatively high.  
This makes sense because decreasing the combustor?s length reduces the passage 
length available for preheat.  When the thermal conductivity of the structure is 
relatively low, heat recirculation is less important and heat loss is more important as 
 
 112 
 
the structure becomes a heat sink for the flame.  A longer structure increases heat loss 
to the environment so decreasing the structure?s length reduces heat loss and raises 
the blow-off limits. 
Figure 5-8 illustrates the effect of combustor height on the stability limits of the 
flame.  It shows that decreasing the height of the combustor raises both stability 
limits.  A smaller gap between the plates enhances thermal coupling and preheating of 
the mixture, thereby increasing the high-velocity blow-off limit.  On the other side, a 
smaller gap also increases the relative importance of heat loss to the environment and 
as a consequence the minimum velocity required to obtain a stable flame, or the low-
velocity blow-off limit, also increases. 
10
-1
10
0
10
1
10
2
0.5
1
1.5
2
2.5
3
k
s
 (W/mK)
U /
 S
L
,fr
Low-Vel BO
High-Vel BO
Low-Vel BO
High-Vel BO
}
}
H=1mm, L=1cm
H=0.5mm, L=1cm
 
Figure 5-8:  Effect of combustor height on the stability limits of the flame 
computed using the analytical model. 
 
 113 
 
5.4.2. Numerical model results 
Figure 5-9 is similar to figure 5-5, except that the curves were generated using the 
1D numerical model described earlier.  The conditions considered are the same as in 
the analytical model (H=1mm, L=1cm,  h
T,e
=50W/m
2
K, k
s
=10 W/mK or 50 W/mK) 
0 0.5 1 1.5 2 2.5 3
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
U / S
L,f r
F
l
am
e pos
i
t
i
on /
 
L
10
50
 
Figure 5-9:  Flame position as a function of inlet velocity for different values 
of the thermal conductivity of the structure computed using the numerical 
model. 
 
The figure shows qualitative behaviors that are similar to those predicted by the 
simple analytical model.  Both the model and the simulation predict similar shaped 
curves that shift in a similar way when the thermal conductivity of the structure is 
increased.  However, there is a discrepancy near the low-velocity blow-off limit 
Slow-burning 
Fast-burning 
 
 114 
 
where the numerical model shows that the flame position increases slightly as U 
decreases and approaches the low-velocity blow-off limit.  This indicates that there 
are two distinct values of the inlet velocity U that are associated with a single flame 
position, as illustrated in figure 5-9.  These two flame regimes will be called the ?fast-
burning? and ?slow-burning? regimes, in order to be consistent with the work of other 
researchers
63
.  This behavior is not observed in the analytical model which predicts a 
single value for the burning velocity at each flame position.  This value corresponds 
to the fast burning regime.   
Since figure 5-9 was generated using the numerical simulation, temperature 
profiles in the channel are available.  Figure 5-10 shows gas and structure 
temperatures non-dimensionalized by the adiabatic flame temperature of the mixture 
(1646K) as a function of non-dimensional downstream distance x/L.   
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
x / L
T /
 T
ad
T
g
T
sFast-burning regime
Slow-burning regime
 
Figure 5-10:  Gas and structure temperature profiles in the fast-burning and 
slow-burning regimes identified above. 
 
 115 
 
 
The figure shows that both gas and structure temperatures are much higher in the 
fast-burning regime than in the slow-burning regime.  This is a because the inlet 
velocity of the fast-burning regime is much higher than that of the slow-burning 
regime.  This, in turn, means that heat loss to the environment is less significant in 
comparison to the overall heat release rate.   
 Figure 5-11 illustrates the fundamental difference between the fast and slow 
burning regimes.  It is a plot of the axial temperature gradient in the structure non-
dimensionalized by an arbitrary reference temperature gradient (T
ad
-T
0
)/L.   
0 0.2 0.4 0.6 0.8 1
-0.1
-0.05
0
0.05
x / L
[d
T
s
/
d
x
]
 / [(
T
ad
-T
0
)/
L
]
Slow-burning regime
Fast-burning regime
 
Figure 5-11: Profile of normalized structural temperature gradient for the 
fast-burning and slow-burning regimes. 
 
 
 116 
 
The figure shows that the temperature gradient in the structure is positive over a 
much larger portion of the channel in the fast burning regime compared to the slow-
burning regime.  This shows that more heat is transferred to the pre-heat zone of the 
flame in the fast burning regime, and this leads to an increased burning velocity.  
Thus, the main difference between the two regimes is the degree of thermal feedback 
upstream through the structure.  This feedback is small to non-existent in the slow 
burning regime and considerably larger in the fast burning regime.  These 
observations are consistent with those made by other researchers who considered 
moving flames as opposed to stationary flames
64
. 
Figure 5-12 compares the influence of the structure thermal conductivity on the 
stability of the flame computed using the analytical model and the numerical 
simulation for H=1mm and L=1cm.  While the stability limits predicted by the 
numerical simulation are much wider than those predicted by the analytical mode, the 
qualitative behaviors are similar.  This indicates that the basic physical processes 
represented in the analytical model are the principal ones governing flame stability at 
the micro-scale.  Furthermore, the discrepancy between the model and the simulation 
occurs in a direction that should be expected:  The distributed nature of the heat 
transfer between the gas and the structure captured in the numerical simulation should 
increase thermal coupling between the gas and the structure and lead to greater pre-
heating of the reactant mixture than would be possible via the single-zone heat 
transfer processes included in the analytical model.  Therefore, one would expect the 
numerical simulation to produce wider stability limits than predicted by the analytical 
model. 
 
 
 117 
 
10
0
10
1
10
2
0
0.5
1
1.5
2
2.5
3
3.5
4
k
s
 (W/mK)
U /
 S
L
,fr
Low-Vel BO
High-Vel BO
Low-Vel BO
High-Vel BO
}
}
Analytical model
Numerical model
 
Figure 5-12: Comparison of the results of the analytical and numerical 
approaches:  effect of thermal conductivity for H=1cm and L=1cm. 
 
Figure 5-13 illustrates the effect of changing the length of the combustor on the 
stability limits of the flame computed using the numerical simulations.  The results 
are also qualitatively similar to the predictions of the analytical model presented in 
figure 5-7.  
 
 118 
 
10
0
10
1
10
2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
k
s
 (W/mK)
U /
 S
L
,fr
Low-Vel BO
High-Vel BO
Low-Vel BO
High-Vel BO
}
}
H=1mm, L=1cm
H=1mm, L=5mm
 
Figure 5-13:  Effect of combustor length on the stability limits of the flame 
computed using the numerical model. 
 
Figure 5-14 illustrates the effect of reducing the height of the combustor on the 
stability limits of the flame computed using the numerical model.  Again, the 
numerical model is in qualitative agreement with the predictions of the analytical 
model presented in figure 5-8 although the low-velocity blow off limit computed 
using the numerical simulation appears to be much less sensitive to the thermal 
conductivity of the structure. 
 
 119 
 
10
0
10
1
10
2
0
1
2
3
4
5
6
k
s
 (W/mK)
U /
 S
L
,fr
Low-Vel BO
High-Vel BO
Low-Vel BO
High-Vel BO
}
}
H=1mm, L=1cm
H=0.5mm, L=1cm
 
Figure 5-14:  Effect of combustor height on the stability limits of the flame 
computed using the numerical model. 
 
5.4.3. Comparison to results of other researchers 
The topic of flame stabilization via structural heat recirculation has also been 
investigated by Norton and Vlachos
65
 recently.  They considered a propane-air flame 
stabilized in a channel and investigated effects of the thermal conductivity of the 
structure and heat loss to the environment on the stability limits of the flame.  Their 
approach was to use a two-dimensional CFD simulation with one-step chemistry and 
fluid-structure interaction.  Their predictions for the flame position as a function of 
inlet velocity for structures with different thermal conductivities are presented in 
figure 5-15. The shapes of the curves and their trends with increasing the thermal 
 
 120 
 
conductivity of the structure are very similar to our predictions made using the 
analytical model  (figure 5-5) and the 1D numerical simulation (figure 5-9).   
 
Figure 5-15:  Flame position as a function of inlet velocity for different values 
of the thermal conductivity of the structure from
15 
 
Their predictions of how thermal conductivity of the structure generally influences 
the stability limits (presented in figure 5-16) are also similar to the predictions of the 
analytical model (figure 5-6) and the 1D numerical simulation (figure 5-12). 
 
 
 121 
 
 
Figure 5-16:  Effect of the thermal conductivity of the structure on the 
stability limits of the flame from
15 
 
This gives us confidence that the basic physical processes represented in the 
simple analytical model are the ones that govern micro-scale combustion.  This also 
demonstrates the vale of simple modeling in developing fundamental understanding 
of complex combustion processes.  In addition, and as a practical matter, it should be 
noted that Figure 5-5 was generated on a desktop computer (Pentium IV processor, 
2.4GHz, 512Mb RAM) in less than 10 seconds, whereas the results of the 2D 
simulations shown in Figure 5-15 required lengthy computations on a PC cluster 
consisting of 58 Pentium IV processors and 58 GB of RAM.  Therefore, considering 
the gross disparity in complexity between the two approaches, the analytical model 
 
 122 
 
predicts flame stability remarkably well and has the advantage of being fast enough to 
be used as a tool to investigate wide parameter spaces as part of the preliminary 
design process. 
 
5.4.4. Experimental validation 
A simple burner consisting of two parallel silicon plates between which a flame is 
stabilized is used to test the predictions of the model and simulation.  The burner and 
diagnostic technique have been developed by Scott Heatwole
66, 67
 and Anand 
Veeraragavan and only a brief description is provided here.  Figure 5-17 is a 3D 
drawing of the silicon-walled micro-burner.  Methane is the fuel, air is the oxidizer, 
and the equivalence ratio is 0.5.  The distance between the parallel plates is 2 mm, the 
thickness of the plates is 500 ?m and their length is 55 mm.   
The silicon walls of the micro-burner are transparent in the infrared region near 
well known absorption features of CH
4
 and CO
2
.  As a result, it is possible to 
interrogate the flow non-intrusively for these species by passing the beam from a 
Fourier Transform Infra-Red (FTIR) spectrometer through the walls of the burner.  
Methane and CO
2
 spectra are collected from which flame position and flame 
temperature are computed using methods reported elsewhere
68,
 
69
.  An infrared 
camera is also used to measure the surface temperature of the parallel plates.  This 
enables the determination of the surface temperature, an estimate of heat loss to the 
environment via a published correlation for natural convection from a vertically 
oriented plate, and provides a visual indication of the position of the flame.  
 
 
 123 
 
 
 
Figure 5-17:  Schematic diagram of the silicon-walled micro-burner. 
 
Figure 5-18 is a series of infrared images of the silicon plate surface temperature 
for different values of the inlet velocity. 
 
 
 
Si Wafers 
Plenum
Adjustable 
Spacing 
Flow Direction 
FTIR Beam
Flame 
 
FLIR 
IR Camera 
 
 124 
 
 
Figure 5-18:  Infrared images of the plates for different values of the inlet 
velocity of the fuel-air mixture.  The reactants are methane and air at ? =0.5. 
 
First, all of the flow velocities are in excess of 1m/s whereas the laminar flame 
speed for Methane-Air at an equivalence ratio of 0.5 is 0.4m/s.  This is a strong 
indication that burning velocity enhancement via heat recirculation through the 
silicon plates is occurring as predicted in Chapter 2.  Second, the images show that 
there is a range of inlet velocities for which a stable flame is obtained between the 
plates.  This demonstrates the shows the enhancement in steady-state flame stability 
discussed earlier in this chapter. 
U=1.1m/s 
U=1.2m/s 
U=1.25m/s U=1.275m/s 
 
 125 
 
Two methods are used to measure the flame position in the experiment.  The first 
uses the FTIR to obtain species concentration profiles where the flame location is 
defined as the position where the methane concentration drops sharply.  The second 
uses the infrared images of the structure temperature where the flame position is 
taken to be the point where the structure temperature is reaches a maximum.   
Figure 5-19 compares the dependence of flame location on inlet velocity 
determined using analytical, numerical, and experimental methods.  In the three 
approaches, the different parameters characterizing the burner (L, H, h
T,e
, and k
s
) have 
been set to the same value such that comparable results could be obtained.  The flame 
position x
flame
 is non-dimensionalized using the length of the combustor (L=55mm) 
and the flame speed is non-dimensionalized using the freely propagating flame speed 
S
L,fr
 which is 40 cm/s for CH
4
 and 90 cm/s for H
2
.  This non-dimensionalization is 
performed so that results from computationally-tractable numerical simulations using 
hydrogen as a fuel can be compared to experimental measurements in methane 
flames.  The geometric parameters in the analytical and numerical model are chosen 
to match the experiment (L=55mm, H=2mm, t=0.5mm).  The equivalence ratio of the 
H
2
 results is 0.5 and the equivalence ratio of the CH
4
 results is 0.8. 
 
 126 
 
0 0.5 1 1.5 2 2.5 3 3.5
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
U / S
L,fr
x
fl
a
m
e
 / 
L
Analytical Model
Numerical Simulation
Experiment (FTIR)
Experiment (IR Camera)
Low -Velocity
Blow -off
High-Velocity
Blow -off
 
Figure 5-19:  Dependence of flame location on inlet velocity predicted using 
analytical, numerical and experimental methods.  
  
The results show that the FTIR and IR camera measurements of flame position are 
consistent with each other and are reasonably close to the predictions of the 1D 
numerical simulation in spite of the fact that the simulations use H
2
 as a fuel whereas 
the experiment uses CH
4
.  This suggests that flame stabilization in this micro-channel 
is governed primarily by heat transfer and less by gas or surface chemistry.  Finally, 
the fact that the flame locations determined using the two experimental techniques 
(IR camera and FTIR spectrometer) agree with each other well suggests that flame 
location can be determined using structure temperature.  This is encouraging as it 
indicates that it will be possible to investigate combustion between plates made of 
materials like ceramics that are not transparent in the infrared. 
 
 127 
 
 
5.5. Transient Stability: Results and Discussion 
5.5.1 Introduction 
The previous section investigated the steady-state stability of flames stabilized in 
micro-channels with conductive walls.  However, another important issue in 
combustion systems is transient stability.  This is especially true in micro-combustors 
because low Reynolds numbers and short residence times may not allow efficient 
mixing.  Poor mixing can cause large fluctuations in the fuel-air composition at the 
entrance of the combustor and potentially lead to large disturbances in the flame 
position.  To investigate this, we sinusoidally perturb the inlet velocity and 
equivalence ratio in the numerical simulation and observe the temporal response of 
the flame.  This is illustrated schematically in figure 5-20.  The perturbation process 
only begins after steady-state has been reached (typically after 100 s). 
 
Figure 5-20: Schematic diagram illustrating the response of a flame 
stabilized in a micro-channel to perturbations in inlet velocity and 
equivalence ratio. 
H 
L
x
flame
(t)
Reactant Products 
U(t) 
? (t) 
 
 128 
 
5.5.2. Sample calculations 
Figure 5-21 shows the transient behavior of a flame subjected to sinusoidal inlet 
velocity perturbations.  The upper plot shows the non-dimensional inlet velocity as a 
function of non-dimensional time where S
L,fr
 is the freely propagating flam speed and 
T is the oscillation period.  For this case, T=0.1s (f=10Hz) and the length of the 
combustor is 1cm.  The lower plot shows the temporal response of non-dimensional 
flame position (L is the length of the passage) as a function of non-dimensional time 
for two values of the distance between the plates.  The figure shows that the 
fluctuations in flame position are much larger at H=5mm than at H=1mm.  Therefore, 
the increased thermal coupling at small passage height seems to enhance transient 
flame stability (just like it enhanced steady-state stability) by suppressing the 
response of the flame to perturbations in inlet velocity  
 
 
 129 
 
0 0.5 1 1.5 2
0
0.5
1
t / T
U /
 S
L
,fr
0 0.5 1 1.5 2
0
0.1
0.2
0.3
t / T
F
l
am
e 
p
o
si
t
i
on /
 L
H=5mm
H=1mm
 
Figure 5-21:  Transient behavior of H
2
-Air flames (F=0.5) subjected to 
perturbations in inlet velocity.  The upper figure shows the velocity 
perturbation and the lower figure shows the response of the flame for two 
different passage heights. 
 
Figure 5-22 shows the effect of sinusoidal oscillations in equivalence ratio where 
the equivalence ratio is defined as: 
.stoic
air
fuel
air
fuel
m
m
m
m
?
?
?
?
?
?
=?                                        (5-23) 
 
 130 
 
0 0.5 1 1.5 2
0.4
0.5
0.6
0.7
t / T
?
0 0.5 1 1.5 2
0
0.1
0.2
0.3
0.4
t / T
Fl
a
m
e 
posi
t
i
on /
 L
H=5mm
H=1mm
 
Figure 5-22:  Transient behavior of H
2
-Air flames subjected to perturbations 
in equivalence ratio.  The upper figure shows the velocity perturbation and 
the lower figure shows the response of the flame for two different passage 
heights. 
 
The figure shows that the flame stabilized in the smaller channel is also more 
stable with respect to oscillations in inlet mixture compositions.  This is consistent 
with the steady modeling and simulation results that show that axial heat transfer 
through the structure broadens the range of stable burning velocities.  The figure also 
shows that there is a 180? phase difference between the oscillation in equivalence 
ratio and the oscillation in the flame position.  This is because as the equivalence ratio 
 
 131 
 
decreases, the burning velocity of the mixture does also.  Since the inlet velocity 
remains constant, this causes the flame position to shift downstream. 
 
5.5.3. Effect of passage height 
While the sample cases shown in the previous section provide some qualitative 
insight regarding the influence of the passage height on the transient stability of 
flames in micro-channels, a more systematic approach needs to be undertaken to 
understand the phenomenon in greater detail.  To this end, the sinusoidal oscillation 
of inlet velocity is defined as: 
()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+= t
T
aUtU
?2
sin1
0
                                            (5-24) 
In equation 5-15, U
0
 is the mean velocity of the flow, T is the oscillation period, 
and a is a non-dimensional parameter reflecting the amplitude of the perturbation.  
Note that |a|>1 U(t) becomes negative.  To overcome this limitation, it is arbitrarily 
decided that  
?
?
?
?
?
?
?
?
?
?
?
?
?
?
+= ?
?
,)
2
sin(1max)(
0
t
T
aUtU                                        (5-25) 
where ?  is a small velocity (1cm/s). 
Figure 5-23 shows how the maximum flame position varies with the amplitude of 
the velocity perturbation for U
0
=0.7m/s and for two values of H. 
 
 132 
 
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0
0.1
0.2
0.3
0.4
0.5
a
M
a
xi
m
u
m
 f
l
am
e 
po
si
t
i
on /
 L
H=5mm
H=1mm
k
s
=20W/mK, T=0.1s, L=1cm, h
T, e
=50W/m
2
K
 
Figure 5-23: Evolution of normalized maximum flame position with 
normalized perturbation amplitude for different passage heights. 
 
Again, the figure shows that flames stabilized in channels with smaller gaps can 
sustain larger inlet velocity perturbations.  This greater stability can be attributed to 
the enhanced thermal coupling that occurs at small scale. 
 
5.5.3. Effect of oscillatory period 
Figure 5-24 shows the evolution of normalized maximum flame position as a 
function of normalized oscillation amplitude for different values of the oscillation 
period. 
 
 
 133 
 
 
0 0.5 1 1.5 2
0
0.05
0.1
0.15
0.2
0.25
a
M
a
xi
m
u
m
 f
l
am
e 
po
si
t
i
on /
 L
T=1s
T=0.1s
T=0.01s
 
Figure 5-24: Evolution of normalized maximum flame position with 
normalized perturbation amplitude for different oscillation periods. 
 
The figure shows that decreasing the period of the oscillation (or increasing its 
driving frequency) enables operation at higher velocity perturbation amplitudes.  This 
is because there is a delay in the response of the flame due to the finite timescale for 
conduction within the structure of the combustor.   
 
5.6. Conclusions 
A simple model has been developed to understand and characterize the steady-state 
stability of flames in micro-channels.  The predictions of the model are in good 
 
 134 
 
qualitative agreement with those of more detailed 1D and 2D numerical simulations, 
as well as with experimental data.  The enhanced stability was shown to depend on 
passage height, length, and the thermal conductivity of the structure.  In general, 
decreasing the passage height and increasing the passage length and thermal 
conductivity of the structure tend to improve steady-state stability, and these 
improvements are caused by axial heat transfer through the combustor structure.  The 
1D numerical model was also used to investigate transient stability and the results 
also showed that smaller passage heights improved stability with respect to 
oscillations in inlet velocity and equivalence ratio.  Taken together, the results 
described in this section demonstrate that another advantage of micro-burners is 
enhanced steady-state and transient stability. 
 
 
 
 
 
 
 135 
 
Chapter 6: Effect of Flame-Structure Thermal 
Interaction on Chemistry  
 
 
6.1. Introduction 
In Chapter 2, enhanced thermal coupling between the reacting gas and the micro-
combustor walls was shown to broaden the flame. In the numerical study, flame 
thickness was determined based on the spatial evolution of gas temperature: 
max
0f
r
dx
dT
TT
?
?
?
?
?
?
?
=?
                                                 (6-1) 
  However, one would also expect the broadening effect to be reflected in the 
spatial evolution of chemical species concentrations.  This is important for 
understanding the effect of thermal coupling on a micro-combustor?s chemical 
efficiency: if species concentration profiles become broader, micro-combustors will 
have to be proportionally longer than their conventional-scale counterparts to achieve 
comparable levels of efficiency.   
In Chapter 3 of this work, the performance of both hydrogen and methane micro-
combustors was investigated.  It was suggested that the behavior of H
2
-air flames and 
CH
4
-air flames was found to be quite different as the combustor?s scale was reduced.  
Methane flames were characterized by much larger temperature overshoots compared 
 
 136 
 
to hydrogen flames, leading to an optimum combustor configuration at a scale much 
larger than for hydrogen.  While the reasons for this discrepancy in behavior remain 
unclear, one hypothesis is that it is due to dissimilarities in the chemical mechanisms 
used to model combustion.  This would mean that the thermal coupling between the 
gas and structure that occurs at micro-scale affects the general characteristics of the 
combustion chemistry as the scale of the combustor is reduced.  Further, it suggests 
that global reaction mechanisms (for example
70, 71
 ) developed for conventional-scale 
flows where heat recirculation is not important may not be appropriate for use at the 
micro-scale where heat recirculation is strong and reaction zones are broadened.  
As a result, this chapter seeks to further our understanding of the influence of the 
thermal interaction between the reacting gas and the structure of the combustor by 
investigating the effect of this interaction on the chemical species profile through the 
flame and on the operation of the chemical mechanism for combustion.  Particular 
attention is paid to the relative importance of the different elementary reactions in the 
mechanism as the scale of the passage, and thus the strength of the interaction 
between the reacting fluid and the combustor structure, is varied.  This is important 
for understanding the effects of structural heat conduction on the progress of the 
individual elementary reactions that compose the overall reaction, and the effect of 
heat transfer on the choice of an appropriate reaction mechanism. 
 
6.2. Alternate definitions of reaction zone thickmess 
One way to define a chemically-based reaction zone thickness is to consider 
radical species concentration distribution.  Figure 6-1 shows the H and O radical 
 
 137 
 
profiles for large and small passage height (H=1cm and H=0.005cm respectively) 
where the passage length is L=1cm. The fuel considered is H
2
 burning in air at an 
equivalence ratio of 0.5 and the numerical code presented in Chapter 2 is used to 
compute the profiles. 
0 0.1 0.2 0.3 0.4 0.5 0.6
0
0.2
0.4
0.6
0.8
1
x / L
Y
O
/Y
O,max
 - H=1cm
Y
H
/Y
H,max
 - H=1cm
Y
O
/Y
O,max
 - H=5 10
-3
cm
Y
H
/Y
H,max
 - H=5 10
-3
cm
 
Figure 6-1: H and O radical mass fraction profiles as a function of non-
dimensionalized streamwise coordinate for L=1cm, ? =0.5 and P=1atm. 
 
The figure clearly shows that the thickness of the radical profiles is larger for the 
situation where H is smaller.  There are several possible approaches for determining 
the flame thickness based on chemical species concentration.  One approach is to 
define the flame thickness as the distance between the points at which the radical 
mass fraction reaches half of its peak value. This is illustrated in figure 6-2. 
 
 138 
 
 
 
Figure 6-2: Definition of the flame thickness based on the evolution of a 
radical species mass fraction. 
 
Another approach is to focus on a reaction product such as H
2
O and define the 
flame thickness as the distance necessary to reach 90% of the equilibrium H
2
O mass 
fraction.  This is illustrated in figure 6-3.  Finally, one could take a similar approach 
with a reactant by, for example, defining the flame thickness as the distance required 
for the H
2
 mass fraction to drop to 10% of its inlet value. 
?
r
x 
Y
max
/2 
Y
max 
 
 139 
 
 
Figure 6-3: Definition of the flame thickness based on the H
2
O mass fraction. 
 
Figure 6-4 compares the evolution of the non-dimensional flame thickness 
determined using the alternative chemically-based metrics discussed above to each 
other and to the flame thickness defined based on the temperature profile (equation 6-
1).  .  The conditions are P=1atm, ? =0.5 and L=1cm.  
?
r
Y
H2O, max 
90% Y
H2O, max 
x 
 
 140 
 
10
-3
10
-1
10
1
1
2
5
10
H / ?
r,fr
?
r
 / 
?
r,
fr
T
Y
O
Y
H
Y
H
2
O
Non-continuum Slip Continuum 
 
Figure 6-4: Evolution of non-dimensional flame thickness as a function of 
non-dimensional passage height for several definitions of the flame thickness. 
The legend indicates the parameter upon which the flame thickness is 
computation is based. 
 
The figure shows that each method of measuring the reaction zone thickness 
shows the same limit-cycle behavior.  This demonstrates that while the presence of 
the flame broadening phenomenon is independent of the measurement method, the 
amplitude of the broadening does depend on the measurement method.  Nevertheless, 
the amplitudes of the species-based metrics do seem to be self-consistent within 30%. 
The fact that the species profiles tend to broaden directly affects chemical 
efficiency.  Chemical efficiency is defined as the ratio of the total heat evolved in the 
 
 141 
 
reaction to the total heat available if the reaction were to go to completion.  In this 
work, we assume that this is equivalent to the fraction of H
2
 that gets consumed 
before the mixture exits the micro-channel: 
i,2
out,2
H
Y
H
Y
1
availableheatMaximum
.rxnbyevolvedheatTotal
ch
?=
=?
                                      (6-2) 
 Figure 6-2 and equation 6-2 suggests that a direct consequence of the broadening 
phenomenon is a reduction in the chemical efficiency.  Figure 6-5 is a plot of the 
chemical efficiency as a function of non dimensional passage height for fixed 
combustor lengths (L=1cm and L=0.5cm).  It shows that decreasing H decreases the 
chemical efficiency as the flame broadens. It also shows that a greater change in 
chemical efficiency occurs if the combustor length is reduced.  
 
 142 
 
10
-3
10
-1
10
1
94
95
96
97
98
99
100
H / ?
r,fr
?
ch
L / ?
r,fr
=18
L / ?
r,fr
=9
Continuum 
Slip Non-continuum 
 
Figure 6-5: Evolution of chemical efficiency as a function of non-
dimensionalized passage height for L=1cm and L=0.5cm.  The fuel is H
2
 
burning in air at an equivalence ratio of 0.5 and a pressure P=1 atm. 
 
6.3. Effect of thermal coupling on the kinetics of a multi-step reaction 
mechanism 
The numerical model presented in Chapter 2 of this work is used to investigate the 
effect of gas-structure interaction on the kinetics of the detailed hydrogen combustion 
mechanism used to model chemistry, i.e. the relative operation of the individual 
elementary reactions that are used to model the overall conversion of the reactants 
 
 143 
 
(H
2
 and air) to products (H
2
O) in the combustion process.  Table 4 in Appendix II 
lists the elementary reactions considered in the analysis as well as their assigned 
reaction number (R#).  Since each reaction has a different activation energy E
a
, each 
reaction has a different sensitivity to changes in temperature.  Therefore, it is 
reasonable to expect that differing levels of heat recirculation could cause the relative 
rates of the elementary reactions composing the overall reaction to change. 
Figure 6-6 shows the rates of progress for the major reactions (those with the 
largest overall reaction rate) - i.e. 1, 2, 3, 9 and 13 - occurring in the atmospheric 
pressure H
2
-air flame as a function of downstream distance x for L=5mm, H=1cm.  
These reactions are listed in table 3: 
Reaction number Reaction Reaction type 
1 H+O2=O+OH Branching 
2 H2+O=H+O Branching 
3 H2+OH=H2O+H Propagation 
9 H+O2+M=HO2+M Propagation 
13 HO2+H=OH+OH Propagation 
 
Table 3:  Major reactions occurring in the flame. 
 
The fuel considered is hydrogen burning in air at an equivalence ratio of 0.5 and a 
pressure p=1 atm.  Figure 6-7 is the same as 6-6 except that the value of H is now 
reduced to 100?m.   
 
 
 
 144 
 
0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
x (cm)
w (
g
m
o
l
/
c
m
3
 s)
R1
R2
R3
R9
R13
 
Figure 6-6: Evolution of rates of progress of the five most important 
reactions as a function of downstream distance for L=5mm and H=1cm. 
 
0.2 0.22 0.24 0.26 0.28 0.3
-0.02
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
x(cm)
w (
g
m
o
l
/
c
m
3
 s)
R1
R2
R3
R9
R13
 
Figure 6-7: Evolution of rates of progress of the five most important 
reactions as a function of downstream distance for L=5mm and H=100?m. 
 
 
 145 
 
Several observations can be made when comparing figures 6-6 and 6-7.  First, the 
rates of progress of the reactions for the smaller passage height are much higher than 
that for the larger passage height.  This is because as H decreases, heat recirculation 
raises the wall temperature which, in turn, raises gas temperature, the rates of the 
individual reactions, and hence the rate of the overall reaction.  This, in turn, results 
in an increase in the burning velocity of the mixture.  Second, the different peaks in 
the curves appear to be wider for the smaller passage height compared to the larger 
passage height (note that the ranges of the x-axis are different in figures 6-6 and 607).  
This is also a consequence of heat recirculation which tends to broaden the reaction 
zone, as discussed in section 6.2.  Lastly, the overall shapes of the rates of progress 
curves appear to be different for the two passage heights.  Taken together, these 
results suggest that the thermal coupling between the gas and the structure of the 
combustor not only has an impact on the burning velocity and flame thickness, but 
also on the elementary reactions that make up the overall reaction.   
To better quantify how the relative importance of the elementary reactions 
composing the overall reaction is affected by the size of the combustor, we define a 
parameter 
i
?
 which represents the total rate at which a particular reaction 
?processes? reactants divided by the rate at which material is processed by the overall 
reaction.  This can be written as the area contained under the curves shown on figures 
6-6 and 6-7 divided by the burning rate: 
 
L
L
0
i
i
S
dx)x(w
?
=?
                                             (6-3) 
 
 146 
 
In equation 6-3, w
i
(x) is the rate of progress of reaction i and S
L
 is the laminar 
burning velocity.  This parameter is representative of the importance of reaction i.   
Figure 6-8 shows 
i
?
 for the five most important reactions identified in figures 6-6 
and 6-7 for two different passage heights: H=1cm and H=100?m.  The figure shows 
that the relative importance of these elementary reactions changes as the size of the 
combustor is reduced:  While 
i
?
 increases slightly for reaction 2 and decreases 
slightly for reaction 1, it decreases significantly for the other reactions.  This shows 
that the reactions that dominate in the detailed mechanism change as the strength of 
the thermal coupling between the gas and the structure increases.   
1 2 3 9 13
0
0.2
0.4
0.6
0.8
1
x 10
-3
R#
?
 (gm
o
l
/
c
m
3
)
H=1c m
H=100?m
 
 
Figure 6-8: Comparison of the relative importance of the five most important 
reactions for two different values of the passage height. 
 
 
 147 
 
The results described above are important because they show that the strong 
interaction between gas and structure that occurs in micro-scale combustors affects 
the relative rates at which individual elementary reactions run.  This, in turn, suggests 
that global reaction mechanisms developed for conventional-scale problems where 
heat recirculation is absent or weak may not be appropriate for use in micro-scale 
combustors.  Therefore, in the absence of validated single-step mechanisms, it may be 
necessary to use detailed reaction mechanisms in order to obtain accurate results. 
 
6.4. Conclusions 
While structural heat conduction affects the general characteristics of a flame 
stabilized in a micro-channel (i.e. burning rate, flame thickness, and stability), the 
results presented in this chapter show that it also influences flame chemistry itself.   
This means that the consequences of flame broadening are not only observed in the 
temperature profile (which becomes less steep), but also in the evolution of chemical 
species (radical species, reactants, products) and of the evolution of the reaction rate 
profiles for each of the different elementary reactions composing the detailed 
mechanism.  Heat recirculation tends to broaden the chemical species profiles, 
suggesting that as the passage height is reduced a proportionally longer combustor is 
necessary to achieve good fuel conversion and therefore chemical efficiency. 
Another important conclusion is that the relative importance of individual 
elementary reactions is also affected by heat recirculation, suggesting that global 
mechanisms developed for large scale combustors where thermal coupling is weak 
might not be applicable to small scale combustors. 
 
 148 
 
While these observations are a first step towards the understanding of the effect of 
gas-structure thermal interaction on the chemical reactions itself, more substantial 
work would be necessary to fully characterize this influence.  In particular, the 
development of single step mechanisms validated for use at the micro-scale would 
facilitate 3D simulations necessary for capturing the full physics of gas-structure 
interaction. 
  
 
 149 
 
Chapter 7:   Conclusions and Perspectives 
 
 
7.1. Contribution and key findings 
This thesis has investigated the effect of structural heat conduction on several 
aspects of microscale combustion and propulsion.  Simple analytical models were 
developed in order to understand the basic physics of combustion in a micro-channel.  
Good qualitative agreement between the results of the analytical models and those of 
a more detailed numerical simulation that included distributed heat transfer between 
gas and structure, heat conducting within the structure, heat loss, and detailed 
chemistry was found, suggesting that the physical processes accounted for in the first-
order models are those governing flame propagation in the presence of a conductive 
structure.  These compact models offer the possibility for first-order design of micro-
combustors by investigating the effect of a wide range of parameters in little 
computational time.  They also provide unique insight regarding the basic physics of 
flame propagation in micro-channels and show that structural heat conduction is a 
dominating effect at small scale.  The latter was shown to affect flame characteristics 
and to greatly influence the performance of micro-combustors and micro-thrusters.  
The key findings of the thesis are summarized below. 
First, the structure of the flame itself was shown to be modified due to heat 
recirculation via the structure of a parallel plate micro-combustor.  Increases of flame 
thickness by an order of magnitude were predicted by a simple analytical model that 
 
 150 
 
extended the theory of Mallard and Le Chatelier.  The analytical results were 
subsequently confirmed by the more detailed numerical simulation.  
Second, burning rate was also shown to increase at small scale due to heat 
recirculation, and substantially raise the power density of micro-combustors as 
compared to their conventional scale counterparts.  However, the price to pay for this 
larger power density is that overall combustion efficiency decreases. 
Third, a similar behavior was found when the power density results were used to 
model the combustion chamber of a micro-thruster.  The possibility for a great 
increase in thrust-to-weight ratio was suggested, at the price of a lower specific 
impulse. 
Fourth, structural heat conduction was shown to greatly influence flame stability.  
In particular, enhanced flame stability was observed analytically, numerically, and 
experimentally.  This behavior was found to be in good qualitative agreement with 
the results of much more detailed 2D numerical simulations performed by other 
researchers. 
Fifth and finally, the thermal interaction between gas and structure was shown to 
influence chemistry at small scale, questioning the validity of global reaction 
mechanisms when thermal coupling is not negligible. 
 
7.2. Future work 
The field of micro-combustion is still developing and this work only investigated a 
specific aspect of the problem, i.e. the effect of structural heat conduction.  Much 
 
 151 
 
more work is still needed in order to find a truly practical application for sub-
millimeter scale combustors. 
Much more work is also needed to fully understand the effect of fluid-structure 
coupling on the characteristic of micro-combustors.   
First, one important area that requires more investigation is catalysis in order to 
take advantage of the large surface to volume ratios characteristic of micro-devices, 
and to model important processes such as radical quenching. 
Second, the results computed using the 1D numerical model should be compared 
to more detailed 2D numerical simulations.  This would allow better modeling of the 
heat transfer between gas and structure by avoiding the necessity to impose the 
Nusselt number. 
Third, pressure loss can be substantial in micro-devices, and the only reliable way 
to predict the magnitude of this loss is by solving the momentum equation for the 
flow field.  Thus, integrating this additional equation in the numerical model would 
be useful and relevant especially for 2D and 3D simulations.   
Fourth, the effect of slip on the thermal feedback process should be investigated in 
greater detail.  Temperature slip has already been studied in chapter 2 of this thesis 
but velocity and species concentration slip could also be important ? especially in 
catalytic systems. 
Fifth, more realistic fuels should be considered in order to provide more sensitive 
estimates of the performance of micro-combustors and micro-thrusters.  For example, 
natural gas or butane could be considered for Earth-based applications, JP10 for 
 
 152 
 
micro-UAV engines, and high pressure liquid hydrogen and oxygen for micro-
rockets. 
Sixth, the effect of thermal coupling on chemistry was only briefly investigated.  
Much more work is necessary in order to understand the effect of heat recirculation 
on detailed chemical processes. 
Finally, and probably most importantly, the results presented in this work are 
mostly analytical and numerical.  Thorough experimental investigation of the 
phenomena predicted by the models is necessary.  In particular, experiments need to 
be developed for measuring flame broadening, the possibility of obtaining high power 
density in micro-combustors, thrust-to-weight ratio augmentation in micro-thrusters, 
steady-state and transient stability, and the effect of thermal feedback on chemical 
conversion. 
 
 153 
 
Appendix I 
The Thermal Theory of Mallard and Le Chatelier 
 
In 1884, two French scientists developed the well-established thermal theory of 
flame propagation
72
 based on the idea that the burning velocity is fixed by the rate at 
which conductive heat transfer from the hot products in the post-flame region pre-
heats the reactants entering the reaction zone. The flame is illustrated in figure A1-1 
and a simplified temperature profile through it is illustrated in figure A1-2.  The 
flame is divided into two zones: 
? A pre-heat zone in which the temperature of the reactants rises from its 
initial value T
0
 to an ignition temperature T
i
.  
? A reaction zone in which the temperature of the gas rises from T
i 
to the 
adiabatic flame temperature T
f
.  
 
 
Figure A1-1: Freely propagating laminar premixed flame. 
S
L
Reactants Products 
?
r,MC 
 
 154 
 
 
Figure A1-2: Temperature profile through the freely propagating flame. 
 
The basic idea behind this simple model is that the energy needed to heat the 
reactants from T
i 
to T
0
 in the preheat zone must be transferred by conduction from the 
reaction zone to the preheat zone. This energy balance is expressed as follows: 
()qTTCm
0ir,pr
&& =?                                             (A1-1) 
where 
r
m&  is the mass flow rate of the reactants, C
p,r
 is the specific heat at constant 
pressure of the reactants and q&  is the rate of conductive heat transfer from the 
reaction zone to the pre-heat zone. 
Assuming that the temperature profile of the gas in the reaction zone is linear gives 
the following expression for the heat transfer from the reaction zone to the preheat 
zone can be expressed as: 
Pre-heat zone 
?
ph 
Reaction Zone 
?
r,MC 
x 
T
f 
T
i
T
0 
T 
q&
 
 155 
 
MC,r
if
rr
TT
Akq
?
?
=&                                              (A1-2) 
In this expression, k
r
 is the thermal conductivity of the reactants, A
r
 is the cross-
sectional area for conduction in the gas-phase, and 
r,MC
 is the flame thickness. 
The mass flow rate of the reactants can be expressed as: 
Lrrr
SAm ?=&                                                (A1-3) 
where S
L
 is the laminar flame speed or the relative speed of the incoming the 
reactants.  The laminar flame speed is related to the reaction rate RR and the flame 
thickness ?
r
 as follows: 
RR
S
MC,r
L
?
=                                                  (A1-4) 
Substituting A2-2 ? A2-4 into A2-1 gives: 
()()
L
rif0ir,pLrr
S
1
RRATTTTCSA ?=??                         (A1-5) 
Solving A1-5 for S
L 
gives: 
RR
C
k
TT
TT
S
r,pr
r
0i
if
L
??
?
=                                     (A1-6) 
This is the well-known expression for the laminar flame speed derived by Mallard 
and Le Chatelier. This expression can be inverted using equation A2-4 to give the 
reaction zone thickness which is the quantity we are interested in here: 
RR
1
C
k
TT
TT
r,pr
r
0i
if
MC,r
?
?
?
?
=                                        (A1-7) 
It will be convenient to note that A1-6 and A1-7 can also be derived using a 
thermo-electrical analogy. The appropriate thermal circuit for heat exchange in the 
flame is depicted in figure A1-3.   
 
 156 
 
?
r,MC 
T
0 
T
i
T
f
R
r
 
 
          
q,I &  
 
 
Figure A1-3: Thermal circuit for a freely-propagating flame  
           
In this analogy, electrical current (I) corresponds to heat flow (q& ), potentials (V) 
correspond to temperatures (T), and resistors relate temperature gradients to heat 
fluxes as follows: 
 
rr
R
T
q
R
V
I
??
=?= &                                                 (A1-8) 
                 Electric       Thermal  
                                                 intensity      heat flux 
 
The thermal resistance to axial conduction through the reactants gas column R
r
 can 
be expressed as: 
rr
MC,r
r
Ak
R
?
=                                                     (A1-9) 
Combining A1-8 and A1-9 gives the conductive heat flow: 
rr
MC,r
if
r
if
Ak
TT
R
TT
q
?
?
=
?
=&                                  (A1-10) 
Note that this is the same expression as A1-2. 
 
 157 
 
Combining A1-1, A1-3 and A1-10, and solving for S
L
 and ?
r
 recovers the same 
expressions for S
L
 and ?
r
 as derived earlier (A1-6 and A1-7). 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
 158 
 
Appendix II 
Detailed Reaction Mechanisms 
 
The mechanisms are written in CHEMKIN format: 
RT
Ea
b
eATk
?
=
 
The units are A: mole-cm-sec-K, E: cal/mole 
II.1. Hydrogen Combustion Mechanism 
Reaction number Reaction A b Ea 
1 H+O2=O+OH 5.10E+16 -0.8 16510 
2 H2+O=H+O 1.80E+10 1 8830 
3 H2+OH=H2O+H 1.20E+09 1.3 3630 
4 OH+OH=H2O+O 6.00E+08 1.3 0 
5 H+OH+M=H2O+M 7.50E+23 -2.6 0 
6 O2+M=O+O+M 1.90E+11 0.5 95560 
7 H2+M=H+H+M 2.20E+12 0.5 92600 
8 H2+O2=OH+OH 1.70E+13 0 47780 
9 H+O2+M=HO2+M 2.10E+18 -1 0 
10 H+O2+O2=HO2+O2 6.70E+19 -1.4 0 
11 H+O2+N2=HO2+N2 6.70E+19 -1.4 0 
12 HO2+H=H2+O2 2.50E+13 0 700 
13 HO2+H=OH+OH 2.50E+14 0 1900 
14 HO2+O=OH+O2 4.80E+13 0 1000 
15 HO2+OH=H2O+O2 5.00E+13 0 1000 
16 HO2+HO2=H2O2+O2 2.00E+12 0 0 
17 H2O2+M=OH+OH+M 1.20E+17 0 45500 
18 H2O2+H=HO2+H2 1.70E+12 0 3750 
19 H2O2+OH=H2O+HO2 1.00E+13 0 1800 
 
Table 4: Hydrogen combustion mechanism
73
 
 
 159 
 
 
II.2. Methane Combustion Mechanism 
Reaction number Reaction A b Ea 
1 H2+CH2(S)=CH3+H 7.23E+13 0 0 
2 H2+O=OH+H 51200 2.7 26.3 
3 H2O+H=H2+O 452000000 1.6 77.1 
4 CH4+O2=CH3+HO2 3.97E+13 0 238 
5 CH4+C=CH+CH3 5E+13 100.5 
6 CH4+H=CH3+H2 13200 3 33.6 
7 CH4+CH=C2H4+H 3.01E+13 0 -1.7 
8 CH4+CH2=CH3+CH3 4.3E+12 0 42 
9 CH4+CH2(S)=CH3+CH3 7E+13 0 0 
10 CH4+C2H=CH3+C2H2 1.81E+12 0 0 
11 CH4+O=CH3+OH 723000000 1.6 35.5 
12 CH4+OH=CH3+H2O 15700000 1.8 11.6 
13 CH4+HO2=CH3+H2O2 9.03E+12 0 103.4 
14 C2H2+C2H2=H2CCCCH+H 2000000000 0 242 
15 C2H2+O2=C2H+HO2 1.2E+13 0 312 
16 H2+C2H=C2H2+H 1.08E+13 9.1 
17 C2H2+H(+M)=C2H3(+M) 8.43E+12 0 10.8 
18 C2H2+CH=C2H+CH2 2.11E+14 0 -0.5 
19 C2H2+CH2=C3H4 1.2E+13 27.7 
20 C2H2+CH2(S)=H2CCCH+H 1.75E+14 0 0 
21 C2H2+C2H=C4H2+H 9.03E+13 0 0 
22 C2H2+O=CH2+CO 2170000 2.1 6.6 
23 C2H2+O=HCCO+H 5060000 2.1 6.6 
24 C2H2+OH=C2H+H2O 6E+13 0 54 
25 C2H2+M=C2H+H+M 1.14E+17 0 447 
26 C2H4+H=C2H3+H2 5.42E+14 0 62.4 
27 C2H4+H(+M)=C2H5(+M) 3970000000 1.3 5.4 
28 C2H4+CH=C3H4+H 1.32E+14 0 -1.4 
29 C2H4+CH2(S)=C3H6 9.64E+13 0 0 
30 C2H4+CH3=CH4+C2H3 4.16E+12 0 46.6 
31 C2H4+O=H+CH2HCO 4740000 1.9 0.8 
32 C2H4+O=CH3+HCO 8130000 1.9 0.8 
33 C2H4+O=CH2CO+H2 680000 1.9 0.8 
34 C2H4+OH=C2H3+H2O 2.05E+13 0 24.9 
35 C2H4+M=C2H2+H2+M 9.97E+16 0 299.3 
36 C2H4+M=C2H3+H+M 7.4E+17 0 404.1 
37 C2H6+H=C2H5+H2 1450000000 1.5 31 
38 C2H6+CH=C2H4+CH3 1.08E+14 0 -1.1 
39 C2H6+CH2(S)=CH3+C2H5 2.4E+14 0 0 
40 C2H6+CH3=C2H5+CH4 0.000000151 6 25.3 
41 C2H6+O=C2H5+OH 1000000000 1.5 24.3 
42 C2H6+OH=C2H5+H2O 7230000 2 3.6 
43 C2H6+HO2=H2O2+C2H5 1.32E+13 0 85.6 
 
 160 
 
44 C4H2+O=C3H2+CO 7.89E+12 0 5.6 
45 C4H2+OH=C3H2+HCO 6.68E+12 0 -1.7 
46 O2+CO=CO2+O 1.26E+13 196.9 
47 O2+CH2O=HCO+HO2 6.02E+13 0 170.1 
48 O2+C=CO+O 1.2E+14 16.7 
49 O2+H+M=HO2+M 2.1E+18 -0.8 0 
50 O2+H+H2O=HO2+H2O 6.89E+15 0 -8.7 
51 O2+H=OH+O 9.76E+13 62.1 
52 O2+CH=CO+OH 1.66E+13 0 0 
53 O2+CH=CO2+ 1.66E+13 
54 O2+CH2=CO2+H2 5.43E+12 0 6.2 
55 O2+CH2=CO2+H+H 5.43E+12 0 6.2 
56 O2+CH2=CO+OH+H 8.15E+12 0 6.2 
57 O2+CH2=CO+H2O 1.48E+12 0 6.2 
58 O2+CH2=CH2O+O 4.2E+12 6.
59 O2+CH2(S)=CO+OH+H 3.13E+13 0 0 
60 O2+CH3=CH2O+OH 3.31E+11 0 37.4 
61 O2+C2H=HCCO+O 9.05E+12 0 0 
62 O2+C2H=CO2+CH 9.05E+12 
63 O2+C2H3=C2H2+HO2 5.42E+12 0 0 
64 O2+C2H5=C2H4+HO2 10200000000 0 -9.2 
65 O2+C3H2=HCO+HCCO 1E+13 0 0 
66 O2+H2CCCH=CH2CO+HCO 30100000000 0 12 
67 O2+HCO=HO2+CO 3.01E+12 0 0 
68 O2+CH3O=CH2O+HO2 21700000000 0 7.3 
69 O2+CH2OH=CH2O+HO2 1.57E+15 -1 0 
70 O2+CH2OH=CH2O+HO2 7.23E+13 0 15 
71 O2+HCCO=CO+CO+OH 1.63E+12 0 3.6 
72 H2O2+H=HO2+H2 1.69E+12 15.7 
73 H2O2+H=OH+H2O 1.02E+13 0 15 
74 H2O2+O=OH+HO2 6.62E+11 0 16.6 
75 H2O2+OH=H2O+HO2 7.83E+12 0 5.6 
76 OH+OH(+M)=H2O2(+M) 7.23E+13 -0.4 0 
77 CO+O+M=CO2+M 1.54E+15 0 12.6 
78 CO+OH=CO2+H 16600000 1.3 -3.2 
79 CO+HO2=CO2+O 1.51E+14 0 99 
80 CO+CH=HCCO 2.77E+11 -7.2
81 CO2+CH=HCO+CO 3.43E+12 0 2.9 
82 CO2+CH2=CH2O+CO 23500000000 0 0 
83 CH2O+H=HCO+H2 126000000 1.6 9.1 
84 CH2O+CH=CH2+HCO 9.64E+13 0 -2.2 
85 CH2O+CH3=CH4+HCO 7.83E-08 6.1 8.2 
86 CH2O+O=HCO+OH 4.16E+11 0.6 11.6 
87 CH2O+OH=HCO+H2O 3430000000 1.2 -1.9 
88 CH2O+HO2=H2O2+HCO 3.01E+12 0 54.7 
89 CH2O+M=HCO+H+M 1.4E+36 -5.5 404.6 
90 CH2O+M=H2+CO+M 3.26E+36 -5.5 404.6 
91 CH2CO+H=CH3+CO 1.81E+13 0 14.1 
92 CH2CO+O=CH2+CO2 1.33E+12 0 5.7 
93 CH2CO+O=CH2O+CO 4.58E+11 0 5.7 
 
 161 
 
94 CH2CO+O=HCO+H+CO 2.52E+11 0 5.7 
95 CH2CO+O=HCO+HCO 2.52E+11 0 5.7 
96 CH2CO+OH=CH3+CO2 2.52E+12 0 0 
97 CH2CO+OH=CH2OH+CO 4.68E+12 0 0 
98 CH2CO+M=CH2+CO+M 6.57E+15 0 241 
99 CH2CO+M=HCCO+H+M 1140000000 0 0 
100 C+CH2=C2H+H 5E+13 
101 C+CH3=C2H2+H 5E+13 0 0 
102 C+OH=CO+ 5E+13 
103 H+H+M=H2+M 1.87E+18 -1 0 
104 H+H+H2=H2+H2 9.79E+16 -0.6
105 H+CH=C+H2 8.43E+12 0 0 
106 H+CH2=CH+H2 6.02E+12 -7.5 
107 H+CH2(S)=CH2+H 2E+14 0 0 
108 H+CH3(+M)=CH4(+M) 1.69E+14 0 0 
109 H+C2H3=C2H2+H2 1.2E+13 0 0 
110 CH3+CH3=C2H5+H 3.01E+13 0 56.5 
111 H+O+M=OH+M 1.18E+19 -1 0 
112 H+OH+M=H2O+M 5.53E+22 -2 0 
113 H+HO2=H2+O2 4.28E+13 0 5.9 
114 H+HO2=OH+OH 1.69E+14 3.7 
115 H+HO2=H2O+O 3.01E+13 0 7.2 
116 H+HCO=CO+H2 9.03E+13 0 
117 H+CH3O=CH2O+H2 1.81E+13 0 0 
118 H+CH2OH=CH3+OH 1.02E+13 0 0 
119 H+CH2OH=CH2O+H2 3.08E+13 0 0 
120 H+HCCO=CH2+CO 1.51E+14 0 0 
121 CH+CH2=C2H2+H 4E+13 0 0 
122 CH+CH3=C2H3+ 3E+13 0 0 
123 CH+C2H3=CH2+C2H2 5E+13 0 0 
124 CH+O=CO+H 3.97E+13 
125 CH+OH=HCO+H 3E+13 0 0 
126 CH+HCCO=C2H2+CO 5E+13 0 0 
127 CH2+CH2=C2H2+H2 1.2E+13 0 3.3 
128 CH2+CH2=C2H2+H+H 1.08E+14 0 3.3 
129 CH2+CH3=C2H4+H 4.22E+13 0 0 
130 CH2+C2H3=C2H2+CH3 1.81E+13 0 0 
131 CH2+O=CO+H+H 7.2E+13 0 0 
132 CH2+O=CO+H2 4.8E+13 
133 CH2+OH=CH2O+H 1.81E+13 0 0 
134 CH2+HCO=CH3+CO 1.81E+13 0 0 
135 CH2+HCCO=C2H3+CO 3E+13 0 0 
136 CH2+HCCO=C2H+CH2O 1E+13 0 8.4 
137 CH2(S)+M=CH2+M 1.51E+13 0 0 
138 CH3+CH3(+M)=C2H6(+M) 3.61E+13 0 0 
139 CH3+O=CH2O+H 8.43E+13 0 0 
140 CH3+OH=CH2(S)+H2O 7.23E+13 0 11.6 
141 CH3+HO2=CH3O+OH 1.8E+13 0 0 
142 CH3+HCO=CH4+CO 1.2E+14 0 0 
143 CH3+M=CH2+H+M 2.91E+16 0 379.1 
 
 162 
 
144 C2H+C2H3=C2H2+C2H2 1.9E+13 0 0 
145 C2H+O=CH+CO 1E+13 
146 C2H+OH=HCCO+H 2E+13 0 0 
147 C2H+OH=CH2+CO 1.81E+13 0 0 
148 C2H3+O=CO+CH3 3E+13 0 0 
149 C2H3+OH=C2H2+H2O 5E+12 0 0 
150 C2H5+O=CH2O+CH3 6.62E+13 0 0 
151 H2CCCH+O=C2H2+CO+H 1.39E+14 0 0 
152 H2CCCH+OH=C3H2+H2O 2E+13 0 0 
153 H2CCCCH+M=C4H2+H+M 1.12E+16 0 194.6 
154 O+O+M=O2+M 5.4E+13 -7.5 
155 O+HO2=O2+OH 3.19E+13 0 0 
156 O+HCO=CO+O 3.01E+13 
157 O+HCO=CO2+H 3.01E+13 0 0 
158 O2+CH3=CH3O+O 4.4E+13 0 131.4 
159 O+CH3O=CH2O+OH 1.81E+12 0 0 
160 O+CH2OH=CH2O+OH 9.03E+13 0 0 
161 O+HCCO=H+CO+CO 9.64E+13 0 0 
162 OH+OH=O+H2O 1510000000 1.1 0.4 
163 OH+HO2=H2O+O2 2.89E+13 0 -2.1 
164 OH+HCO=H2O+CO 1.02E+14 0 0 
165 OH+CH3O=CH2O+H2O 1.81E+13 0 0 
166 OH+CH2OH=CH2O+H2O 2.41E+13 0 0 
167 OH+HCCO=HCO+HCO 1E+13 0 0 
168 OH+HCCO=CH2O+CO 1E+13 0 0 
169 HO2+HO2=H2O2+O2 4.22E+14 0 50.1 
170 HO2+HO2=H2O2+O2 1.32E+11 0 -6.8 
171 HCO+HCO=CH2O+CO 3.01E+13 0 0 
172 HCO+M=H+CO+M 4.49E+14 0 65.9 
173 CH3O+M=CH2O+H+M 1.55E+14 0 56.5 
174 CH2OH+M=CH2O+H+M 1.26E+16 0 125.6 
175 HCCO+HCCO=C2H2+CO+CO 1E+13 0 0 
 
Table 5: The Leeds reaction mechanism for methane combustion
74
 
 
 163 
 
Bibliography 
 
1
 Juergen Mueller, ?Micropropulsion for small spacecraft?, Progress in Aeronautics 
and Astronautics, Volume 187, AIAA, 2000 
2
 Mehra, A, ?Development of a High Power Density Combustion System for a Silicon 
Micro Gas Turbine Engine?, Ph.D. Thesis, MIT, Feb 2000. 
3 
D. G. Shepherd, ?Aerospace Propulsion?, American Elsevier Publishing Company, 
Inc., p. 24. 
4
 Mehra, A, ?Development of a High Power Density Combustion System for a Silicon 
Micro Gas Turbine Engine?, Ph.D. Thesis, MIT, Feb 2000. 
5
 Zeldovich, Ya B., Zh. Eksp. Theor. Fiz., 11:159-168 (1941) 
6 
D. B. Spalding, Proc. Royal Soc. A, 240, 1220 (1957) 83-100. 
7
 Hardesty, D. and Weinberg, F., Combustion Sci. and Tech. 8 (1974) 201-214. 
8
 T. Takeno, K. and Sato, K Combustion Sci. and Tech. 20 (1979) 73-84. 
9
 Peterson, R. B., Microscale Thermophysical Engineering, 2, pp 121-131, 1998 
10
 Aichlmayr, H. T., Kittelson, D.B., and Zachariah, M. R., Chemical Engineering 
Science, v. 57, n. 19, Oct 21 2002, p. 4161-4171 
11
 Daou, J. and Matalon, M., Combustion and Flame 128 (2002) pp. 321-339. 
12
 Daou, J. and Matalon, M., ?Flame Propagation in Channels: Differential Diffusion 
Effects?, 3
rd
 Joint Meeting of the U.S. Sections of The Combustion Institute, March 
2003. 
13
 Ronney, P.D., Combustion and Flame 135 (2003) pp. 421-439. 
14
 Y. Ju, B. Xu, Proc. Combust. Inst., vol. 30, 2005, pp. 2445-2453. 
 
 164 
 
15
 D. G. Norton, D. G. Vlachos, Combustion and Flame, 138, pp 97-107 (2004). 
16
 S. W. Churchill, Chem. Eng. Technol., 12, pp 249-254 (1989). 
17
 Alexeenko AA, Levin DA, Fedosov DA, Gimelshein SF, Collins RJ , Journal of 
Propulsion and Power 21 (1) 95-101 Jan-Feb 2005 
18
 Ketsdever AD, Clabough MT, Gimelshein SF, Alexeenko A, AIAA Journal 43 (3): 
633-641 March 2005 
19
 Epstein, A. and Senturia, S. D., Science 1997, 276 1211 
20
 Mehra, A, ?Development of a High Power Density Combustion System for a 
Silicon Micro Gas Turbine Engine?, Ph.D. Thesis, MIT, Feb 2000. 
21
 Spadaccini, C. M, Zhang, X., Cadou, C. P., Miki, N., and Waitz, I. A., 
?Development of a Catalytic Silicon Microcombustor for Hydrocarbon-fueled Power 
MEMS?, Proceedings of the IEEE Micro Electro Mechanical Systems (MEMS), 
2002, p. 228-231 
22
 Jeongmin Ahn, Craig Eastwood, Lars Sitzki and Paul D. Ronney, Proc. Combust. 
Inst., vol. 30, 2005 
23
 Miesse C, Masel RI, Short M, Shannon MA,  Combustion Theory And Modelling 9 
(1): 77-92 Feb 2005  
24
 A.P. London, A.A. Ayon, A.H. Epstein, S.M. Spearing, T. Harrison and al, Sensors 
and Actuators A 92 pp351-357, 2001. 
25
 Hitt D.L, Zakrzwski, C.M. and Thomas M.K., 2001, Smart Mater. Struct., vol. 10, 
1163-1175.  
26
 Ketsdever AD, Clabough MT, Gimelshein SF, Alexeenko A, AIAA Journal 43 (3): 
633-641 Mar 2005 
 
 165 
 
27
 Leach, T. T., ?Modeling and Simulation of the Micro-Combustion between Two 
Parallel Plates?, M.S. Thesis, University of Maryland, 2003. 
28
 T. Leach, C. Cadou, G. Jackson, Combustion Theory and Modelling, Vol. 10, No. 
1, February 2006, 85--103. 
29
 Mallard, E. and Le Chatelier, H. L., Ann. Mines 4 (1883) 379. 
30
 D. B. Spalding, Proc. Royal Soc. A, 240, 1220 (1957) 83-100. 
31
 Hardesty, D. and Weinberg, F., Combustion Sci. and Tech. 8 (1974) 201-214. 
32
 T. Takeno, K. and Sato, K Combustion Sci. and Tech. 20 (1979) 73-84. 
33
 Coffee, T.P., Kotlar, A.J., and Miller, M.S, Combustion and Flame 54 (1983) 155-
169. 
34
 Zhu, H. and Jackson, G.S., ?Transient Modeling for Assessing Catalytic Combustor 
Performance in Small Gas Turbine Applications?, 2001-GT-0520, Proceedings of 
ASME Turbo Expo 2001: 46
th
 ASME International Gas Turbine & Aeroengine 
Technical Congress, New Orleans, LA, June 4-7, 2001. 
35
 Miller, J. A., and Bowman, C. T., Progress in Energy and Combustion Science 15 
(1989) 287-338. 
36
 Deuflhard, P., Harier, E., and Zugck, J., Numer. Math. 51 (1987) 501-516. 
37
 Miller, J.A., Mitchell, R.E., Smooke, M.D., and Kee, R.J., Proc. Comb. Inst.  19 
(1982) 181. 
38
 Law, C.K, ?A Compilation of Experimental Data on Laminar Burning Velocities? 
in: N. Peters and B. Rogg, (Eds) Reduced Kinetic Mechanisms for Applications in 
Combustion Systems, Springer-Verlag, New York, 1993, pp. 15-26. 
39
 Law, C.K. and Sung, C.J., Progress in Energy and Combustion Science 26 (2000) 
459-505. 
 
 166 
 
40
 Glassman, I., ?Combustion?, Third edition, Academic Press, 1996. 
41
 Karniadakis, G.E. and Beskok, A., ?Micro-Flows: Fundamentals and Simulation?, 
Springer, New York, pp. 15. 
42
 Karniadakis, G.E. and Beskok, A., ?Micro-Flows: Fundamentals and Simulation, 
Springer?, New York, pp. 15. 
43
 Mohamed Gad-el-Hak, ASME Journal of Fluids Engineering, 121 (1999). 
44
 N. G. Hadjiconstantinou, O. Simek, ASME Journal of Heat Transfer, 124 (2), 356-
364. 
45
 R. L. Panton, ?Incompressible flow?, Second edition, John Wiley and Sons, Inc. 
46
Ronney, P.D., Combustion and Flame 135 (2003) pp. 421-439. 
47
 Mills, A. F., ?Basic Heat and Mass Transfer?, Irwin, pp. 241, 1995. 
48
 Lide, D.R (ed.), ?CRC Handbook of Chemistry and Physics?, 81th ed. 
49
 Miller, J. A., and Bowman, C. T., Progress in Energy and Combustion Science 15 
(1989) 287-338. 
50
 T. Takeno, K. and Sato, K Combustion Sci. and Tech. 20 (1979) 73-84 
51
 A. Carlos Fernandez-Pello, Proc. Combust. Inst. 29 (2002) 
52
 W. A. Strauss, R. Edse, Proc. Combust. Inst. 7 (1959) 377. 
53
 K.J. Hughes, T. Tur?nyi, A. Clague, M.J.Pilling, Int. J. Chem. Kinet., 2001, 33, 
513-538 
54
 John D. Anderson, ?Modern Compressible Flow?, Third Edition, McGraw Hill. 
55
 www.boeing.com/defense-space/space/propul/SSME.html 
56
 J. Buckmaster, T. Takeno, Combustion Science and Technology, 25, pp 153-158 
(1981). 
 
 167 
 
57
 S. W. Churchill, Chem. Eng. Technol., 12, pp 249-254 (1989). 
58
 A. F. Mills, ?Basic Heat and Mass Transfer?, Irwin. 
59
 A. F. Mills, ?Basic Heat and Mass Transfer?, Irwin. 
60
 S. R. Turns, ?An Introduction to Combustion?, 2
nd
 Edition, McGraw and Hill. 
61
 T. P. Coffee, A. J. Kotlar, M. S. Miller, Combust. And Flame 54 pp 155-169 
(1983). 
62
 Glassman, I., ?Combustion?, Third edition, Academic Press, 1996,  p. 174. 
63
 Y. Ju, B. Xu, Proc. Combust. Inst., vol. 30, 2005, pp. 2445-2453. 
64
 Y. Ju, B. Xu, Proc. Combust. Inst., vol. 30, 2005, pp. 2445-2453. 
65
 D. G. Norton, D. G. Vlachos, Combustion and Flame, 138, pp 97-107 (2004). 
66
 Scott Heatwole, ?In Situ Infrared Diagnostics for a Silicon Walled Micro-Scale 
Combustion Reactor?, MS Thesis, University of Maryland, 2004. 
67
 Heatwole S, Cadou CP, Buckley SG, Combustion Science And Technology 177 (8): 
1449-1461 Aug 2005 
68
 Scott Heatwole, ?In Situ Infrared Diagnostics for a Silicon Walled Micro-Scale 
Combustion Reactor?, MS Thesis, University of Maryland, 2004. 
69
 Heatwole S, Cadou CP, Buckley SG, Combustion Science And Technology 177 (8): 
1449-1461 Aug 2005 
70
 T. Coffee, A. Kotlar, M. Miller, Combust. Flame, 58 (1984) 59-67. 
71
 C.K. Westbrook, F.L. Dryer, Combust. Sci. Technol. 27 (1981) 31?43. 
72
 Mallard, E. and Le Chatelier, H. L., Ann. Mines 4 (1883) 379. 
73
 Miller, J. A., and Bowman, C. T., Progress in Energy and Combustion Science 15 
(1989) 287-338. 
 
 168 
 
74
 K.J. Hughes, T. Tur?nyi, A. Clague, M.J.Pilling, Int. J. Chem. Kinet., 2001, 33, 
513-538