DISS. ETH NO. 21763
Direct Numerical Simulation of Catalytic
Ignition
A thesis submitted to attain the degree of
DOCTOR OF SCIENCES of ETH ZURICH
(Dr. sc. ETH Zurich)
presented by
Andrea Brambilla
M.Sc. Energy Engineering, Politecnico di Milano
born on 14.06.1982
citizen of Italy
accepted on the recommendation of
Prof. Dr. Konstantinos Boulouchos, examiner
PD Dr. Ioannis Mantzaras, co-examiner
Dr. Christos Frouzakis, co-examiner
Prof. Dr. Gianpiero Groppi, co-examiner
2014
Abstract
Combustion of hydrogen and syngas in micro- or mesoscale channels has
recently attracted increasing interest for the development of both portable
power generation devices and large-scale “zero-emissions” power plants. Although gas-phase combustion can be stabilized even in small confinements
by adopting technical measures like annealed walls and excess enthalpy combustion, flame/wall interactions, mainly in terms of radical quenching and
heat losses, may result in many unstable combustion modes. On the other
hand, in small geometrical confinements catalytic combustion is more favorable compared to pure gas-phase combustion due to the increase in surfaceto-volume ratio and the suppression of most intrinsic flame instabilities in
the presence of a catalyst. Within this context, Direct Numerical Simulation assumes a great relevance for fundamental studies by resolving all the
involved spatiotemporal scales. In the present work, a Direct Numerical
Simulation low-Mach-number reactive flow solver including detailed transport, and gas-phase and surface chemistry, was used to study two different
topics relevant to micro-/mesoscale combustion.
In the first part, instabilities occurring in lean syngas/air homogeneous
combustion in an inert mesoscale channel were investigated. Experimental
images obtained by collecting the OH* chemiluminescence signal in a channel flow reactor were used to validate the code against steady flames while
an observed oscillatory flame was the starting point for a parametric nui
merical study. The combustion modes were mapped as a function of wall
temperature and CO to H2 molar ratio and new instabilities, not previously
reported in the literature, were found. Finally, recent experimental results
with the fast OH-LIF technique revealed steady and oscillatory modes in
agreement with the computed ones.
In the second part, the numerical solver was further developed by implementing conjugate heat transfer that was successfully validated against
analytical and numerical results available in literature. The startup of lean
hydrogen/air hetero-/homogeneous combustion in a mesoscale channel was
studied by resolving for the first time all relevant spatiotemporal scales.
Fast transient processes like light-off and gas-phase ignition, with a duration ranging from fractions of milliseconds to a few milliseconds, were studied in detail. Direct Numerical Simulation was also used for benchmarking
a quasisteady-state code and the limitations of its underlying assumptions
were identified.
Throughout the work, the Computational Singular Perturbation (CSP)
analysis has been extensively used as a diagnostic tool to obtain insights
into the physics behind the observed phenomena. Some important results
achieved are the identification of competition for OH radical between H2
and CO as a main mechanism responsible for weak combustion in lean
syngas/air homogeneous combustion, and the recognition of OH desorption
from the catalytic surface as a key step for triggering gas-phase ignition in
lean hydrogen/air hetero-/homogeneous combustion.
ii
Sommario
La combustione di idrogeno e syngas in microcanali e mesocanali sta attraendo sempre piu’ interesse a livello industriale sia per lo sviluppo di
dispositivi portatili di generazione di potenza che nell’ambito di centrali
termoelettriche a ridotte emissioni inquinanti. Sebbene la combustione in
fase gassosa possa essere stabilizzata anche in canali di piccole dimensioni
(da frazioni di millimetro a pochi millimetri) adottando opportuni accorgimenti tecnici quali l’inertizzazione delle pareti ed il preriscaldamento dei
reagenti, l’interazione fiamma/parete rende la combustione instabile a causa
delle perdite termiche e della ricombinazione dei radicali. Per canali di piccole dimensioni la combustione catalitica e’ preferibile rispetto alla combustione in fase gassosa dati l’elevata superficie per unita’ di volume e l’effetto
stabilizzante del catalizzatore. All’interno del contesto appena introdotto,
la “Direct Numerical Simulation” riveste notevole importanza a livello di
ricerca di base per la sua capacita’ di studiare tutte le scale spaziali e temporali coinvolte. Per il presente studio, un codice numerico per “Direct
Numerical Simulation” di fluidi reagenti a bassi numeri di Mach che include
cinetica chimica dettagliata (sia in fase gassosa che eterogenea) e trasporto
di massa, e’ stato utilizzato per studiare due differenti problematiche relative alla combustione in microcanali e mesocanali.
La prima parte della tesi e’ dedicata allo studio di instabilita’ nella combustione in fase gassosa di miscele magre di syngas e aria in mesocanali. Il
iii
codice numerico e’ stato dapprima validato con immagini sperimentali di fiamme stazionarie ottenute rilevando la chemiluminescenza del radicale OH
eccitato. Successivamente, partendo da condizioni sperimentali alle quali
una fiamma instabile e’ stata osservata, si e’ effettuato uno studio parametrico con lo scopo di mappare i diversi tipi di instabilita’ in funzione della
temperatura di parete e della composizione del syngas in termini del rapporto CO:H2 . Nuove instabilita’ mai riportate nella letteratura scientifica
sono state osservate. Infine, risultati sperimentali ottenuti recentemente
con OH-LIF ad alta frequenza hanno rivelato l’esistenza di fiamme stabili
ed oscillatorie in accordo con i risultati numerici.
Nella seconda parte, il codice di calcolo e’ stato ulteriormente sviluppato con l’aggiunta dello scambio termico con la parete solida (scambio termico coniugato) e della conduzione termica multidimensionale all’interno
della parete che sono stati validati con dati analitici e numerici disponibili
nella letteratura scientifica. Il transitorio di avviamento della combustione
di miscele magre di idrogeno e aria all’interno di mesocanali catalitici e’
stato quindi studiato. Per la prima volta in questo tipo di studi e’ stato
possibile simulare tutti i processi coinvolti e le relative scale spaziali e temporali. Fenomeni transitori come l’ignizione catalitica ed in fase gassosa, la
cui durata varia da frazioni di millisecondo a pochi millisecondi, sono stati
studiati in dettaglio. Il codice numerico per “Direct Numerical Simulation”
e’ stato anche usato per verificare le previsioni e le limitazioni di un codice
di calcolo quasi stazionario.
I dati ottenuti dalle simulazioni sono stati spesso analizzati con la
tecnica chiamata “Computational Singular Perturbation” (CSP) per rivelare i processi chimici e fisici coinvolti, e alcuni importanti risultati sono
stati ottenuti. Nell’ambito della combustione in fase gassosa di miscele
magre di syngas e aria, la competizione per il radicale OH tra idrogeno e
monossido di carbonio e’ stata ritenuta responsabile per il fenomeno definito
“weak combustion“, mentre nell’ambito della combustione catalitica il desorbimento del radicale OH dalla parete e’ stato identificato come processo
chiave nell’ignizione in fase gassosa.
iv
Acknowledgements
First, I would like to thank all the people who gave me technical (but not
only) support and in particular the LAV team headed by Professor Konstantinos Boulouchos for providing a pleasant working place, my office mate
Martin Schmitt, my supervisors Dr. Yiannis Mantzaras and Dr. Christos
Frouzakis for their guidance, the experimental team at CFG/PSI for the
experimental data, Prof. Ananias Tomboulides and Dr. Stefan Kerkemeier
for helping with code development, Prof. Gianpiero Groppi for accepting
to be my co-examiner, Paul Scherrer Institute and Swiss National Science
Foundation for financial support, and the team at Brutus cluster of ETHZ
for providing the computational resources.
My sincere gratitude goes to all those who are co-authors of the great
experience I had in Z¨
urich, sharing with me joy but also frustration moments. I would like to thank the colleagues who were also good friends
out of ETH, all the colleagues of the Italian-Greek-Persian-Indian-... lunch
team with whom I enjoyed the lunch break regardless of the food served in
the canteen, the Italian connection in Z¨
urich, my flatmate, and the Z¨
urichMilano weekend commuters with whom I spent most of my Sunday evenings
watching movies and discussing during the epic trips with Cisalpino train.
I gratefully acknowledge all the friends and people that I met in Z¨
urich,
whose names are not listed because the list would be longer than the thesis
itself, but each of whom gave an important contribution.
v
Last but not least I would like to thank the people who stood on my
side even being geographically far away from me and (hopefully!) suffering a
bit for my departure. A very special thank to my girlfriend, who encouraged
me continuously, to my family and all my friends in Italy.
vi
Contents
Contents
vii
List of Tables
1
List of Figures
2
1 Introduction
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . .
5
5
10
2 Numerical Approach
2.1 Governing equations . . . . . . . . . . . . . .
2.1.1 Gas-phase . . . . . . . . . . . . . . . .
2.1.2 Catalytic wall . . . . . . . . . . . . . .
2.2 Numerical method . . . . . . . . . . . . . . .
2.3 Code validation: Conjugate Heat Transfer . .
2.3.1 Validation of the steady-state solution
2.3.2 Validation of the transient solution . .
2.4 Computational Singular Perturbation Analysis
2.4.1 Methodology . . . . . . . . . . . . . .
13
13
13
18
21
23
24
25
28
29
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3 Flame dynamics in lean CO/H2 /air combustion
33
3.1 Review on flame dynamics . . . . . . . . . . . . . . . . . . . 33
vii
Contents
3.2
3.3
3.4
Computational domain, boundary conditions and chemistry
Results and discussion . . . . . . . . . . . . . . . . . . . . .
3.3.1 Steady experimental vs DNS results . . . . . . . . . .
3.3.2 Effect of wall temperature . . . . . . . . . . . . . . .
3.3.3 Effect of the CO to H2 ratio . . . . . . . . . . . . . .
3.3.4 CSP analysis of the V-shaped and weak flames . . . .
3.3.5 Experimental results with fast OH-LIF . . . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
40
41
45
65
68
74
87
4 Hydrogen catalytic combustion in Pt-coated channels
91
4.1 Review on numerical models for catalytic combustion . . . . 91
4.2 Computational domain, boundary conditions and chemistry
95
4.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . 98
4.3.1 CSP analysis of a transient SPSR . . . . . . . . . . . 99
4.3.2 Detailed Numerical Simulations . . . . . . . . . . . . 103
4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5 Conclusions and future work
141
5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
Bibliography
viii
149
List of Tables
3.1
3.2
3.3
3.4
Effective Lewis number . . . . . . . . . .
Summary of non-ignited and weak flames
Explosive mode analysis V-shaped flame
Explosive mode analysis weak flame . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
58
62
70
72
4.1
4.2
4.3
4.4
4.5
4.6
4.7
Wall thermophysical properties . . . . . . . .
CSP analysis SPSR at Tr = 1200 K . . . . . .
CSP analysis SPSR: influence initial coverage
Numerical conditions . . . . . . . . . . . . . .
Ignition times . . . . . . . . . . . . . . . . . .
Hydrogen conversion . . . . . . . . . . . . . .
CSP analysis of gas-phase ignition . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
97
101
102
103
113
114
120
1
List of Figures
2
2.1
2.2
Validation of steady conjugate heat transfer. . . . . . . . . . . .
Validation of transient conjugate heat transfer. . . . . . . . . .
25
27
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17
3.18
3.19
Domain and boundary conditions . . . . . . . . . . . . . . .
First comparison DNS vs steady experimental flame . . . . .
Second comparison DNS vs steady experimental flame . . .
Oscillatory flame with OH* chemiluminescence . . . . . . .
iHRR and conversion oscillatory ignition (964 K≤ TW ≤ 969
Summary of oscillatory ignition . . . . . . . . . . . . . . . .
Frequency of oscillatory ignition . . . . . . . . . . . . . . . .
Conversion at TW = 994 K . . . . . . . . . . . . . . . . . . .
2-D OH maps of oscillatory ignition at TW = 974 K . . . . .
Tmax and heat transfer rate at TW = 974 K . . . . . . . . . .
Propagation speed and stretch at TW = 974K . . . . . . . .
Mass fractions 2-D maps at TW = 974 K . . . . . . . . . . .
Reaction rates 2-D maps at TW = 974 K . . . . . . . . . . .
Random ignition spots . . . . . . . . . . . . . . . . . . . . .
Weak flames . . . . . . . . . . . . . . . . . . . . . . . . . . .
V-shaped flames . . . . . . . . . . . . . . . . . . . . . . . .
Numerical stability map . . . . . . . . . . . . . . . . . . . .
Flame tip opening . . . . . . . . . . . . . . . . . . . . . . .
Explosive modes . . . . . . . . . . . . . . . . . . . . . . . .
38
43
44
46
49
50
51
52
53
54
56
59
60
61
63
64
65
68
69
. .
. .
. .
. .
K)
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
. .
List of Figures
3.20
3.21
3.22
3.23
3.24
3.25
3.26
3.27
3.28
Reaction rates R3 and R25 . . . . . . . . . . .
Test rig with fast OH-LIF . . . . . . . . . . .
Experimental stability map with fast OH-LIF
V-shaped flames with fast OH-LIF . . . . . .
Asymmetric flames with fast OH-LIF . . . . .
First oscillatory flame with fast OH-LIF . . .
Measured wall temperature oscillatory flames
Spatially integrated OH oscillatory flames . .
Second oscillatory flame with fast OH-LIF . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
73
76
79
80
82
84
84
86
86
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18
4.19
4.20
4.21
4.22
Domain and boundary conditions . . . . . . . . . . . . . .
CSP analysis transient SPSR . . . . . . . . . . . . . . . .
Rates of progress during light-off . . . . . . . . . . . . . .
Light-off Case 1: iSHRR, Stefan velocity, TW , coverages .
Hydrogen profile evolution during light-off . . . . . . . . .
Hydrogen conversion rate during light-off . . . . . . . . . .
Wall temperature evolution during heat-up . . . . . . . . .
Coverages evolution during heat-up . . . . . . . . . . . . .
2-D maps gas-phase radicals during heat-up . . . . . . . .
Explosive mode gas-phase ignition . . . . . . . . . . . . . .
Gas-phase ignition Case 3: iHRR, iSHRR and conversion
Gas-phase ignition Case 3: coverages . . . . . . . . . . . .
Gas-phase ignition Case 3: 2-D OH maps . . . . . . . . . .
Gas-phase ignition Case 6: iHRR, iSHRR and conversion
Gas-phase ignition Case 6: 2-D OH maps . . . . . . . . . .
Gas-phase ignition Case 7: 2-D OH maps . . . . . . . . . .
Gas-phase ignition Case 7: Stefan velocity, Pt(s) and O(s)
DNS vs QSS predictions: light-off . . . . . . . . . . . . . .
DNS vs QSS predictions: gas-phase ignition . . . . . . . .
QSS failure: 2-D H2 maps . . . . . . . . . . . . . . . . . .
QSS failure: H(s) and O(s) . . . . . . . . . . . . . . . . . .
QSS failure: TW . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
96
100
105
107
110
111
115
118
119
120
122
123
124
125
126
128
129
130
131
134
135
136
3
Chapter 1
Introduction
1.1
Motivation
Recent advances in micro-fabrication technologies have sustained a growing
trend in the miniaturization of electromechanical devices, which nowadays
rely heavily on batteries [1, 2]. One of the main disadvantages of batteries, especially at the millimeter/centimeter-scale range, is their low energy
density. Even the most advanced lithium-ion batteries have an energy density of merely 0.2 kWh/kg [3] that cannot compete with the typical energy
densities of hydrocarbon fuels, e.g. about 15 kWh/kg for methane and
12 kWh/kg for diesel. The comparison clearly shows the attractiveness
of hydrocarbon-based power generation systems for powering microdevices.
Combustion appears to be a strong candidate for the conversion of chemical
energy. A miniature combustion device with a mere 3 % efficiency would
favorably compete with state-of-the-art batteries [1] simply because the fuel
is easily replaceable and thus the device does not suffer from long recharging
5
1. Introduction
times. Chigier et al. [4] estimated that a microcombustor can achieve energy
density up to 10 kWh/kg, while micro-engines may have much higher energy density than conventional engines [2]. Furthermore, there are specific
applications where mechanical power and/or heat are required and combustion devices can provide them directly. Several research projects have been
carried out in order to develop combustion-based micro power generation
systems, including Swiss-roll-type microcombustors coupled to thermoelectric devices [5, 6], microscale Wankel-type rotary engines [1, 7], micro gas
turbines [8, 9], and microthrusters [10, 11].
Pure gas-phase combustion at the micro- and mesoscale, defined respectively as sizes being smaller or larger than 1 mm [2], poses severe challenges
in terms of stability because with the increase of surface-to-volume ratio
wall heat losses and radical quenching play a larger role and proper technical measures have to be taken in order to stabilize a flame. Pure gas-phase
combustion of methane has been demonstrated at the sub-millimeter scale
in reactors with appropriately annealed walls [12]. Practical configurations
such as the mesoscale Swiss-roll-type burner, where the burned gas energy
is recycled to the unburned gas mixture, have been shown to stabilize nearstoichiometric propane/air flames at channel widths of 3.5 mm [13] by using
the concept of “excess-enthalpy” combustion introduced in the pioneering
works of Weinberg et al. [14] and Takeno et al. [15]. On the other side, high
surface-to-volume ratios favor catalytic (heterogeneous) combustion, which
becomes the preferred route compared to pure gaseous (homogeneous) combustion. Catalytic reactions broaden the classical flammability limits and
6
1.1. Motivation
lead to stable combustion even in the presence of high heat losses [13,16]. Recent works [17, 18] have demonstrated the potential of catalytically coated
walls in suppressing intrinsic flame instabilities in micro- and mesoscale
channels. Applications of catalytic combustion at the micro-/mesoscale
also include stationary gas turbines used in large-scale power generation,
whereby combustion is attained in honeycomb reactors comprising a multitude of catalytically coated channels with individual hydraulic diameters
ranging from less than 1 mm to about 2 mm [19]. Relevant combustion
methodologies comprise the fuel-lean catalytically stabilized thermal combustion (CST) [20, 21], the fuel-rich catalytic/fuel-lean gaseous combustion [22–24], and the fuel-lean inverse catalytically stabilized combustion
(i-CST) [25], with the potential to achieve single-digit NOx emissions (≤ 3
ppmv) [20].
Within the context of growing environmental and fuel-supply security
concerns, combustion of hydrogen and hydrogen-rich fuels like synthesis gas
(syngas) appears to be particularly attractive for both large-scale power
plants and microreactors for portable power generation [26]. Hydrogen-rich
syngases can be produced in large power plants through gasification of various solid fuels like coal and biomass [27], refinery residues [28] and municipal solid wastes [29]. A potential application of these fuels is in Integrated
Gasification Combined Cycle (IGCC) approaches, which use a variety of fuel
feedstocks for power generation. When precombustion capture and storage
of CO2 are applied [30, 31], IGCC provides a viable way towards a carbonfree power generation. Combined production of power and hydrogen can
7
1. Introduction
also be achieved in some IGCC approaches [28] and in technologies based on
natural gas decarbonization [32]. Hydrogen-containing fuels are also investigated in postcombustion CO2 capture concepts for natural-gas-fired power
plants [33, 34]. Finally, hydrogen and hydrogen-rich fuels are of interest in
microreactors for portable power generation units [35] where hydrocarbon
fuels can be processed in suitable on-board micro-reformers [36–39].
Motivated by the aforementioned interest to develop “zero-emissions”
micro-/mesoscale combustors for both large-scale and portable power generation, the present doctoral thesis undertakes a numerical study of hetero/homogeneous combustion of hydrogen and syngas in mesoscale channels.
The themes that are tackled are two: in the first part, instabilities occurring in lean syngas combustion in inert mesoscale channels are investigated,
while in the second part a fully transient direct numerical simulation code
for catalytic combustion is extended to account for conjugate heat transfer
and is used to study dynamic processes (light-off and gas-phase ignition) in
catalytic combustion of lean hydrogen in mesoscale channels. The scientific
relevance and the novelty behind these two studies is summarized below.
When combustion occurs at a small scale, the increase in surface-tovolume ratio and the reduction of the characteristic thermal inertia time of
the combustor structure will significantly enrich the combustion phenomena
via reduced mixing and enhanced wall heat loss, and wall-flame thermal and
kinetic couplings, which lead to various new phenomena such as flameless
combustion, weak flame propagation, flame street formation, kinetic extinction, pulsating and spinning instabilities [2]. A review of flame dynamics
8
1.1. Motivation
in micro-/mesoscale channels is provided in Chapter 3. The dynamics of
syngas fuels in spatially multidimensional confined systems has not been
investigated yet. The present study aims at extending the knowledge about
flame dynamics in small confinements by undertaking a detailed numerical
investigation of lean syngas/air flame dynamics in a mesoscale inert (noncatalytic) 2-D channel with a height of 7 mm. The code is first validated
against steady experimental flames where the OH* chemiluminescence signal is detected, while parametric numerical studies are carried out by varying the wall temperature and the CO:H2 ratio. The fast OH-LIF technique
is finally used to assess the capability of the numerical code in reproducing
both steady and dynamic combustion modes.
The development of advanced numerical catalytic combustion models poses severe challenges due to the different physicochemical processes
with disparate spatio-/temporal scales that must be accounted for: surface
and gas-phase chemistry, convective and diffusive heat and mass transport
within the bulk flow (interphase transport), intraphase transport within the
catalyst washcoat, heat conduction in the channel wall and radiative heat
transfer. The corresponding time scales range from fractions of milliseconds
(surface and gaseous chemical times) to seconds (solid axial heat conduction time). Steady and transient models are required for the description
of an individual catalytic channel during steady operation and light-off, respectively. Given the long characteristic solid wall heat-up time compared
to the chemical, convective and diffusive time scales of a catalytic reactor, the quasisteady-state (QSS) assumption for the gas-phase and surface
9
1. Introduction
chemistry is often invoked. A review of available numerical models is provided in Chapter 4. In this thesis a fully transient (without invoking the
quasiteady-state assumption) 2-D numerical study of fuel-lean hydrogen/air
hetero-/homogeneous combustion in a planar mesoscale (2 mm height) Ptcoated channel is undertaken, using detailed gas-phase and catalytic chemistry, and 2-D heat conduction in the solid wall. Fast processes such as
light-off, homogeneous ignition and hetero-/homogeneous combustion interactions during dynamic operation, with time scales ranging from fractions
of milliseconds to a few milliseconds, were studied for the first time given
the model’s capability to resolve all relevant time scales.
1.2
Thesis Structure
The thesis is organized as follows. In Chapter 2, the mathematical model
and the employed numerical methodology are briefly described, while the
code extension for conjugate heat transfer is validated against analytical and
numerical data available in literature. In addition, the Computational Singular Perturbation (CSP) analysis, a computational tool used extensively
for post-processing purposes, is introduced. Chapter 3 presents results of
the detailed numerical study on combustion dynamics in fuel-lean premixed
syngas/air mixtures in a mesoscale inert channel. The steady and oscillatory
modes revealed by the numerical study are compared against experimental
results obtained both by detecting the OH* chemiluminescence signal and
with the fast OH-LIF technique. Chapter 4 is dedicated to the results
10
1.2. Thesis Structure
of the fully transient detailed numerical study on fuel-lean hydrogen/air
hetero-/homogeneous combustion in a Pt-coated mesoscale channel. Finally, Chapter 5 summarizes the main conclusions and delineates possible
directions for future studies and code developments.
11
Chapter 2
Numerical Approach
This chapter describes the mathematical model, the discretization method
and the Direct Numerical Simulation (DNS) solver used for the detailed
numerical simulations presented in Chapter 3 and Chapter 4. Validation
of the numerical code is also provided. Finally, the basic concepts of CSP
analysis, a computational tool that was used to extract physical insights
from the numerical results, are introduced and explained.
2.1
2.1.1
Governing equations
Gas-phase
Several combustion applications are characterized by low Mach number
flows (M a 1) with compressibility effects due to strongly exothermic
chemical reactions. Both compressible and incompressible schemes are not
suitable to these conditions, unless they are modified [40]. The numerical
13
2. Numerical Approach
schemes developed for solving low-Mach-number compressible flows try to
retain compressibility effects due to heat release and changes in composition,
while acoustic waves are filtered out. In this way a significant increase of
the integration time step can be achieved since the fast acoustic time scales
are not resolved. The mathematical formulation employed for the present
study stems from the asymptotic analysis of Majda and Sethian [41] which
leads to a decomposition of the pressure
p(x, t) = p0 (t) + p1 (x, t),
(2.1)
where p1 /p0 = O(M a2 ), while the hydrodynamic pressure, p1 , is decoupled
from the thermodynamic pressure computed with the equation of state.
The low-Mach-number governing equations for a reactive gaseous mixture in an open system, ignoring body forces and radiative heat transfer,
are:
Continuity
∂ρ
+ ∇ · (ρu) = 0
∂t
Momentum
ρ
∂u
+ u · ∇u
∂t
= −∇p1 + ∇ · (µS)
S = ∇ u + (∇ u)T −
14
(2.2)
2
(∇ · u) I
3
(2.3)
(2.4)
2.1. Governing equations
Energy
ρcp
∂T
+ u · ∇T
∂t
= ∇ · (λ∇T ) −
Ng
X
hi ω˙i − ρ
i=1
cp =
Ng
X
Ng
X
!
cp,i Yi Vi
· ∇T (2.5)
i=1
cp,i Yi
(2.6)
i=1
Gas-phase species
ρ
∂Yi
+ u · ∇Yi
∂t
= −∇ · ρYi Vi + ω˙i i = 1, ..., Ng
(2.7)
with Ng the total number of gaseous species. In Eqs. (2.2) to (2.7), hi , ω˙ i , Yi ,
V i , Wi , cp,i are the enthalpy, chemical mass production term, mass fraction,
diffusion velocity vector, molecular weight, and heat capacity of species i,
respectively, while λ is the thermal conductivity, µ the dynamic viscosity, ρ
the density, and cp the heat capacity of the mixture at constant pressure. I
is the identity matrix.
For a detailed gas-phase reaction mechanism with Nrg elementary reactions and Ng species χi
Ng
X
i=1
νij0 χi
Ng
X
νij00 χi ,
j = 1, ..., Nrg
(2.8)
i=1
the rate of progress of reaction j is function of the product of the molar
concentrations of the species [Xi ] = ρYi /Wi according to the law of mass
15
2. Numerical Approach
action
qj =
kjf
Ng
Y
0
νij
[Xi ]
−
kjr
i=1
Ng
Y
00
[Xi ]νij ,
(2.9)
i=1
where νij0 , νij00 are the stoichiometric coefficients of species i in the forward
and reverse direction, respectively. The forward rate constant kjf of reaction
j is given by an Arrhenius-type expression
kjf
= Aj T
βj
exp
−Ej
RT
,
(2.10)
where the pre-exponential factor Aj , the temperature exponent βj , and
the activation energy Ej are specified in the detailed reaction mechanism;
the forward and reverse reaction rate constants are related through the
equilibrium constant Kc,j (T ) = kjf /kjr . The production/destruction rate of
species i is the sum of the rates of progress qj of all reactions j = 1, . . . , Nrg
involving the ith species:
ω˙ i = Wi
Nrg
X
νij qj ,
(2.11)
j=1
with νij = νij00 − νij0 .
The diffusion velocities Vi of gas-phase species are given, according
to the kinetic theory of gases, by the Stefan-Maxwell equation [42] which
requires the solution of the (Ng × Ng ) multicomponent system at every
grid point. In order to reduce the computational cost, the HirschfelderCurtiss approximation [43] is used, which accounts for the first term in the
convergent series expansion and is the best first-order approximation to the
16
2.1. Governing equations
full matrix diffusion as shown in [44]. V i is given by
V i = V˜i + V c ,
(2.12)
where the first term in the right-hand side of Eq. (2.12), by neglecting Soret
and Dufour effects and by using the mixture-averaged diffusion model [45],
becomes:
1 − Yi
,
j6=i Xj /Dji
Di = P
V˜i = − (Di /Xi ) ∇Xi ,
(2.13)
with Dji , Di and Xi being the ith species binary diffusion coefficients,
mixture-averaged diffusivity and mole fraction, respectively. A correction
velocity,
Vc =−
Ng
X
Yi V˜i
(2.14)
i=1
needs to be introduced [46] in order to ensure global mass conservation, i.e.
the expressions
Ng
X
Yi ≡ 1,
i=1
Ng
X
Yi V i ≡ 0,
(2.15)
i=1
are satisfied. The gaseous mixture is considered to be ideal and the equation
of state is given by the ideal gas law:
Equation of state (ideal gas law)
p0 =
ρRT
¯
W
(2.16)
17
2. Numerical Approach
where
¯ =
W
Ng
X
Yi
Wi
i=1
!−1
(2.17)
is the mean molecular weight of the mixture, and R the universal gas constant.
In the low-Mach-number formulation, the continuity equation (Eq. (2.2))
is replaced by Eq. (2.18):
Ng
X
1 Dρ
1 DT
1 DYi
¯
∇·v =−
=
+W
=: Qt
ρ Dt
T Dt
Wi Dt
i=1
(2.18)
where Qt is the so-called “thermal divergence”. By substituting energy (Eq.
(2.5)) and species (Eq. (2.7)) equations, Eq. (2.18) becomes:
+
2.1.2
1
ρcp T
Ng
¯
1XW
∇·u=
−∇ · ρYi Vi + ω˙ i +
ρ i=1 Wi
"
!
#
Ng
Ng
X
X
∇ · (λ∇T ) −
hi ω˙ i − ρ
cp,i Yi Vi · ∇T =: Qt
i=1
(2.19)
i=1
Catalytic wall
Within the solid walls, the energy equation under the assumption of constant themophysical properties becomes:
ρs cs
∂T
= λs ∇ 2 T
∂t
(2.20)
where ρs , cs and λs are the density, specific heat capacity and thermal
conductivity of the solid, respectively. At the catalytic coating, the surface
18
2.1. Governing equations
species coverage (θNg +k ) equations are solved
∂θNg +k
s˙ N +k
= σNg +k g
k = 1, ..., Ns
∂t
Γ
(2.21)
with Γ the surface site density of the catalyst, σNg +k the site occupancy
of species Ng + k, s˙ Ng +k the catalytic molar production rate of the Ng +
kth surface species and Ns the number of surface species. The production/destruction of gas-phase species at the catalytically active wall via
adsorption/desorption reactions introduces a flux boundary condition at
the gas-wall interface:
n · [ρYi (u + V i )]+ = s˙i Wi i = 1, ..., Ng
(2.22)
where (+) indicates the properties just above the gas-wall interface, n is
the unit outward-pointing vector normal to the surface and u is the Stefan
velocity given by
1X
n · u|+ =
s˙i Wi i = 1, ..., Ng
ρ i=1
Ng
(2.23)
The interfacial energy condition is:
n · λ∇T |+ −
Ng
X
i=1
X
Ng +Ns
n · [ρYi (u + V i )hi ]+ =
s˙k Wk hk + n · λs ∇T |−(2.24)
k=Ng +1
By substituting boundary condition of Eq. (2.22) into the flux term on the
left-hand side of Eq. (2.24), a more compact form of the interfacial energy
19
2. Numerical Approach
balance can be obtained:
X
Ng +Ns
n · λ∇T |+ =
s˙i Wi hi + n · λs ∇T |−
(2.25)
i=1
The surface molar reaction rates s˙i are computed from the law of mass action
in the same way as the gas-phase reaction rates (see Eq. (2.8) to Eq. (2.11))
with the surface species molar concentrations computed as [Xk ] = Γ θk /σk .
Some species coverages can modify the Arrhenius expression given in Eq.
(2.10); in this case a modified expression is derived, which is function of
three coverage parameters ηkj , µkj and εkj for species k and reaction j [47]:
kjf
= Aj T
βj
exp
−Ej
RT
NY
g +Ns
µ
10ηkj θk θk kj exp[
k=Ng +1
−εkj θk
]
RT
(2.26)
where the product in Eq. (2.26) considers only surface species responsible for the coverage modification. For some simple surface reactions the
surface reaction rate is specified as a function of a “sticking coefficient” γj
(probability that a process takes place when a collision occurs). The sticking coefficient γj can be converted into the usual mass-action kinetic rate
constant with the relation [47]:
kjf
γj
= m
Γ
r
RT
2πWk
(2.27)
where Wk is the molecular weight of the gas-phase species and m is the sum
of all the stoichiometric coefficients of reactants that are surface species. It
is assumed that the sticking coefficient is much less than one, so that the
20
2.2. Numerical method
molecular motion close to the solid surface is random and the presence of
the solid surface does not affect the collision frequency of gas-phase species
with the surface itself [47]. When the sticking coefficient is close to one,
species which are carried close to the surface due to their random motion
have a high probability of staying there [47]. This phenomena causes a nonMaxwellian velocity distribution and a correction factor provided by Motz
and Wise [48] is applied, resulting in a modified expression of the forward
rate constant [47]:
kjf
2.2
γj
1
=
1 − γj /2 Γ m
r
RT
2πWk
(2.28)
Numerical method
The governing equations presented in Section 2.1 were discretized in space
using the spectral element method (SEM). The SEM is a high-order weighted
residual technique developed by Patera and coworkers [49] in the 80’s that
couples the tensor product efficiency of global spectral methods with the
geometric flexibility of the finite element method (FEM) [40]. The SEM has
usually a much higher approximation order (typically N = 5÷15) compared
to the FEM, and it yields exponentially convergent dispersion error [40].
The SEM features excellent transport properties (minimal numerical diffusion and dispersion) for a significantly larger fraction of the resolved modes
than is possible with the higher order finite difference and finite volume
methods [50]. The starting point of SEM discretization is the Galerkin formulation where the set of partial differential equations is converted into its
21
2. Numerical Approach
weak form. The weak formulation is then discretized by replacing integration with quadrature and by constructing the approximation subspace using
polynomials of order N . In our case, the computational domain is split into
quadrilateral conforming elements and the solution and geometry are expressed as sums of N th order tensor product Lagrange polynomials based
on the Gauss-Lobatto-Legendre quadrature points. On a parallel computer
architecture, each processor accommodates a subset of elements where only
a nearest neighbor communication at partition boundaries is needed [40].
A complete reference on the subject can be found in the books of Gottlieb
and Orszag [51], Deville et al. [50] and Canuto et al. [52].
The discretized equations are solved with a parallel code based on
the incompressible flow solver nek5000 [53] by using a high-order splitting
scheme for low-Mach-number reactive flows based on the formulation of
Tomboulides et al. [54]. Since detailed kinetic mechanisms are usually quite
stiff, with time scales being orders of magnitude faster compared to the fluid
motion, it appears reasonable to split the low-Mach-number system of equations (Eq. (2.2) to (2.7)) into two subsystems, namely the hydrodynamic
subsystem (continuity and momentum equations) and the thermochemistry
subsystem (energy and species equations). The two subsystems are coupled
through velocity in the species and energy equations, while thermophysical properties appearing in the flow equations depend on temperature and
mixture composition. Furthermore the thermal divergence constrains the
evolution of the flow field. The result is a large system of coupled non-linear
partial differential equations. The algorithm decouples the two subsystems.
22
2.3. Code validation: Conjugate Heat Transfer
First, the nonlinear thermochemistry subsystem is solved implicitly using a
stiff integrator and an extrapolated velocity field. The spatial discretization
of the thermochemistry subsystem yields an ordinary differential equations
(ODE) inital value problem which is integrated in time with the variable
qth-order (q = 1, ..., 5), variable-step multistep integrator CVODE [55],
based on stiffly stable backward differentiation formulas (BDF). The thermodynamic state of the system is then known. During the next substep
the hydrodynamic subsystem is advanced in time using a projection type
velocity correction scheme introduced by Orszag et al. [56]. A semi-explicit
integration approach is used for the hydrodynamic subsystem.
Thermodynamic properties and reaction rates, of both gas-phase and
surface species, are evaluated using CHEMKIN [47, 57], while transport
properties are computed using CHEMKIN’s transport library [45].
2.3
Code validation: Conjugate Heat Transfer
Details about the validation of the incompressible and the low-Mach-number
non-reactive flow solvers can be found in the PhD thesis of Kerkemeier
[40]. Furthermore, the low-Mach-number reactive flow solver with gas-phase
chemistry was successfully validated in [40] against the results of a laminar
planar premixed lean hydrogen/air flame computed with the well established PREMIX code of CHEMKIN [58]. Finally, the catalytic model implemented during the completed PhD thesis of Pizza [59] has been validated
in [17] against steady catalytic combustion results from a well established
23
2. Numerical Approach
code developed in [60] which successfully simulated experiments in lean
methane [61], lean hydrogen [62], and more recently in lean CO/H2 [63]
and in rich hydrogen [64] over Pt. For the present work the capability of
the low-Mach-number catalytic solver was extended by inclusion of conjugate heat transfer. Validation of conjugate heat transfer under steady and
transient operation is reported in the next two sections.
2.3.1
Validation of the steady-state solution
Validation of steady-state conjugate heat transfer solutions was carried out
by comparing the computed interfacial gas-solid temperature TW with analytical results from [65]. The test case was a 2-D channel with solid walls
and constant thermophysical properties for both fluid and solid. At the
inflow, fully developed laminar velocity profile and constant temperature
TIN were specified. At the outflow, zero Neumann boundary conditions
were imposed. The outer wall temperature, Tb > TIN , was fixed, while the
vertical front solid wall faces were adiabatic.
Figure 2.1 shows the nondimensional interfacial gas-solid temperature,
ΘW = (TW − TIN )/(Tb − TIN ), as a function of the nondimensional streamwise coordinate x/h, where h is the channel height, for solid to fluid thermal
conductivity ratios λs /λ equal to 100 (Fig. 2.1(a)) and 1000 (Fig. 2.1(b)).
A solid to fluid thermal conductivity ratio of practical interest is within
the range delimited by the two selected values. The agreement between
simulated (solid line) and analytical (symbols) results is excellent.
24
2.3. Code validation: Conjugate Heat Transfer
Figure 2.1: Comparison between simulated (solid line) and analytical (symbols) results for the conjugate heat transfer problem in [65] and for solid to
fluid thermal conductivity ratios λs /λ of (a) 100 and (b) 1000. The abscissa
is the streamwise coordinate nondimensionalized with the channel height,
while the ordinate is the nondimensional interfacial gas-solid temperature
ΘW = (TW − TIN )/(TIN − Tb ).
2.3.2
Validation of the transient solution
Validation of transient conjugate heat transfer solutions was carried out
by comparing the computed interfacial gas-solid temperature TW with the
finite difference solution presented in [66]. The test case was a 2-D channel
with constant thermophysical properties for both fluid and solid. At the
inflow, fully developed laminar velocity profile and constant temperature
TIN were specified. At the outflow, zero Neumann boundary conditions
were imposed. The outside boundaries of the solid walls were adiabatic.
The initial temperature of the whole domain was TIN and the initial velocity
25
2. Numerical Approach
was equal to the inflow velocity profile. The channel walls were heated up
through a source term qW in the energy equation that was sinusoidal in the
axial coordinate and could be either constant (Eq. (2.29)) or exponential
(Eq. (2.30)) in time:
qW
= f (X) = (1 + CR sin(πKX))
qg
(2.29)
qW
= f (X, F ) = (eDF − 1)sin(πKX)
qg
(2.30)
where qg is a reference heat generation rate per unit volume that can
be chosen arbitrarily, while the nondimensional axial coordinate is X =
(αx)/[um (h/2)2 ], with α being the thermal diffusivity of the fluid and um
the average fluid velocity, and the nondimensional time is F = αt/(h/2)2 ,
with t being the physical time. The values of the parameters in Eq. (2.29)
and Eq. (2.30) are CR = 1.0, K = 1.0 and D = 4.0. The first generation
term (Eq. (2.29)) is a model for a steady-state generation function, while
the second one (Eq. (2.30)) can be used to model either an initial startup
or a change in the steady-state power level [66].
Figure 2.2 reports the nondimensional interfacial gas-solid temperature,
ΦW = λ(TW − TIN )/[qg δ(h/2)] as a function of the nodimensional axial
coordinate X. The source term was sinusoidal in space and exponential
in time (Eq. (2.30)) for Fig. 2.2(a) while it was sinusoidal in space and
constant in time (Eq. (2.29)) for Fig. 2.2(b). The temperature in Fig.
2.2(a) is plotted at different nondimensional times F , and for a solid-to-
26
2.3. Code validation: Conjugate Heat Transfer
Figure 2.2: Comparison between computed (solid lines) and finite difference [66] (symbols) nondimensional gas-solid interfacial temperature
ΦW = λ(TW − TIN )/[qg δs (h/2)] for the conjugate heat transfer problem
in [66]. The comparison in (a) refers to a solid to fluid thermal capacity
ratio B = (ρs cs δs )/[ρcp (h/2)] equal to 1.0 and to different nondimensional
time instants, F = αt/(h/2)2 , while the comparison in (b) refers to F = 0.2
and to different values of B. The nondimensional axial coordinate is defined
as X = (αx)/[um (h/2)2 ]. The heat generation source term was sinusoidal in
the streamwise coordinate and exponential in time (Eq. (2.30)) in (a), and
sinusoidal in the streamwise coordinate and constant in time (Eq. (2.29))
in (b).
fluid thermal capacity ratio B = (ρs cs δs )/[ρcp (h/2)] equal to 1, while the
temperature in Fig. 2.2(b) is plotted at F = 0.2 and for different values of
B. The comparison between the results of our simulations (solid lines) and
those of the finite difference code (symbols) in Fig. 2.2 shows an excellent
agreement.
27
2. Numerical Approach
2.4
Computational Singular Perturbation Analysis
The Computational Singular Perturbation (CSP) method [67,68], originally
developed for the reduction of detailed reaction mechanisms, has been proposed as a tool to analyze different types of flames. Applications include the
interaction of a 2-D counter-rotating vortex-pair with a premixed methaneair flame [69], a methane-air edge flame [70] and n-heptane/air triple flames
in partially premixed mixing layers [71]. Chemical Explosive Mode Analysis
(CEMA), a variant of CSP, was recently proposed to investigate the modes
with positive eigenvalues [72], which are characterized by an exponential
increase in time, and was used to detect the location of premixed flame
fronts and to identify auto-ignition as dominant stabilization mechanism
for turbulent lifted hydrogen [72] and ethylene [73] jet flames. To the best
of the author’s knowledge, CSP analysis has never been applied to study
catalytic combustion. In the next section, the CSP methodology will be
introduced for the general case where both gas-phase and surface reactions
occur, while the case where only homogeneous chemistry is present can be
derived as a particular case by removing the heterogeneous chemistry part.
CSP analyis has been used to analyze both gas-phase syngas combustion in
Chapter 3 and hetero-/homogeneous hydrogen combustion in Chapter 4.
28
2.4. Computational Singular Perturbation Analysis
2.4.1
Methodology
The conservation equations for energy (Eq. (2.5)) and species (Eq. (2.7))
can be written in the CSP notation as
∂w
= L(w) + g(w) = f
∂t
(2.31)
where w is the vector of temperature, gas-phase species mass fractions and
surface species coverages (T, Y1 , . . . , YNg , θNg +1 , . . . , θNg +Ns ), L represents
the convective-diffusive terms, and g is the chemical source term
g(w) =
Nrg
X
k
S k R (w) +
k=1
Nrs
X
S j Rj (w)
(2.32)
j=1
where Nrg and Nrs are the number of gas-phase and surface reactions, respectively, while Rk and Rj are the net rates of progress of the k th gaseous
reaction and of the j th surface reaction, respectively. S k and S j are the
stoichiometric vectors of the k th gaseous reaction and of the j th surface
reaction, respectively, defined as:


PNg νik Wi hi
− i=1 ρcp 


ν1k W1




ρ


.


..




 νNg k WNg

Sk = 

ρ






0




.


..




0
(2.33)
29
2. Numerical Approach


PNg +Ns A νij Wi hi
V
ρcp 
− i=1


ν
W
A 1j 1




V
ρ


..


.






A νNg j WNg
Sj = 

V
ρ




ν(Ng +1)j
 σN +1

g


Γ


..


.




ν(Ng +Ns )j
σNg +Ns
Γ
(2.34)
where A is the catalytic surface area and V is the volume, while νik and νij
are the net stoichiometric coefficients of species i in gaseous reaction k and
surface reaction j, respectively.
The right-hand side of Eq. (2.31) can be decomposed into fast r and
N +Ns +1
g
slow s modes by using the right {ai }i=1
N +Ns +1
g
and left {bi }i=1
eigenvec-
tors of the Jacobian matrix of the source term g, which provide the firstorder approximation of the CSP basis vectors/co-vectors [74]. The modal
decomposition of Eq. (2.31) becomes
f = f r (w) + f s (w) =
M
X
r=1
X
Ng +Ns +1
r
ar h +
as hs
(2.35)
s=M +1
with M the number of fast modes and hn = bn · f the mode amplitudes.
The different modes identified by the eigenvalues and eigenvectors of
the source term Jacobian, can be analyzed by computing the CSP pointer
30
2.4. Computational Singular Perturbation Analysis
Qki of species i with respect to mode k
ai bki
Qki = PNg +Nks +1
i=1
|aik bki |
(2.36)
which measures the projection of the vector of species i in composition
space on the eigenvector ak [71] and indicates how parallel the axis of the
ith species is to the corresponding mode. Pairs of complex conjugate modes
are converted into pairs of real modes [75]
hn = bnR · f ,
hn+1 = bnI · f
(2.37)
where bnR and bnI are the real and imaginary part of the eigenvector of the
nth mode, respectively. The CSP pointers associated with a pair of complex
conjugate modes are
2ai bniR
Qni = PNg +NnR
,
s +1
|ain bni |
i=1
2ainI bniI
Qn+1
=
P
i
Ng +Ns +1 i n
|an bi |
i=1
(2.38)
CSP also introduced the participation indices of the chemical reactions
k with respect to mode j
Pkj
(bj · S k )Rk
= PNrg +Nrs
k0 =1
|(bj · S k0 )Rk0 |
(2.39)
For the gas phase, participation indices can be defined also for a transport
process l
bjl Ll
P
j l0
j
k0
l0 |bl0 L | +
k0 |(b · S k0 )R |
Plj = P
(2.40)
31
2. Numerical Approach
where Ll is the diffusive/convective term in the right-hand side of energy
and gas-phase species equations (Eq. (2.5) and Eq. (2.7)).
The eigenvectors are chosen in the direction of increasing temperature
since it was found that this convention leads to participation index signs in
agreement with the signs of the sensitivity coefficients [71]. The participation indices determine the contribution of the chemical reactions and of the
transport processes to the mode amplitudes.
32
Chapter 3
Flame dynamics in lean
CO/H2/air combustion 1
3.1
Review on flame dynamics
The oxidation mechanisms of hydrogen and carbon monoxide are of fundamental importance, as they are subsets of all hydrocarbon reaction mechanisms [76]. It is well known that the oxidation of CO requires traces of
hydrogen-containing species. In the presence of a small amount of hydrogen or water, OH radicals are formed leading to the oxidation of CO and
to the production of H radicals, which feed into the chain-branching mechanism [76].
1
Sections 3.1, 3.2 and 3.3.1 to 3.3.4 are published in “A. Brambilla, C. E. Frouzakis, J.
Mantzaras, R. Bombach, K. Boulouchos. Flame dynamics in lean premixed CO/H2 /air
combustion in a mesoscale channel, Combust. Flame, 161:1268-1281, 2014” while Section 3.3.5 was accepted for oral presentation at the 35th International Symposium on
Combustion as “A. Brambilla, M. Schultze, C. E. Frouzakis, J. Mantzaras, R. Bombach,
K. Boulouchos. An experimental and numerical investigation of premixed syngas combustion dynamics in mesoscale channels with controlled wall temperature profiles, 35th
International Symposium on Combustion, San Francisco, August 3-8 2014”.
33
3. Flame dynamics in lean CO/H2 /air combustion
Syngas oxidation exhibits complex dynamic behaviors even in simple
homogeneous systems such as closed vessels or continuous stirred-tank reactors (CSTR). Experimental studies on the oxidation of carbon monoxide
with different amounts of hydrogen in both closed [77, 78] and open systems [79] revealed the existence of four reaction modes: slow oxidation
without light emission (steady dark reaction), steady and oscillatory reactions accompanied by chemiluminescence (steady and oscillatory glow, respectively), and a fourth mode characterized by oscillatory ignition with
large, transient temperature excess and complete fuel consumption (oscillatory ignition). Oscillatory ignition appeared only in the presence of more
than 1500 ppmv of H2 [79]. The oscillatory behavior was also revealed by
theoretical studies [80] and CSTR simulations [81]. Other studies on syngas
combustion in open systems have also reported complex dynamic phenomena. In the experimental investigation of Johnson et al. [82], combustion
of CO with various amounts of hydrogen in a CSTR showed bistability,
oscillatory behavior and a period-doubling cascade to chaos. Simple and
complex oscillations were also reported for pure hydrogen combustion in a
CSTR [83–86].
In practical channel-flow reactors, in addition to chemical kinetics, diffusive and advective transport as well as interactions with the walls, mainly
in terms of heat losses and termination reactions [12], become important.
These effects are exemplified at the micro- and meso-scales, which are characterized by high surface-to-volume ratios. Premixed flames in channels
experience hydrodynamic instabilities [87, 88], caused by the density jump
34
3.1. Review on flame dynamics
across the flame front, and diffusional-thermal instabilities [89], driven by
the imbalance between heat and mass transport of the deficient reactant.
Heat losses are an additional destabilizing factor which narrows down the
combustion stability envelop [90].
Numerical studies of flame dynamics in micro- and meso-scale channels
have also revealed a wealth of combustion modes. Pizza et al. investigated
lean premixed hydrogen/air flames in micro- [91] and meso-scale [92] planar channels at atmospheric pressure and a prescribed wall temperature
by means of 2-D simulations using detailed chemistry and transport. The
impact of channel height and inlet velocity on flame dynamics was assessed
and stability maps were constructed. In micro-scale channels, a rich variety
of dynamics was observed with increasing channel height. For a channel
height h = 1 mm, four different flame modes were observed with increasing inlet velocity [91]: repetitive ignition/extinction, V-shaped stationary
flames, asymmetric stationary flames, oscillating flames and again asymmetric stationary flames. Oscillating flames were suppressed for channels
narrower than h = 0.8 mm, while asymmetric flames vanished for h < 0.7
mm. Simulations in mesoscale channels with h = 2, 4 and 7 mm [92]
revealed mild combustion, repetitive ignition/extinction, “closed” (singlebranch) and “open” (double-branch) symmetric flames and oscillating flames
exhibiting harmonic, mixed mode or chaotic behavior [93]. A simpler thermodiffusive model with a single-step irreversible reaction and unity Lewis
number [94] reproduced some of these modes, indicating that they were not
necessarily of diffusive or kinetic origin but rather a result of flame-wall in-
35
3. Flame dynamics in lean CO/H2 /air combustion
teractions. An increase of the wall temperature was shown to have a strong
impact, dampening the dynamics [94]. The effect of Lewis number and flow
rate on the linear stability of symmetric and asymmetric flames propagating
in adiabatic channels was addressed in [95] using a thermodiffusive model.
Experimentally, Dogwiler et al. [61] observed asymmetric flames in lean
methane/air combustion in an inert (non-catalytic) channel, which were anchored on either the upper or the lower wall. The sensitivity of the asymmetric flame to external perturbations resulted in transitions between the
two flame structures at random time intervals. Maruta et al. [96] performed
experiments with lean premixed methane/air and propane/air mixtures in a
2 mm diameter tube. For a propane/air mixture with an average velocity of
30 cm/s and equivalence ratio ϕ = 0.5, they reported different flame types
depending on the wall temperature: stable flames, pulsating flames and a
combination of repetitive ignition/extinction with pulsating flames. For the
lean methane/air mixture, stable flames were identified at both high and
low velocities (a vigorously burning flame and a mild combustion mode, respectively), while dynamic modes were reported at intermediate velocities
in the form of repetitive ignition/extinction and a combination of repetitive
ignition/extinction with pulsating flames. Jackson et al. [97] studied numerically a problem very similar to that of [96], showing instabilities characterized by regular oscillations or by a violent process of ignition and extinction.
The repetitive ignition/extinction mode was also experimentally observed
by Richecoeur and Kyritsis [98] for a non-premixed methane/oxygen flame
in a 4 mm diameter curved duct.
36
3.1. Review on flame dynamics
Heat loss effects on flames propagating in channels were studied by
Cui et al. [99]. They investigated pulsating instabilities for Le > 1 flames
using a thermodiffusive model and concluded that heat losses to the walls
promoted flame oscillations by reducing the critical Lewis number for flame
stability, so that oscillations could be observed at near-extinction conditions.
Daou and Matalon [100] studied numerically the impact of conductive heat
losses on flames propagating in channels and reported total flame extinction
in narrow channels and partial flame extinction in wider channels, with the
flame extinguishing only at the walls in the latter case. Using a thermodiffusive model, Chakraborty et al. [101] studied the interplay between non-unity
Lewis number and heat losses in laminar premixed flames propagating in a
channel. They observed simultaneous tip-opening and dead (extinguished)
space near the wall for Lewis numbers significantly below unity and for
high heat losses, which in turn gave rise to a multicellular flame with a
“funnel-like” shape.
The dynamics of syngas fuels in spatially multidimensional systems
has not been investigated yet. The present study undertakes a detailed
numerical investigation of lean syngas/air flame dynamics in a meso-scale
channel with a height of 7 mm. Parametric studies are carried out by varying the wall temperature and the CO:H2 ratio. The numerical results are
compared against experiments based on both the OH* chemiluminescence
signal detection and the fast OH-LIF technique. The chapter is organized
as follows. The computational domain, boundary conditions and chemical
reaction mechanism are presented in Section 3.2. The numerical results
37
3. Flame dynamics in lean CO/H2 /air combustion
follow in Section 3.3. Comparisons between computed steady flames and
OH* chemiluminescence images is given in Section 3.3.1 while in Sections
3.3.2.1 to 3.3.2.5 the various combustion modes computed by varying the
wall temperature and for a CO:H2 molar ratio of rCO/H2 = 7 are elaborated
upon. The effect of varying rCO/H2 is discussed in Section 3.3.3 and finally
the Computational Singular Perturbation analysis is presented in Section
3.3.4. Recent experimental results obtained with the fast OH-LIF technique
are compared against the numerical combustion modes in Section 3.3.5.
3.2
Computational domain, boundary conditions
and chemistry
A schematic of the computational domain with the boundary conditions is
shown in Fig. 3.1.
Figure 3.1: Schematic of the computational domain (all lengths in mm).
For the species, flux boundary conditions are imposed at the inlet and the
wall is non-reactive (zero normal fluxes).
38
3.2. Computational domain, boundary conditions and chemistry
It consists of a 2-D channel of height h = 7 mm and length L depending
on the flame position and extent. At the inflow, a uniform velocity profile
UIN , constant inlet temperature TIN and flux boundary conditions for the
species [102]
ρYi (u + V i ) = ρuYi,ref
(3.1)
were imposed, where Yi,ref was the reference inlet composition. Along the
channel walls, no slip for both velocity components and zero flux for all
species (i.e. chemically inert walls) were applied. Although radicals could
recombine on the wall, such reactions, which depend strongly on the specific
material and surface treatment method, were not included, similarly to
previous studies which successfully simulated many of the experimentally
observed flame dynamics [91–93]. In this respect, a quenchless wall was
simulated, mimicking the nearly quenchless materials reported in [12].
The wall temperature TW was prescribed. When experimental conditions were simulated, the wall temperature profiles were obtained from
thermocouple measurements via cubic spline interpolation. In the parametric numerical studies, the wall temperature ramped up in the first 3.5 mm of
the domain (a length equal to the channel half-height) via a hyperbolic tangent function from the inlet temperature to the final value, which was kept
constant farther downstream. A similar wall temperature profile was used
in previous studies ( [91] and [92]). At the outlet, zero Neumann boundary
conditions (outflow in Fig. 3.1) were imposed. In all cases, the pressure
was atmospheric.
39
3. Flame dynamics in lean CO/H2 /air combustion
The resolution requirements were determined by calculating the thermal thickness [103] of a 1-D flame using the PREMIX code of CHEMKIN
[58]:
δT =
Tad − Tu
max(dT /dx)
(3.2)
with Tu the unburnt flame temperature, Tad the adiabatic flame temperature
and max(dT /dx) the maximum temperature gradient across the flame. The
grid and the polynomial order were chosen in order to have at least 10
points within the thermal flame thickness; more details will be provided in
Section 3.3. Grid independence of the solution was verified by repeating
the simulations with different polynomial orders.
In the following, non-dimensional quantities will be reported in terms
of the inlet velocity, a reference temperature Tref = 994 K that was representative of the wall temperature, the channel height, and the thermophysical
properties evaluated at Tref and the inlet composition.
The CO/H2 /O2 elementary scheme of Li et al. [104] was used for the
chemistry description with 32 non-duplicate elementary reversible reactions
and 12 chemical species.
3.3
Results and discussion
The predictive capability of the model and of the reaction mechanism was
firstly investigated against experiments at steady combustion based on the
40
3.3. Results and discussion
OH* chemiluminescence signal detection. Parametric simulations assessing
the impact of wall temperature and CO:H2 ratio on the observed flame dynamics were subsequently carried out. Finally, recent experimental results
with the fast OH-LIF technique are reported and compared against the
numerical simulations.
3.3.1
Steady experimental vs DNS results
Experiments were performed in an optically accessible channel reactor employed in earlier works [105, 106], which is here only briefly described while
the details are provided in Section 3.3.5.1 where the focus is on the experiments. The L = 300 mm long and W = 104 mm wide channel consists of two 9-mm thick Si[SiC] horizontal plates, separated by two vertical
quartz windows with a height of h = 7 mm. The large aspect ratio (W/h
≈ 15) yielded a nearly 2-D flow as elaborated in previous steady experiments [62]. Twelve thermocouples per plate, embedded 0.9 mm beneath
the channel surfaces, monitored the temperatures along the vertical symmetry axis, and the measured upper and lower wall surface temperatures were
used as boundary conditions in the simulations. Experimental uncertainties
on the ignition distance have been analyzed and found to be 2% in [62]. The
excited OH* chemiluminescence signal was collected from one side window
via an intensified CCD camera and a bandpass filter centered at 308 nm.
Channel lengths of 120 mm were recorded, sufficient to cover the combustion zone in the present experiments. The OH* chemiluminescence images
were subsequently compared with the computed heat release rate profiles,
41
3. Flame dynamics in lean CO/H2 /air combustion
HRR = −
PNg
i=1
hi ω˙i . Experiments in premixed counterflow methane-air
flames by Panoutsos et al. [107] demonstrated that the excited OH* can be
used as a heat release rate marker.
The first examined case corresponded to a steady flame at equivalence ratio ϕ = 0.42 (based on both fuel components), CO:H2 molar ratio rCO/H2 = 3, inlet velocity UIN = 3.21 m/s and inlet temperature
TIN = 335 K. The Reynolds number based on the inflow mixture viscosity, velocity and channel height was ReIN = 1130. The chemiluminescence
measurements yielded a flame anchored close to the channel inlet, such that
a computational domain with a length L = 105 mm was sufficient for the
simulations. The domain was discretized with 75 equally spaced elements
in the streamwise direction and 10 elements in the transverse direction, the
latter refined towards the wall. A polynomial order of 15 was chosen, resulting in 16 interpolating points per element and spatial direction. The
1-D flame thickness computed using PREMIX [58] was δT = 0.95 mm, the
laminar flame speed SL = 17.2 cm/s and the adiabatic flame temperature
Tad = 1578 K, such that the thermal flame thickness was resolved with ten
grid points.
The measured upper and lower wall temperatures are plotted in Fig.
3.2(a), while the measured OH* chemiluminescence and the computed heat
release rate distributions are shown in Figs. 3.2(b) and 3.2(c), respectively.
Predictions were in good agreement with the measurements in terms of the
flame anchoring position as well as the ensuing flame shape.
Comparison was also made with a second steady flame, which was
42
3.3. Results and discussion
Figure 3.2: (a) Measured wall temperature along the upper (black solid line)
and lower (red dashed line) walls, (b) measured OH* chemiluminescence and
(c) computed heat release rate for the steady flame at ϕ = 0.42, rCO/H2 = 3,
UIN = 3.21 m/s, TIN = 335 K.
anchored farther downstream from the inlet. The simulated domain had
a length of L = 175 mm, discretized with 125 equally-spaced elements
in the streamwise direction and 10 elements in the wall-normal direction,
refined towards the walls. At ϕ = 0.6, rCO/H2 = 36, UIN = 4.01 m/s, and
TIN = 323 K (ReIN = 1550), the corresponding 1-D flame had a thickness of
δT = 1.14 mm, a laminar flame speed SL = 18.4 cm/s and an adiabatic flame
temperature Tad = 1980 K. For a polynomial order of N = 15, the number of
grid points within the flame thickness was 13. The experimental asymmetric
flame structure in Fig. 3.3(b) captured also numerically in Fig. 3.3(c) was
attributed to the pronounced asymmetry between the upper and lower wall
43
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.3: (a) Measured wall temperature along the upper (black solid
line) and lower (red dashed line) channel walls, (b) measured OH* chemiluminescence and (c) computed heat release rate for the steady flame at
ϕ = 0.6, rCO/H2 = 36, UIN = 4.01 m/s, TIN = 323 K.
temperatures (Fig. 3.3(a)).
The comparisons between experimental and numerical results in Figs.
3.2 and 3.3 indicate that in both cases the model reproduced well the ignition locations as well as the flame shapes. The flame asymmetry was also
well reproduced by the model and the difference between measured and predicted ignition locations in Fig. 3.3 was 6 mm on the lower wall and 3 mm
on the upper one.
44
3.3. Results and discussion
3.3.2
Effect of wall temperature
Images of a dynamic mode identified experimentally for ϕ = 0.42, rCO/H2 =
7, UIN = 1.66 m/s, and TIN = 338 K (ReIN = 585) are shown in Fig. 3.4.
The experimental images, which integrate the OH* chemiluminescence across
the lateral direction, show an oscillatory flame with multiple branches oscillating out of phase. However, these images could not be used for comparing
the numerical oscillatory flames described in the following sections for different reasons. First, this would have necessitated a high-speed camera in
conjunction with a fast planar OH laser induced fluorescence (LIF) system,
so as to isolate potential 3-D effects of the dynamic modes, as it was done for
the experimental study reported in Section 3.3.5. In addition, such dynamics would also require 3-D simulations using temperature profiles measured
over the channel width, which were not available with the current experimental setup. In order to numerically investigate the processes leading to
the oscillatory behavior, a detailed numerical study was performed using
the same operating conditions and varying the wall temperature to cover
the measured range.
As mentioned in Section 3.2, the wall temperature (TW ) was ramped
up in the first 3.5 mm of the computational domain with a hyperbolic
tangent profile from TIN to a constant value in the range TW =954 K to
1274 K. The corresponding 1-D premixed flame thickness for the conditions
of Fig. 3.4 was δT = 1.16 mm, and the laminar flame speed and adiabatic
flame temperature were SL = 13.8 cm/s and Tad = 1594 K, respectively.
The computational domain length L = 105 mm was discretized with 75
45
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.4: (a) Measured wall temperature along the upper (black solid
line) and lower (red dashed line) channel walls, (b) to (d) measured OH*
chemiluminescence for an oscillatory flame at t = 64 ms, 76 ms and 87 ms,
respectively (ϕ = 0.42, rCO/H2 = 7, UIN = 1.66 m/s, TIN = 338 K).
equally spaced elements in the streamwise direction and 10 elements in the
transverse direction, refined towards the walls, and a polynomial order of
12 was chosen. When the flame extended to the outflow, simulations were
repeated with longer domains (up to L = 140 mm) to ensure a minimal effect
of the outflow boundary conditions. Apart for a slightly higher normalized
integral of the heat release rate (iHRR) over the whole domain
iHRR = −
Z X
Ng
ω˙i hi dV
(3.3)
V i=1
which was due to residual conversion in the exhaust gases, the results were
not affected by the domain length.
46
3.3. Results and discussion
Four combustion modes were identified: (i) no ignition for TW ≤ 960 K,
(ii) dynamic modes, consisting of oscillatory ignition with periodic behavior
for 964 K ≤ TW ≤ 1004 K and random strongly burning ignition spots
alternating with periods of weak combustion for 1006 K ≤ TW ≤ 1008 K,
(iii) weakly burning steady flames for 1010 K ≤ TW ≤ 1034 K, and finally
(iv) strongly burning steady V-shaped flames for 1039 K ≤ TW ≤ 1274 K.
The observed modes are discussed in the following.
3.3.2.1
No ignition: TW ≤ 960 K
The main characteristics of the non-ignited case at TW = 954 K are summarized in the top row of Table 3.2 (see Section 3.3.2.4). Insignificant amounts
of H2 and CO were converted (0.26% and 0.04%, respectively), with the conversion defined as
R
Conversion = 1 −
R
ρYi (u + V i ) · ndS
OU T
−
ρYi (u + V i ) · ndS
(3.4)
IN
n being the unit outward pointing normal vector.
Small amounts of radicals were produced near the walls, and their
mass fraction increased monotonically with increasing streamwise distance
but never became high enough to promote runaway for a channel length
of L = 105 mm. The temperature within the channel was nowhere higher
than TW and the iHRR reached a steady value of 2.3x10−3 , three orders
of magnitude lower than the iHRR obtained for the ignited cases. Similar
observations were made at TW = 960 K where the iHRR at steady state
47
3. Flame dynamics in lean CO/H2 /air combustion
attained a maximum value of 0.01.
For a longer domain with L = 140 mm, at TW = 960 K the higher residence time available for radical buildup resulted in the oscillatory ignition
described in the next section. Therefore, the presence of the non-ignited
case depended not only on TW but also on the domain length.
3.3.2.2
Oscillatory ignition: 964 K ≤ TW ≤ 1004 K
For wall temperatures in the range 964 K ≤ TW ≤ 1004 K, runaway within
the channel length L = 105 mm was obtained, characterized by periodic
ignition-extinction phenomena.
At the lowest considered wall temperature of TW = 964 K, the flame
was periodically blown out of the domain. The corresponding iHRR time
history (Fig. 3.5(a)) showed that the blow-off intervals (where iHRR dropped down to 0.01) were followed by periods of radical buildup along the hot
walls, which about 23 ms later resulted in re-ignition close to the outflow.
The periods of low heat release rate were characterized by incomplete fuel
conversion (Fig. 3.5(c)).
For TW ≥ 969 K, a flame existed within the channel during the whole
simulated time interval, as it could be inferred from the iHRR signals of
Figs. 3.5(b), 3.6 and 3.9(a), which never dropped down to the values of the
non-ignited case. Hydrogen and CO were almost fully consumed, although
at TW = 969 K the conversion of CO could still drop to about 60% for short
time intervals (Fig. 3.5(d)).
48
3.3. Results and discussion
(a)
(b)
(c)
(d)
Figure 3.5: iHRR for (a) TW = 964 K and (b) TW = 969 K; conversion of H2
(solid line) and CO (dashed line) for (c) TW = 964 K and (d) TW = 969 K.
Characteristic instantaneous flame structures and the time histories of
the iHRR with increasing wall temperature are shown in Fig. 3.6, while
the oscillation frequencies extracted from the fast Fourier transform of the
iHRR signal and the maximum temperatures reached during the ignitionextinction cycles are plotted in Fig. 3.7. The frequency of the ignition
events increased exponentially with increasing TW due to the faster radical
buildup close to the hot walls, while the maximum temperature in the domain reached its peak value for TW = 974 K. Two modes were observed
within the range of wall temperatures of the oscillatory behavior: an asymmetric mode characterized by one asymmetric flame propagating along the
upper or the lower wall during each ignition-extinction cycle (Figs. 3.6(a1),
3.6(b1), 3.6(e1) and Figs. 3.9(b) to 3.9(f)), which will be discussed in de49
3. Flame dynamics in lean CO/H2 /air combustion
(a)
(b)
Figure 3.6: Instantaneous distributions of YOH and time histories of the
iHRR at (a1),(a2) TW = 984 K (asymmetric oscillatory flame), (b1),(b2)
TW = 994 K (asymmetric oscillatory flame), (c1),(c2) TW = 994 K (symmetric oscillatory flame), (d1),(d2) TW = 1000 K (symmetric oscillatory
flame) and (e1),(e2) TW = 1004 K (asymmetric oscillatory flame).
tail below for the TW = 974 K case, and a symmetric mode where ignition
kernels appeared concurrently and at the same axial location on both walls
and resulted in two propagating flames (Figs. 3.6(c1) and 3.6(d1)). The
50
3.3. Results and discussion
Figure 3.7: (a) Frequency of the iHRR and (b) maximum temperature in
the computational domain achieved within an ignition-extinction cycle for
different TW . Solid and dashed lines represent asymmetric and symmetric
combustion modes, respectively.
symmetric mode firstly appeared at TW = 994 K and for this temperature
it coexisted with the asymmetric mode. The symmetric mode was characterized by fuel conversion (Fig. 3.8(b)) lower than the asymmetric one
(Fig. 3.8(a)) due to fuel leakage close to the channel midplane. The average iHRR of the asymmetric mode was about 3.7 (Figs. 3.6(a2), 3.6(b2)
and 3.6(f2)), which corresponded to complete fuel (H2 and CO) conversion,
while the symmetric mode was characterized by an average iHRR of 3.2
(Figs. 3.6(c2) and 3.6(d2)).
Different stages of the asymmetric ignition-extinction process for TW =
974 K are presented in Fig. 3.9. At time t1 a kernel ignited on the lower
51
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.8: Conversion of H2 (solid line) and CO (dashed line) for TW =
994 K and two coexisting modes: (a) asymmetric flame, (b) symmetric
flame.
wall at x ≈ 26 mm (Fig. 3.9(b)), but the peaks of temperature and radical
concentrations were located farther downstream, within a flame that had
ignited on the upper wall during the previous ignition-extinction cycle. The
iHRR at t1 was close to the minimum value (Fig. 3.9(a)). The flame propagated both upstream and downstream from the new kernel (Fig. 3.9(c)),
consuming rapidly the unburnt fuel trapped between the ignition spot and
the flame of the previous cycle. The iHRR at t2 approached the peak of
the entire cycle and the radicals reached their maximum concentrations at
the tip of the downstream propagating front (Fig. 3.9(c)). The flame of the
previous cycle was pushed downstream by the expanding gases from the
new propagating flame, and was weakened due to upstream fuel depletion.
Once the unburnt fuel was consumed, the downstream propagation
could not be sustained any longer, the iHRR started to decrease (t3 in
52
3.3. Results and discussion
Figure 3.9: Oscillatory flame at TW = 974 K. (a) time history of iHRR;
(b) to (f) YOH distributions at the time instants t1 to t5 marked in (a). The
white line in Figs. (d) and (e) is the YH2 = 10−4 isoline used to define the
flame front for the plots of Fig. 3.11.
Fig. 3.9(a)) and Tmax dropped sharply to the minimum value (Fig. 3.10(a)
shortly after t2 ). An unstable asymmetric flame was formed (Fig. 3.9(d))
with its lower and upper ends propagating in opposite directions. The flame
shrunk till a flame structure (which was the mirror image of the one at t1 )
53
3. Flame dynamics in lean CO/H2 /air combustion
formed with a new ignition kernel appearing on the upper wall (Fig. 3.9(f)).
Figure 3.10 shows the time series of the maximum temperature and of
the integral over the upper (solid line) and lower (dashed line) walls of the
local nondimensional heat transfer rate,
Z
λ∇T · n dx
iHT R =
(3.5)
L
Figure 3.10: Time histories of (a) the maximum temperature inside the
channel and (b) the non-dimensional integral heat transfer rate from the
fluid to the upper (solid line) and lower (dashed line) walls for the oscillatory
flame at TW = 974 K; negative values signify that heat is transferred from
the fluid to the wall. The circles indicate the five time instants for the
flames in Fig. 3.9.
At t1 , the maximum temperature was located within the weakening
54
3.3. Results and discussion
flame from the previous ignition-extinction cycle. Shortly after, the temperature peak was inside the newly-ignited flame and Tmax increased sharply
to reach the peak of the entire cycle. Regarding iHT R (Fig. 3.10(b)), the
maximum of the heat losses along the lower wall was attained due to the
propagation of the flame from the new kernel close to it, while on the upper
wall the extinction of the previously ignited flame rendered the heat losses
positive (i.e. heat was transferred from the wall into the channel).
Figure 3.11(a) shows the absolute propagation speed along the flame
fronts at times t3 and t4 defined as
Sa = Sd + u · n
(3.6)
with Sd the displacement speed, the front speed relative to the flow velocity
u [103],
Sd = −
1 DΘ
,
|∇Θ| Dt
(3.7)
and n the normal to the front vector, pointing towards the fresh mixture
n=
∇Θ
|∇Θ|
(3.8)
Here, Θ is the hydrogen mass fraction isovalue YH2 = 10−4 , defining the
flame front.
Both fronts were characterized by absolute propagation speeds of opposite sign at their ends: the upper end propagated with a positive velocity
towards the unburnt gases while the lower end regressed with a negative
55
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.11: (a) Local absolute propagation speed, (b) stretch, and (c)
local heat release rate along the isoline YH2 = 10−4 for the flame fronts of
Figs. 3.9(d) (solid lines) and 3.9(e) (dashed lines); TW = 974 K.
velocity towards the exhaust gases. In the central part of the flame, the
displacement speed almost balanced the flow velocity resulting in a nearly
zero value of the absolute propagation speed. As a result of this distribution
of propagation speed along the flame front, the flame length decreased from
56
3.3. Results and discussion
about 6.5 cm (Fig. 3.9(d)) to less than 4 cm (Fig. 3.9(e)).
The flame front propagated faster at the upper end for two reasons.
Firstly, the upper-end was located farther downstream, so that the fresh
mixture just upstream of the front had a higher temperature compared to
the lower-end due to preheating from the hot wall. Another reason was the
interplay of the two fuel components with the flame stretch defined as [103]
K=
1 DA
= a t + Sd κ
A Dt
(3.9)
Here, A is the flame area, κ the front curvature
κ=∇·n
(3.10)
and at the total aerodynamic strain rate of the front
at = (δij − ni nj )
∂ui
∂xj
(3.11)
The stretch for the flames in Figs. 3.9(d) and 3.9(e) was negative (positive)
at the upper (lower) end (see Fig. 3.11(b)). Regarding the Lewis numbers,
LeH2 = 0.291 while LeCO = 1.116. Had the fuel been only H2 , the burning
rate would have intensified at the lower end and weakened at the upper
end (the opposite would have happened if the fuel was only CO). Since
rCO/H2 = 7, the fuel was mainly CO and the effective Lewis number based
on the inflow composition and computed according to Law et al. [108] was
slightly larger than unity (see Table 3.1). Since most hydrogen was depleted
57
3. Flame dynamics in lean CO/H2 /air combustion
before reaching the upper flame tip (compare Figs. 3.12(a) and 3.12(c)), the
unburnt mixture at the upper end had a rCO/H2 higher than the one at the
inflow and hence the effective Lewis number was greater than the one based
on the inflow composition. The HRR along the two fronts was larger at the
negatively-stretched upper end compared to the positively-stretched lower
end (see Fig. 3.11(c)). At t3 , the local heat release rate HRR was stronger in
the neighborhood of both walls compared to the almost unstretched central
part (Fig. 3.11(c), solid line) presumably due to the effect of preheating
from the hot walls.
rCO/H2
0
1
5
7
10
LeH2
0.366
0.327
0.295
0.291
0.287
LeCO
1.274
1.135
1.116
1.101
Leef f
0.366
0.836
1.012
1.025
1.036
Tad [K]
1471
1542
1591
1594
1604
SL [cm/s]
36.1
22.7
15.2
13.8
12.7
Table 3.1: Lewis numbers of H2 and CO, effective Lewis number according to [108], adiabatic flame temperature and flame speed based on the
unburned mixture computed with PREMIX [58] for different rCO/H2 and
ϕ = 0.42.
For the flame at t3 , Figs. 3.12 and 3.13 show the mass fractions and
the non-dimensional mass reaction rates of H2 , H2 O, CO and CO2 . Water
was mainly produced at the lower flame tip (Figs. 3.12(b) and 3.13(b))
resulting in an uneven distribution in the exhaust gases. The consumption
rate of hydrogen was stronger at the lower flame tip (Fig. 3.13(a)), while
CO was consumed more at the upper end (Fig. 3.13(c)). Close to the upper
tip, there was a region where CO reacted with water to produce hydrogen
58
3.3. Results and discussion
Figure 3.12: Distribution of the mass fractions of (a) H2 , (b) H2 O, (c) CO
and (d) CO2 for the flame of Fig. 3.9(d); TW = 974 K.
(Fig. 3.13(a) and 3.13(b)).
During t3 ≤ t ≤ t4 , the upper wall was in extended contact with
fresh mixture (Fig. 3.9(d) and 3.9(e)), and radical buildup resulted in the
formation of a new ignition kernel (Fig. 3.9(f)). Just before this kernel
ignited, the iHRR attained its minimum value, while at the same time the
iHT R summed on both walls was about -0.6 (Fig. 3.10(b)), suggesting that
its influence on the flame weakening was considerable (being about half of
the iHRR in magnitude).
As TW increased (Fig. 3.6), the propagation mechanism of the asymmetric flames did not change significantly, however: (i) the ignition points
shifted closer to the inflow, (ii) the maximum length of a single flame decreased, (iii) for values above TW = 994 K the flame was no longer able
to propagate upstream due to the higher frequency of the ignition events,
59
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.13: Non-dimensional reaction rates in mass of (a) H2 , (b) H2 O,
(c) CO and (d) CO2 for the flame of Fig. 3.9(d); TW = 974 K. The white
vertical lines in the color scales of (a) and (b) mark the zero values.
resulting in downstream thermal expansion and upstream fuel depletion.
3.3.2.3
Random ignition spots: 1004 K ≤ TW ≤ 1008 K
Further increase of the wall temperature resulted in a combustion mode
with multiple ignition kernels formed at random time instants and locations on the lower and/or upper wall. Since the same dynamic behavior
was observed in domains of L =105 mm as well as 140 mm, the observed
flame structures were independent of the channel length. Periods of a seemingly steady weak combustion mode characterized by a low heat release
rate similar to the weak flame described in the next section were observed.
Chemical activity was constrained in narrow layers close to the walls where
radical concentrations varied in a non-monotonic fashion, eventually leading
to ignition kernel formation.
60
3.3. Results and discussion
Figure 3.14: Time history of (a) iHRR and (b) maximum temperature for
TW = 1008 K; (c), (d) and (e) YOH distributions at the time instants t1 , t2 ,
t3 marked in (a) and (b).
Figure 3.14 shows the time history of iHRR and maximum temperature as well as YOH distributions at three time instants for TW = 1008 K.
The focus is on three time instants: the last interval of the non-periodic
evolution, weak combustion and reignition. Just before reignition (t3 in
Fig. 3.14), the maximum attained temperature was Tmax = 1010 K, only
2 K above the wall temperature. After reignition, Tmax exceeded 1600 K
(Fig. 3.14(b)). The flame at t1 (Fig. 3.14(c)) evolved from two kernels, one
on the upper and another on the lower wall, while at t2 (Fig. 3.14(d)) one
61
3. Flame dynamics in lean CO/H2 /air combustion
ignition kernel on the upper wall and three on the lower wall can be seen.
Just before the first reignition (Fig. 3.14(e)) the layer of radicals close to
the wall attained a wavy shape, leading to the creation of new kernels (one
readily visible at x ≈ 110 mm). Furthermore, it could be observed that the
second reignition in Fig. 3.14(a) and 3.14(b) at t = 0.92 s was stronger than
the first one, indicating that the perturbations responsible for this behavior
were not damped with time.
3.3.2.4
Weak flames: 1010 K ≤ TW ≤ 1034 K
For TW ≥ 1010 K only steady flames were formed and two regimes with
very different combustion types could be identified: weakly burning flames
for 1010 K ≤ TW ≤ 1034 K and strongly-burning V-shaped flames at higher
wall temperature (the latter will be elaborated upon in Section 3.3.2.5).
TW
K
954
1010
1012
1034
1034
L
mm
105
105
105
105
210
Tmax
K
954.0
1012.4
1015.1
1065.2
1065.2
Tmax location
mm
wall
33.64 , 0.06
34.71 , 6.93
47.47 , 0.35
47.47 , 0.35
H2 conv.
%
0.26
63.69
64.47
78.71
93.33
CO conv.
%
0.04
18.07
18.69
30.59
41.87
iHRR
0.002
0.865
0.889
1.342
1.778
Table 3.2: Main characteristics of the weak flames: maximum temperature
within the domain (Tmax ), x and y coordinates of the maximum temperature
point (channel length is L =105 mm, height h =7 mm and position [0,0] is
the inlet corner on the lower wall, see fig. 3.1), conversion of H2 and CO
and non-dimensional integral heat release rate iHRR. The first row refers
to the non-ignited case at TW = 954 K for comparison.
Weak flames were characterized by incomplete fuel conversion and max62
3.3. Results and discussion
imum gas temperatures only a few degrees above TW (see Table 3.2), with
reactions and heat release occurring in narrow layers close to the walls.
These flames were different from the non-ignited case presented in Section
3.3.2.1 (see first row of Table 3.2) since the maximum temperature was
higher than TW while the fuel conversion and the iHRR were orders of
magnitude higher than those of the non-ignited case.
Figure 3.15: Weak combustion case for TW = 1034 K: (a) Nondimensional
mass flow rate of H2 (solid line) and CO (dashed line) through the channel
normalized by the corresponding inflow values; (b) reaction rate (RR) of
H2 (solid line) and CO (dashed line) averaged over the channel height; (c)
YOH distribution.
Figure 3.15 depicts a case of weak combustion at the highest TW =
63
3. Flame dynamics in lean CO/H2 /air combustion
1034 K where this mode was observed; fuel conversion, reaction rates and
production of radicals peaked close to the inlet and they decreased gradually with increasing streamwise distance. The simulation was repeated by
doubling the channel length to 210 mm, but the increase in fuel conversion
remained modest (Table 3.2), indicating that this combustion mode did not
depend on the channel length, while for the non-ignited case the domain
length affected conversion. The weak regime presented in this paper differed
from the one of previous studies ( [92] and [96]), where it was only observed
at very low inflow velocities (≤ 0.8 cm/s in [92] and ≤ 3.0 cm/s in [96]).
Figure 3.16: Distribution of YOH of V-shaped flames for (a) TW = 1074 K
and (b) TW = 1274 K.
3.3.2.5
V-shaped flames: 1044 K ≤ TW ≤ 1274 K
For sufficiently high TW , stable vigorous combustion and complete fuel conversion were attained with the non-dimensional iHRR constantly close to
3.7 (complete fuel conversion). V-shaped flames without tip opening were
64
3.3. Results and discussion
established within the range 1044 K ≤ TW ≤ 1274 K (Fig. 3.16). By increasing the wall temperature from 1044 K to 1274 K the maximum temperature
in the domain increased from 1403 K to 1522 K, values lower than the adiabatic flame temperature of 1594 K, while for dynamic flames this value was
usually exceeded (see Fig. 3.7(b)). An increase in TW enhanced the local
flame speed, which in turn led to a reduction in the flame area (Fig. 3.16).
3.3.3
Effect of the CO to H2 ratio
A limited parametric study was carried out by varying the CO to H2 molar
ratio in the range 1 ≤ rCO/H2 ≤ 20, together with the wall temperature. The
observed transitions and flame modes in the TW versus rCO/H2 parameter
space are summarized in Fig. 3.17. Because of the high computational cost,
the numerical study is not as detailed as for rCO/H2 = 7.
Figure 3.17: Flame stability map in the TW – rCO/H2 parameter space.
65
3. Flame dynamics in lean CO/H2 /air combustion
3.3.3.1
Observed dynamics
For rCO/H2 = 20, the range 984 K ≤ TW ≤ 1274 K was investigated. At
TW = 984 K and 994 K, the mixture did not ignite. The radical concentration along the wall increased with increasing streamwise distance,
but the heat release and the fuel conversion were negligible. Starting at
TW = 1004 K, weak combustion was observed where the radical concentrations peaked close to the inflow and then decreased downstream, although the maximum temperature did not exceed the wall temperature.
At TW = 1034 K, the maximum temperature Tmax = 1034.29 K barely exceeded the wall temperature, while for TW = 1044 K the maximum temperature was 1046.61 K. By further increasing TW , the difference between wall
and maximum domain temperature increased, resulting in higher fuel conversion and heat release rate. Transition from weak combustion to strongly
burning V-shaped flames was smooth and finally, at TW = 1174 K the first
strongly-burning V-shaped flame fully contained within the domain was established with a maximum temperature of 1461 K. No dynamic combustion
mode was found for rCO/H2 = 20.
For rCO/H2 = 10, five combustion regimes were observed: no ignition
for TW = 969 K, oscillatory ignition for 974 K ≤ TW ≤ 994 K with frequency increasing with TW , random ignition spots between TW = 1004 K
and 1006 K. At TW = 1007 K the weak combustion regime started and
transitioned to strongly burning V-shaped flames at TW = 1074 K.
For rCO/H2 = 6, the mixture did not ignite at TW = 960 K, and the
oscillatory ignition regime spanned over the wall temperature range 965 to
66
3.3. Results and discussion
1006 K. The oscillation frequency increased again with TW , exhibiting a
symmetric mode of flame propagation between 985 K and 1003 K. Random
ignition spots and chaotic behavior were found between 1007 K and 1010 K,
while for TW = 1015 K and 1020 K weak combustion was attained.
For rCO/H2 ≤ 5 the weak combustion mode vanished. The stronglyburning V-shaped flame was observed for rCO/H2 = 5 at TW = 1025 K,
while at TW = 1020 K random ignition spots were formed; the oscillatory
ignition regime appeared for TW ≤ 1010 K.
Finally, the random ignition spots disappeared for rCO/H2 < 4.
3.3.3.2
Flame tip opening
The effect of rCO/H2 on the flame structure at the tip of the V-shaped
flame was investigated at TW = 1274 K for CO to H2 ratio between 0 (pure
hydrogen) and 10. The heat release rate distributions in Fig. 3.18 show that
the flame length increased with increasing rCO/H2 ratio due to the decreased
laminar flame speed (Table 3.1). The combined effect of the subunity Lewis
number (Table 3.1) and the negative flame curvature resulted in reactant
defocusing prevailing over the focusing of heat. As a consequence, the
burning rate at the tip weakened and the flame opened (Fig. 3.18(c) and
3.18(d)) for rCO/H2 = 1 and 0.
As already discussed, the effective Lewis number increased at higher
rCO/H2 ratios and the intensified heat release at the tip resulted in a closed
V-shaped flame (Fig. 3.18(a) and (b)) for rCO/H2 = 10 and 5. It should be
67
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.18: Nondimensional heat release rate (HRR) profiles for V-shaped
flames at TW = 1274 K and (a) rCO/H2 = 10, (b) rCO/H2 = 5, (c) rCO/H2 = 1
and (d) rCO/H2 = 0 (pure hydrogen).
noted, however, that in all cases no fuel leaked through the flame tip.
3.3.4
CSP analysis of the V-shaped and weak flames
The CSP analysis was applied to the flames with rCO/H2 = 7 presented in the
previous sections. In particular, the explosive modes (those characterized by
a positive real part of the eigenvalue and therefore by an exponential growth
in time) of V-shaped and weak flames were studied and compared in order
to identify the physics responsible for the absence of strong combustion in
the weak flames.
68
3.3. Results and discussion
3.3.4.1
V-shaped flame at TW = 1074 K and rCO/H2 = 7
The explosive mode analysis was first tested on the steady V-shaped vigorously burning flame at TW = 1074 K and rCO/H2 = 7. For a premixed
hydrogen/air flame, Lu et al. [72] found that an explosive mode exists only
in the preheat zone, involving mainly the thermal runaway, while radicals
(the most important being H) played only a secondary role. Radical explosion became important for auto-ignition.
Figure 3.19: Distribution of the log10 of the real part of the explosive mode
eigenvalue for (a) a V-shaped flame at TW = 1074 K and (b) a weak flame
at TW = 1034 K.
Figure 3.19(a) shows the explosive region for the flame of Fig. 3.16(a),
while Table 3.3 summarizes the results of the CSP analysis. Two regions
with different behavior could be distinguished at the flame anchoring points
(point A) and the preheat zone (point B) in Fig. 3.19(a). The peak in the
positive real part of the explosive mode was located at the flame anchoring
region where radical formation played a dominant role. CSP indicated that
at the flame anchoring points the radicals H, O and OH had major contri-
69
3. Flame dynamics in lean CO/H2 /air combustion
butions to the explosive behavior, while the temperature pointer was one
order of magnitude smaller. The preheat zone was dominated by thermal
runaway, with temperature being the most important variable followed by
O, H and OH.
TW [K]
1074
Point
A
λexp
r
178.67
x, y [mm]
4.2, 0
1074
B
7.17
72.8, 3.5
Q
0.373 H
0.345 O
0.149 OH
0.083 T
0.642 T
0.084 O
0.040 H
0.035 OH
P
-0.656 LT
-0.113 LH
0.033 R1
0.023 R2
-0.270 CT
-0.087 R9
0.124 R25
0.124 LT
Table 3.3: Real part of the eigenvalue of the most explosive mode in nondimensional units (λexp
r ), together with the CSP pointers (Q) and participation indices (P ) for selected points of the flame shown in Fig. 3.19(a). L, C
and R stand respectively for the diffusive, convective and reactive part in
the right hand side of the governing equations for temperature and species
(Eq. (2.5) and Eq. (2.7). The subscripts of the diffusive and convective
terms refer either to a species or to temperature while for the reactive term
they refer to the reaction number according to [104]: R1 : H + O2 = O +
OH, R2 : O + H2 = H + OH, R9 : H + O2 (+M) = HO2 (+M) and R25 :
CO + OH = CO2 + H.
The two zones displayed very different characteristics in terms of the
participation indices. At the anchoring points, the diffusive processes, especially of temperature (heat losses) and H radical, had a negative impact
since they opposed the radical build up process, while chain branching reactions involving H, O and OH radicals were the major positive contributors.
In the preheat zone, the advective term of the energy equation together
with reaction R9 (which consumed the very reactive H radical in favor of
70
3.3. Results and discussion
the more stable HO2 ) had a negative contribution, while diffusion of heat
and the strongly exothermic reaction R25 (see Table 3.3) had a positive contribution. The contributions identified with the participation indices are in
agreement with those reported in [71].
At the preheat zone boundary, a second explosive mode was found on
the burned side, which was slower than the one dominating the preheat region; however at the boundary it grew and formed a complex conjugate pair
with the first explosive mode. Within a narrow region, the real part of the
complex conjugate eigenvalues decreased and became negative, suppressing
the explosive behavior (as also observed in [71]).
3.3.4.2
Analysis of the weak flame at TW = 1034 K and rCO/H2 = 7
The explosive region of the weak flame at TW = 1034 K (Fig. 3.15) is
reported in Fig. 3.19(b), while the main results of the CSP analysis are
summarized in Table 3.4. At point A, the same characteristics as for the
anchoring point of the V-shaped flame were found. The positive real part
of the eigenvalues peaked in this region, which was dominated by radical
explosion with chain branching reactions having a positive contribution, and
diffusive processes of temperature and H radical having a negative one. The
behavior was more sensitive to the heat conduction and less sensitive to the
H radical diffusion compared to that of the V-shaped flame (compare Tables
3.3 and 3.4). The chain branching reaction R11 also appears in Table 3.4,
although its positive contribution was quite small. This reaction was also
identified in [71] as having a positive contribution to the explosive mode.
71
3. Flame dynamics in lean CO/H2 /air combustion
TW (K)
1034
Point
A
λexp
r
184.52
x, y [mm]
4.2, 0
Q
0.429 H
0.394 O
0.164 OH
0.009 T
1034
B
2.86
47.6, 0
0.686 T
0.033 O
0.014 OH
0.011 H
P
-0.79 LT
-0.05 LH
0.036 R1
0.019 R2
0.012 R11
-0.015 LH2
-0.005 R9
0.956 LT
0.004 R25
Table 3.4: Real part of the eigenvalue of the most explosive mode in nondimensional units (λexp
r ), together with the CSP pointers (Q) and participation indices (P ) for selected points of the flame shown in Fig. 3.19(b). L, C
and R stand respectively for the diffusive, convective and reactive part of
the right hand side of the governing equations for temperature and species
(Eq. (2.5) and Eq. (2.7). The subscripts of the diffusive and convective
terms refer either to a species or to temperature while for the reactive term
they refer to the reaction number according to [104]: R1 : H + O2 = O +
OH, R2 : O + H2 = H + OH, R11 : HO2 + H = OH + OH, R9 : H + O2
(+M) = HO2 (+M) and R25 : CO + OH = CO2 + H.
At point B (Fig. 3.19(b)), located farther downstream on the wall, temperature played the most important role. Participation indices provided a
possible explanation for the weak reactivity of this combustion mode. The
main positive contribution to the explosive mode was due to heat transfer
from the wall (Table 3.4), while the heat-releasing reaction R25 had a participation index two orders of magnitude smaller than in the preheat zone of
the V-shaped flame (Tables 3.3 and 3.4). The diffusive term in the H2 equation had the main negative contribution. The picture emerging from the
participation indices suggests that the explosive mode could no longer be
sustained by the exothermicity of R25 , and reactions could only take place
close to the hot walls supplying heat to the gas by conduction. Diffusion of
72
3.3. Results and discussion
H2 opposed the explosive mode, and this was attributed to the competition
for the OH radical between reactions R3 : H2 + OH = H2 O + H, and R25 :
CO + OH = CO2 + H, the most important reactions for the consumption
of H2 and CO (see Figs. 3.20(a) and 3.20(b)).
Figure 3.20: Reaction rate (RR) of the two most important reactions for
the depletion of (a) H2 and (b) CO, averaged over the channel height. (c)
Product of the kinetic constants of reactions R3 (solid line) and R25 (dashed
line) with the inflow concentrations of H2 and CO respectively, as a function
of temperature. The chemical reactions are numbered according to [104]:
R2 : O + H2 = H + OH, R3 : H2 + OH = H2 O + H, R22 : CO + O (+ M)
= CO2 (+ M), R25 : CO + OH = CO2 + H.
73
3. Flame dynamics in lean CO/H2 /air combustion
The product of the R3 and R25 kinetic constants with the respective
inflow concentrations of H2 and CO plotted in Fig. 3.20(c) as a function
of temperature indicated that, close to the wall temperature TW where
weak flames were observed, the two curves crossed, and the rate of the heat
releasing reaction R25 was overtaken by the chain propagating reaction consuming hydrogen. This could also explain why at higher H2 concentration,
i.e. rCO/H2 ≤ 5, the weak combustion mode was not observed since the
crossing point appears at a much lower temperature and H2 chemistry always dominates within the investigated TW range. On the other side, for
higher CO concentrations the intersection shifts to higher temperatures and
CO chemistry dominates, extending the range where weak combustion can
be observed (Fig. 3.17). In the neighborhood of the crossing points, chaotic
ignition spots were observed, which were suppressed either at high H2 or
high CO concentrations.
For the weak flame no pair of complex conjugate eigenvalues were found
at the boundary of the explosive region, suggesting that this kind of behavior
was typical of ignited cases like the V-shaped flame where combustion was
self sustained.
3.3.5
Experimental results with fast OH-LIF
A recent experimental campaign investigated lean premixed syngas/air combustion in a 7 mm height channel by adopting the fast OH-LIF technique.
Parametric studies were performed by varying the wall temperature, the
total equivalence ratio, and the CO:H2 ratio. The main objective was to
74
3.3. Results and discussion
assess the capability of the numerical code in reproducing the different combustion modes, especially the dynamic ones which could not be compared
against the OH* chemiluminescence images for the reasons explained in
Section 3.3.2. Experimental flames were compared against 2-D transient
simulations. The syngas flame modes were also mapped over the studied
parameter range and the impact of CO:H2 ratio on flame stability was investigated. The numerical setup is the one described in Section 3.2 while the
experimental setup will be described in Sections 3.3.5.1 and 3.3.5.2. Results
follow in Section 3.3.5.3.
3.3.5.1
Channel reactor
The test rig and OH-LIF setup are depicted in Fig. 3.21. The rectangular
reactor comprised two 9-mm-thick horizontal Si[SiC] ceramic plates (300
mm long (-x), 104 mm wide (-z) and placed 7 mm apart (-y)), and two
3-mm-thick vertical quartz windows. Wall temperatures were measured by
S-type thermocouples, embedded 0.9 mm beneath the inner Si[SiC] channel
surfaces by eroding 8.1 mm deep and 1 mm diameter holes from the outer
Si[SiC] surfaces. Twenty-two thermocouples per plate were used, arranged
in five axial rows, z = 0, ± 22.5 and ± 45 mm, with z = 0 denoting the x-y
symmetry plane (Fig. 3.21(a)). To control the wall temperature, a combined
heating/cooling arrangement was employed. The front vertical Si[SiC] faces
were in contact to a water-cooled metal frame (see Fig. 3.21(b)) so as to
moderate the reactor entry temperatures and mitigate flame flashback of
the highly reactive hydrogen-containing mixtures. On the other hand, two
75
3. Flame dynamics in lean CO/H2 /air combustion
(a)
00000000000
00 00 00 00 00 00 00 00 00 00 00
00000000000
x y
z
0000000000
00000000000000000000000000000000000000000000
00 00 00 00 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
00000000000000000000000000000000000000000000
00 00 00 00 00 00 00 00 00 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
0 0 0 0 0 0 0 0 0 0 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0y
00 00 00 00 00 00 00 00 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
7
0 0 0 0 0 0 0 0 0 0 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
00 00 00 00 00 00 00 00 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0x
00
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0000000000000000000000000000000000000000000000000000
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 0 0 0 0 0 0 0 0 0 0 0
z
00000000000000000000000000000000000000000
(b)
O2,N2 TCB TCA 0 0 0 Thermocouples(TC)
0 0 0 00 00 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
00 00 00 00 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0Heater
Static
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0coils
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00
000
mixers 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
000
y
7
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 00 00 00 00
00 00 00
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Grids 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0300
00000000000000000000000000000000
x
00 00 00 00 00 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
H2, CO Cooling
water
Figure 3.21: (a) Test rig with OH-LIF setup, (b) reactor detail. All distances in mm.
resistive heater coils were placed above the Si[SiC] plates over the length
100 < x < 300 mm (with x = 0 denoting the start of the channel reactor) to
counteract external heat losses. By adjusting the coil power, ramped axial
wall temperature profiles with varying magnitudes could be achieved.
The reactor was mounted inside a cylindrical steel-tank that controlled
the pressure; the present studies were conducted at atmospheric pressure.
Two 350-mm-long and 35-mm-thick quartz windows at both sides of the
76
3.3. Results and discussion
tank allowed for optical access in the reactor (Fig. 3.21(a)), while two additional quartz windows at the rear tank flange and the reactor exhaust, respectively, facilitated streamwise optical access for the LIF excitation beam.
H2 , CO, O2 and N2 were supplied by pressurized gas bottles and their
flows were controlled by four dedicated Brooks mass-flow meters; the O2 /N2
composition was that of air. The gases flowed through two sequential static
mixers and then through a rectangular flow-straightening section that had
the same cross section as the reactor (104 x 7 mm2 ) and was equipped with
cross-flow grids to provide uniform flow at the reactor entry (Fig. 3.21(b)).
The inlet temperature at the reactor entry was monitored by a sheathed
S-type thermocouple (TCA in Fig. 3.21(b)), the gas temperature after the
static mixers by another sheathed K-thermocouple (TCB ) and, finally, the
upper wall temperatures along the flow-straightening section by three Ktype thermocouples located at x = -42, -62, and –94 mm.
3.3.5.2
Laser diagnostics
The 526.5 nm second harmonic of a high-repetition rate pulsed Nd:YLF
laser (Quantronix Darwin-Duo 527-80-M, pulse duration 120 ns, repetition
rate up to 10 kHz) pumped a Radiant-Dyes Narrow Scan High-Rep laser
using Rhodamine 6G dye. The output radiation was frequency-doubled to
generate the 285 nm excitation beam, which had a low pulse energy (∼0.5
mJ) to avoid saturation of the A(ν 0 = 1) ← X (ν 00 = 0) OH transition.
Repetition rates of 750 Hz and 1000 Hz were used, sufficient to temporally
resolve the observed flame dynamics. A diverging cylindrical lens (f = -100
77
3. Flame dynamics in lean CO/H2 /air combustion
mm) and a converging spherical lens (f = 600 mm) transformed the 285
nm beam into vertical light sheet that propagated counterflow along the
z = 0 symmetry plane (Fig. 3.21(a)). Counterflow propagation offset from
the symmetry plane, in the range -30 < z < 30 mm, was also possible.
A fast CMOS camera (LaVision HighSpeedStar6 with 1024 x 1024 pixels,
1024 x 512 used) with an image intensifier (HighSpeedIRO) recorded the
OH fluorescence from both (1-1) and (0-0) transitions (308 and 314 nm,
respectively) at 90◦ through the reactor and tank side-windows. Channel
areas of 100 x 7 mm2 were recorded onto a 1024 x 74 pixel area of the CMOS
chip. Over the 100 mm recorded length, the light sheet thickness varied from
0.3 to 0.4 mm. The camera was traversed axially to the particular position
of each investigated flame, while at each position a movie consisting of 400
to 2000 images was recorded.
3.3.5.3
Results and discussion
Experiments were conducted at three lean equivalence ratios (ϕ = 0.25, 0.30
and 0.42, based on both fuels) and CO:H2 volumetric ratios from 1:1 to 20:1.
Experiments at the highest ϕ = 0.42 were not possible for CO:H2 < 3:1, due
to flame flashback, while CO:H2 = 1:1 was achievable only for the lowest ϕ
= 0.25. Wall temperatures, controlled by adjusting the heating coil powers,
ranged from 550 to 1320 K and inlet velocities changed modestly from 1.54
to 1.96 m/s. Surface temperatures measured by the central thermocouples
(z = 0) differed by as much as 16 K from the corresponding off-center values
(z = ± 22.5, ± 45 mm). All observed combustion modes are summarized in
78
3.3. Results and discussion
CO:H2 volumetric ratio
20
10
8
6
4
2
os cillatory
mode s
s ta tio na ry
mode s
1
600
700
800
900
1000
1100
1200
Average wall temperature [K]
Figure 3.22: Experimentally-deduced stability maps. Circles: stationary
modes, triangles: oscillatory modes. Filled symbols: modes with nearly
symmetric shapes, open symbols: asymmetric modes. Red symbols: ϕ =
0.42, blue symbols: ϕ = 0.30, green symbols: ϕ = 0.25. Dashed lines
demarcate regions of stationary (right) and oscillatory (left) modes, for
each investigated ϕ.
Fig. 3.22, parameterized in terms of the equivalence ratio, CO:H2 ratio, and
wall temperature, where the latter was the average value of the upper-wall
and lower-wall thermocouples along z = 0. Obviously, the non-uniform wall
temperatures in each experimental case did not always render the mean
temperature as a meaningful parameter for the discussion of the ensuing
dynamics. Nonetheless, Fig. 3.22 identified stationary combustion modes
at the higher temperatures and higher CO:H2 ratios and oscillatory modes
at lower temperatures and lower CO:H2 ratios. Detailed description of the
observed modes is provided next.
79
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.23: Stationary V-shaped flames: (a) ϕ = 0.42, CO:H2 = 3:1, UIN
= 1.71 m/s, TIN = 335 K, (b) ϕ = 0.30, CO:H2 = 7:1, UIN = 1.70 m/s,
TIN = 324 K. (a1, b1) measured (symbols) and fitted (lines) streamwise
wall temperature profiles, (a2, b2) OH-LIF images, (a3, b3) simulations
(color bars: predicted OH in ppm-mass). Arrows marked xig denote ignition
positions.
Steady V-shaped flames
At the higher wall temperatures and higher CO:H2 ratios, stationary Vshaped flames were formed (filled circles in Fig. 3.22). Figure 3.23 provides
LIF-measured and predicted OH maps, as well as streamwise profiles of
upper-wall and lower-wall temperatures at z = 0, for two selected cases.
The flames were modestly asymmetric, due to small temperature differences
between the walls (typically by as much as 25 K), but they always formed
two connected branches extending to the walls. The agreement between
80
3.3. Results and discussion
measurements and predictions in Fig. 3.23 was particularly good in terms of
flame shapes and ignition positions (the latter, indicated by arrows marked
xig , were defined as the far-upstream positions whereby OH rose to 5% of
its maximum value in the channel).
Ignition for both cases in Fig. 3.23 occurred at the slightly hotter
lower wall, close to the entry in Fig. 3.23(a3) and at the reactor rear in
Fig. 3.23(b3) (xig ≈ 4 and 174 mm, respectively).
Steady asymmetric flames
Stationary asymmetric flames (open circles in Fig. 3.22) are reported in
Fig. 3.24. OH-LIF images and predictions in Fig. 3.24(a2, a3) indicated
formation of a “lower asymmetric flame”. Asymmetric flames exhibited a
single flame branch extending from the upper to the lower wall and should
not be confused to the V-shaped flames with modest asymmetry shown in
Fig. 3.23. To clarify this issue, additional simulations were performed with
a steady elliptic code (code details have been provided elsewhere [25, 61]),
which can also find the unstable solutions.
In Fig. 3.24(a4) a V-shaped flame and in Fig. 3.24(a5) an “upper asymmetric flame” were obtained with the steady code (in addition to the lower
asymmetric flame, which is not shown), depending on the starting guess.
We have shown in earlier works [91,92] that steady asymmetric flames stem
from bifurcations of V-shaped steady flames, such that Fig. 3.24(a4) was
indeed an unstable solution. While both lower and upper asymmetric flame
solutions have been numerically shown to occur over the same parameter
81
3. Flame dynamics in lean CO/H2 /air combustion
Figure 3.24: Stationary asymmetric flames: (a) ϕ = 0.42, CO:H2 = 3:1, UIN
= 1.69 m/s, TIN = 334 K, (b) ϕ = 0.30, CO:H2 = 15:1, UIN = 1.70m/s,
TIN = 340 K. (a1, b1) measured (symbols) and fitted (lines) streamwise
wall temperature profiles, (a2, b2) OH-LIF images, (a3, b3) simulations
(color bars: predicted OH in ppm-mass). Images (a4, a5) provide additional
solutions with a steady solver. Arrows marked xig denote ignition positions.
space in channels with symmetric wall boundary conditions (same wall temperatures) [91, 92], it is not clear whether the upper asymmetric solution
in Fig. 3.24(a5) is realizable. The unequal wall temperatures in the experiment may render the upper asymmetric solution unstable, or may require
significant external perturbations to transit from the lower to the upper
82
3.3. Results and discussion
asymmetric solution. Stability of the obtained solutions could only be assessed with the transient code; however, when using as initial condition for
the transient simulations the upper-asymmetric solution from the steady
code in Fig. 3.24(a5), the flame reverted to the lower-asymmetric solution
in Fig. 3.24(a3). Hence it appears that the only stable solution is the lower
asymmetric in Fig. 3.24(a2, a3).
Figure 3.24(a2, a3) (lower asymmetric flame) and Fig. 3.24(b2, b3) (upper asymmetric flame) indicate again very good agreement between OH-LIF
measurements and predictions. Ignition positions (xig ) were well-captured
by the model both for the upstream-igniting flame in Fig. 3.24(a2, a3) and
the downstream-igniting flame in Fig. 3.24(b2, b3).
Oscillatory flames
Theoretical analysis based on a simpler thermodiffusive model has demonstrated that a decrease in wall temperature promoted oscillatory combustion
modes in channels [94]. In the present experiment, various types of oscillatory modes were also observed at the lower range of wall temperatures
and predominantly at the lower CO:H2 ratios (triangles in Fig. 3.22). The
measured oscillatory flame dynamics were the same, irrespective of the LIF
light-sheet lateral position (-30 < z < 30 mm). Finally, most oscillatory
modes involved asymmetric flame propagation (open triangles in Fig. 3.22)
while nearly symmetric flame propagation was attained in very few cases
(solid triangles in Fig. 3.22).
Figure 3.25(a1-a4) depicts OH-LIF images at four times during the os83
3. Flame dynamics in lean CO/H2 /air combustion
(a1)
(b1)
(a2)
(b2)
(a3)
(b3)
(a4)
(b4)
50
70
90 110
x [mm]
130 150 20
0 2000
40
60
80
x [mm]
100 120
Figure 3.25: Oscillatory flame, ϕ =0.42, CO:H2 = 7:1, UIN = 1.66 m/s, TIN
= 338 K. (a1-a4) OH-LIF images at times 0, 28, 46, and 78 ms, respectively,
(b1-b4) predicted OH (color bar: OH in ppm-mass) for constant TW = 969
K and same inlet properties, at times 0, 16, 26, and 41 ms, respectively.
Figure 3.26: Measured wall temperature profiles for the oscillatory flames
of: (a) Fig. 3.25, (b) Fig. 3.28.
cillatory period of a flame with ϕ = 0.42, CO:H2 = 7:1, UIN = 1.66 m/s
and TIN = 338 K. The measured wall temperature profiles are provided in
Fig. 3.26(a); during the flame oscillations, no change occurred in the wall
temperatures due to the large thermal inertia of the 9-mm-thick Si[SiC]
plates. While the numerical model accurately captured the foregoing stationary combustion modes, the oscillatory modes proved very sensitive to
the specific boundary conditions. This may not be surprising, as oscillatory
flames have been known to be sensitive to gas radiation losses [109] (not con84
3.3. Results and discussion
sidered in the model), to specific kinetic rates (discussed in CSTR syngas
simulations [81]) that might not be accurately reproduced by the reaction
mechanism, and possibly to 3-D effects. Moreover, the experiments could
be subject to small external perturbations, which could trigger the oscillatory modes and which could not be included in the model. Nonetheless,
qualitative agreement between measurements and simulations was obtained
at somewhat different operating conditions.
Simulations using the measured temperature profiles in Fig. 3.26(a)
did not yield an oscillatory flame but an asymmetric one; this was also the
case when using axial wall temperature profiles acquired by the offset-inz thermocouples (Fig. 3.21(a)). However, simulations at a constant wall
temperature TW = 969 K (all other conditions being the same) shown in
Fig. 3.25(b1-b4) compared very favorably against the experiments. Measurements and predictions in Fig. 3.25 pertained to slightly different domain
lengths (50-150 mm and 20-120 mm, respectively). The simulations correctly reproduced all flame structures and propagation stages in Fig. 3.25:
(a1-b1) ignition at the lower wall with the ignited flame of the previous
cycle (visible downstream) being blown out of the domain; (a2-b2) establishment of a lower asymmetric flame; (a3-b3) shrinkage of the asymmetric
flame towards a more symmetric structure; (a4-b4) ignition of a new flame
on the opposite, upper wall. The process repeated itself in a periodic way.
Times histories of integrated predicted OH mass fraction and the integrated
OH-LIF signal (Fig. 3.27) indicated periodic solutions. The calculated oscillation frequencies were 22 Hz and 10.5 Hz for the simulation and experiment,
85
3. Flame dynamics in lean CO/H2 /air combustion
OH (a.u.)
OHx10-3
respectively.
8
(a)
4
0
12
8 (b)
4
0
0.0
0.1 Tim e (s ) 0.3
0.4
Figure 3.27: Time series of integrated OH in the domains of Fig. 3.25: (a)
predictions, (b) OH-LIF.
Figure 3.28: Oscillatory flame, ϕ =0.42, CO:H2 = 7:1, UIN = 1.60 m/s, TIN
= 326 K. (a1-a4) OH-LIF images at times 0, 92, 210, and 318 ms, respectively, (b1-b4) predicted OH (color bar: ppm-mass) for constant TW = 964
K and same inlet properties, at times 0, 53, 116, and 169 ms, respectively.
Another oscillatory flame, exhibiting very broad axial flame movement,
is illustrated in Fig. 3.28. The experimental conditions were ϕ = 0.42,
CO:H2 = 7:1, UIN = 1.60 m/s and TIN = 326 K, while the measured wall
temperature profiles are shown in Fig. 3.26(b). Simulations of a flame with
similar behavior are shown in Fig. 3.28(b1-b4), obtained for a case with the
86
3.4. Conclusions
same inlet conditions but a constant wall temperature TW = 964 K. The
flame ignited at the end of the channel and either an upper-asymmetric
flame (Fig. 3.28(a1, b1) and Fig. 3.28(a2, b2)), or a lower-asymmetric flame
(Fig. 3.28(a3, b3) and Fig. 3.28(a4, b4)) propagated upstream. Reverting
between upper- and lower-asymmetric flame propagation was perfectly repeatable in both experiments and predictions. This mode resembled oscillatory ignition/extinction, as the downstream flame propagation could reach
the channel outlet in certain cases. Deduced oscillation frequencies from the
OH time histories were 8 Hz and 4 Hz for the simulation and experiment,
respectively.
3.4
Conclusions
Premixed combustion of CO/H2 /air mixtures in a 7 mm height channel
with chemically inert walls was studied numerically. The predictive capability of the computational tool was validated against OH* chemiluminescence
measurements under conditions at which steady flames were obtained. Experiments further detected dynamic phenomena for a mixture composition
characterized by ϕ = 0.42 and CO to H2 molar ratio rCO/H2 = 7, which was
adopted as the starting point for a parametric numerical study to map the
combustion modes as a function of TW and rCO/H2 .
For ϕ = 0.42 and rCO/H2 = 7, four combustion modes were identified with increasing wall temperature: (i) no ignition (ii) dynamic modes,
including oscillatory ignition with periodic behavior and random strongly
87
3. Flame dynamics in lean CO/H2 /air combustion
burning ignition spots alternating with periods of weak combustion, (iii)
weakly burning steady flames, and finally (iv) vigorously burning steady
V-shaped flames.
By increasing TW within the oscillatory ignition regime, the frequency
of the ignition events increased exponentially. Above TW = 1004 K another
dynamic mode appeared which has not been previously reported in the
literature. It is characterized by ignition spots appearing at apparently
random locations and time instants close to the walls and evolving into
strongly burning flames. Periods of intense burning alternated with periods
of weak combustion. As the wall temperature exceeded TW = 1008 K,
combustion became steady and two stationary modes were observed: (i)
weak combustion characterized by low chemical activity confined close to
the walls, maximum temperature barely exceeding the wall temperature and
incomplete fuel conversion for 1010 K ≤ TW ≤ 1034 K, and (ii) a V-shaped
premixed flame at higher wall temperatures.
The weak flame regime reported here differs from the one of previous
studies [92, 96] where it was only observed at very low inflow velocities and
was mainly driven by flame/wall interactions. In the present case it appears
to be mainly kinetically driven.
The Computational Singular Perturbation (CSP) technique was employed to gain insight into the physicochemical processes responsible for
the weak combustion mode. The rates of the two main reactions responsible for depleting the fuel, respectively H2 + OH = H2 O + H and CO +
OH = CO2 + H, crossed at a temperature close to the range of the weak
88
3.4. Conclusions
combustion mode. CSP results provided evidence that a plausible reason
for weak combustion was the competition for the OH radical between the
two reactions.
Variation of the rCO/H2 ratio within the range 1 ≤ rCO/H2 ≤ 20 showed
that at high CO concentrations the weak combustion mode range extended
over a wider TW range, while the dynamic modes disappeared. At high H2
concentrations the weak mode vanished and the oscillatory mode extended
over a wider TW range, but the chaotic mode with random ignition spots
was suppressed.
Finally, the fast OH-LIF technique was adopted in a recent experimental campaign whose main target was to assess the predictive capability of
the numerical code in reproducing dynamic modes. The numerical results
compared favorably to steady experimental flames under the same operating conditions, while dynamic modes were correctly predicted by imposing
a constant wall temperature profile instead of the measured one. Possible
reasons for the discrepancy might be some approximations used in the numerical code like the absence of radiative heat losses, the adequacy of the
detailed reaction mechanism in reproducing correctly specific kinetic rates
under the given conditions, and possibly 3-D effects, phenomena to which
the oscillatory modes have been shown to be very sensitive.
89
Chapter 4
Hydrogen catalytic combustion in
Pt-coated channels 1
4.1
Review on numerical models for catalytic
combustion
The development of advanced numerical catalytic combustion models poses
severe challenges due to the different physicochemical processes with disparate spatio-/temporal scales: surface and gas-phase chemistry, convective
and diffusive heat and mass transport within the bulk flow (interphase transport), intraphase transport within the catalyst washcoat, heat conduction in
the channel wall and radiative heat transfer. The corresponding time scales
range from fractions of milliseconds (surface and gaseous chemical times) to
1
The content of the present chapter is published in “A. Brambilla, C. E. Frouzakis, J.
Mantzaras, A. Tomboulides, S. Kerkemeier, K. Boulouchos. Detailed transient numerical
simulation of H2 /air hetero-/homogeneous combustion in platinum-coated channels with
conjugate heat transfer, Combust. Flame, doi:10.1016/j.combustflame.2014.04.003, in
press”.
91
4. Hydrogen catalytic combustion in Pt-coated channels
seconds (solid axial heat conduction time). Steady and transient models are
required for the description of an individual catalytic channel during steady
operation and light-off, respectively. Furthermore, a spatial dimensionality
of at least two is necessary to correctly describe interphase transport and
particularly homogeneous combustion. The reason is that gaseous combustion, which is in many cases detrimental to the catalyst integrity and hence
of interest for safety reasons, is largely controlled by the reactant species and
temperature boundary layer profiles [110, 111]. Homogeneous combustion
is often neglected in catalytic combustion modeling even though its influence, especially on steady-state combustion characteristics and at elevated
pressures, may be significant [16, 112, 113]. Flow laminarization induced
by the intense heating from the catalytic surfaces renders laminar modeling appropriate at inlet Reynolds numbers where an isothermal flow would
have otherwise transitioned to turbulence [114]. Three-dimensional effects
do play a role in commercial straight-channel honeycomb reactors [115,116]
despite the rounding of the channel corners during coating [117] that in turn
suppresses strong secondary flows.
One-dimensional channel flow models with lumped heat/mass interphase transport coefficients and simplified chemistry have been extensively
used in the past (e.g. [118–121]). The models are appealing for their simplicity, however unique lumped heat and mass transport coefficients (Nusselt
and Sherwood numbers) for developing channel flows relevant to practical
catalytic combustion systems are not warranted in the presence of hetero/homogeneous chemistry [120,122]. A further simplification in 1-D channel-
92
4.1. Review on numerical models for catalytic combustion
flow models is the classical plug-flow reactor model, which is valid either for
very low or very high Reynolds numbers [123].
Two-dimensional models reproduce realistically the fluid mechanical
transport and, when coupled with detailed hetero-/homogeneous chemistry,
provide a powerful tool for both reactor design and fundamental research.
For moderate to high Reynolds numbers, the boundary-layer models provide accurate results at an affordable computational cost [123]. Elliptic
(full Navier-Stokes) 2-D models with detailed hetero-/homogeneous chemistry [19, 60, 62, 105, 114, 123–130] have become common research tools for
steady simulations, but to the best of the author’s knowledge research
codes based on such detailed models have not been applied in transient
simulations. Transient catalytic channel simulations with detailed chemistry were reported in [131] using a commercial CFD software although this
study did not address topics like light-off, gas-phase ignition and hetero/homogeneous chemistry coupling.
Given the long characteristic solid wall heat-up time compared to the
chemical, convective and diffusive time scales of a catalytic reactor, the
quasisteady-state (QSS) assumption for the gas-phase and surface chemistry
is often invoked. It is nonetheless understood that the QSS assumption may
not hold at all spatial locations and/or all times during transient reactor
operation [112, 132]. Two-dimensional transient CFD models for a single
channel invoking the QSS approximation and using one-step surface chemistry were reported in [132–134]. Recently, a QSS transient 2-D channel
model including detailed surface and gas-phase kinetics, heat conduction
93
4. Hydrogen catalytic combustion in Pt-coated channels
in the solid wall and surface radiation heat transfer, was able to reproduce
ignition and extinction experiments in catalytic partial oxidation (CPO)
of methane over a rhodium-coated reactor [135]; the same code was further applied to study the light-off and hetero-/homogeneous combustion of
syngas mixtures over platinum [136] and of lean methane/air mixtures in Ptcoated microreactors [112]. The quasisteady-state approximation imposes
in many cases a 1-D discretization of the solid wall [112,135], otherwise the
transverse heat conduction in the solid wall would become faster than the
gas-phase equilibration process and thus invalidate the QSS approximation.
Fully transient 2-D models with inclusion of intraphase washcoat diffusion
but simplified catalytic chemistry were presented in [137]. Fundamental
studies have also been carried out using transient codes with detailed chemistry for 1-D stagnation flows (SF) [138,139], which are characterized by an
affordable computational cost, although the SF geometry is not particularly
relevant to practical systems.
The present study undertakes a fully transient (without invoking the
quasiteady-state assumption) 2-D numerical study of fuel-lean hydrogen/air
hetero-/homogeneous combustion in a planar mesoscale (2 mm height) Ptcoated channel, using detailed gas-phase and catalytic chemistry, along with
2-D heat conduction in the solid wall. Fast processes such as light-off, homogeneous ignition and hetero-/homogeneous combustion interactions during
dynamic operation, with time scales ranging from fractions of milliseconds
to a few milliseconds, were studied for the first time given the model’s capability to resolve all relevant time scales. The code was further used to assess
94
4.2. Computational domain, boundary conditions and chemistry
the regimes of validity for the QSS assumption by comparing the present
predictions with those of the QSS code used in [135, 136]. Implications of
such comparisons are important by providing guidelines for proper use of
the simpler QSS models.
The chapter is organized as follows. The computational domain, boundary conditions and chemical reaction mechanisms are presented in Section
4.2 while the results follow in Section 4.3. In Section 4.3.1, the Computational Singular Perturbation method is applied to a transient surface perfectly stirred reactor (SPSR) in order to discuss the stiffness of the surface
reaction mechanism. Section 4.3.2 presents the results of 2-D transient simulations. The physics of light-off, wall heat-up and gas-phase ignition are
discussed in Section 4.3.2.1 to 4.3.2.3, while the comparison between fully
transient and QSS predictions is elaborated in Section 4.3.2.4.
4.2
Computational domain, boundary conditions
and chemistry
The gas-phase computational domain was a 2-D channel of height h = 2 mm
and length L = 10 mm (see Fig. 4.1). Symmetry was not invoked in
order to allow for potential asymmetric flame structures, as reported in
previous studies [17, 18]. At the inflow, uniform velocity profile, constant
inlet temperature and flux boundary conditions for the species (Eq. (3.1))
were specified. The first 1 mm channel length was considered chemically
inert (non-catalytic) while the remaining Lc = 9 mm were coated with Pt.
95
4. Hydrogen catalytic combustion in Pt-coated channels
A similar configuration of the channel walls was adopted in [112] and [140]
to avoid numerical complications associated with the catalyst leading edge
coinciding with the inlet boundary. Along the inert part of the channel
walls, no slip for both velocity components and zero flux for all species
(i.e. chemically inert walls) were applied. Along the Pt-coated part of the
channel walls the normal component of the fluid velocity was set equal to
the Stefan velocity (Eq. (2.23)) while the tangential component was set
equal to zero. The flux boundary conditions of Eq. (2.22) were applied for
gas-phase species. At the outlet, zero Neumann boundary conditions were
imposed. The inlet pressure was in all cases atmospheric.
L = 10 mm
s
s cs
s= 25 m
Figure 4.1: Schematic of the computational domain.
The solid walls had a thickness of δs = 25 µm and the materials employed were either cordierite (ceramic) or FeCr-alloy (metallic), with their
respective properties summarized in Table 4.1 [141]. The outer horizontal
channel walls were treated as adiabatic because of the presence of nearby
identical channels in a honeycomb reactor. The vertical front solid wall
faces towards the inlet and the outlet were also adiabatic (see Fig 4.1). The
96
4.2. Computational domain, boundary conditions and chemistry
catalyst has been treated as an extremely thin layer deposited on the solid
wall with the catalytically active area equal to the geometric area. Thus,
internal mass transfer limitation has not been considered.
Material
cordierite
FeCr-alloy
λs
[W/mK]
2
16
ρs
[kg/m3 ]
2600
7200
cs
[J/kgK]
1464
615
ρs cs
[kJ/kgK]
3806
4428
Table 4.1: Properties for cordierite and FeCr-alloy materials.
The fluid domain was discretized with 50 equally-spaced elements in
the streamwise direction (x) and 20 elements in the transverse direction
(y), the latter refined towards the wall to resolve strong gradients at the
gas-solid catalytic interface. The solid walls were discretized such that at
the gas-solid interface the elements had the same size as the corresponding
elements on the fluid side. In total, 50 elements in the streamwise direction
and 2 elements in the transverse direction were used for each wall. Thus, the
element size in the x-direction was uniform and equal to 0.2 mm while in the
y-direction it was varying from 0.2 mm at the channel axis to 0.01 mm at
the wall/fluid interface. Ninth-order interpolating polynomials per spatial
direction and element were employed. Grid-independence of the solution
was checked by varying the polynomial order. Due to the high stiffness of
the surface reaction mechanism, as will be discussed in Section 4.3.1, the
integration time step was 4 ns; it was further reduced down to 1 ns during
fast transient events like light-off and gas-phase ignition.
The oxidation of hydrogen on Pt was described by the detailed (mean-
97
4. Hydrogen catalytic combustion in Pt-coated channels
field) scheme of Deutschmann et al. [126] with 14 elementary reactions,
6 gas-phase species and 5 surface species. The surface site density was
Γ =2.7x10−9 mol/cm2 . The elementary gas-phase scheme of Li et al. [142]
was employed, which comprised 21 reversible reactions and 8 gaseous species
with accompanying gas-phase thermodynamic data. The capacity of this
scheme to reproduce measured gas-phase ignition characteristics under both
fuel-lean and fuel-rich H2 /air hetero-/homogeneous combustion has been
demonstrated in [64, 143]. Thermodynamic properties and reaction rates
were evaluated using CHEMKIN [47, 57], while transport properties were
computed using CHEMKIN’s transport library [45].
In the following, non-dimensional quantities will be reported in terms
of inlet velocity and composition, a reference temperature of 650 K and a
reference length equal to the channel height.
4.3
Results and discussion
The first part of the present section is dedicated to an analysis of the hetero/homogeneous chemical reactions by using the Computational Singular Perturbation (CSP) method [67, 68]. To the best of the author’s knowledge
CSP analysis has never been applied before on catalytic chemistry. The
different modes and the corresponding time scales were identified with an
in-house transient surface perfectly stirred reactor (SPSR) code based on
CHEMKIN [47, 57], and the key issue of numerical stiffness is further discussed. The second and main part of this section is devoted to the results of
98
4.3. Results and discussion
the fully transient detailed numerical simulations of a 2-D catalytic channel.
The physicochemical processes governing the start-up of a catalytic channel
are discussed in detail and fully transient numerical results are compared
with those of a code relying on the quasi-steady state assumption.
4.3.1
CSP analysis of a transient SPSR
CSP analysis was initially performed on a transient SPSR. The operating
conditions were representative of those used for the detailed numerical simulations presented later, although a 2-D catalytic channel cannot obviously
be reduced to a 0-D SPSR. Lean hydrogen/air mixtures with equivalence
ratio ϕ = 0.3, inlet temperature (and initial temperature) TIN = 650 K and
atmospheric pressure were investigated. The catalytic surface-to-volume
ratio of the SPSR was A/V = 9 cm−1 , same as that of the 2-D channel
with h = 2 mm and L = 10 mm (1 mm inert, 9 mm Pt coated, see Fig.
4.1). The SPSR residence time was τ = 10−3 s, which was equal to the
flow residence time in the 2-D channel of length L = 10 mm based on an
inflow velocity UIN = 10 m/s. The initial SPSR gas-phase composition was
equal to the inflow composition and the surface was initially covered only
with O(s). The elementary gas-phase mechanism of Li et al. [142] and the
detailed surface mechanism of Deutschmann et al. [126] discussed in Section
4.2 were used.
Figure 4.2(a) shows the SPSR temperature as a function of time while
Fig. 4.2(b) shows the time history of the absolute value of the real part
of the eigenvalue associated to the fastest mode, which is indicative of the
99
4. Hydrogen catalytic combustion in Pt-coated channels
Figure 4.2: Predictions in a transient SPSR with ϕ = 0.3, TIN = 650 K,
atmospheric pressure, A/V = 9 cm−1 , τ = 10−3 s. Time history of (a) temperature and (b) absolute value of the real part of the eigenvalue associated
to the fastest mode.
stiffness of the reaction mechanism. Throughout the simulation the fastest
time scale was about 10−12 s, which is comparable to that reported in laminar n-heptane premixed flames [144]. The latter was obtained using the
reaction mechanism of Curran et al. [145] with 560 species and 2538 reactions, where the smallest active time scale attained its lowest value 10−12 s
in the reaction zone. Table 4.2 reports the results of the CSP analysis for
the four fastest modes with an SPSR temperature Tr = 1200 K, a temperature at which both gas-phase and surface chemistry were active. The table
clearly shows the presence of three fast time scales associated exclusively to
surface species, namely H(s), H2 O(s) and OH(s) according to CSP pointer
values Q. The main reactions associated with the three fastest modes according to their participation indices P are adsorption of hydrogen and the
reversible reaction between H(s) and O(s) to form OH(s) for the first mode,
100
4.3. Results and discussion
and adsorption/desorption of water for the second and third modes. Gasphase species appeared only in the fourth fastest mode, whose time scale
was separated by three orders of magnitude from that of the third mode.
The other two surface species of the reaction mechanism, i.e. Pt(s) and
O(s), became important in the sixth fastest mode (not reported in Table
4.2), whose eigenvalue had a real part of -8.3x105 . Throughout the entire
simulation, the surface species H(s) had a CSP pointer very close to unity
for the fastest mode, apart from the short time period between 1.1x10−5 s
and 1.8x10−5 s, which corresponded to the rapid fall of |λr | in Fig. 4.2(b)
and it was dominated by O(s) and H2 O(s).
λr [s−1 ]
-1.13x1012
-2.19x1011
-1.37x109
-2.19x106
Q
9.95x10−1 H(s)
4.59x10−3 OH(s)
9.98x10−1 H2 O(s)
1.50x10−3 OH(s)
9.89x10−1 OH(s)
4.61x10−3 H(s)
7.44x10−1 H
2.51x10−1 OH
P
-4.84x10−1 H2 +2Pt(s)=>2H(s)
4.60x10−1 H(s)+O(s)=OH(s)+Pt(s)
-5.00x10−1 H2 O(s)=>H2 O+Pt(s)
4.69x10−1 H2 O+Pt(s)=>H2 O(s)
-3.97x10−1 H2 O(s)=>H2 O+Pt(s)
3.73x10−1 H2 O+Pt(s)=>H2 O(s)
3.19x10−1 H2 +OH=H2 O+H
-2.15x10−1 OH(s)=>OH+Pt(s)
Table 4.2: CSP analysis performed in a transient SPSR with reactor temperature Tr = 1200 K (see Fig. 4.2). The operating conditions are: lean
H2 /air with ϕ = 0.3, TIN = 650 K, A/V = 9 cm−1 , τ = 10−3 s. The first
column is the real part of the eigenvalues of the four fastest modes, while
the second and third columns are the most important CSP pointers (Q)
and participation indices (P ) of each mode, respectively.
The initial surface coverage had strong impact on stiffness during the
initial time steps. Table 4.3 presents CSP results of the fastest time scale
in the SPSR at the second integration time step (t = 2x10−6 s) and for
101
4. Hydrogen catalytic combustion in Pt-coated channels
different initial surface coverages. The case affected by the highest stiffness
was the one with an initial surface coverage of 100 % O(s), while the fastest
time scale of the other two cases, 100 % H(s) and 100 % Pt(s), respectively,
was about two orders of magnitude slower. By comparing the participation
indices for the three different cases, it could be clearly noticed that the
fastest process, i.e. the one mainly responsible for the high stiffness of the
system, was hydrogen adsorption.
θinitial
O(s)
λr [s−1 ]
-1.19x1012
Q
1.00 H(s)
H(s)
-6.08x1010
1.00 O(s)
Pt(s)
-5.65x1010
1.00 O(s)
P
5.00x10 H2 +2Pt(s)=>2H(s)
-5.00x10−1 H(s)+O(s)=OH(s)+Pt(s)
-5.00x10−1 H(s)+O(s)=OH(s)+Pt(s)
2.80x10−1 O2 +2Pt(s)=>2O(s)
-5.00x10−1 H(s)+O(s)=OH(s)+Pt(s)
2.79x10−1 O2 +2Pt(s)=>2O(s)
−1
Table 4.3: CSP analysis of the fastest time scale in a transient SPSR at
the second integration time step (2x10−6 s) and for different initial surface
coverages. Operating conditions as in Table 4.2. The second column is the
real part of the eigenvalues of the four fastest modes, while the third and
fourth columns are the most important CSP pointers (Q) and participation
indices (P ) of each mode, respectively.
In order to study the influence of the heterogeneous reaction mechanism on the resulting stiffness, the foregoing SPSR problem was solved by
replacing the detailed surface mechanism of Deutschmann et al. [126] with
a one-step surface mechanism based on the early work of Schefer [146]. This
one-step reaction mechanism has been used recently in a three-dimensional
direct numerical simulation of turbulent channel flow catalytic combustion [147]. The one-step reaction mechanism reduced dramatically the stiffness of the problem and the eigenvalue associated to the fastest time scale
102
4.3. Results and discussion
had a real part of -1.26x106 s−1 , i.e. six orders of magnitude slower than
the fastest time scale obtained with the detailed mechanism.
4.3.2
Detailed Numerical Simulations
In the next sections results of fully transient 2-D detailed numerical simulations are presented and discussed. Simulations were performed under the
operating conditions summarized in Table 4.4. Light-off, solid wall heat-up
and gas-phase ignition were investigated, the underlying physicochemical
processes were identified and CSP analysis was used as diagnostic tool. The
validity of the QSS assumption under different operating conditions and its
ability to capture key phenomena in catalytic combustion was assessed by
comparing the DNS results against those of a QSS code [135, 136].
Case
1
2
3
4
5
6
7
Material
cordierite
cordierite
cordierite
FeCr-alloy
FeCr-alloy
FeCr-alloy
FeCr-alloy
ϕ (H2 /air)
0.3
0.3
0.3
0.3
0.3
0.3
0.3
UIN [m/s]
10
10
10
10
10
10
5
TIN [K]
450
550
650
450
550
650
650
Table 4.4: Numerical conditions.
4.3.2.1
Light-off
Transient simulations were started from an initial temperature equal to the
inlet temperature (see Table 4.4) for both the fluid and the solid domain.
103
4. Hydrogen catalytic combustion in Pt-coated channels
The catalytic surface was initially covered with 100% O(s), which is realistic
in practical systems when only air flows in the channel before the fuel is
injected. Simulations with the transient SPSR at an inlet air composition
resulted in a surface O(s) coverage ranging from 100 % to 99.8 % for inlet
temperatures between 450 K and 650 K, respectively. In order to facilitate
the comparison with the QSS code, at t = 0 s the air was instantaneously
replaced with a H2 /air mixture having equivalence ratio ϕ = 0.3 (the QSS
code could not simulate the process of hydrogen injection). The physics of
the light-off process will be analyzed in detail for one selected case (Case 1 in
Table 4.4) and the main differences with the other cases will be subsequently
outlined.
Figure 4.3 shows the integral over the catalytic surface (Scatalytic ) of the
nondimensional rate of progress (iROP ) of the main reactions involved in
light-off for Case 1 (cordierite, TIN = 450 K, UIN = 10 m/s), defined as:
Z
iROPi =
qi dS
(4.1)
Scatalytic
with qi being the rate of progress of the ith reaction. The reaction numbering in Fig. 4.3 is identified in the figure caption. The reaction mainly
responsible for the delay in light-off was oxygen desorption (R6 ), due to its
high activation energy, as it was initially necessary to create free Pt sites
where hydrogen could adsorb (R1 ). Once adsorbed, H(s) reacted rapidly
104
4.3. Results and discussion
Figure 4.3: Integral over the catalytic surface of the nondimensional rate of
progress of the main reactions involved in the light-off of Case 1 in Table 4.4
(cordierite, TIN = 450 K, initial surface coverage 100% O(s)). The reaction
numbers on the right refer to the following reactions of Deutschmann’s
surface mechanism [126]: R1 : H2 + 2Pt(s) => 2H(s), R6 : 2O(s) => O2 +
2Pt(s), R8 : H2 O + Pt(s) => H2 O(s), R9 : H2 O(s) => H2 O + Pt(s), R12 :
H(s) + O(s) = OH(s) + Pt(s), R13 : H(s) + OH(s) = H2 O(s) + Pt(s), R14 :
OH(s) + OH(s) = H2 O(s) + O(s). The y-axes are in logarithmic scale.
with O(s) to produce OH(s) and free Pt sites according to the reaction:
R12 : H(s) + O(s) = OH(s) + Pt(s).
This reaction sequence cannot be responsible for runaway because it de-
105
4. Hydrogen catalytic combustion in Pt-coated channels
pends on the slow oxygen desorption reaction for the creation of new Pt
sites and moreover hydrogen adsorption followed by R12 cannot increase
the number of free Pt sites: when two Pt sites are taken by dissociative
adsorption of hydrogen, the same amount of Pt sites is freed again after
completion of reaction R12 . Thermal feedback is also not responsible for
the runaway because the nondimensional integral surface heat release rate
(iSHRR),
Z
iSHRR = −
Scatalytic
X
Ng +Ns
s˙i hi dS
(4.2)
i=1
was almost negligible before light-off (see Fig. 4.4(a)) and, given the high
heat capacity of the wall, the entire light-off process was almost isothermal
(see Fig. 4.4(c)).
The mechanism responsible for the catalytic runaway was identified by
using CSP analysis. An arbitrary point on the catalytic channel surface was
considered and studied in detail before the light-off process had commenced.
The results were independent on the location of the chosen point since the
conditions before light-off were the same at any position on the catalytic
surface. An explosive mode was found, whose corresponding eigenvalue
had a positive real part of about 3x105 s−1 . The main species of this mode,
identified by the largest CSP pointers, were Pt(s) and OH(s), while the
main reactions with a positive participation index were:
R1 : H2 + 2Pt(s) => 2H(s)
106
4.3. Results and discussion
OH(s) O(s)
H2O(s)
Pt(s)
H(s)
Figure 4.4: Light-off for Case 1 in Table 4.4 (cordierite, TIN = 450 K, initial
surface coverage 100% O(s)): (a) Integral over the catalytic surface of the
nondimensional heat release rate due to gas-phase species (solid line) and
surface species (dashed line); (b) Stefan velocity, (c) surface temperature
and (d) surface coverages plotted at a reference point (x = 10 mm) on the
catalytic surface of the upper wall.
R9 : H2 O(s) => H2 O + Pt(s)
R14 : OH(s) + OH(s) = H2 O(s) + O(s)
The results of CSP analysis, supported by the rates of progress shown in
Fig. 4.3, can be interpreted as follows. The buildup of OH(s) due to the
combination of hydrogen adsorption and R12 , which conserved the number
of existing free Pt sites as previously said, was accompanied by a fast in-
107
4. Hydrogen catalytic combustion in Pt-coated channels
crease in the rate of reaction R14 . Production of H2 O(s) from R14 was the
key step leading to runaway since this species desorbed rapidly (R9 ) leaving
one new free Pt site. Water adsorption (R8 ) before light-off was negligible,
therefore the new free Pt sites could be occupied by hydrogen adsorption.
Reaction R13 ,
R13 : H(s) + OH(s) = H2 O(s) + Pt(s)
was initially unimportant given the low coverage of H(s), which was mainly
consumed through reaction R12 , but its contribution became considerable
just before the runaway due to the increased hydrogen adsorption. Given
the fast water desorption rate, reaction R13 together with R9 had the potential to create two new Pt sites for every H(s) converted, thus leading
to runaway. Oxygen desorption became negligible after light-off and the
species O(s) was consumed almost entirely through reaction R12 .
Figure 4.4 summarizes the light-off process for Case 1 of Table 4.4.
Time histories are plotted in Fig. 4.4 at a reference catalytic point located
at the channel outlet. Figure 4.4(a) shows the iSHRR which was split into
heat release rate due to gas-phase species (solid line) and due to surface
species (dashed line). Surface heat release in catalytic combustion is almost
entirely due to gas-phase species (conversion of hydrogen and oxygen into
water) but during fast transient events, e.g. light-off, the heat released by
surface species can be noticeable due to a fast rearrangement of surface
coverages (for our specific initial conditions iSHRR became negative). For
108
4.3. Results and discussion
the same reason, the Stefan velocity became relevant only during light-off as
can be seen in Fig. 4.4(b), where positive values denote a Stefan flux towards
the wall. The Stefan velocity was initially positive due to a net deposition
of H atoms on the surface, which were mainly converted into OH(s) and
H2 O(s) (see Fig. 4.4(d)), while at the end of the light-off process it became
slightly negative due to a net water desorption. Figure 4.4(c) shows a small
surface temperature increase (∼ 2.5 K) during light-off, which confirms the
aforementioned negligible impact of the thermal feedback owing to the high
heat capacity of the wall. Finally Fig. 4.4(d) shows the time evolution of
surface coverage: it is seen that the rise in OH(s) was accompanied by a
rise in Pt sites. The observed behavior of surface species agrees with the
previously proposed runaway mechanism: the initial buildup of OH(s) on
the surface was necessary to activate reactions R14 and R13 which, combined
with water desorption R9 , released new Pt sites and hence accelerated the
light-off process.
The transition from the kinetically controlled to the mass transport
limited regime can be seen in Fig. 4.5 and Fig. 4.6(b), which provide the
time evolution of the transverse hydrogen mass fraction profile and hydrogen conversion rate, respectively. Once the runaway started, conversion of
hydrogen increased very sharply in the kinetically controlled regime due to
the fast surface chemical kinetics (even at 450 K) and to the high hydrogen
concentration at the walls. At a certain time diffusion became limiting for
the catalytic reactions, leading to a drop in hydrogen conversion (Fig. 4.6)
and the mass transport limited regime was reached.
109
4. Hydrogen catalytic combustion in Pt-coated channels
Figure 4.5: Time evolution during light-off of the transverse hydrogen mass
fraction profiles for Case 1 in Table 4.4 (y = 0 and 2 mm are the channel
walls while y = 1 mm is the channel midplane). The plotted time instants
are: t = 0.32 ms solid line, t = 0.34 ms dashed line, t = 0.40 ms dashed
dotted line, t = 0.80 ms dashed double dotted line, t = 1.20 ms dotted line.
The profiles are plotted at the channel outlet x = 10 mm.
Figure 4.6 shows nondimensional hydrogen conversion rate (H2 CR) during light-off for the three cases with cordierite walls (Cases 1 to 3 in Table
4.4). Hydrogen nondimensional conversion has been computed in two ways:
difference between inlet and outlet mass flow rate (Fig. 4.6(a))
Z
H2 CR = −
Z
−
ρYH2 (u + V H2 ) · ndS
(4.3)
IN
ρYH2 (u + V H2 ) · ndS
OU T
wit n being the outward pointing normal vector, and integral over the cat-
110
4.3. Results and discussion
Figure 4.6: Computed hydrogen conversion rate during light-off for the three
cases with cordierite wall (Cases 1 to 3 in Table 4.4): (a) difference between
inlet and outlet hydrogen mass flow rate in nondimensional units and (b)
hydrogen conversion rate computed from the integral over the catalytic
surface of the nondimensional hydrogen reaction rate. Solid line refers to
Case 1 with TIN = 450 K, dashed line refers to Case 2 with TIN = 550 K
and dotted lines refer to the Case 3 with TIN = 650 K. For comparison,
the nondimensional inlet mass flow rate of hydrogen is 0.0126, 0.0103 and
0.0088 for Cases 1 to 3 in Table 4.4, respectively.
alytic surface of the hydrogen reaction rate (Fig. 4.6(b))
Z
H2 CR =
sH˙ 2 WH2 dS
(4.4)
Scatalytic
The two definitions of hydrogen conversion rate converge to the same value
only at the end of light-off, when a new quasi-steady state controlled by mass
transport limitation was achieved. Both methods give the same results in
terms of the light-off times reported in Table 4.5, which are identified by
111
4. Hydrogen catalytic combustion in Pt-coated channels
a sudden increase in the H2 conversion rate (see Fig. 4.6(a) and 4.6(b)),
although the peaks have different characteristics for the two methods due to
different transient effects. The peaks in Fig. 4.6(b), computed by using the
surface reaction rate, were achieved due to the aforementioned acceleration
of chemical kinetics and to the high hydrogen concentration at the catalytic
walls at earlier times, since initially there were no spatial gradients. The
subsequent drop in hydrogen conversion rate was due to mass transport
limitation. The peaks in Fig. 4.6(a), computed by using inlet and outlet
mass flow rates, were due to the sudden fall in the outlet mass flow rate
caused by the Stefan flux (see Fig. 4.4(b)). At its peak the Stefan velocity
was about 11 cm/s and when integrated over the total catalytic length of 18
mm yielded a mass flow rate towards the catalytic surface equal to about
10 % of the inlet/outlet mass flow rates. Therefore, the presence of Stefan
velocity had a considerable effect on the flow field, although its impact
was limited to a very short time during light-off. It is worth noting that
such effects of Stefan velocity have not been reported in practical channel
geometries, as most spatially multidimensional simulations invoked the QSS
assumption, where by definition the Stefan velocity is equal to zero.
Light-off was achieved at 0.33 ms for Case 1 at TIN = 450 K, while it
was almost instantaneous for the cases with TIN = 650 K and 550 K. When
the simulations were initialized with a surface coverage either 100% H(s) or
100% Pt(s), light-off was almost instantaneous (i.e. occurring during the
first integration time step of 4 ns) for any initial temperature in Table 4.4.
The reason was that the only chemical reaction that could inhibit light-off
112
4.3. Results and discussion
Case
1
2
3
4
5
6
7
tlight−of f,DN S
ms
0.33
0.03
0.00*
0.33
0.03
0.00*
0.00*
tgas−phase ign,DN S
ms
351.60
262.72
192.49
not computed**
not computed**
258.20
311.32
tgas−phase ign,QSS
ms
350÷375
250÷275
180÷200
475÷500
350÷375
260÷280
300÷320
Table 4.5: Ignition times for the conditions specified in Table 4.4. Both
DNS and QSS times are reported for gas-phase ignition. For the QSS code
a time interval is given due to the limitation in the time step size which was
20 ms for the case at TIN = 650 K and 25 ms for the other cases.
* Cases 3, 6 and 7 had light-off at the first integration time step of 4 ns.
** Cases 4 and 5 have not been integrated down to the time of gas-phase
ignition due to very high computational cost.
at the considered initial temperatures was oxygen desorption, as previously
discussed. After light-off, the hydrogen conversion rate was the same for
the three different temperatures because it was dictated by mass transport
limitations.
The cases with FeCr-alloy wall are not reported in Fig. 4.6 since they
do not differ from their respective cordierite cases with the same inlet temperature. Due to the high thermal capacity of the wall, thermal feedback
did not play an important role in light-off, therefore the wall material had
almost no effect.
Finally, Case 7 with the lowest inlet velocity UIN = 5 m/s experienced
the same light-off behavior as Case 6 with UIN = 10 m/s and same inlet
temperature (TIN = 650 K) but the mass transport limited regime was
characterized by a higher hydrogen conversion relative to the inlet mass,
113
4. Hydrogen catalytic combustion in Pt-coated channels
Case
1
2
3
4
5
6
7
mH2 ,converted
H2 conversion
surface
(t = 20 ms)
(t = 20 ms)
4.1x10−3
4.1x10−3
4.1x10−3
4.1x10−3
4.1x10−3
4.1x10−3
2.9x10−3
%
32.7
39.9
46.9
32.7
39.9
46.9
65.2
H2 conversion
gas-phase
(t > tgas−phase
ign. )
%
63.4 (t=400 ms)
87.1 (t=300 ms)
94.1 (t=200 ms)
not computed
not computed
98.4 (t=280 ms)
99.6 (t=315 ms)
H2 conversion
surface
(t > tgas−phase
ign. )
%
4.3 (t=400 ms)
4.1 (t=300 ms)
5.4 (t=200 ms)
not computed
not computed
1.5 (t=280 ms)
0.4 (t=315 ms)
Table 4.6: DNS-computed hydrogen conversion for the conditions specified
in Table 4.4. All mass flow rates have been nondimensionalized with UIN
= 10 m/s to make them comparable.
although the total amount of converted hydrogen was lower (see Table 4.6).
Columns 2 and 3 of Table 4.6 summarizes hydrogen conversion wellafter light-off (t = 20 ms) for all the cases. For the cases with UIN = 10
m/s, the total amount of hydrogen converted was the same although the
conversion increased with increasing inlet temperature due to a lower inlet
density and thus a lower inlet hydrogen mass flow rate. At later times,
between light-off and gas-phase ignition, the hydrogen conversion increased
only slightly due to the small influence of the increasing wall temperature
on the mass transport.
4.3.2.2
Wall heat-up
The time period between light-off and gas-phase ignition was characterized
by wall heat-up due to the surface heat release, while the heat released
by gas-phase chemistry was negligible. Stefan velocity and surface heat
114
4.3. Results and discussion
release rate due to adsorbed species were minimal given the absence of any
fast transient event: the changes occurring on the catalytic surface were
mainly driven by the slow heat-up of the solid wall. The temperature and
radical distributions established during this time period were responsible
for the later gas-phase ignition. The present section focuses mainly on
the comparison between cordierite and FeCr-alloy, since the wall material
played a key role during heat-up, unlike the light-off period: heat capacity
affected the rise of temperature, while thermal conductivity was responsible
for redistribution of thermal energy within the solid.
Figure 4.7: Time evolution of interfacial gas/solid temperature profiles between light-off and gas-phase ignition for two cases with UIN = 10 m/s and
TIN = 650 K: (a) cordierite (Case 3 in Table 4.4) and (b) FeCr-alloy (Case
6 in Table 4.4). Solid lines are DNS results and symbols are QSS results.
The dashed line represents the value of the adiabatic flame temperature.
115
4. Hydrogen catalytic combustion in Pt-coated channels
Figure 4.7 shows a comparison between cordierite (Fig. 4.7(a)) and
FeCr-alloy (Fig. 4.7(b)) for Case 3 and Case 6 with UIN = 10 m/s and
TIN = 650 K. The differences in time evolution of the gas/solid interfacial
temperature profiles were exclusively due to the different properties of the
wall materials. Under mass transport limitation, hydrogen conversion was
controlled by transport of the reactants at the catalytic surface, which was
only slightly influenced by wall temperature. Thus, the total amount of
converted hydrogen (see Table 4.6) and consequently the surface heat release
rate was the same for both cases. Wall heat-up was faster for cordierite
compared to FeCr-alloy because of the slightly lower heat capacity ρs cs (see
Table 4.1). The temperature peak was sharper for cordierite due to its
lower thermal conductivity and could thus promote gas-phase ignition by
creating a hot spot, as also reported in [112]. The peak temperature for
the cordierite case was superadiabatic at t = 180 ms (see Fig. 4.7(a)) while
for FeCr-alloy superadiabatic temperature was achieved only well after gasphase ignition. Superadiabatic temperature will be discussed in Section
4.3.2.3. For FeCr-alloy, heat conduction from the catalytic section to the
upstream inert section was considerable (see Fig. 4.7(b)) and it could create
favorable conditions for gas-phase ignition within the non-catalytic zone
where hydrogen was not catalytically depleted, as it will be discussed in the
next section.
The different properties of the two wall materials played a key role
when considering the mechanical integrity of a catalytic reactor that could
be compromised by the thermally-induced mechanical stresses. For both
116
4.3. Results and discussion
considered cases, the thermal spatial gradients in the solid were maximum
close to the beginning of the catalytic section, where the surface heat release
rate suddenly jumped from zero to the maximum value attained within
the channel. The maximum thermal spatial gradients were achieved just
prior to gas-phase ignition. Regarding the transverse thermal gradient, the
maximum value in the entire gas/solid domain was attained on the gas side
of the catalytic surface: this was due to the interfacial energy conservation
in Eq. (2.24), where the thermal conductivity of the gas was much lower
than that of the fluid. At t = 180 ms cordierite had a maximum transverse
thermal gradient equal to 5094 K/mm on the gas side and 835 K/mm on
the solid side, while for FeCr-alloy the maximum values obtained at t = 240
ms were 2823 K/mm and 125 K/mm, respectively. The maximum value of
the axial thermal gradient was 3763 K/mm for cordierite and 826 K/mm for
FeCr-alloy. As expected, due to its lower thermal conductivity, cordierite
favored larger solid thermal stresses compared to FeCr-alloy.
Figure 4.8 shows the time evolution of surface coverage during the
heat-up phase for Case 3. With increasing wall temperature, more free Pt
sites became available. Figure 4.9 shows gas-phase radicals just before gasphase ignition (t = 180 ms). The most abundant species were H2 O2 and
the radical HO2 , which were produced even at a relatively low temperature
and were precursors of the more active radicals H, O and OH. These two
radicals were produced only in the gas-phase since they did not participate
in the surface reaction mechanism. Also H and O radicals were produced
only in the gas-phase as their desorption was recombinative. Due to a
117
4. Hydrogen catalytic combustion in Pt-coated channels
O(s)
Pt(s)
OH(s)
H2O(s)
H(s)
Figure 4.8: Time evolution of surface coverage at two times between lightoff and gas-phase ignition for Case 3 in Table 4.4 (cordierite, UIN = 10 m/s
and TIN = 650 K). Solid lines: t = 20 ms, dashed lines: t = 180 ms.
net desorption, the peak in OH radical was maximum at the catalytic wall
unlike the other radicals whose peaks were away from the catalytic wall (see
Fig. 4.9).
Regarding the influence of the operating conditions, Case 7 with reduced inlet velocity (UIN = 5 m/s) converted a lower amount of hydrogen
as previously discussed (see Table 4.6), therefore the surface heat release
rate was smaller compared to the cases with UIN = 10 m/s, causing a slower
wall heat-up. Finally, a lower inlet temperature had the obvious effect of
slowing down the wall heat-up.
118
y [mm]
4.3. Results and discussion
2
1 H
0
2
1 O
0
2
1 OH
0
2
1 HO2
0
2
1 H 2O 2
0
2
0
1.75E-06
0.00E+00
1.81E-05
0.00E+00
2.87E-05
0.00E+00
1.26E-04
0.00E+00
3.23E-05
0.00E+00
6
4
x [mm]
8
10
Figure 4.9: Two-dimensional maps of gas-phase radical mass fractions for
Case 3 with cordierite wall, UIN = 10 m/s and TIN = 650 K just before
gas-phase ignition (t = 180 ms).
4.3.2.3
Gas-phase ignition
Gas-phase ignition in catalytic channels requires higher wall temperatures
compared to an inert channel with same operating conditions due to the
near wall catalytic depletion of the limiting reactant. Surface chemistry
has also a modest inhibiting effect by adsorbing H and O radicals [62]. On
the other side surface chemistry could promote gas-phase ignition in two
ways: through wall heat-up due to surface heat release and through OH
radical desorption. In order to identify the crucial processes responsible for
gas-phase ignition, CSP analysis was used.
CSP analysis was performed to study the explosive modes, i.e. those
characterized by a positive real part of the corresponding eigenvalue and
hence by an exponential increase in time. Explosive modes are typical of
119
y [mm]
4. Hydrogen catalytic combustion in Pt-coated channels
2
1
0
0
2
6
4
8
10
x [mm]
0.00E+00
1.70E+00
Figure 4.10: Two-dimensional maps of the log10 of the real part of the
eigenvalue of the explosive mode for the case in Fig. 4.9. The explosive
mode is defined as the mode whose eigenvalue has a positive real part.
λr [s−1 ]
4.68x10+1
Q
5.34x10−1 H
2.98x10−1 O
1.67x10−1 OH
P
3.42x10 OH diffusion
7.19x10−2 H2 +OH=H2 O+H
5.04x10−2 H+O2 =O+OH
-4.28x10−1 H diffusion
-3.14x10−2 O diffusion
-1.53x10−2 OH convection
−1
Table 4.7: Results of the CSP analysis for the conditions in Fig. 4.9. The
first column is the real part of the eigenvalue of the explosive mode, the
second column contains the CSP pointers (Q) and the third column contains
the participation indices (P ).
runaway processes like gas-phase ignition. Fig. 4.10 shows the explosive
region for the conditions in Fig. 4.9, which refers to Case 3 of Table 4.4
at a time instant just prior to gas-phase ignition. The explosive region extended, as expected, only over a zone with higher temperature and higher
radical concentrations. Table 4.7 reports the result of CSP analysis refer-
120
4.3. Results and discussion
ring to a point where the real part of the explosive eigenvalue was close
to the maximum value obtained within the explosive region. CSP results
clearly identified OH diffusion from the catalytic wall towards the explosive
region as the main process responsible for the gas-phase runaway, while gasphase chain branching reactions had a much smaller effect, at least at this
time instant. Later, during the runaway, chain branching reactions became
dominant. Processes inhibiting the runaway were those that removed the
most reactive radicals from the explosive region, i.e. diffusion of H and O
radicals, which were adsorbed by the catalytic walls, and the convective
transport of OH radical. The participation index of temperature diffusion
did not appear among the three largest positive participation indices, therefore its contribution to the explosive mode was considerably weaker than
the contribution of OH diffusion.
Figures 4.11 to 4.13 depict the gas-phase ignition for Case 3 in Table
4.4. Figure 4.11(a) shows the increase in time of iHRR caused by gasphase ignition and the subsequent decrease in time of the iSHRR due to
the competition for hydrogen consumption between the gas-phase and the
catalytic pathways; gas-phase combustion prevented hydrogen from reaching the catalytic surface. Figure 4.11(b) shows the axial profiles of the
gas-phase and catalytic hydrogen conversion rates; gas-phase conversion
has been integrated over the channel half-height to facilitate comparisons
with the surface conversion on one wall. In a narrow region close to the beginning of the catalytic section, surface conversion dominated. Gas-phase
ignition could not be established within the inert section due to the low ther-
121
4. Hydrogen catalytic combustion in Pt-coated channels
195 ms
192.5 ms
190 ms
Figure 4.11: Gas-phase ignition for Case 3 in Table 4.4 (cordierite, TIN =
650 K, UIN = 10 m/s): (a) time evolution of iHRR (solid line) and iSHRR
(dashed line) and (b) time evolution of the spatial distribution of gas-phase
(solid lines) and catalytic (dashed lines) hydrogen conversion rates. The
curves in (b) refer to t = 190 ms, t = 192.5 ms and t = 195 ms.
mal conductivity of cordierite (see Fig. 4.7(a)). Farther downstream, the
presence of a region dominated by gas-phase chemistry reduced catalytic
conversion to almost zero. The time response of the surface coverages following gas-phase ignition is reported in Fig. 4.12. Surface coverages did
not change within the region dominated by catalytic chemistry, whereas
downstream the depletion of hydrogen by gas-phase combustion caused an
increase in the surface coverage of the excess reactant, which was oxygen,
and a decrease of all other surface species.
Stefan velocity (not plotted) was negligible during gas-phase ignition,
achieving a maximum value of about 0.5 mm/s within the zone dominated
by gas-phase combustion. Fig. 4.13 illustrates the gas-phase ignition sequence. Initially, the peak in OH, which was previously at the wall (see Fig.
122
4.3. Results and discussion
Figure 4.12: Time evolution of surface coverages during gas-phase ignition
for Case 3 in Table 4.4 (cordierite, TIN = 650 K, UIN = 10 m/s). The
curves refer to t = 190 ms, t = 192.5 ms and t = 195 ms and the arrows
indicate the direction of increasing time.
4.9) due to OH desorption, shifted to the gas-phase (Fig. 4.13(a)). The
runaway process started and a flame propagated downstream staying close
to the wall (Fig. 4.13(b)) where conditions were more favorable to gasphase combustion, i.e. the temperature and radical concentrations were
higher while the flow velocity was lower. Finally, an open (two separate
branches) flame was established (Fig. 4.13(c)). Due to the high DNS computational cost, the steady state could not be computed and simulations
were terminated at t = 200 ms. However, steady state was computed with
the QSS code which showed that, due to heat conduction towards the inert
section, the flame moved upstream and eventually the gas-phase ignition
123
4. Hydrogen catalytic combustion in Pt-coated channels
(a)
2
1.37E-04
y [mm]
1
0
2
(b)
0.00E+00
4.91E-04
(c)
0.00E+00
1.60E-03
1
0
2
1
0
0
2
4
6
x [mm]
8
10
0.00E+00
Figure 4.13: Two-dimensional maps of OH radical mass fractions during
gas-phase ignition for Case 3 in Table 4.4 (cordierite, TIN = 650 K, UIN =
10 m/s): (a) t = 190 ms, (b) t = 192.5 ms and (c) t = 195 ms.
point shifted within the inert zone. The other two cases with cordierite
wall but lower inlet temperatures (Case 2 and Case 3) showed a similar
gas-phase ignition process, although at steady state the gas-phase ignition
point could not reach the inert zone. For the three cordierite cases the
maximum temperature at steady state was superadiabatic. Superadiabatic
temperature is theoretically achieved under mass transport limited conversion and adiabatic operation in fuel lean catalytic combustion of hydrogen,
due to the lower than unity Lewis number and to the nature of the catalytic
diffusion reaction zone. Superadiabatic temperature can also be achieved
in practical non adiabatic combustion systems [148].
The FeCr-alloy wall cases were characterized by a higher temperature of
the upstream inert section before gas-phase ignition when compared to the
cordierite cases (see Fig. 4.7) because of their higher thermal conductivity
124
4.3. Results and discussion
(see Table 4.1). For this reason, the FeCr-alloy cases were more likely to
ignite within the inert section where the hydrogen concentration at the
wall was higher. The two cases with FeCr-alloy wall and TIN = 450 and
550 K could not be computed until the occurrence of gas-phase ignition
due to the relatively long physical time required (see Table 4.5) and to the
high computational cost of the DNS catalytic simulations, thus emphasis is
placed on the other two cases (Cases 6 and 7 in Table 4.4).
Figure 4.14: Gas-phase ignition for Case 6 in Table 4.4 (FeCr-alloy, TIN =
650 K, UIN = 10 m/s): (a) time evolution of iHRR (solid line) and iSHRR
(dashed line); time evolution of the axial profiles of hydrogen conversion
rates (b) gas-phase and (c) catalytic. The curves in (b) and (c) refer to t =
250 ms, t = 258 ms, t = 267.2 ms and t = 270 ms and the arrows indicate
the direction of increasing time.
Figure 4.14 shows the integral heat release rate and hydrogen conversion rate during gas-phase ignition for Case 6 in Table 4.4 (FeCr-alloy, TIN
= 650 K, UIN = 10 m/s). The initial time evolution of the iHRR and
125
y [mm]
4. Hydrogen catalytic combustion in Pt-coated channels
2
1
0
2
1
0
2
1
0
2
1
0
0
(a)
2
2.38E-05
(b)
0.00E+00
2.32E-04
(c)
0.00E+00
1.67E-03
(d)
0.00E+00
2.42E-03
4
6
x [mm]
8
10
0.00E+00
Figure 4.15: Two-dimensional maps of OH radical mass fractions during
gas-phase ignition for Case 6 in Table 4.4 (FeCr-alloy, TIN = 650 K, UIN =
10 m/s): (a) t = 250 ms, (b) t = 258 ms, (c) t = 267.2 ms and (d) t = 270
ms.
iSHRR in Fig. 4.14(a) was similar to that in Fig. 4.11(a), which referred
to a case with same operating conditions but cordierite wall, with the iHRR
increasing and the iSHRR subsequently decreasing due to competition for
hydrogen between gas-phase and surface mechanisms. At t = 267 ms the
FeCr-alloy case had a sharp peak in iHRR due to the occurrence of gasphase ignition within the inert section. Figures 4.14(b) and 4.14(c) show
the spatial distribution of gas-phase and catalytic hydrogen conversion rate,
respectively. Initially, gas-phase ignition occurred downstream of the catalytic section, while surface chemistry was still dominating the upstream
region. Later, a second ignition point appeared within the inert section and
was connected with the downstream gas-phase combustion zone. Catalytic
conversion decreased close to the beginning of the catalytic section unlike
the cordierite case, since the surface was effectively shielded by gas-phase
126
4.3. Results and discussion
combustion. The gas-phase hydrogen conversion rate had a minimum at
the location where the catalytic conversion rate was maximum (see Fig.
4.14(b,c)). Hydrogen leaked through the flame in a narrow region close to
the beginning of the catalytic section and it was converted on the surface.
Time evolution of surface coverages is not reported but it was similar to the
one of the cordierite case (Fig. 4.12) where the zone affected by gas-phase
combustion underwent an increase in oxygen coverage.
Figure 4.15 shows the gas-phase ignition sequence for Case 6 in Table 4.4. Gas-phase combustion started within the catalytic section (Fig.
4.15(a)), the flame propagated downstream close to the wall (Fig. 4.15(b)),
two ignition spots appeared within the inert section (Fig. 4.15(c)) and finally a nearly stationary flame was established (Fig. 4.15(d)), where the
ignition points were located in the inert section and gas-phase combustion
was weakened at the beginning of the catalytic section as previously discussed.
The other computed case with FeCr-alloy wall (Case 7 in Table 4.4,
TIN = 650 K, UIN = 5 m/s) achieved gas-phase ignition directly within
the inert zone. Figure 4.16 shows the ignition sequence: while gas-phase
ignition started within the catalytic section, two ignition points appeared
in the inert section (Fig. 4.16(a)) and developed quickly into a flame that
propagated downstream (Fig. 4.16(b)); at a certain time a “flame ball”
detached from the main flame (Fig. 4.16(c)) and consumed the unburnt
fuel left downstream, while a nearly stationary flame with an open tip, due
to the smaller than unity Lewis number, established (Fig. 4.16(d)). Ste-
127
y [mm]
4. Hydrogen catalytic combustion in Pt-coated channels
2
1
0
2
1
0
2
1
0
2
1
0
0
(a)
2
6.20E-06
(b)
0.00E+00
2.90E-03
(c)
0.00E+00
2.52E-03
(d)
0.00E+00
2.55E-03
4
6
x [mm]
8
10
0.00E+00
Figure 4.16: Two-dimensional maps of OH radical mass fractions during
gas-phase ignition for Case 7 in Table 4.4 (FeCr-alloy, TIN = 650 K, UIN =
5 m/s): (a) t = 310.84 ms, (b) t = 311.4 ms, (c) t = 311.72 ms and (d) t =
312.12 ms.
fan velocity and surface coverages adapted almost instantaneously to the
new conditions determined by the downstream propagating flame as seen
in Fig. 4.17(a), and Figs. 4.17(b,c), respectively. The steady state solution
computed with the QSS code showed that for all the FeCr-alloy cases the ignition points were located within the inert section and the wall temperature
was superadiabatic due to a residual catalytic hydrogen conversion.
Table 4.6 reports the hydrogen catalytic and gaseous conversions: gasphase combustion increased significantly hydrogen consumption while the
residual catalytic conversion after gas-phase ignition was very small and
confined in a narrow region close to the beginning of the catalytic section.
128
4.3. Results and discussion
Figure 4.17: Gas-phase ignition for Case 7 in Table 4.4 (FeCr-alloy, TIN =
650 K, UIN = 5 m/s): (a) time evolution of the spatial distribution of Stefan
velocity; time evolution of the spatial distribution of (b) Pt(s) coverage and
(c) O(s) coverage. The curves refer to t = 311.16 ms (solid line), t = 311.40
ms (dashed line), t = 311.60 ms (dashed-dotted line) and t = 311.80 ms
(dotted line).
4.3.2.4
Comparison between DNS and QSS code results
DNS results were compared with results from the QSS code developed
in [135, 136]. Within the QSS assumption, gas-phase processes are considered to be fast enough to equilibrate within an integration time step of
the transient solid energy equation. Therefore, while the fully transient energy equation is integrated for the solid domain, gas-phase equations are
solved with a steady solver at every time step. The minimum time step for
the QSS code was limited by the characteristic time for radial gas-phase
diffusion (tg,r ) and it was computed as tg,r = (h/2)2 /αg , where αg is the
129
4. Hydrogen catalytic combustion in Pt-coated channels
thermal diffusivity of the gas-phase. Within the temperature range 450-650
K the characteristic time tg,r varied between 16 ms and 8 ms, respectively,
therefore a time step of 20 ms for the cases with TIN = 650 K and a time
step of 25 ms for the other cases were chosen, sufficiently long for gas-phase
equilibration. The solid wall was discretized only in the axial direction since
the characteristic time for transverse conduction was much faster (1.19 ms
for cordierite and 0.17 ms for FeCr-alloy) than tg,r and hence it could not
be resolved within the gas-phase quasisteady assumption.
O(s)
Pt(s)
OH(s)
H2O(s)
H(s)
Figure 4.18: Comparison between DNS and QSS code for Case 1 in Table
4.4 at t = 25 ms: (a) comparison of the gas/solid interfacial wall temperature, and (b) comparison of the surface coverages. Solid lines are DNS and
symbols are QSS predictions.
A first obvious remark is that the light-off and gas-phase ignition processes could not be studied with the QSS code under the conditions in Table
130
y [mm]
4.3. Results and discussion
2
1
0
2
1
0
0
t = 375 ms
t = 375 ms
2
0.00E+00
2.10E-04
6.32E-04 0.00E+00
0.00E+00
y [mm]
QSS
t = 350 ms
DNS
t = 350 ms
4
6
x [mm]
8
10 0
2
1.16E-03 0.00E+00
4
6
x [mm]
8
10
8.10E-04
Figure 4.19: Comparison between DNS and QSS code predictions of 2-D
OH mass fractions for Case 1. The solutions of the two codes are compared
immediately before and immediately after the gas-phase ignition time.
4.4 since they were faster than the QSS time resolution: light-off duration
was less than one millisecond for all cases (see Table 4.4), while gas-phase
ignition duration ranged from a few milliseconds (see Fig. 4.13 and Fig.
4.16) to a time period comparable to the QSS time step (see Fig. 4.15).
Comparison between QSS and DNS results was carried out for Case 1, for
which slow chemistry effects (due to lower TIN ) could potentially invalidate
the QSS assumption. Figure 4.18 shows comparison between DNS and QSS
code at t = 25 ms, which is the first time step available with the QSS code,
in terms of interfacial gas/solid temperature (Fig. 4.18(a)) and surface coverage (Fig. 4.18(b)). The agreement is quite good, with the maximum
difference between the two temperature profiles being about 30 K at the
beginning of the catalytic zone. The discrepancy was attributed to the
presence of transverse gradients in the solid that could not be resolved with
131
4. Hydrogen catalytic combustion in Pt-coated channels
the 1-D discretization of the QSS code. Figure 4.19 shows the comparison
of 2-D OH maps between QSS and DNS codes for Case 1, immediately
before and after gas-phase ignition. The QSS code slightly underpredicted
OH mass fractions due to a lower peak temperature owing to the 1-D discretization of the wall. Nonetheless the time interval within which gas-phase
ignition occurred was correctly predicted by the QSS code and the solution
converged to the DNS one after gas-phase ignition.
The QSS code never failed in describing the correct physics at any
condition considered in Table 4.4, because surface chemistry was still fast
enough for the quasisteady-state assumption to hold. A further comparison
between QSS and DNS code is given in Fig. 4.7, providing the time evolution of surface temperature during wall heat-up for Cases 3 and 6 with
the highest inlet temperature (TIN = 650 K) and different wall materials.
During the wall heat-up phase there were no substantial differences between
DNS and QSS code. Table 4.5 shows gas-phase ignition times for the DNS
and QSS code predictions. The ignition time for DNS has been defined
as the time instant when the time derivative of the iHRR was maximum
while for the QSS code a visual inspection of the solution sufficed. The
agreement between the ignition times computed with the two codes was
again very good.
The initial surface coverage did not impact the comparison. All DNS
simulations were repeated with an initial coverage of either 100% H(s) or
100% Pt(s): in these cases light-off was almost instantaneous, as previously
mentioned, but just after 2 ms the differences with the solution obtained
132
4.3. Results and discussion
with an initial coverage of 100% O(s) were negligible.
At the initial temperatures of Table 4.4 surface chemistry was still
too fast to invalidate the quasisteady-state assumption: the only limiting
reaction was oxygen desorption during light-off and, once some free Pt sites
were available for hydrogen adsorption, light-off could be achieved leading
quickly to the mass transport limited regime. In order to observe the failure
of the QSS code, a lower initial temperature was required.
A test case was run using the same conditions as Case 3 in Table
4.4 (cordierite, UIN = 10 m/s, TIN = 650 K) but starting from an initial
temperature of 400 K. In order to make the comparison with the QSS
assumption fairer, at t = 0 s the fluid at 400 K inside the channel was
instantaneously replaced with fluid at 650 K since the QSS code could
not simulate the process of injection. The initial surface coverage was 100
% H(s), unlike the other simulated cases, because the simulation started
below the light-off temperature, as will be subsequently discussed. An initial
surface coverage of 100 % H(s) was a realistic condition when the catalyst
was in contact with the H2 /air mixture below light-off temperature: Fassihi
et al. [149] and Rinnemo et al. [150] reported that the catalytic surface
was initially covered with hydrogen when pH2 /(pH2 + pO2 ) > 0.1 − 0.2 and
in our case already pH2 /(pH2 + pO2 ) = 0.375. The numerical experiment
had similarities with the catalytic ignition studies in [139, 149–151] where,
starting from a low initial temperature, either a platinum wire or a foil were
heated through an external source until self-acceleration of the reaction rate
and therefore light-off was achieved. In our case the temperature of the
133
4. Hydrogen catalytic combustion in Pt-coated channels
catalytic surface was increased from the initial 400 K through heat transfer
with the hotter flow at TIN = 650 K. The increase of the wall temperature
was necessary to achieve light-off unlike the other cases previously described
y [mm]
where the temperature was already high enough for light-off.
2
1
0
2
1
0
2
1
0
2
1
0
0
DNS t = 17 ms
DNS t = 18 ms
DNS t = 20 ms
QSS t = 20 ms
2
0.00E+00
6
4
x [mm]
8
10
8.73E-03
Figure 4.20: Two-dimensional maps of H2 mass fraction for a case with
cordierite wall, UIN = 10 m/s, TIN = 650 K and initial temperature equal
to 400 K.
Figure 4.20 shows hydrogen mass fraction distributions obtained with
the DNS code during the light-off sequence and the comparison with the
QSS code after 20 ms. At t = 17 ms the hydrogen conversion was still
kinetically limited everywhere. The hydrogen concentration at the wall
decreased with increasing axial distance due to conversion and, since the
light-off temperature increased with rising hydrogen concentration due to
the well known self-poisoning effect [149, 150], conditions were more favorable close to the channel outlet for light-off to happen. Figure 4.21 shows
134
4.3. Results and discussion
Figure 4.21: Time evolution of (a) hydrogen coverage and (b) oxygen coverage for a case with cordierite wall, UIN = 10 m/s, TIN = 650 K and initial
temperature equal to 400 K. Solid lines are DNS solutions and the dashed
line is the QSS solution after the first time step.
that the surface coverage at t = 17 ms was still hydrogen dominated although at the channel outlet the hydrogen coverage had dropped to 85 %
(note the logarithmic scale). Using a simplified model Fassihi et al. [149]
predicted that catalytic ignition occurred when the hydrogen coverage was
95 % for pH2 /(pH2 + pO2 ) > 0.17, but of course they had different heat
transfer characteristics and surface mechanism compared to our case.
At t = 18 ms light-off occurred close to the channel outlet and then
propagated upstream as can be seen in Fig. 4.20 where the light-off distance
is identified by a sharp drop in the hydrogen mass fraction at the wall. The
surface coverage switched from hydrogen-dominated to oxygen-dominated
135
4. Hydrogen catalytic combustion in Pt-coated channels
Figure 4.22: Gas/solid interfacial temperature profile for a case with
cordierite wall, UIN = 10 m/s, TIN = 650 K and initial temperature equal
to 400 K. Solid lines are DNS solutions and the dashed line is the QSS
solution after the first time step.
(see Fig. 4.21), in agreement with [149, 150]. The first solution available
from the QSS code was at t = 20 ms, which could be compared with that
from the DNS code at the same time. Figure 4.20 shows that the hydrogen
mass fraction at the outlet was higher for the DNS code compared to the
QSS code. Furthermore, the QSS code predicted that the surface reactions
were already mass transport limited everywhere while according to the DNS
code there was a border between kinetically limited and mass transport limited zones at x ≈ 2 mm. Figure 4.21 shows that according to the QSS code
the surface was already oxygen dominated everywhere at t = 20 ms, while
in the DNS results there was still part of the surface which was hydrogen
dominated. Finally, Fig. 4.22 highlights a large difference in wall tempera136
4.4. Conclusions
ture between the two codes: since the QSS code predicted light-off within
its first time step, the surface had a higher temperature due to the bigger
chemical heat release rate.
The results presented in this section highlight the importance of using
DNS for fundamental studies of light-off and gas-phase ignition, and as a
benchmark to assess the validity of the QSS assumption.
4.4
Conclusions
Transient detailed 2-D numerical simulations of fuel-lean hydrogen hetero/homogeneous combustion in a mesoscale channel were performed by means
of a spectral element code already including detailed transport, gas-phase
and surface chemistry, whose capability was further augmented by implementing and validating 2-D conjugate heat transfer. The numerical code,
which was free from any assumption, allowed for the first time the resolution
of all the spatiotemporal scales in a practical 2-D catalytic channel geometry. A parametric numerical study was carried out, where the influence of
wall material (cordierite and FeCr-alloy), inlet velocity (UIN = 5, 10 m/s)
and inlet temperature (TIN = 450, 550, 650 K) on the startup process, from
light-off until gas-phase ignition, was investigated.
Light-off was almost instantaneous for any operating condition due to
the fast heterogeneous chemistry. Even at the lower initial temperature, i.e.
450 K, catalytic ignition occured after only 0.33 ms, while for cases with an
initial temperature of 650 K it happened at the first reported integration
137
4. Hydrogen catalytic combustion in Pt-coated channels
time step (4 ns). Flow velocity and wall material had no influence since
thermal feedback did not play any role given the high heat capacity of
the solid. CSP analysis allowed the identification of the reaction sequence
responsible for the catalytic runaway: after an initial buildup of OH(s),
water production and desorption accelerated, leading to an increase of free
Pt sites available for hydrogen adsorption. Stefan velocity achieved a peak
value of about 11 cm/s and had a strong impact on the flow field; such a
result has been observed for the first time in a practical channel geometry.
The wall heat-up phase was mainly influenced by the material properties, in particular the thermal conductivity since heat capacity was very similar for the two materials. Given the lower thermal conductivity, cordierite
had sharper temperature peaks that could promote earlier gas-phase ignition within the catalytic section, while FeCr-alloy had higher heat conduction towards the inert section where conditions for gas-phase ignition
were more favorable. The lower thermal conductivity of cordierite was also
responsible for higher thermal gradients leading to more severe thermallyinduced mechanical stresses.
CSP analysis identified OH desorption from the catalytic surface as the
main process leading to gas-phase ignition, while heat conduction from the
hot wall played only a secondary role. Following downstream propagation
of the flame, a fast rearrangement of surface coverages could be observed,
with the time step being reduced to 1 ns due to the high stiffness. The total
duration of gas-phase ignition ranged from about 1.5 ms for the case at TIN
= 650 K and UIN = 5 m/s, to about 17 ms for the case at TIN = 450 K
138
4.4. Conclusions
and UIN = 10 m/s.
Predictions of the DNS code were compared with those of a code relying
on the quasisteady-state (QSS) assumption in order to define the regime of
validity of the latter. The time step of the QSS code was limited by the
time required to equilibrate the gas phase and it was 20 ms for cases at TIN
= 650 K, and 25 ms for all other cases. Under the operating conditions used
for DNS simulations, predictions of the QSS code were in good agreement
since chemistry was fast enough to reach equilibrium within one time step
of the QSS code. Nonetheless, the QSS code could not address the study
of fast processes like light-off and gas-phase ignition whose duration was
shorter than the QSS integration time step.
An additional test case was run with TIN = 650 K but a lower initial
temperature for the solid equal to 400 K in order to bring the QSS assumption to its limit. In this case comparison of the DNS and QSS predictions
after one QSS time step, i.e. 20 ms, showed a considerable discrepancy with
the QSS code predicting a surface temperature more than 100 K higher than
the DNS temperature. Furthermore, while DNS predicted the achievement
of light-off only on part of the catalytic surface after 20 ms, the QSS code
predicted that light-off was already achieved everywhere.
The reason for the high computational cost of transient DNS simulations was investigated in a transient SPSR by means of CSP analysis. The
hydrogen surface reaction mechanism exhibited an extremely high stiffness
with the fastest time scale equal to 10−12 s pertaining to hydrogen adsorption and the reaction O(s) + H(s) => OH(s) + Pt(s), which had the same
139
4. Hydrogen catalytic combustion in Pt-coated channels
order of magnitude as the fastest time scale of n-heptane detailed gas-phase
reaction mechanism with 560 species and 2538 reactions. For comparison,
the fastest time scale for a one-step hydrogen surface reaction mechanism
was 10−6 s.
In conclusion, despite the much higher computational cost of DNS
compared to simplified models, fundamental studies cannot avoid relying
on DNS due to the need to resolve all the relevant spatiotemporal scales.
Codes relying on the QSS assumption are computationally affordable and
robust tools for design purposes, but they have to be benchmarked against
DNS results in order to define the range of operating conditions under which
their predictions are reliable.
140
Chapter 5
Conclusions and future work
5.1
Conclusions
A Direct Numerical Simulation solver for low-Mach-number reactive flows
including detailed transport, and elementary gas-phase and heterogeneous
chemistry, was further developed by implementing conjugate heat transfer,
and was successfully validated against literature analytical and numerical
results. The code allows time-accurate simulation of hetero-/homogeneous
combustion by resolving for the first time all the relevant spatiotemporal
scales and represents an important step in fundamental catalytic combustion. The main goal of this thesis was to investigate numerically hetero/homogeneous combustion of lean hydrogen and syngas/air mixtures in
mesoscale channels, which assumes great relevance in the development of
“zero-emissions” power generation at both small and large scales. The major
findings and conclusions of this work are summarized below.
In the first part, numerical simulations of gas-phase combustion in lean
141
5. Conclusions and future work
syngas/air mixtures in an inert mesoscale channel have been performed
with the aim to map combustion instabilities. The starting point of this
study was an experimental campaign conducted with OH* chemiluminescence in a 7-mm height channel, which showed the existence of many stable
and oscillatory combustion modes. The code was firstly validated against
steady experimental flames whose OH* chemiluminescence signal was detected. Starting from the operating conditions where an oscillatory flame
was found experimentally, a parametric numerical study was carried out
wherein the wall temperature was varied within the range of measured wall
temperature.
New unsteady combustion modes have been numerically revealed: a
first periodic mode named “oscillatory ignition”, and a second non-periodic
mode named “random ignition spots”. Weak combustion has been observed
for the first time at relatively high flow velocities. In terms of physical understanding, two were the major results of the study. For the “oscillatory
ignition” mode, the interplay between two fuels, i.e. hydrogen and carbon
monoxide, with Lewis number smaller and larger than unity, respectively,
was found to be a possible explanation for the observed instability along
with the heat exchange with the wall. Regarding the weak combustion
mode, CSP analysis has been used as diagnostic tool and the results supported the hypothesis that competition for OH radical between H2 and CO
could be the driving force for weak combustion. Finally, a stability map
in the wall temperature/CO:H2 ratio space has been obtained. The results
of a recent experimental campaign adopting fast OH-LIF compared favor-
142
5.1. Conclusions
ably against the computed combustion modes. Steady flames showed an
excellent agreement while dynamic modes, in particular oscillatory ignition,
were correctly predicted but imposing a different wall temperature profile.
Possible reasons for the discrepancy might be approximations in the numerical code, in particular the fact that radiative heat transfer was neglected,
the adequacy of the detailed reaction mechanism in reproducing key reaction rates, and 3-D effects, phenomena to which the dynamic modes are
particularly sensitive.
In the second part, the DNS solver was used to investigate hetero/homogeneous combustion of lean hydrogen/air mixtures in a 2-mm height
channel, including heat conduction in the solid wall. Given the absence of
any limitation in the spatiotemporal resolution, the code could be used for
fundamental studies on transient processes, e.g light-off and gas-phase ignition, and for benchmarking codes relying on the quasisteady-state (QSS)
assumption. A parametric numerical study was carried out, where the influence of inlet velocity, inlet temperature and wall material was investigated.
For all the simulated cases, light-off was much shorter than 1 ms,
therefore it could not be studied with the QSS code. The key reactions
responsible for catalytic ignition were identified, involving initially oxygen
desorption and OH(s) buildup, while the runaway was triggered by water
production and desorption. Stefan velocity achieved a peak value of 11
cm/s and had a large impact on the flow field, albeit for a very short time
period. This effect has been observed for the first time in multidimensional
systems. During the wall heat-up phase, the wall material played a con-
143
5. Conclusions and future work
trolling role and the DNS code allowed computing the thermal gradients in
the solid responsible for thermally induced mechanical stresses, while the
1-D wall discretization of the QSS code did not allow such an estimation.
Finally, gas-phase ignition was achieved and the DNS code allowed to study
the flame propagation process and hetero-/homogeneous pathways coupling.
Surface species adapted almost instantaneously to the conditions dictated
by the propagating flame due to the very fast surface chemistry. The QSS
code could not be used to study gas-phase ignition since the latter lasted
from a few milliseconds to a time period comparable to the QSS code time
step. CSP analysis identified the net OH desorptive flux from the catalytic
surface as the main process responsible for runaway in gas-phase ignition.
For all the simulated cases the QSS code, although not being able to
resolve the fast transient processes, showed a very good agreement with
the DNS code. A new test case at a lower initial temperature, i.e. 400 K,
was run and under these conditions the QSS code predicted an earlier lightoff and therefore a higher wall temperature due to the invalidation of the
fast surface chemistry assumption. The results exemplified the usefulness
of DNS for fundamental studies and for assessing the range of validity of
codes based on simplified assumptions like the QSS one.
5.2
Future Work
Although the DNS code available at the Aerothermochemistry and Combustion Systems Laboratory of ETH Zurich already represents a state-of-the-
144
5.2. Future Work
art tool in combustion research, further developments have been identified
including both the implementation of new physical processes and the improvement of the computational performance.
Regarding the improvement of the physical model, the implementation
of radiative heat transfer, at least between the solid walls, could significantly augment the code predicting capability given its relevance at both
steady and transient reactor operation [60, 152]. A preliminary study to
implement radiative heat transfer has been undertaken and the main identified problem consisted in the increase of the computational cost. This was
because radiative heat transfer necessitates iterations within a single time
step and communication between different CPUs on a parallel architecture,
as the interacting surface grid points are distributed among many different
processors.
For future studies, the real challenge would be to simulate turbulent
combustion in catalytic channels with detailed gas-phase and surface chemistry. The increased pressure requirements of new generation turbines have
resulted in incoming Reynolds numbers of up to 30,000 based on the hydraulic diameter of an individual channel of the honeycomb reactor [19, 20].
However, turbulent catalytic combustion modeling has not yet received
proper attention. Up to now, the DNS code could simulate turbulent catalytic combustion only with one-step surface chemistry [147] due to the
excessively high computational cost of detailed surface chemistry owing to
the higher stiffness as discussed in Section 4.3.1. With the aid of CSP
analysis, the present work identified the physical reasons for the stiffness of
145
5. Conclusions and future work
the detailed surface reaction mechanism, with the fastest time scale associated to hydrogen adsorption. Further code development is required but a
starting point has already been identified and the first changes have been
implemented into the DNS code.
A first idea was to replace the current implicit integrator used for the
thermochemistry part, i.e. CVODE [55], with a new integrator which is
particularly suited for multirate problems, i.e. ARKode [153]. ARKode
is a time integration package for stiff, nonstiff and multi-rate systems of
ordinary differential equations, which is based on adaptive-step additive
Runge-Kutta methods. The right-hand side of the equations can be partitioned into two components, a first one containing the “slow” time scales
to be integrated explicitly, and a second one containing the “fast” time
scales to be integrated implicitly. The partition of the right-hand side into
“fast” and “slow” time scales can be guided by CSP analysis. The multirate formulation seems to be particularly suited for hetero-/homogeneous
combustion problems with conjugate heat transfer, where the time scales
due to chemistry are much faster than those associated to diffusive (in the
gas and solid) and convective processes. At the moment, ARKode package
has been coupled to the DNS code and successfully validated against the
results of the simulations presented in Chapter 4, but all the attempts to efficiently split the right-hand side of the equations resulted in an even higher
computational cost. Therefore additional work is needed.
A second approach was to couple the DNS code with a differential
algebraic solver in order to implement a quasisteady-state assumption that
146
5.2. Future Work
was less restrictive than the one presented in the thesis, where only the
solid energy equation was solved as transient. CSP analysis can identify
the candidate species to be treated as steady, while transient equations are
solved for the other species and temperature. Table 4.2 shows that the three
fastest time scales are associated exclusively to H(s), H2 O(s) and OH(s), and
there is a three orders of magnitude gap between the third and the fourth
fastest time scale, which means that treating the three surface species as
steady within an integration time step would be a reasonable assumption.
The idea is to find a trade-off between the advantages of QSS codes in terms
of lower computational costs, with the absence of time step restrictions of
DNS codes.
This formulation was first tested on an in-house quasisteady-state SPSR
and preliminary results showed a negligible difference between the fully transient solution and the quasisteady-state solution, both in the description of
the transient and of the steady state. The same approach was tested on
the DNS code by replacing the CVODE integrator with the differentialalgebraic solver IDA [154]. The approach has been successfully validated
against the results presented in Chapter 4 but again the computational cost
was higher. An important remark is that IDA had a much higher ratio of
linear to nonlinear number of iterations compared to CVODE and ARKode,
which suggests that a preconditioner would be particularly beneficial in this
case.
147
Bibliography
[1]
A. C. Fernandez-Pello. Micropower generation using combustion: issues and approaches. Proc. Combust. Inst., 29:883–899, 2002.
[2]
Y. Ju and K. Maruta. Microscale combustion: Technology development and fundamental research. Prog. Energ. Combust., 37:669–715,
2011.
[3]
K. Kang, Y. S. Meng, J. Breger, C. P. Grey, and G. Ceder. Electrodes
with high power and high capacity for rechargeable lithium batteries.
Science, 311:977–980, 2006.
[4]
N. Chigier and T. Gemci. A review of micro propulsion technology.
41st aerospace sciences meeting and exhibit, 2003.
[5]
K. Maruta, K. Takeda, L. Sitzki, K. Borer, P. D. Ronney, S. Wussow, and O. Deutschmann. Catalytic combustion in microchannel for
MEMS power generation. Third Asia-Pacific conference on combustion, 2001.
[6]
J. Vican, B. F. Gajdeczko, F. L. Dryer, D. L. Milius, I. A. Aksay,
and R. A. Yetter. Development of a microreactor as a thermal source
for microelectromechanical systems power generation. Proc. Combust.
Inst., 29:909–916, 2002.
[7]
K. Fu, A. Knobloch, F. Martinez, D. Walther, A. C. Fernandez-Pello,
A. Pisano, D. Liepmann, K. Miyaska, and K. Maruta. Design and experimental results of small-scale rotary engines. Proc. IMECE-ASME,
2001.
149
Bibliography
[8]
A. H. Epstein. Millimeter-scale, micro-electro-mechanical systems gas
turbine engines. J. Eng. Gas Turb. Power Trans. ASME, 126:205–226,
2004.
[9]
K. Isomura, S Tanaka, S-i Togo, and M. Esashi. Development of highspeed micro- gas bearings for three-dimensional micro-turbo machines.
J. Micromech. Microeng., 15:222–227, 2005.
[10] D. Lewis, S. Janson, R. Cohen, and E. Antonsson. Digital micropropulsion. Sensor Actuat. A-Phys., 80:143–154, 2000.
[11] A. London, A. Ayon, A. H. Epstein, S. Spearing, T. Harrison, Y. Peles,
and et al. Microfabrication of a high pressure bipropellant rocket
engine. Sensor Actuat. A-Phys., 92:351–357, 2001.
[12] C. Miesse, R. Masel, M. Short, and M. Shannon. Experimental
observations of methane-oxygen diffusion flame structure in a submillimeter microburner. Combust. Theor. Model., 9(1):77–92, 2005.
[13] J. Ahn, C. Eastwood, L. Sitzki, and P. D. Ronney. Gas-phase and
catalytic combustion in heat -recirculating burners. Proc. Combust.
Inst., 30:2463–2472, 2005.
[14] F. J. Weinberg and S. A. Lloyd. A burner for mixtures of very low
heat content. Nature, 251:47–49, 1974.
[15] T. Takeno and K. Sato. An excess enthalpy flame theory. Combust.
Sci. Technol., 20:73–84, 1979.
[16] S. Karagiannidis, J. Mantzaras, G. Jackson, and K. Boulouchos.
Hetero-/homogeneous combustion and stability maps in methanefueled catalytic microreactors. Proc. Combust. Inst., 31:3309–3317,
2007.
[17] G. Pizza, J. Mantzaras, C. E. Frouzakis, A. G. Tomboulides, and
K. Boulouchos. Suppression of combustion instabilities of premixed
hydrogen/air flames in microchannels using heterogeneous reactions.
Proc. Combust. Inst., 32:3051–3058, 2009.
[18] G. Pizza, J. Mantzaras, and C. E. Frouzakis. Flame dynamics in
catalytic and non-catalytic mesoscale microreactors. Catal. Today,
155:123–130, 2010.
150
Bibliography
[19] R. Carroni, T. Griffin, J. Mantzaras, and M. Reinke. High-pressure
experiments and modeling of methane/air catalytic combustion for
power-generation applications. Catal. Today, 83:157–170, 2003.
[20] K. W. Beebe, K. D. Cairns, V. K. Pareek, S. G. Nickolas, J. C. Schlatter, and T. Tsuchiya. Development of catalytic combustion technology
for single-digit emissions from industrial gas turbines. Catal. Today,
59:95–115, 2000.
[21] R. Carroni, V. Schmidt, and T. Griffin. Catalytic combustion for
power generation. Catal. Today, 75:287–295, 2002.
[22] L. L. Smith, H. Karim, M. J. Castaldi, S. Etemad, and W. C. Pfefferle.
Rich-catalytic lean-burn combustion for fuel-flexible operation with
ultra-low emissions. Catal. Today, 117:438–446, 2006.
[23] T. Griffin, D. Winkler, M. Wolf, C. Appel, and J. Mantzaras. Staged
catalytic combustion method for the advanced zero emissions gas turbine power plant. ASME, pages 2004–54101, 2004.
[24] A. Schneider, J. Mantzaras, and P. Jahnson. Experimental and numerical investigation of the catalytic partial oxidation of CH4 /O2 mixtures diluted with H2 O and CO2 in a short contact time reactor. Chem.
Eng. Sci., 61:4634–4646, 2006.
[25] Y. Ghermay, J. Mantzaras, and R. Bombach. Effects of hydrogen preconversion on the homogeneous ignition of fuel-lean H2 /O2 /N2 /CO2
mixtures over platinum at moderate pressures. Combust. Flame,
157:1942–1958, 2010.
[26] S. Dunn. Hydrogen futures: toward a sustainable energy system. Int.
J. Hydrogen Energ., 27:4933–4943, 2002.
[27] K. Li, R. Zhang, and J. Bi. Experimental study on syngas production by co-gasification of coal and biomass in a fluidized bed. Int. J.
Hydrogen Energ., 35:2722–2726, 2010.
[28] R. Domenichini, M. Gallio, and A. Lazzaretto. Combined production of hydrogen and power from heavy oil gasification: Pinch analysis, thermodynamic and economic evaluations. Energy, 35:2184–2193,
2010.
151
Bibliography
[29] M. He, Z. Hu, B. Xiao, J. Li, X. Guo, S. Luo, F. Yang, Y. Feng,
G. Yang, and S. Liu. Hydrogen-rich gas from catalytic steam gasification of municipal solid waste (MSW): Influence of catalyst and
temperature on yield and product composition. Int. J. Hydrogen Energ., 34:195–203, 2009.
[30] M. J. Castaldi and J. P. Doocher. Investigation into a catalytically
controlled reaction gasifier (CCRG) for coal to hydrogen. Int. J. Hydrogen Energ., 32:4170–4179, 2007.
[31] N. V. Gnanapragasam, B. V. Reddy, and M. A. Rosen. Hydrogen production from coal gasification for effective downstream CO2 capture.
Int. J. Hydrogen Energ., 35:4933–4943, 2010.
[32] P. Chiesa, G. Lozza, A. Malandrino, M. Romano, and V. Piccolo.
Three-reactors chemical looping process for hydrogen production. Int.
J. Hydrogen Energ., 33:2233–2245, 2008.
[33] C. C. Cosmos. Evaluation of power generation schemes based on
hydrogen-fuelled combined cycle with carbon capture and storage
(CCS). Int. J. Hydrogen Energ., 36:3726–3738, 2011.
[34] L. Tock and F. Marechal. H2 processes with CO2 mitigation: thermoeconomic modeling and process integration. Int. J. Hydrogen Energ.,
37:11785–11795, 2012.
[35] D. G. Norton, E. D. Wetzel, and D. G. Vlachos. Fabrication of singlechannel catalytic microburners: effect of confinement on the oxidation
of hydrogen/air mixtures. Ind. Eng. Chem. Res., 43:4833–4840, 2004.
[36] J Thormann, L. Maier, P. Pfeifer, U. Kunz, O. Deutschmann, and
K. Scubert. Steam reforming of hexadecane over Rh/CeO2 catalyst
in microchannels: experimental and numerical investigation. Int. J.
Hydrogen Energ., 34:5108–5120, 2009.
[37] E. Lopez, A. Irigoyen, T. Trifonov, A. Rodriguez, and J. A. Llorca.
A million-channel reformer on a fingertip: moving down the scale in
hydrogen production. Int. J. Hydrogen Energ., 35:3472–3479, 2010.
[38] A. J. Santis-Alvarez, M. Nabavi, N. Hild, D. Poulikakos, and W. J.
Stark. A fast hybrid start-up process for thermally self-sustained catalytic n-butane reforming in micro-sofc power plants. Energ. Environ.
Sci., 4:3041–3050, 2011.
152
Bibliography
[39] U. Izquierdo, V. L. Barrio, J. F. Cambra, J. Requies, M. B. Guemez,
P. L. Arias, and et al. Hydrogen production from methane and natural gas steam reforming in conventional and microreactor reaction
systems. Int. J. Hydrogen Energ., 37:7026–7033, 2012.
[40] S. G. Kerkemeier. Direct numerical simulation of combustion on petascale platforms: application to turbulent non-premixed hydrogen autoignition. PhD thesis, ETH Zurich, 2010.
[41] A. Majda and J. Sethian. The derivation and numerical solution of the
equations for zero Mach number combustion. Combust. Sci. Technol.,
42:185–205, 1985.
[42] R. B. Bird, W. E. Stewart, and E. N. Lightfoot. Transport phenomena.
John Wiley and Sons, Inc., 1960.
[43] J. O. Hirschfelder and C. F. Curtiss. Theory of propagation of flames.
Part I: General equations. Symposium on Combustion and Flame,
and Explosion Phenomena, 3:121–127, 1949.
[44] V. Giovangigli. Convergent iterative methods for multicomponent
diffusion. Impact of computing in science and engineering, 3:244–276.,
1991.
[45] R. J. Kee, G. Dixon-Lewis, J. Warnatz, M. E. Coltrin, and J. A. Miller.
A Fortran computer code package for the evaluation of gas-phase multicomponent transport properties, Report No. SAND86-8246. Technical report, Sandia National Laboratories, 1996.
[46] T. P. Coffee and J. M. Heimerl. Transport algorithms for premixed
laminar, steady-state flames. Combust. Flame, 43:273–289, 1982.
[47] M. E. Coltrin, R. J. Kee, F. M. Rupley, and E. Meeks. Surface
Chemkin: A fortran package for analyzing heterogeneous chemical
kinetics at a solid surface/gas-phase interface, Report No. SAND968217. Technical report, Sandia National Laboratories, 1996.
[48] H. Motz and H. Wise. Diffusion and heterogeneous reaction. III. atom
recombination at a catalytic boundary. J. Chem. Phys., 31:1893–1894,
1960.
[49] A. T. Patera. A spectral element method for fluid dynamics: Laminar
flow in a channel expansion. J. Comput. Phys., 54:468–488, 1984.
153
Bibliography
[50] M.O. Deville, P.F. Fischer, and E.H. Mund. High-order methods for incompressible fluid flow. Cambridge University Press, New-York, 2002.
[51] D. Gottlieb and S. Orzag. Numerical analysis of spectral methods:
theory and applications. SIAM, Philadelphia, 26, 1977.
[52] C. Canuto, A. Hussaini, A. Quarteroni, and T. A. Zang. Spectral
methods in fluid dynamics. Berlin; New York: Springer-Verlag, 1988.
[53] P. F. Fischer, J. W. Lottes, and S. G. Kerkemeier. Nek5000 web page,
2008. http://nek5000.mcs.anl.gov.
[54] A. G. Tomboulides, J. C. Y. Lee, and S. A. Orszag. Numerical simulation of low Mach number reactive flows. J. Sci. Comput., 12:139–167,
1997.
[55] G. D. Byrne and A. C. Hindmarsh. PVODE, an ODE solver for
parallel computers. Int. J. High Perform. C., 13:354–365, 1999.
[56] S. Orzag, M. Israeli, and M.O. Deville. Boundary conditions for incompressible flows. J. Sci. Comput., 1:75–111, 1986.
[57] R. J. Kee, F. M. Rupley, E. Meeks, and J. A. Miller. Chemkin: A Fortran chemical kinetics package for the analysis of gas-phase chemical
kinetics, Report No. SAND96-8216. Technical report, Sandia National
Laboratories, 1996.
[58] F. M. Rupley, R. J. Kee, and J. A. Miller. Premix: A Fortran program
for modeling steady laminar one-dimensional premixed flames, Report
No. SAND85-8240. Technical report, Sandia National Laboratories,
1995.
[59] G. Pizza. Numerical simulation of premixed flame dynamics in catalytic of and non-catalytic ducts. PhD thesis, ETH Zurich, 2009.
[60] U. Dogwiler, P. Benz, and J. Mantzaras. Two-dimensional modeling
for catalytically stabilized combustion of a lean methane-air mixture
with elementary homogeneous and heterogeneous chemical reactions.
Combust. Flame, 116:243–258, 1999.
[61] U. Dogwiler, J. Mantzaras, P. Benz, B. Kaeppeli, R. Bombach, and
A. Arnold. Homogeneous ignition of methane-air mixtures over platinum: comparison of measurements and detailed numerical predictions. Proc. Combust. Inst., 27:2275–2282, 1998.
154
Bibliography
[62] C. Appel, J. Mantzaras, R. Schaeren, R. Bombach, A. Inauen,
B. Kaeppeli, B. Hemmerling, and A. Stampanoni. An experimental
and numerical investigation of homogeneous ignition in catalytically
stabilized combustion of hydrogen/air mixtures over platinum. Combust. Flame, 128:340–368, 2002.
[63] Y. Ghermay, J. Mantzaras, and R. Bombach.
Experimental
and numerical investigation of hetero-/homogeneous combustion of
CO/H2 /O2 /N2 mixtures over platinum at pressures up to 5 bar. Proc.
Combust. Inst., 33:1827–1835, 2011.
[64] M. Schultze, J. Mantzaras, R. Bombach, and K. Boulouchos. An
experimental and numerical investigation of the hetero-/homogeneous
combustion of fuel-rich hydrogen/air mixtures over platinum. Proc.
Combust. Inst., 34:2269–2277, 2013.
[65] A. Pozzi and M. Lupo. The coupling of conduction with forced convection in Graetz problems. J. Heat Transf., 112:323–328, 1990.
[66] J. Sucec. Unsteady forced convection with sinusoidal duct wall generation: the conjugate heat transfer problem. Int. J. Heat Mass Tran.,
45:1631–1642, 2002.
[67] S. H. Lam. Using CSP to understand complex chemical kinetics. Combust. Sci. Technol., 89:375–404, 1993.
[68] S. H. Lam and D. A. Goussis. The CSP method for simplifying kinetics. Int. J. Chem. Kinet., 26:461–486, 1994.
[69] M. Valorani, H. N. Najm, and D. A. Goussis. CSP analysis of a transient flame-vortex interaction: time scales and manifolds. Combust.
Flame, 134:35–53, 2003.
[70] H. N. Najm, M. Valorani, D. A. Goussis, and J. Prager. Analysis of
methane-air edge flame structure. Combust. Theor. Model., 14(2):257–
294, 2010.
[71] J. Prager, H. N. Najm, M. Valorani, and D. A. Goussis. Structure of
n-heptane/air triple flames in partially-premixed mixing layers. Combust. Flame, 158:2128–2144, 2011.
[72] T. F. Lu, C. S. Yoo, J. H. Chen, and C. K. Law. Three-dimensional
direct numerical simulation of a turbulent lifted hydrogen jet flame
155
Bibliography
in heated coflow: a chemical explosive mode analysis. J. Fluid Mech.,
652:45–64, 2010.
[73] Z. Luo, C. S. Yoo, E. S. Richardson, J. H. Chen, C. K. Law, and T. F.
Lu. Chemical explosive mode analysis for a turbulent lifted ethylene
jet flame in highly-heated coflow. Combust. Flame, 159:265–274, 2012.
[74] M. Valorani, D. A. Goussis, F. Creta, and H. N. Najm. High order
corrections in the approximation of inertial manifolds and the construction of simplified problems with the CSP method. J. Comput.
Phys., 209:754–786, 2005.
[75] T. Lu, Y. Ju, and C. K. Law. Complex CSP chemistry reduction and
anlysis. Combust. Flame, 126:1445–1455, 2001.
[76] I. Glassman and R. A. Yetter. Combustion. Academic Press, Fourth
Edition, 2008.
[77] J. R. Bond, P. Gray, and J. F. Griffiths. Oscillations, glow and ignition
in carbon monoxide oxidation I. Glow and ignition in a closed reaction
vessel and the effect of added hydrogen . Proc. R. Soc. Lond. A,
375:43–64, 1981.
[78] J. R. Bond, P. Gray, J. F. Griffiths, and S. K. Scott. Oscillations, glow
and ignition in carbon monoxide oxidation II. Oscillations in the gasphase reaction in a closed system. Proc. R. Soc. Lond. A, 381:293–314,
1982.
[79] P. Gray, J. F. Griffiths, and S. K. Scott. Oscillations, glow and ignition in carbon monoxide oxidation in an open system I. Experimental
studies of the ignition diagram and the effects of added hydrogen.
Proc. R. Soc. Lond. A, 397:21–44, 1985.
[80] P. Gray, J. F. Griffiths, and S. K. Scott. Oscillations, glow and ignition
in carbon monoxide oxidation in an open system II. Theory of the
oscillatory ignition limit in the c.s.t.r. Proc. R. Soc. Lond. A, 402:187–
204, 1985.
[81] J. F. Griffiths and A. F. Sykes. Numerical interpretation of oscillatory
glow and ignition during carbon monoxide oxidation in a well-stirred
flow reactor. J. Chem. Soc. Faraday T., 85:3059–3069, 1989.
156
Bibliography
[82] B. R. Johnson, J. F. Griffiths, and S. K. Scott. Oscillations and chaos
in CO+O2 combustion. Chaos, 1(4):387–395, 1991.
[83] D. L. Baulch, J. F. Griffiths, A. J. Pappin, and A. F. Sykes. Stationarystate and oscillatory combustion of hydrogen in a well-stirred flow
reactor. Combust. Flame, 73:163–185, 1988.
[84] B. R. Johnson, J. F. Griffiths, and S. K. Scott. Characterisation of
oscillations in the H2 + O2 reaction in a continuous-flow reactor. J.
Chem. Soc. Faraday T., 87(4):523–533, 1991.
[85] B. R. Johnson, S. K. Scott, and A. S. Tomlin. Modelling complex
oscillations for the H2 + O2 reaction in an open system. J. Chem.
Soc. Faraday T., 87(16):2539–2548, 1991.
[86] B. R. Johnson and S. K. Scott. Complex and non-periodic oscillations
in the H2 + O2 reaction. J. Chem. Soc. Faraday T., 93(17):2997–3004,
1997.
[87] G. Darrieus. Propagation d’un front de flamme, Essai de Theorie
des vitesses arnoles de Deflagration par developpment spontanne de
la Turbulence. Sixth International Congress of Applied Mathematics,
1946.
[88] L. Landau. On the theory of slow combustion. Acta Physicochim.
URS., 19:77–85, 1945.
[89] G. I. Barenblatt, Y. B. Zeldovich, and A. G. Istratov. On the diffusionthermal stability of a laminar flame. J. Appl. Mech. Tech., 4:21–26,
1962.
[90] G. Joulin and P. Clavin. Linear stability analysis of nonadiabatic
flames - diffusional-thermal model. Combust. Flame, 35:139–153,
1979.
[91] G. Pizza, C. E. Frouzakis, J. Mantzaras, A. G. Tomboulides, and
K. Boulouchos. Dynamics of premixed hydrogen/air flames in microchannels. Combust. Flame, 152:433–450, 2008.
[92] G. Pizza, C. E. Frouzakis, J. Mantzaras, A. G. Tomboulides, and
K. Boulouchos. Dynamics of premixed hydrogen/air flames in
mesoscale channels. Combust. Flame, 155:2–20, 2008.
157
Bibliography
[93] G. Pizza, C. E. Frouzakis, and J. Mantzaras. Chaotic dynamics in
premixed hydrogen/air channel flow combustion. Combust. Theor.
Model., 16(2):275–299, 2012.
[94] V. N. Kurdyumov, G. Pizza, C. E. Frouzakis, and J. Mantzaras. Dynamics of premixed flames in a narrow channel with a step-wise wall
temperature. Combust. Flame, 156:2190–2200, 2009.
[95] V. N. Kurdyumov. Lewis number effect on the propagation of premixed flames in narrow adiabatic channels: Symmetric and nonsymmetric flames and their linear stability analysis. Combust. Flame,
158:1307–1317, 2011.
[96] K. Maruta, T. Kataoka, N. I. Kim, S. Minaev, and R. Fursenko. Characteristics of combustion in a narrow channel with a temperature gradient. Proc. Combust. Inst., 30:2429–2436, 2005.
[97] T. L. Jackson, J. Buckmaster, Z. Lu, D. C. Kyritsis, and L. Massa.
Flames in narrow circular tubes. Proc. Combust. Inst., 31:955–962,
2007.
[98] F. Richecoeur and D. C. Kyritsis. Experimental study of flame stabilization in low Reynolds and Dean number flows in curved mesoscale
ducts. Proc. Combust. Inst., 30:2419–2427, 2005.
[99] C. Cui, M. Matalon, and T. L. Jackson. Pulsating mode of flame
propagation in two-dimensional channels. AIAA J., 43:1284–1292,
2005.
[100] J. Daou and M. Matalon. Influence of conductive heat losses on
the propagation of premixed flames in channels. Combust. Flame,
128:321–339, 2002.
[101] S. Chakraborty, A. Mukhopadhyay, and S. Sen. Interaction of Lewis
number and heat loss effects for a laminar premixed flame propagating
in a channel. Int. J. Therm Sci., 47:84–92, 2008.
[102] A. E Lutz, R. J. Kee, J. F. Grcar, and F. M. Rupley. Oppdiff: A
Fortran program for computing opposed-flow diffusion flames, Report
No. SAND96-8243. Technical report, Sandia National Laboratories,
1997.
158
Bibliography
[103] T. Poinsot and D. Veynante. Theoretical and numerical combustion.
R.T. Edwards, Inc., 2005.
[104] J. Li, Z. Zhao, A. Kazakov, M. Chaos, F. L. Dryer, and J. J. Scire. A
comprehensive kinetic mechanism for CO, CH2 O, and CH3 OH combustion. Int. J. Chem. Kinet., 39(3):109–136, 2007.
[105] M. Reinke, J. Mantzaras, R. Schaeren, R. Bombach, A. Inauen, and
S. Schenker. High-pressure catalytic combustion of methane over platinum: in-situ experiments and detailed numerical predictions. Combust. Flame, 136:217–240, 2004.
[106] J. Mantzaras, R. Bombach, and R. Schaeren. Hetero-/homogeneous
combustion of hydrogen/air mixtures over platinum at pressures up
to 10 bar. Proc. Combust. Inst., 32:1947–1955, 2009.
[107] C. S. Panoutsos, Y. Hardalupas, and A. M. K. P. Taylor. Numerical evaluation of equivalence ratio measurement using OH∗ and
CH∗ chemiluminescence in premixed and non-premixed methane–air
flames. Combust. Flame, 156:273–291, 2009.
[108] C. K. Law, G. Jomaas, and J. K. Bechtold. Cellular instabilities of
expanding hydrogen/propane spherical flames at elevated pressures:
theory and experiment. Proc. Combust. Inst., 30:159–167, 2005.
[109] J. K. Bechtold, C. Cui, and M. Matalon. The role of radiative heat
losses in self-extinguishing and self-wrinkling flames. Proc. Combust.
Inst., 30:177–184, 2005.
[110] J. Mantzaras and P. Benz. An asymptotic and numerical investigation of homogeneous ignition in catalytically stabilized channel flow
combustion. Combust. Flame, 119:455–472, 1999.
[111] J. Mantzaras and C. Appel. Effects of finite rate heterogenous kinetics on homogeneous ignition in catalytically stabilized channel flow
combustion. Combust. Flame, 130:336–351, 2002.
[112] S. Karagiannidis and J. Mantzaras. Numerical investigation on the
start-up of methane-fueled catalytic microreactors. Combust. Flame,
157:1400–1413, 2010.
[113] M. Reinke, J. Mantzaras, R. Bombach, S. Schenker, and A. Inauen.
Gas phase chemistry in catalytic combustion of methane/air mixtures
159
Bibliography
over platinum at pressures of 1 bar to 16 bar. Combust. Flame,
141:448–468, 2005.
[114] C. Appel, J. Mantzaras, R. Schaeren, R. Bombach, and A. Inauen.
Turbulent catalytically stabilized combustion of hydrogen/air mixtures in entry channel flows. Combust. Flame, 140:70–92, 2005.
[115] P. Canu and S. Vecchi. CFD simulation of reactive flows: catalytic
combustion in a monolith. AIChe J., (48):2921–2935, 2002.
[116] A. Di Benedetto, V. Di Sarli, and G. Russo. Effect of geometry on
the thermal behavior of catalytic micro-combustors. Catal. Today,
155:116–122, 2010.
[117] R. E. Hayes and S. T. Kolaczkowski. Mass and heat-transfer effects
in catalytic monolith reactors. Chem. Eng. Sci., 49:3578–3599, 1994.
[118] A. E. Cerkanowicz, R. B. Cole, and J. G. Stevens. Catalytic combustion modeling - Comparison with experimental data. J. Eng. Power-T
ASME, 99:593–600, 1977.
[119] S. T. Lee and R. Aris. On the effects of radiative heat transfer in
monoliths. Chem. Eng. Sci., 37:827–837, 1977.
[120] G. Groppi, A. Belloli, E. Tronconi, and P. Forzatti. A comparison
of lumped and distributed models for monolith catalytic combustors.
Chem. Eng. Sci., 50:2705–2715, 1995.
[121] S. T. Kolaczkowski and S. Serbetcioglu. Development of combustion
catalysts for monolith reactors: a consideration of transport limitations. Appl. Catal. A-Gen., 138:199–214, 1996.
[122] R. E. Hayes and S. T. Kolaczkowski. A study of Nusselt and Sherwood
numbers in a monolith reactor. Catal. Today, 47:295–303, 1999.
[123] L. L. Raja, R. J. Kee, O. Deutschmann, J. Warnatz, and L. D.
Schmidt. A critical evaluation of Navier-Stokes, boundary layer, and
plug-flow models of the flow and chemistry in a catalytic-combustion
monolith. Catal. Today, 59:47–60, 2000.
[124] J. Mantzaras, C. Appel, and P. Benz. Catalytic combustion of
methane/air mixtures over platinum: homogeneous ignition distances
in channel flow configurations. Proc. Combust. Inst., 28:1349–1357,
2000.
160
Bibliography
[125] O. Deutschmann and L. D. Schmidt. Modeling the partial oxidation
of methane in a short-contact-time reactor. AIChe J., 44:2465–2477,
1998.
[126] O. Deutschmann, L. I. Maier, U. Riedel, A. H. Stroemman, and R. W.
Dibble. Hydrogen assisted catalytic combustion of methane on platinum. Catal. Today, 59:141–150, 2000.
[127] M. Reinke, J. Mantzaras, R. Schaeren, R. Bombach, and W. Kreutner.
Homogeneous ignition in high-pressure combustion of methane/air
over platinum: comparison of measurements and detailed numerical
predictions. Proc. Combust. Inst., 29:1021–1029, 2002.
[128] M. Reinke, J. Mantzaras, R. Schaeren, R. Bombach, A. Inauen, and
S. Schenker. Homogeneous ignition of CH4 /air and H2 O- and CO2 diluted mixtures over platinum: an experimental and numerical investigation at pressures up to 16 bar . Proc. Combust. Inst., 30:2519–2527,
2004.
[129] C. Appel, J. Mantzaras, R. Schaeren, R. Bombach, A. Inauen, N. Tylli,
M. Wolf, T. Griffin, D. Winkler, and R. Carroni. Partial catalytic oxidation of methane to synthesis gas over rhodium: in situ Raman measurements and detailed simulations. Proc. Combust. Inst., 30:2509–
2517, 2004.
[130] L. I. Maier, M. Hartmann, S. Tischer, and O. Deutschmann. Interaction of heterogeneous and homogeneous kinetics with mass and
heat transfer in catalytic reforming of logistic fuels. Combust. Flame,
158:796–808, 2011.
[131] S. Mazumder and M. Grimm. Numerical investigation of wall heat
conduction effects on catalytic combustion in split and continuous
monolith tubes. Comput. Chem. Eng., 32:552–560, 2008.
[132] N. Sinha, C. Bruno, and F. V. Bracco. Two-dimensional, transient
catalytic combustion of CO-air on platinum. Physicochem. Hydrodyn.,
6:373–391, 1985.
[133] J. S. Tien. Transient catalytic combustor model. Combust. Sci. Technol., 26:65–75, 1981.
[134] R. E. Hayes, S. T. Kolaczkowski, W. J. Thomas, and J. Titiloye.
Transient experiments and modeling of the catalytic combustion of
161
Bibliography
methane in a monolith reactor. Ind. Eng. Chem. Res., 35:406–414,
1996.
[135] A. Schneider, J. Mantzaras, and S. Eriksson. Ignition and extinction
in catalytic partial oxidation of methane-oxygen mixtures with large
H2 O and CO2 dilution. Combust. Sci. Technol., 180:89–126, 2008.
[136] J. Mantzaras. Catalytic combustion of syngas. Combust. Sci. Technol.,
180:1137–1168, 2008.
[137] R. Wanker, H. T. Raupenstrauch, and G. Staudinger. A fully distributed model for the simulation of a catalytic combustor. Chem.
Eng. Sci., 55:4709–4718, 2000.
[138] O. Deutschmann, F. Behrendt, and J. Warnatz. Formal treatment
of catalytic combustion and catalytic conversion of methane. Catal.
Today, 46:155–163, 1998.
[139] L. L. Raja, R. J. Kee, and L. R. Petzold. Simulation of the transient,
compressible, gas-dynamic behavior of catalytic combustion ignition
in stagnation flows. Proc. Combust. Inst., 27:2249–2257, 1998.
[140] K. Maruta, K. Takeda, J. Ahn, K. Borer, L. Sitzki, P. D. Ronney,
and O. Deutschmann. Extinction limits of catalytic combustion in
microchannels. Proc. Combust. Inst., 29:957–963, 2002.
[141] Y. S. Touloukian, R. W. Powel, C. Y. Ho, and P. G. Klemens. Thermophysical properties of matter. Plenum, New York, 1970.
[142] J. Li, Z. Zhao, A. Kazakov, and F. L. Dryer. An updated comprehensive kinetic model of hydrogen combustion. Int. J. Chem. Kinet.,
36:566–575, 2004.
[143] Y. Ghermay, J. Mantzaras, R. Bombach, and K. Boulouchos. Homogeneous combustion of fuel lean H2 /O2 /N2 mixtures over platinum
at elevated pressures and preheats. Combust. Flame, 158:1491–1506,
2011.
[144] J. Prager, H. N. Najm, M. Valorani, and D. A. Goussis. Skeletal mechanism generation with CSP and validation for premixed n-heptane
flames. Proc. Combust. Inst., 32:509–517, 2009.
162
Bibliography
[145] H. J. Curran, P. Gaffuri, W. J. Pitz, and C. K. Westbrook. A comprehensive modeling study of n-heptane oxidation. Combust. Flame,
114:149–177, 1998.
[146] R. W. Schefer. Catalyzed combustion of H2 /air mixtures in a flat plate
boundary layer: II. Numerical model . Combust. Flame, 45:171–190,
1982.
[147] F. Lucci, C. E. Frouzakis, and J. Mantzaras. Three-dimensional direct
numerical simulation of turbulent channel flow catalytic combustion
of hydrogen over platinum. Proc. Combust. Inst., 34:2295–2302, 2013.
[148] M. Schultze and J. Mantzaras. Hetero-/homogeneous combustion of
hydrogen/air mixtures over platinum: fuel-lean versus fuel-rich combustion modes. Int. J. Hydrogen Energ., 38:10654–10670, 2013.
[149] M. Fassihi, V. P. Zhdanov, M. Rinnemo, K. E. Keck, and B. Kasemo.
A theoretical and experimental study of catalytic ignition in the
hydorgen-oxygen reaction on platinum. J. Catal., 141:438–452, 1993.
[150] M. Rinnemo, O. Deutschmann, F. Behrendt, and B. Kasemo. Experimental and numerical investigation of the catalytic ignition of mixtures of hydrogen and oxygen on platinum. Combust. Flame, 111:312–
326, 1997.
[151] O. Deutschmann, R. Schmidt, F. Behrendt, and J. Warnatz. Numerical modeling of catalytic ignition. Proc. Combust. Inst., 26:1747–1754,
1996.
[152] A. L. Boehman. Radiation heat transfer in catalytic monoliths. AIChe
J., 44:24745–2755, 1998.
[153] D. R. Reynolds.
ARKode Example Documentation, Southern
Methodist University Center for Scientific Computation technical report. Technical report, 2013.
[154] A. C. Hindmarsh. The PVODE and IDA algorithms, LLNL technical
report UCRL-ID-141558. Technical report, 2000.
163

1999 冬学期 流体工学特論 -非定常流れの解析と計測-

The state and outlook of LPG in Japan

EU Directive - Karcher KäN 050.032 56 KB

Combustion and Energy - Flame Selection System

The load history shown in the Fig. 1 is repeatedly applied to

CO - Tekna

Promotion Flyer DocPad Final - Fuchs International Technology

Download our Arreis FR Brochure

PCB TERMINAL BLOCKS