@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Mechanical Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Huang , Jian"@en ; dcterms:issued "2010-01-16T17:02:17Z"@en, "2006"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description "This thesis focuses on the study of natural gas combustion under engine relevant conditions. The work begins with the development of a detailed chemical kinetic mechanism that represents the ignition characteristics of methane with various minor additives over a wider range of operating conditions than previously existing mechanisms. The mechanism includes a NOx submechanism selected from the literature that yields good agreement with experimental data in various methane/air combustion systems. The excessive computational load associated with detailed chemistry is alleviated using a trajectory generated low-dimensional manifold (TGLDM) method. The TGLDMs generated in this work provide a satisfactory approximation of calculation using detailed chemistry in various methane/air reaction systems with a significant reduction of the computational cost. An innovative combustion model for simulating turbulent diffusion flames is presented at the end of this thesis. The model employs the Conditional Source-term Estimation method for the closure of the chemical source term. It obtains production/consumption rates of reaction scalars through TGLDMs generated with the new reaction mechanism. The model was used to simulate ignition and combustion of transient turbulent methane jets under engine-relevant conditions; it has achieved encouraging results in comparison with the experimental data from this work as well as the literature."@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/18215?expand=metadata"@en ; skos:note "NATURAL GAS COMBUSTION UNDER ENGINE-RELEVANT CONDITIONS by J I A N H U A N G B.Sc. Shanghai Maritime University, 1996 M.A.Sc. The University of British Columbia, 2002 A THESIS SUBMITTED I N P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E F A C U L T Y OF G R A D U A T E STUDIES (Mechanical Engineering) T H E UNIVERSITY OF BRITISH C O L U M B I A March 2006 © Jian Huang, 2006 Abstract This thesis focuses on the study of natural gas combustion under eng-ine relevant conditions. The work begins with the development of a de-tailed chemical kinetic mechanism that represents the ignition charac-teristics of methane with various minor additives over a wider range of operating conditions than previously existing mechanisms. The mech-anism includes a N O x submechanism selected from the literature that yields good agreement with experimental data in various methane/air combustion systems. The excessive computational load associated with detailed chemi-stry is alleviated using a trajectory generated low-dimensional mani-fold (TGLDM) method. The T G L D M s generated in this work provide a satisfactory approximation of calculation using detailed chemistry in various methane/air reaction systems with a significant reduction of the computational cost. A n innovative combustion model for simulating turbulent diffu-sion flames is presented at the end of this thesis. The model employs the Conditional Source-term Estimation method for the closure of the chemical source term. It obtains production/consumption rates of re-action scalars through T G L D M s generated with the new reaction me-chanism. The model was used to simulate ignition and combustion of transient turbulent methane jets under engine-relevant conditions; it has achieved encouraging results in comparison with the experimen-tal data from this work as wel l as the literature. i i Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Nomenclature and Acronyms x Acknowledgement xii Dedication xiv Chapter 1. Introduction and Thesis Outline 1 1.1 Introduction 1 1.2 Outline 2 Chapter 2. Natural Gas Combustion - A Review 5 2.1 Introduction 5 2.2 Reaction Kinetics for Natural Gas Combustion 5 2.2.1 Experimental Studies of Methane Ignition 6 2.2.2 Experimental Studies of Natural Gas Ignition 8 2.2.3 Ignition Chemistry of Methane 10 2.2.4 Ignition Chemistry of Natural Gas 16 2.2.5 Ignition Studies of Methane with Hydrogen 17 2.2.6 Species Concentration Profile 19 2.2.7 Laminar Premixed Flame 21 2.3 NOx Chemistry 24 2.3.1 Thermal NO 25 2.3.2 Prompt NO 27 2.3.3 NO Formation from Fuel-Bound Nitrogen 28 2.3.4 Formation of N0 2 30 2.3.5 In-flame NOx Reduction Strategies 30 2.4 Analysis and Reduction of Detailed Chemistry 33 2.4.1 Sensitivity Analysis 34 2.4.2 Reaction Flow Analysis 35 2.4.3 Chemical Time Scale and Eigenvalue Analysis 37 iii Table of Contents 2.5 Modeling of Turbulent Reactive Flow 52 2.5.1 Laminar Flamelet Model 54 2.5.2 Conditional Moment Closure 59 2.5.3 Conditional Source Term Estimation 63 2.6 Summary 65 Chapter 3. Ignition of Methane with Higher Alkane Additives 66 3.1 Introduction 66 3.2 Experimental Apparatus 67 3.3 Kinetic Model 71 3.4 Results and Discussion 79 3.4.1 Experimental and Mode l Results 79 3.4.2 Kinetic Analysis—Methane/Ethane Mixtures 86 3.4.3 Kinetic Analysis—Methane/Propane Mixtures 92 3.4.4 Analytical Mode l 93 3.5 Conclusions 112 Chapter 4. Ignition of Methane with Hydrogen Addition 113 4.1 Introduction 113 4.2 Experiments 115 4.3 Kinetic Mode l 116 4.4 Results and Discussion 123 4.5 Conclusions 141 Chapter 5. Evaluation of Detailed Nitrogen Oxides Mechanisms in Methane Flames 143 5.1 Introduction 143 5.2 Experimental Data 144 5.3 NOx Mechanisms 145 5.4 Results and Discussion 146 5.5 Conclusions 153 Chapter 6. Trajectory Generated Low Dimensional Manifold in C H 4 - Air Combustion Systems 154 6.1 Introduction 154 6.2 Construction of T G L D M 160 6.3 Results and Discussion 161 6.3.1 Unstrained Premixed Laminar Flame 163 6.3.2 Laminar Flamelet 166 6.3.3 Perfectly Stirred Reactor 171 6.4 Conclusions 175 iv Table of Contents Chapter 7. Simulation of Turbulent Reactive Methane Jets Un-der Engine Relevant Conditions 176 7.1 Introduction 176 7.2 CSE-TGLDM Combustion Model 177 7.2.1 Conditional Source-term Estimation 177 7.2.2 Trajectory Generated Low-Dimensional Manifold . 180 7.3 Experiments • • 181 7.4 Model Formulation and Validation 185 7.4.1 CFD Model Formulation 185 7.4.2 Combustion Model Formulation 186 7.5 Results and Discussion 189 7.5.1 Non-reactive Jet Validation 189 7.5.2 Shock-tube Ignition 194 7.5.3 NO Prediction 198 7.5.4 Combustion Bomb Ignition 208 7.5.5 Computational Time 212 7.6 Conclusions 213 Chapter 8. Conclusions and Future Work 214 8.1 Summary and Conclusions 214 8.2 Future Work 222 Bibliography 224 Appendix A. Derivations of Key Equations in CSP 241 Appendix B. Eigen Decomposition and Schur Decomposition 243 Appendix C. Solution of Ricatti's Differential Equation 245 Appendix D. Flux Corrected Transport 248 Appendix E. Numerical Schlieren 250 v List of Tables 2.1 Experimental conditions and empirical coefficients for methane ignition from the literature 7 2.2 Mechanisms for natural gas combustion in the literature 14 3.1 Modified and extended reactions in the model 73 3.2 Molar fractions of main constituents in the test mixtures (0= 1 for all mixtures, balanced by air (79%N 2, 21%02)) 80 3.3 Experimental conditions and ignition delay results 81 3.4 The skeletal mechanism used in the analytical model 101 4.1 Molar fractions of main constituents in the test mixtures (0= 1 for all mixtures, balanced by air (79%N 2, 21%02)) 116 4.2 Modified and extended reactions in the new model (based on GRI-Mech 1.2) 119 4.3 Experimental conditions and ignition delay results 125 4.4 Correlation coefficients for Eq. 4.2 127 5.1 Selected NOx mechanisms from the literature 145 5.2 Comparison of key reaction rates at 1600 and 2000K 152 7.1 Experimental conditions of this work 184 7.2 Operating conditions for the experiment of Naber et al 211 vi List of Figures 2.1 Main reaction paths during methane ignition 12 2.2 Main reaction paths of prompt NO in methane flame 29 2.3 Sensitivity of ignition delay in methane/air mixture with respect to elementary reaction rates 36 2.4 Histogram of eigenvalues of a combustion system 42 2.5 Convergence of reaction trajectories onto an ILDM 43 2.6 Comparison of ILDM and TGLDM for a hypothetical reaction system 50 2.7 Coordinate transformation in laminar flamelet model 55 3.1 Schematic of the shock tube facility 67 3.2 Determination of experimental time 69 3.3 Ignition delay in methane/ethane mixtures 83 3.4 Ignition delay in methane/propane mixtures 84 3.5 Ignition delay in methane/ethane/propane mixtures 85 3.6 Sensitivity of ignition delay with respect to rates of elementary reactions for the 3.7% C 2 H 6 /96.3%CH 4 mixture 87 3.7 Cumulative contribution to OH formation 90 3.8 Main reaction path for methane oxidation with the presence of minor ethane and propane additives 94 3.9 Definition of ignition delay in the analytical model 99 3.10 Comparison of ignition delay calculated from the model with that from the detailed chemistry 103 3.11 Species concentrations predicted by the analytical model com-pared with the calculation from the detailed kinetic model for the stoichiometric methane/air mixture 104 3.12 Ignition delay with ethane/propane additives 109 3.13 Species concentrations predicted by the analytical model com-pared with the calculation from the detailed kinetic model for mixture No.3 110 3.14 Species concentrations predicted by the analytical model com-pared with the calculation from the detailed kinetic model for mixture No.6 I l l 4.1 Calculated ignition delay of hydrogen/oxygen/argon mixtures with experimental data from the literature 122 4.2 Predicted methane ignition delay using the new mechanism 124 4.3 Measured and calculated ignition delay of test mixtures — 126 4.4 Comparison of measured ignition delay with formula 4.2 .. 128 vn List of Figures 4.5 Normalized sensitivity of ignition delay with respect to individ-ual reaction rate for methane and methane/hydrogen mixtures at 40 bar. 130 4.6 H concentration during the induction period 134 4.7 Integrated contribution of major H-generation reactions — 135 4.8 Integrated contribution of major H-consumption reactions . 136 4.9 Comparison of rates of progress of recombination reactions (sum of R33-R38) and the chain branching reaction (R39) 137 4.10 Integral reaction flow analysis showing main reaction pathway for mixture #2 at 1300k and 40 bar 139 4.11 Calculated ignition limit of hydrogen for a given ignition delay in this work 140 5.1 NOx prediction in simulated PSR 147 5.2 N 2 0 prediction in simulated PSR 148 5.3 NO prediction in a laminar counterflow diffusion flame 149 5.4 NO prediction in a laminar counterflow diffusion flame 150 6.1 Comparison of TGLDM with ILDM 159 6.2 Delaunay-triangulated TGLDM in YCo2 ~ YH2O plane 162 6.3 Comparison of calculated flame structure of a one-dimensional premixed laminar methane/air flame 165 6.4 Laminar flame velocity of a stoichiometric methane/air mixture calculated using TGLDM method 166 6.5 Comparison of calculated temperature and species profile of a steady state methane/air flamelet at 30 bar in the mixture fraction space. Xst = 100s-1 169 6.6 Comparison of calculated temperature and species profile of a steady state methane/air flamelet at 30 bar in the mixture fraction space. Xst = 300s-1 170 6.7 Calculated temperature and species mass fraction history in a PSR 173 6.8 Calculated temperature and species mass fraction history during a quenching process in a PSR 174 7.1 Schematic of the shock tube facility 182 7.2 Structure of CSE-TGLDM method in the simulation 188 7.3 Comparison of a non-reactive methane jet, P ratio=3 191 7.4 Comparison of a non-reactive methane jet, P ratio=5 192 7.5 Comparison of simulated and experimental Schlieren pictures for a methane jet at 200/is from the start of injection 193 7.6 Definition of ignition using the scalar history 195 viii List of Figures 7.7 Comparison of ignition delay calculated using the C S E - T G L D M combustion model with experimental measurements 196 7.8 Comparison of axial locations of initial ignition kernel from the model with those from the experiment 197 7.9 Profiles of reaction scalars, mixture fraction and its variance in the computational domain shortly after ignition 199 7.10 Profiles of reaction scalars, mixture fraction and its variance in the computational domain before the injection stops 200 7.11 Profiles of reaction scalars, mixture fraction and its variance in the computational domain shortly after the injection stops . 201 7.12 Profiles of reaction scalars, mixture fraction and its variance in the computational domain after the injection stops 202 7.13 Profiles of reaction scalars, mixture fraction and its variance in the computational domain during ignition 203 7.14 Profiles of reaction scalars, mixture fraction and its variance in the computational domain shortly after ignition 204 7.15 Correction of N O mass fraction 208 7.16 Comparison of predicted NOx mass fraction and experimental measurements 209 7.17 Simulation results for the experiment of Nabef et al 211 7.18 Comparison of C P U time on C S E - T G L D M module with that on C F D module in a typical simulation 212 ix Nomenclature and Acronyms Symbol Meaning Units D Diff usi vity m2.s~l E Activation Energy cal.mol*1 h Specific Enthalpy J.kg'1 I Identity Matrix J Jacobian Matrix s\"1 k Reaction Rate p Pressure Pa or bar P Probability Density Function P Projection Matrix R Universal Gas Constant J.mol'1 V Eigen(Schur) Basis Matrix S Sensitivity T Temperature K Y Mass Fraction Z Mixture Fraction (Y\\Z — rj) Conditional Average of Y with Con-dition Z — 7] (Z) Mean of Mixture Fraction (Z\"2} Variance of Mixture Fraction CMC Conditional Moment Closure CSE Conditional Source-term Estimation EGR Exhaust Gas Recirculation FCT Flux Corrected Transport ILDM Intrinsically Low Dimensional Man-ifold TGLDM Trajectory-Generated Low Dimensio-nal Manifold x Nomenclature and Acronyms Symbol Meaning Units e Dissipation Rate of Turbulent Kinetic J.kg~l.s~l Energy A Thermal Diffusivity 9 — 1 m .s A Diagonal Matrix of Eigenvalues s-1 P Density kg.m~3 4> Equivalence Ratio

1300 -0.02 -1.20 1.26E-14 32.7 [29] 40-260 <1300 -0.38 -1.31 4.99E-14 19.0 The global activation energy, E, indicates the sensitivity of ignition delay with respect to changing temperature. It can be seen from Ta-ble 2.1 that at relatively high temperature, the experimentally obtained value for E is around 50 kca l /mol , while it reduces significantly to 19 kca l /mol at temperatures below 1300 K. The reduction in the ac-tivation energy with reducing temperature implies that the reactions which are rate-limiting in methane system are different at different temperatures; it also highlights the limitation of the above empirical coefficients, which should not be used beyond the experimental ranges within which they were obtained. While most of the earlier studies focused on ignition at high tem-perature and low pressure, ignition delay data at elevated pressures and moderate temperatures have become more available in the liter-ature recently [28-30]. Petersen et al. [29,30] conducted shock tube ^nits (s(cm 3/mol) x + y) 7 Chapter 2. Natural Gas Combustion - A Review experiments on ignition of methane/air and methane/oxygen/argon mixtures at pressures from 40 to 260 atm and temperatures from 1040 to 1500 K. .The objective of their study was to understand the methane ignition mechanism for ram-propulsion applications so that they cov-ered the equivalence ratios (4>) in the fuel-lean (cp = 0.4) and.fuel-rich (cp > 3.0) regions. Later, Huang et al. [28] reported shock-tube igni-tion results for undiluted methane/air mixtures at pressures from 16 to 40 bar and temperatures from 1000 to 1350 K. The equivalence ratios ranged from slightly lean (cp = 0.7) to slightly rich (

= 1.0 and 0 = 1 . 1 with a value around 39 cm/s; it drops to around 9cm/sat = 0.6 and l l c m / s a t 0 = 1.5 [58]. The classical laminar premixed flame propagation theory based on Zel'dovich's analysis [63] suggests that the laminar flame velocity VL depends on the diffusivity D of mass and energy, and the characteristic time of reaction r . The formula (2.4) 21 Chapter 2. Natural Gas Combustion - A Review indicates a balance between the rate of species transport as a physical process and the rate of chemical reactions. Mathematically vL is the eigenvalue of a set of partial differential equations (PDEs) that govern the steady-state mass and energy conservations of a reaction system. According to the theory of gas diffusion [57], species diffusivity has an inverse first-order dependence on pressure. The pressure dependence of the chemical time scale is related to reaction order n by Thus the laminar flame velocity has an overall pressure dependence of order n/2 — 1. For methane/natural gas combustion, the main reac-tions are mostly second-order so that one can expect the laminar flame velocity to be roughly proportional to P ° 5 . For a single-step reaction with an Arrhenius rate expression, the temperature dependence of the laminar burning velocity would then be where E is the activation energy and Tb is the temperature of the burned gas. Warnatz [62] compared calculated laminar flame velocities using the above correlations for methane/air mixtures with measurements from the literature and reported good agreement over a wide range of temperature and pressure. Using advanced digital computing technology, the contemporary method of solving for the laminar flame velocity has changed signifi-cantly from the classical approaches introduced above. Solutions can now be obtained by solving energy and species mass conservation equations with detailed chemistry [57,62]. For a one-dimensional pre-mixed diffusion flame, the governing equation of species mass conser-r oc P ,1-n (2.5) (2.6) 22 Chapter 2. Natural Gas Combustion - A Review vation is given by where p is the density, Y is the mass fraction, t is time, u is the bulk convection velocity, V is the diffusion velocity, and Co is the rate change of mass fraction due to chemical reactions. The subscript i indicates the ith species. The mixture averaged diffusion velocity Vir which in-corporates both mass diffusion and thermal diffusion velocities, can be obtained from the Fickian formula [57] Vi = ~ A V Xi - v T , (2.8) Xi pYiT where D is the mixture-averaged diffusion coefficient, X is the mole fraction, and DT is the thermal diffusivity at temperature T. The energy conservation equation is given by 1=1 2=1 where Cp is the specific heat at constant pressure, A is the thermal con-ductivity, and Ah°f is enthalpy of formation. The boundary conditions are Yi = YUti, T = TU at x = 0, (2.10) and dYi dT ~ = 0, — = 0 at x oo. (2.11) ox ox The subscript u indicates unburned mixture. The above n+l (n species conservation and one energy) partial differential equations (PDEs) can be solved numerically, and the laminar flame velocity can be obtained from the steady-state solution of the convection velocity u. Detailed species and temperature profiles across a laminar flame, i.e. the struc-ture of flame front, can also be extracted from the solution of the above PDEs. 23 Chapter 2. Natural Gas Combustion - A Review Bernstein et al. [64] performed numerical calculations for the spatial profiles of key intermediate species (0, H , OH, C H , CO, HCO, and C H 3 ) in a premixed methane/air flame at low pressures. The model results are in good agreement with their experimental measurements obtained using two laser-based diagnostic methods. Heard et al. [65] measured the O H , C H , and NO radicals in slightly rich methane/air flames burn-ing at 30, 70, and 120 torr using the Laser-Induced Fluorescence (LIF). Absolute NO and O H concentrations were determined using separate calibration experiments. They performed a simulation of the flame us-ing a numerical model similar to that introduced above. The predicted profiles of the radical concentrations as a function of height above the burner agree well with the experimental results. 2.3 NOx Chemistry Nitrogen oxides (NOx) are the main pollutant emissions from combus-tion processes. Severe environmental and health issues can be caused by an elevated level of NOx in the atmosphere. Nitrogen oxides are critical components of photochemical smog; Nitrogen dioxides cause damage to the mechanisms that protect the human respiratory tract and can increase a person's susceptibility to, and the severity of, res-piratory infections and asthma. Long-term exposure to high levels of nitrogen dioxide can cause chronic lung disease. Because internal com-bustion engines and gas turbines are a significant source of NOx emis-sions, stringent regulations have been imposed by government agen-cies worldwide to control NOx emissions from these devices. Conse-quently, the formation of NOx and its reduction have become and re-main to be a major focus of combustion studies by researchers from both industry and academia. 24 Chapter 2. Natural Gas Combustion - A Review NO and NO2 are the main nitrogen oxides generated in IC engines. Experimental data on NO and N 0 2 formation in methane/air mixtures in stirred reactors or in laminar premixed flames are readily available in the literature [66-69]. Detailed discussion of NOx formation mech-anisms can be found in combustion books [57,62]. In general, for the formation of NO, three major mechanisms have been identified: the thermal NO (Zel'dovich mechanism), the prompt NO (Fenimore mec-hanism) and the fuel-bound nitrogen mechanism. 2.3.1 Thermal NO The Zel 'dovich [63] mechanism for NO formation involves three major steps5 0 + N 2 ^ N O + N , (R21) 0 2 + N ^ NO + O, (R22) and N + OH ^ NO + H . (R23) The rate of reaction R21, the rate-limiting step in the Zel'dovich me-chanism, has a very high activation energy: around 318 kJ /mol [70]. As a result, the Zel 'dovich mechanism is very sensitive to tempera-ture change. That is why the mechanism is also called thermal NO route. As a rule-of-thumb, the contribution from the Zel 'dovich me-chanism becomes significant only when the temperature is above 1800 K [57]. In internal combustion engines, the rate of reaction R21 is sig-nificantly lower than the rate of fuel oxidation and change of engine conditions, which makes the formation of NO a kinetically-controlled 5This is often called the \"extended\" Zel'dovich mechanism. Reactions 21 and 22 were initially suggested by Zel'dovich [63]. Reaction R23 was added by Lavoie et al. [63]. 25 Chapter 2. Natural Gas Combustion - A Review slow process [71]. The rate expression of the formation of NO from the Zel 'dovich mechanism is ^ 1 = kf21[0}[N2] + kf22[N][02] + H3[N][OH] (2.12) - kb21[NO][N] - kb22[NO}[0) - kb23[NO][H], where k denotes reaction rates; the square bracket denotes molar con-centrations; superscript f denotes forward reactions; superscript b de-notes backward reactions. Simplification of the rate expression for NO can be achieved by applying the quasi-steady-state assumption to the nitrogen atom, which leads to g g l = 2 4 [ o i r a i - i w i v w f t i w ) . ( 2 . 1 3 ) dt 'l + k\\l[NO]/(k!n[02\\ + kUOH\\) where K = k21k22/(k2lk22). In an IC engine, the amount of N O formed in the postflame gas is usually much higher than that in the flame front due to the short residence time of the latter [71]. For thermal NOx, it is thus reasonalbe to approximate the concentration of O, 02, OH, H and N 2 using the equilibrium values at the local pressure and temper-ature. The validity of the equilibrium assumption has been proven by Kuo [57] and Warnatz [62]. Warnatz reported that the equilibrium as-sumption yields a reasonable approximation to detailed chemistry for temperatures above 1700 K [62]. For conditions when the concentration of NO is significantly lower than it equilibrium concentration, the frac-tional term on the r.h.s of Eq. 2.13 reduces to unity. Correspondingly, the rate of NO formation reduces to d ^ = 2k{1[0]e[N2]e. (2.14) where subscript e denotes equilibrium properties. 26 Chapter 2. Natural Gas Combustion - A Review 2.3.2 Prompt NO Fenimore [72,73] measured the NO concentration in a flat premixed hydrocarbon flame; he found that the NO concentration close to the unburned mixture side is significantly higher than that predicted by the Zel 'dovich mechanism. He thus postulated a new mechanism for the prompt formation of NO in a flame with the presence of hydrocar-bon radicals. First, nitrogen-containing radicals are formed, via: C H + N 2 =± H C N + N , (R24) C H 2 + N 2 =± H C N + NH, (R25) C H 2 + N 2 =i H 2 C N + N , (R26) C + N 2 =± C N + N . (R27) According to Glarborg et al. [74], reaction R24 is the most important initiation step in the prompt mechanism. The activation energy of this reaction reported in the literature ranges from 57 k j / m o l [75,76] to 87 k j / m o l [77], which is significantly smaller than that of R21 in the thermal mechanism. This is the main reason for the dominance of the prompt NO in low temperature zones in hydrocarbon flames. In lean and slightly rich flames, the hydrocyanic acid (HCN) and nitrogen radical formed in the above reactions are quickly converted to NO via H C N + O =i N C O + H, (R28) N C O + H =± N H + CO, (R29) N H + H =± N + H 2 , (R30) N + OH ^ NO + H. (R31) Mil ler and Bowman [78] reported that for equivalence ratios greater than 1.2, the rate of conversion from H C N to NO decreases due to: 27 Chapter 2. Natural Gas Combustion - A Review 1. recycling of NO to H C N by C, C H and C H 2 radicals, and 2. the shift of direction in reaction N + NO =^ N 2 + 0 : The most important hydrocarbon radical in the formation of prompt NO is C H , whose.concentration is established through C H 3 + X =i C H 2 + H X , (R32) C H 2 + X =^ C H + H X , • (R33) C H + X =i C + H X , (R34) were X denotes H or OH radicals. Figure 2.2 shows the main reaction path for the formation of prompt NO in methane flames. For methane/air premixed flames at ambient pressure, Mil ler and Bowman [78] showed that C H and total fixed ni -trogen (TFN) both peak at 4> = 1-4. Beyond this equivalence ratio, the concentration of C H decreases due to the depletion of H and OH radi-cals, which leads to a reduction in prompt NO. The prompt mechanism contributes most significantly to NO formation at T < 2000K, while the thermal mechanism dominates for T > 2500K [57]. 2.3.3 NO Formation from Fuel-Bound Nitrogen Nitrogen compounds contained in natural gas may be released in a combustion process to form NOx. Studies have shown that in the gas phase, fuel-bound nitrogen (FBN) converts quickly into hydrogen cy-anide (HCN) and ammonia (NH 3) [57], which then proceed along their respective reaction paths to form NO. In general, F B N is not a major issue in natural gas combustion since most natural gas contains very low concentrations of chemically bound nitrogen; it is much more im-portant for coal combustion since the F B N mass fraction can be as high as a few percent in coal [62]. 28 Chapter 2. Natural Gas Combustion - A Review H, OH H, OH H, OH CH3 • CH2 • C H • C Figure 2.2: M a i n reaction paths of prompt NO in methane flame. The paths marked by dashed lines are more important in the fuel-rich re-gion. 29 Chapter 2. Natural Gas Combustion - A Review 2.3.4 Formation of N0 2 In typical engine exhaust gas, the majority of NOx emissions is in the form of NO; however, for conditions where the flame temperature is low (such as in engines with significant fractions of recirculated ex-haust gas), a considerable fraction of NO can be converted to N0 2 thr-ough NO + H02 ^ N0 2 + OH, (R35) N0 2 + H —NO + OH, (R36) N02 + 0=iNO + 02. (R37) Reaction R35 is the dominating reaction in this process because of the high concentration of H0 2 in low temperature flames. Reaction R35 can be enhanced by the migration of H radicals from high-temperature regions to low-temperature ones in a. diffusion flame. The hydrogen radicals can be readily combined with molecular oxygen to form H0 2 in H + 0 2 + M =± H0 2 + M. (R38) Reaction R38 is a third-order reaction, which implies that the conver-sion from NO to N02 has a higher rate at higher pressures, which makes this mechanism more relevant to combustion in IC engines. 2.3.5 In-flame NOx Reduction Strategies Various in-flame NOx reduction technologies have been developed in the past decades for internal combustion engines and gas turbines to address increasing environmental concerns. Technologies such as lean combustion, staged combustion and exhaust gas recirculation (EGR) [79] have been widely adopted in the design of modern combustion de-vices. For example, in direct-injection natural gas engines, recent stud-30 Chapter 2. Natural Gas Combustion - A Review ies by McTaggart-Cowan et al. [5,6] show that up to an 80% reduction of NOx in the engine exhaust can be achieved with a high fraction of EGR compared with baseline results under identical load conditions. While these techniques are mainly based on the concept of thermal NOx reduction by lowering the flame temperature, their kinetic impact should also be considered. Li and Williams [36,80] studied the effects of H 2 0, C0 2 and N 2 addition in a two-stage methane/air counterflow flame. Water was found to be the most effective agent for reducing flame temperature due to its high specific heat, which led to the great-est reduction in thermal NOx. Although C 0 2 and N 2 have similar spe-cific heats, C0 2 is more effective than N 2 in suppressing NOx forma-tion. This can be attributed to the kinetic interaction between C 0 2 and CH through reaction CH + C 0 2 =i CHO + CO. (R39) Similarly, kinetic benefits can be obtained with water addition due to reaction • CH + H 2 0 = i CH 2 0 + H. (R40) The high concentration of C 0 2 or H 2 0 shifts the directions of R39 and R40, which reduces the rate of CH production, and consequently, the prompt NO formation. NOx reburning is another effective technology for reducing NOx emission. Reburning is a chemically complex process in which NO is reduced in a hydrocarbon flame. Experimental studies of NOx re-burning by doping reaction mixtures with NO, N0 2 and N 2 0 in plug-flow reactors or premixed flames have been reported by numerous re-searchers [68,81-83]. The consumption of NO in a reburning process is found to be strongly related to temperature, equivalence ratio and 31 Chapter 2. Natural Gas Combustion - A Review residence time. Fuel-rich mixtures were found to favor NO reduction in most studies. Mil ler and Bowman [78] postulated a mechanism for NOx reburn-ing in hydrocarbon flames. It is believed that during reburning, nitric oxide is converted into hydrogen cyanide (HCN) and cyano (CN) by hydrocarbon radicals through C H 2 + NO ^ H C N O + H, (R41) H C N O + H ^ H C N + OH, (R42) C H + N O ^ H C N + 0, (R43) C + NO ^ C N + O. (R44) Their numerical study shows that the concentrations of C H 2 and C H radicals peak in slightly fuel-rich flames so that the reduction efficiency is the highest in this region. This is consistent with the experimental observations introduced above. For flame temperatures between 800 and 1500 K , Glarborg et al. [84] suggested that the most important reactions in the reburning mechan-ism are those related to H C C O + NO and C H 3 + NO. Specifically reac-tions C H 3 + NO(+M) ^ CH 3 NO(+M) , (R45) C H 3 + NO ^ H 2 C N + OH, (R46) and C H 3 + NO ^ H C N + H 2 0 (R47) are dominant reactions for NO reburning in C H 4 / a i r flames. Reaction H C C O + NO =± H C N O + CO (R48) 32 Chapter 2. Natural Gas Combustion - A Review dominates NO reburning in natural gas/C2 flames because of the rel-atively high concentration of acetylene (C2H2) in these flames. Acety-lene is a main precursor for the formation of H C C O radicals. Hydrogen cyanide formed in the above reactions is partially recirculated back to N2 via reactions O + HCN =^ NCO + H, (R49) NCO + H=iNH + CO, (R50) NH + H =i N + H2, (R51) N + NO=±N 2 + H. (R52) Although the NOx-formation mechanisms in hydrocarbon flames have been studied intensively, the uncertainties in the rate constants of some important elementary reactions remain high (e.g. HCCO + NO reaction [84]). Much research work is still necessary for us to fully un-derstand the formation of NOx in practical combustion systems. 2.4 Analysis and Reduction of Detailed Chemistry As shown in the reference cited in Table 2.2, a typical detailed reaction mechanism for natural gas combustion may contain dozens of species and hundreds of reactions. Depending on the combustion system, the thermal and kinetic contributions from some reactions are more sig-nificant than others. Numerical analyses help to identify the most sig-nificant reactions. They also serves as valuable tools for the reduction of detailed reaction mechanisms, whose implementation may be oth-erwise beyond the reach of the available computational resource. Sen-sitivity analysis, reaction flow analysis and eigenvalue-based analysis are three categories of analysis methods which are extensively used in 33 Chapter 2. Natural Gas Combustion - A Review studying reaction mechanisms. Sensitivity analyses identify the rate-limiting reaction steps. Reaction flow analyses investigate the main re-action paths and integrated contribution of elementary reactions in any time period of interest. Eigenvalue/eigenvector analyses determine the characteristic time scales and directions of the chemical reactions [62]. 2.4.1 Sensitivity Analysis For a reaction system with defined initial and boundary conditions, the rate at which the system changes its state is a function of the species concentrations and the rates of elementary reactions involved. Mathe-matically, for a system with n species and m reactions, the above state-ment can be expressed as H = f(c1...cn,k1...km), (2.15) where £ denotes a state variable of the reaction system (e.g. tempera-ture, pressure, individual species mass fraction, etc.), c is the species concentration and k is the reaction rate. Some of the elementary re-actions have nearly no effect on the solution of £, while others are more significant in determining the state of the system [62]. The rate-determining or rate-limiting reaction steps can be identified by exam-ining the dependence of the solution of £ on the parameter ki by means of two definitions. One, the absolute sensitivity is *'=!!•• • < 2- 1 6 ) in which the subscript j denotes the jth reaction. The second definition is the relative sensitivity ^ ^ _ d l n £ 3 ~ fdkj'dhxkj- [ ] 34 Chapter 2. Natural Gas Combustion - A Review The differential equation governing the change of Sj can be derived by differentiating Eq. 2.15 with respect to kj d di d di dt dkj dkj dt dSf _df (^L^\\ (2.18) dt dkj 4—f \\dci ^ i v , . i=i j in which the superscript i denotes the zth species. If one is interested in the sensitivity of concentration of species m with respect to the rate constant kj, Eq. 2.18 reduces to dS, 3 (2.19) dt VdkJ ^ WdcJ •> i=i 1 Equation 2.18 or 2.19 can be integrated numerically to obtain the value of the sensitivity as the reaction system evolves. A n approximate evaluation of Eq. 2.18, which has the merit of sim-plicity in application, is given by b> ~ J dkj ~ J (n-l/n)kj ' U - 2 U ) where n is a multiplication factor, which is typically given a value of two in the literature [22,29]. Equation 2.20 is sometimes referred as a \"brute-force\" sensitivity analysis. Figure 2.3 presents an example of brute-force sensitivity analysis [28]. It shows that at the higher tem-perature, the recombination of two methyl radicals to form an ethane molecule is a main rate-limiting step in the induction period; while at the lower temperature, reactions that are related to the formation and consumption of O H radicals are more important in determining the ignition delay. 2.4.2 Reaction Flow Analysis Sensitivity analyses provide important information regarding the rate-limiting steps of a reaction system; however, the information is usually 35 Chapter 2. Natural Gas Combustion - A Review CH302-lCH4c->CH302H-iCH3 (R1&5) CH302-tCH3<->CH30-tCH30 (R184) CH3+O2<->CH302 (R179) CH3+CH20<»HCO-tCH4 (R161) 2CH3(+M)<->C2H6(+M) (R158) H02-tCH4<«>CH3+H2O2 (R157) CH3-tO2<->0H4CH20 (RIM) CH3+O2<->0+CH30 (RI JJ) KO2+CH2O<»>HC0+H202 (Rt21) K02-(CH3<->OH+CH30 (RI 19) H02+CH3<->02-lCH4 (RI 1ST) 2H02<->02+H202 (RI 16) OH-1CH20<->HCO+H20 (RI01) OH+CH4<->CH3+H20 (R98) OH+H02<->02+H20 (R87) 20H(+M)<->H202(+M) (R8J) H+CH4<->CH3+H2 (RJ3) H-H32<->CHOH (R3S) 02+CH20<->H02-tHCO (R32) -0.10 0.00 O.10 020 Normalized sensitivity Figure 2.3: Sensitivity of ignition delay in methane/air mixture with respect to elementary reaction rates insufficient to determine the main reaction path the system takes to evolve from one state to another. Such information can be obtained from reaction flow analyses. For a reaction process involving n species and m elementary reactions, the rate of change of species mass fraction at any moment t can be written in matrix form as G i = \\W = SijFj> ( 2 - 2 1 ) where y is the rate array of mass fractions of n species; i = l...n; S is a n x m matrix of stoichiometric coefficients; j = l...m, and F is the rate array of m reactions. The relative contribution, 1. (2.32) The CSP method solves the stiff O D E system of chemical reactions us-ing the following procedure: 1. Integrate the governing equation 2.21 using sufficiently small time steps. 2. Compute the Jacobian matrix numerically (this is where the name of the method comes from). Decompose the Jacobian matrix to 40 Chapter 2. Natural Gas Combustion - A Review find the eigenvalues and corresponding eigenvectors. Determine the dimension of the fast subspace and its basis vectors. 3. Solve Eq. 2.29 for f. 4. Repeat Steps 1 to 3 until G ^ s t in the fast subspace becomes triv-ial. 5. Bypass reactions related to the fast mode and integrate Eq. 2.21 using large time steps; update the Jacobian matrix and its eigen-values at each time step. 6. Whenever the dimension of the fast subspace changes, return to step 1. The CSP method resolves the stiffness problem in the governing ODEs of a reaction system by separating the time scales. In this way, the efficiency of integration is greatly increased. Because the reactions in the fast subspace are eliminated only when the projections in G™ast become insignificant, the accuracy of the detailed reaction mechanism is well preserved. L u et al. [89] compared model results using CSP with direct numerical integration for premixed C H 4 / a i r and H 2 / a i r flames. They reported very good agreement in terms of flame velocity and flame structure. Intrinsically Low-Dimensional Manifold The eigenvalues of the Jacobian matrix give the characteristic time scale of a reaction system. The number of large negative eigenvalues deter-mines the dimension of the fast subspace. For a typical combustion system, Maas and Pope [90] have shown that most of the eigenvalues of the Jacobian matrix are large and negative, which implies that the 41 Chapter 2. Natural Gas Combustion - A Review 25 Figure 2.4: Histogram of eigenvalues of natural gas - air combustion system at 2100 K , = 1. The system contains 71 species, 5 elements. There are five zero eigenvalues corresponding to the element mass con-servations that are not shown in the graph. rate of the system is controlled by a small number of slow reactions. Figure 2.4 shows an example of the histogram of eigenvalues in a ho-mogeneous, constant pressure natural gas - air reaction system. It can be seen that for the conditions investigated, a majority of eigen modes have very short time scales (r < le — 4s). In other words, a constrained reaction system tends to converge to a low-dimensional manifold (at-tractive subspace) as the reactions in the fast subspace exhaust rapidly. A constrained system is one that conserves elements and has defined boundary conditions, which lead to a unique equilibrium state, e.g. a constant-pressure, adiabatic reaction system. Figure 2.5 shows the con-42 Chapter 2. Natural Gas Combustion - A Review t=1 E-4 s t=1 E-3 s Figure 2.5: Convergence of reaction trajectories onto a one-dimensional I L D M . The reaction system is a stoichiometirc H 2 — CO — air system containing 13 species [90].The manifold is projected onto 4>HIO — on the manifold is a function of two parameters u and v: d<& _ d&(u, v) du ^ <9 is the species mass fraction on the manifold. The values of <1>, u and v are tabulated. In a reactive flow, the chemical reactions are perturbed by mass transport in the flow field through diffusion and convection. The perturbation alters the direction of reactions from the tangent plane on the manifold. Thus the value of d$>/dt is usually dif-ferent from that of dY/dt. It is necessary to project the perturbation into the manifold space. According to Pope and Maas [94], parameter u evolves following du — = P U - ( S + F) , (2.41) where S is the chemical source term and F is the perturbation; P u is a projection vector for u, which is given by ( * „ • * „ ) ( * „ • * „ ) - ( * u - * „ ) 2 ' K' ) where v is dQ/dv, and P is a specified n x n projection matrix. Similar equations apply to the evolution of v. Ideally, the pro-jection matrix should be based on the eigenvectors of the Jacobian ma-trix [94]; however, perpendicular projection in the composition space 47 Chapter 2. Natural Gas Combustion - A Review is sometimes used to simplify computation. This means that the pro-jection matrix is reduced to an identity matrix. A n important feature of the T G L D M is that the reaction vector S is in the tangent plane of the manifold so that the selection of a projection matrix only affects the perturbation term F but not the chemical source term S. The difference between the I L D M and the T G L D M with tailored ini -tial states can be examined by considering the following hypothetical reaction system: S1 h S2 % S3, • (2.43) with the initial conditions ^Sl,t=0 — 1) Xs2,t=0 — ^53,t=0 = 0. (2.44) Here S2 is an intermediate species, and ki S3. (2.45) The Jacobian matrix of this reaction system is given by J = (2.46) The three right eigenvalues of J are —k2, —ki and 0. The corresponding eigenvector matrix is written as / 0 V V 1 ^2— 1 0 (2.47) The matrix inverse of V gives V = 1 0 0 V 1 1 1 (2.48) 48 Chapter 2. Natural Gas Combustion - A Review The first row of V designates the projection vector to the fast sub-space. Based on eigenvectors corresponding to fast eigen modes, a one-dimensional I D L M is given by dY32 h dYSl dt k2 — ki dt whose solution is given by == 0, (2.49) Ys^jr^-Ys,, (2.50) k2 — ki which is the same as one would get from the quasi-steady-state analy-sis assuming d[S2]/dt = 0. The projection of the constrain due to mass conservation in the Ys1 — Ys2 plane is a straight line intersecting both axes at unity. The final I L D M subject to this constraint is shown in Fig. 2.6. For demonstration, the ratio of k2/k1 is set at 10 to empha-size the difference between fast and slow processes. It can be seen that, along the section of the I L D M defined by the linear constraint, the di-rection of the projected reaction vector is given by dYS2 = k2 kl + k2 8YSl hYSl h 1 • ; which is apparently not in the tangent plane of the I L D M . Thus a back projection is required to return the solution to the manifold. For the T G L D M , the initial state can be obtained based on the species involved in the global reaction. The linear system given by Eq. 6.5 is used to obtained the manifold generator [94] that defines the initial state of the trajectories. In order to solve Eq. 6.5, the number of major species to be involved in the global reactions should be equal to nc+np, ie. the total number of elements plus the dimension of the manifold. For the current simple system with a one-dimensional manifold, the initial state is defined by : ' ; ) ( £ ' ) - ( ! ) • 49 Chapter 2. Natural Gas Combustion - A Review Figure 2.6: Comparison of 1-D I L D M and T G L D M calculated for a hy-pothetical reaction system. The ratio of reaction rate between the fast and slow processes is 10. The manifolds are projected into Ysx - Ys2 plane. The thin lines represent trajectories subject to perturbations of various time scales. 50 Chapter 2. Natural Gas Combustion - A Review In forming Eq. 2.52, the a priori knowledge that YSl = 1 at t = 0 has been used. In reactive flow calculations, such information is usually available from initial conditions. The resulting T G L D M , which rep-resents the trajectory of an unperturbed reaction system, is shown in Fig. 2.6. When the system is perturbed by some physical process, such as mass transport by diffusion and convection, the Jacobian matrix takes the form / - f c i + ei 0 0 \\ J = I kx -k2 + e2 0 , (2.53) V - e i k2 - e2 0 / where e\\ and e2 represent arbitrary time scales associated with the per-turbations. We assume the perturbations occur in the same subspace as defined in Eq. 6.5. The assumption is consistent with that of identical diffusivity and unity Lewis number adopted by many turbulent com-bustion models [96,97]. The three eigenvalues of this perturbed system are 0, — k2 + e2, —k2-\\-e\\. The closed form solution is given by YSl = exp[-(k1-e1)t] (2.54) YS2 = ^ - ^ — ) { l - e X p [ - ^ - e 2 ) ] } YS3 = l - Y S l - Y S 3 . Under conditions when ei where the chemical reaction cannot sus-tain the heat loss due to turbulent mixing, a quenching of the flame occurs. In turbulent flame calculations, it is preferable to represent the profile of x with a single parameter. Following Law and Chung [102], Peters [96] proposed using the scalar dissipation at the stoichiometric surface, i.e. Xst = x(^at) as the representative parameter. Assuming flamelets of the mixing layer type are predominant in turbulent diffu-sion flames, Peters chose to use the counter-flow geometry to describe 56 Chapter 2. Natural Gas Combustion - A Review the Z dependence of the x profile. From the analytical solution, it can be shown that X{Z) = -exp{-2[er(c-1(2Z)}2}, (2.61) where erfc - 1 is the inverse of the complementary error function. The functional dependence of x{Z) on Xst can be derived from Eq. 2.61 X(Z) = Xstf(Z)/f(Zst), (2.62) where f(Z) is the exponential term in Eq. 2.61. In R A N S turbulence models, the mean value of scalar dissipation can be related to turbulent fluctuation through [96,103] X = cxjZ^, (2.63) k where e is the dissipation rate of turbulent kinetic energy; k is the tur-bulent kinetic energy; Z\"2 is the variance of the mixture fraction; cx is a scaling coefficient with a standard value of 2.0 [103]. The ~ indi-cates that the variable is a Favre-averaged quantity. The representing parameter, Xst in Eq. 2.62 is often equated its Favre-average, which can be calculated from the integral equation Once the x profile in the flow filed is solved, the mean value of a species mass fraction can be computed from the joint PDF of Z and x Yi= / Yi(Z;X)P(Z,x)dZdX. (2.65) Jo Jo Further simplification can be achieved if Z and x a r e assumed to be statically independent. In that case, Eq. 2.65 is reduced to poo rl Yt= / Yi(Z;x)P(Z)P(x)dZdX. (2.66) Jo Jo 57 Chapter 2. Natural Gas Combustion - A Review In the most direct implementation of the laminar flamelet model, the steady-state flamelet equations with various profiles of the scalar dissi-pation are solved numerically. The solutions are then tabulated to form a flamelet library with Z and x as the table dimensions. The mean mass fractions can then be obtained by solving the conservation equations of various moments of Z and x, and substituting the resulting PDF and conditional reaction rates from the pre-computed flamelet library into Eq. 2.66. The steady-state flamelet model relaxes the fast chemistry assump-tion significantly; however, experimental and theoretical studies have shown that the flamelet structure cannot respond instantaneously to changes in scalar dissipation [104-106]. Hence, the steady-state flam-elets are not suitable for modeling processes where the chemical time scale is comparable to or longer than the flow time scale. To address this issue, various unsteady flamelet models have been developed and tested with different levels of success. Coelho and Peters [107] studied a piloted methane/air diffusion flame using an Eulerian particle flamelet model. The unsteady calcu-lations were performed in the post-processing stage by transporting fluid particles that trace the temporal evolution of the scalar dissipa-tion rate and solve the unsteady flamelet equation with the varying x value. The results show an improvement in the predicted species con-centration profile compared with the steady-state model. A similar ap-proach, called the representative interactive flamelet model(RIF) pro-posed by Barths et al. [108], solves the unsteady flamelet equations in-teractively with the solution of the flow field. This method has been im-plemented in simulating combustion and pollutant formation in diesel engines [108-111] where the transient process dominates. Rao and Rut-58 Chapter 2. Natural Gas Combustion - A Review land [112] proposed a flamelet time scale model, which features lower computational cost compared with the RIF model. The model takes a first-order expansion of the steady-state flamelet solution. A chemical time scale determined from the Jacobian matrix is used to compute the rate of change of species mass fraction from the steady-state solution. Although the laminar flamelet model and its various derivatives are being used extensively in modeling turbulent combustion, it is impor-tant to realize their inherent limitations. The underlying assumption of flamelet models is that the turbulent flame is an ensemble of laminar flamelets. For this assumption to be valid, the structure of the flame front must remain locally laminar. In other words, the thickness of the flame must be thinner than the smallest length scale of turbulence - the Kolmogorov length scale. It is now generally accepted that the flamelet assumption is only valid in the region of large turbulent Damkohler number, Da, which is defined by Da = ^ = ±P-, (2.67) where rv is the macroscopic time scale of the flow field; TL is the charac-teristic time scale of the laminar flame; lp is the largest turbulent eddy length scale; v ' is the turbulence intensity ; lL is the laminar flame thickness, and vi is the laminar burning velocity. When the physi-cal time scale approaches the chemical time scale, or the value of Da is small, the suitability of the flamelet assumption becomes question-able [62,113]. 2.5.2 Conditional Moment Closure The conditional moment closure (CMC) model was proposed indepen-dently by Klimenko [114] and Bilger [115,116], and was described in 59 Chapter 2. Natural Gas Combustion - A Review details in their joint review [97]. Although the final form of the CMC equations are unified, the mathematical methods and model assump-tions adopted by Klimenko and Bilger in their derivations are quite different. Klimenko started his derivation from the transport equation of a two-dimensional joint PDF, P(Z) !M¥L+div((pv\\Z)P) (2.68) The above equation follows the convention used by Klimenko [97] where Z is the sample space variable for Y, Z = (Z1: Z2), W\\ is the chemical source term; the expression (a\\c) is short for (a\\b = c) which is the conditional expectation of a conditioned on the variable b being equal to c. Equation 2.68 is multiplied by Z and integrated over all Z to get dm™P^Kd,v(mmr,)P(v)) (2.69) = {P\\V)(W\\V)P(JI) + ^ - , ' where Jr . 2{pM{D(VY • V«|„)P(,) - ^ y » f W . (2.70) Here N = D(V£)2 is the dissipation rate of the conserved scalar; the variable £ = 77 = Z2. Closure of the last term on the r.h.s, J y , which is the reaction scalar flux in conserved scalar space, was achieved thro-ugh a diffusion approximation with the form JY = A(Y\\ri) + B ^ \\ l ! l . (2.71) ot] 60 Chapter 2. Natural Gas Combustion - A Review Here A and B are the drift coefficient and diffusion turbulent coeffi-cient respectively. This closure assumption leads to the final form of the basic CMC equation which governs the evolution of the conditional values of reaction scalars; (2.72) = {W\\v), . • • where v\" = v — (v\\ri) is the velocity fluctuation about the conditional mean; similarly, Y\" is the fluctuation of species mass fraction. The phys-ical meaning of the second and third terms on the l.h.s. of Eq. 2.72 are the convection of the conditional value of the reaction scalar; the forth term on the l.h.s. represents the effect of turbulent diffusion on the conditional expectation; the term on the r.h.s. is the conditional source term. In Bilger's derivation of the CMC equation [115,116], the condi-tional value of the reaction scalar is decomposed into its mean and fluctuation, which are substituted into the transport equation. After taking the conditional average, the equation becomes = {p\\v){W\\v) + eQ + eY with eQ = (div(pDV(Y\\r))) + pDV^ • V^j^-\\n), (2.74) eY = -(P-QJ- + pv • V F \" - div{DpVY\")\\'|| + A||(Y|r/)* - *\"*ll}, (2-84) where Q is the original coefficient matrix for the discrete integral equa-tions; the superscripts t and t — dt is the times at which the scalars 64 Chapter 2. Natural Gas Combustion - A Review are evaluated; A is a weighting coefficient specified by the modeler. In Grout's [131] implementation, A was chosen to add just enough a pri-ori information to produce a well-behaved solution. The regularization term A|| (Y^) 4 — (Y|77)*-d*|| limits the change of conditional average be-tween two consecutive time steps and acts to stabilize the solution. A main advantage of the CSE method is that the computational cost is substantially lower than that of CMC. Meanwhile, it does not involve constraining assumptions such as those employed by the lam-inar flamelet model, and is thus applicable to a wide range of turbulent non-premixed flames. 2.6 Summary This chapter has reviewed experimental and kinetic studies on natural gas combustion. The chemical kinetic mechanisms that govern the oxi-dation of pure methane and methane with additions of higher alkanes and hydrogen have been discussed. The need for reliable experimen-tal data on natural gas ignition under engine-relevant conditions has been established. The main mechanisms for nitric oxide formation in hydrocarbon flames and several useful techniques for in-flame NOx re-duction have been introduced. Numerical techniques and mathemat-ical methods that are used frequently in analyzing and reducing de-tailed reaction mechanisms have been discussed. Finally, the closure problem in turbulent combustion modeling has been introduced, and three popular methods for closing the chemical source term using two-parameter representation of the PDF of reactive scalar have been dis-cussed. 65 C h a p t e r 3 Ignition of Methane with Higher Alkane Additives7 3.1 Introduction Ethane and propane are two non-methane alkanes often present in nat-ural gas; in combustion studies, they are often added to methane to represented typical natural gas [15]. The variation of ethane and pr-opane concentration in natural gas can significantly change the igni-tion characteristics of the base fuel; this is particularly relevant to the performance of H C C I engines [17,132] as well as to forced-ignition natural-gas engines and gas turbines in the sense of controlling au-toignition [16,18]. The issue is especially important for mobile applica-tions since natural gas composition can vary significantly from place to place. Previous researchers have conducted a large number of ex-perimental and numerical studies to understand the ignition behavior in homogeneous methane/ethane or methane/propane mixtures [10, 15,17-25,43]; however, most of these studies—particularly the experi-mental ones—were focused on systems at near atmospheric pressures and very high temperatures. Our knowledge of natural gas ignition 7 A version of this chapter has been accepted for publication. Huang. J. and Bushe, W. K. (2005) Experimental and kinetic study of autoignition in meth-ane/ethane/air and methane/propane/air mixtures under engine relevant conditions, Combust. Flame, In Press 66 Chapter 3. Ignition of Methane with Higher Alkane Additives Figure 3.1: Schematic of the shock tube facility used in premixed ex-periments under conditions relevant to practical combustion devices (such as IC engines or gas turbines) is still insufficient. Thus, the first objective of the work described in this chapter was to obtain reliable experimen-tal data that extend the ignition database of various methane/ethane, methane/propane and simulated natural gas mixtures burning in air under engine-relevant conditions. The second was to obtain a better understanding of the kinetic interaction between methane and the two higher alkanes when undergoing ignition in air. 3.2 Experimental Apparatus The experiments were conducted in a shock tube facility previously de-scribed in detail [28]. A schematic of the shock tube facility is shown in Fig. 3.1. The stainless-steel shock tube has a circular cross section with an inner diameter of 59 mm. The length of the driver section is 3.18m, 67 Chapter 3. Ignition of Methane with Higher Alkane Additives and length of the driven section is 4.25m. Five PCB dynamic pressure transducers were flush mounted along the driven section to measure the incident shock velocity. The ignition is indicated by the rise of the pressure signal from the transducer mounted at the end plate of the driven section. A double-diaphragm technique was used to guaran-tee the rupture of diaphragms at the desired pressure ratio. The pres-sure and temperature behind the reflected shock were calculated us-ing standard normal shock relations [133] incorporating temperature-dependent thermodynamic properties and the ideal gas state equation. Peterson et al. [29] have shown that the difference between the experi-mental conditions calculated based on the ideal gas state equation and those based on a real gas state equation is not significant in the pressure range covered by this study. In order to measure the long ignition delay times associated with relatively low temperatures, tailored interface conditions [133] were obtained by tuning the thermodynamic properties of the driver gas us-ing different mixtures of air and helium. The result is that only a Mach wave is generated at the contact surface when the reflected shock wave passes across; for an ideal shock reflection, the conditions in the exper-imental region remain unaffected until the arrival of the rarefaction wave [133]. In this study, undiluted air-fuel mixtures were used as the driven gas for the best simulation of conditions in IC engines. Since the driven gas is not monoatomic, interactions between the reflected shock wave and the boundary layer could lead to reflected-shock bi-furcations. It has been reported that streams of cold driven gas emerg-ing from the bifurcated shock can arrive at the end plate along the wal l of the shock tube and eventually surround a core of hot gas [134]. The effect of this cooling flow in the current experiments was examined 68 Chapter 3. Ignition oi Methane with Higher Alkane Additives 0.01 0.008 0.006 0.004 0.002 Calculated Maximum Experimental Time Calculated Cooling Flow Penetration Time Current Experimental Range 2.6 2.8 3 Mach Number 3.2 3.4 Figure 3.2: Determination of experimental time by comparing calcu-lated time for tailored interface conditions and the time for cooling flow to reach the end plate due to the reflected-shock/boundary-layer interaction. The dotted line indicates the upper limit of Mach number of the incident shock in the current study. numerically using models proposed by Mark [135] and Davis [136]. In Fig. 3.2, a comparison is made between the typical time scale for the cooling gas to reach the end plate and the experimental time scale de-termined by the arrival of the rarefaction fan. It can be seen that for the current shock-tube configuration and experimental conditions, the cooling flow from the bifurcated shock foot cannot affect the experi-mental conditions before the arrival of the rarefaction wave. The contact surface instability as a result of the reflected shock pas-sage is another phenomenon, which could lead to a premature termi-nation of the experiment. This is because an unstable contact surface can promote mixing between the cold driver gas and hot driven gas. 69 Chapter 3. Ignition of Methane with Higher Alkane Additives Markstein [137] applied the interface stability theory of Taylor [138] to the shocked contact surface to show that the perturbation velocity at the contact surface can be expressed by the formula u(x, y) = A U k A ^ ^ e ^ ± x + ^ (3.1) 92 + P3 Here, AU is the difference in velocities of driver gas immediately be-hind the incident and reflected shocks; k — 2ix/\\ where A is the wave-length; A is the amplitude of the perturbation; and p2 and p 3 are den-sities of driven and driver gas behind the incident shock respectively. The sign of this perturbation velocity depends on the sign of p 2 — Pz-Using this formula, it can be shown that the direction of the perturba-tion under the current experimental conditions is towards the driver gas side so that its effect on the test mixture is small. To minimize the uncertainty in the experimental temperature, the actual experimental time was limited below 50 percent of the valid tailored interface pe-riod, which is determined from the pressure trace measured at the end plate. The experimental uncertainty was calculated using a standard er-ror analysis procedure [139] based on the uncertainty in measured in-cident shock velocities and correcting for attenuation in both incident and reflected shock velocities. The root-mean-square of the tempera-ture uncertainty behind the reflected shock is 14 K , which agrees with the uncertainty calculated from measured reflected-shock Mach num-bers using a method suggested by Skinner [140]. Research-grade or ultra-high purity gases (methane-99.97%, ethane-99.999%, propane-99.99%) were used for testing mixtures, which were prepared in a stainless steel vessel using the partial pressure method. The uncertainty of ethane/propane concentration in the final mixture is less than 5% of the target concentration. 70 Chapter 3. Ignition of Methane with Higher Alkane Additives 3.3 Kinetic Model For the autoignition of CH 4 -a i r mixtures in high pressure and moder-ate temperature, recent kinetic studies reported by Petersen et al. [29] and Huang et al. [28] show an increased significance of methylperoxy (CH3O2) chemistry for temperatures below HOOK, which causes the global activation energy to change. The proposed kinetic mechanisms in the above two studies show good agreement with the experimental ignition delay for a wide range of initial conditions. GRI mechanism (vl.2) was taken as the base mechanism as in our previous work [28]. Due to the importance of alkylperoxy chemistry in the low-medium temperature range, reactions related to the formation of methylperoxy, ethylperoxy and propylperoxy were added. The reaction rates were taken from the work of Tsang et al. [141]. Unfortunately, there has been no experimental measurement reported for these reactions at condi-tions close to those of this study and the rates reported in the literature show a relatively large uncertainty. The rate coefficients of some re-actions were thus optimized based on the experimental results in this work. Hunter et al. [44] developed a detailed kinetic mechanism for the oxidation of ethane at elevated pressures (3-10 bar) and moderate tem-peratures (916-966K) based on experimental data from their flow re-actor. A n important reaction path identified in their analysis includes the formation and oxidation of ethylperoxy ( C 2 H 5 0 2 ) . The extended C2 reactions in the new mechanism were mostly taken from Hunter et al. [44]. The C3 sub-mechanism was based on the work of Frenklach et al. [22]. Reactions between iso/normal propylperoxy and major alkane species were taken from Curran et al. [142]. The final mechanism con-71 Chapter 3. Ignition of Methane with Higher Alkane Additives tains 55 species and 278 elementary reactions. The extended reactions and their rate coefficients are listed in Table 3.1. 72 Table 3.1: Modified and extended reactions in the mo-del (based on GRI-Mech 1.2) N o Reaction A B E(cal) Ref. 157 H 0 2 + C H 4 =i C H 3 + H 2 0 2 4.48E+13 0.0 24629.0 [38] 174 C 2 H 5 + 02 =i H 0 2 + C 2 H 4 2.90E+11 0.0 6856.0 [44] 177 C H 3 + OH =i C H 2 0 + H 2 8.00E+12 0.0 0.0 [29] 178 C H 3 + 0 2 ^ C H 3 0 2 8.52E+58 -15.0 17018.0 [141] 179 C H 3 0 + H 0 2 =i C H 2 0 + H 2 0 2 1.20E+13 0.0 0.0 [29] 180 C H 3 0 + C H 3 ^ C H 2 0 + C H 4 2.41E+13 0.0 0.0 [29] 181 C H 3 0 2 H ^ C H 3 0 + OH 4.00E+15 0.0 42920.0 [143] 182 C H 3 0 2 + C H 3 =i C H 3 0 + C H 3 0 3.00E+13 0.0 -1200.0 [29] 183 C H 3 0 2 + C H 4 =i C H 3 0 2 H + C H 3 5.40E+11 0.0 18580.0 [141], a 184 C H 3 0 2 + H 2 0 2 C H 3 0 2 H + H 0 2 2.40E+12 0.0 9942.0 [29] 185 C H 3 0 2 + C H 2 0 ^ C H 3 0 2 H •+ HCO 2.00E+12 0.0 11663.0 [29] 186 C H 3 0 2 + H 0 2 =± C H 3 0 2 H • + 0 2 4.60E+10 0.0 -2600.0 [29] 187 C H 3 0 2 + C H 3 0 2 =± 0 2 + C H 3 0 + C H 3 0 3.70E+11 0.0 2200.0 [44] 188 C H 3 0 2 + C H 3 C H O =i C H 3 0 2 H + C H 3 C O 2.80E+12 0.0 13592.0 [44] 189 C 2 H5 + 0 2 =^ C 2 H 5 0 2 1.10E+47 -10.6 14830.0 [44] 190 C 2 H 5 + 0 2 ^ C 2 H 5 0 + 0 1.10E+13 -0.2 27926.0 [44] 191 C 2 H 5 + 0 2 =± C H 3 C H O + OH 1.60E+14 -1.2 10392.0 [44] 192 C 2 H 5 + H 0 2 ^ C 2 H 5 0 + OH 1.15E+13 0.0 0.0 [44] Continued on next page Table 3.1 - continued from previous page N o Reaction A B E(cal) Ref. 193 C2H5 + CH3 ;=± C 2 H 4 + CH4 7.94E+11 0.0 0.0 [44] 194 C 2 H 5 + C H 3 0 2 ^ C H 3 O + C 2 H 5 0 2.40E+13 0.0 0.0 [44] 195 C 2 H 5 + C 2 H 5 ^ C 2 H 6 + C 2 H 4 1.45E+12 0.0 0.0 [44] 196 C 2 H 4 + CH3O2 ^ C 2 H 4 0 + CH3O 2.82E+12 0.0 17104.0 [44] 197 C2H4O ^ C H 3 + HCO 3.40E+13 0.0 57788.0 [44] 198 C2H4O ^ C H 3 C H O 5.48E+11 0.0 52341.0 [44] 199 C2H4O + H ^ C 2 H 3 0 + H 2 7.90E+13 0.0 9795.0 [44] 200 C 2 H 4 0 + OH ^ C 2 H 3 0 + H 2 0 1.78E+13 0.0 3607.0 [44] 201 C 2 H 5 0 ^ C H 2 0 + C H 3 1.00E+15 0.0 21600.0 [44] 202 C 2 H 5 0 ^ C H 3 C H O + H 2.51E+14 0.0 23387.0 [44] 203 C 2 H 5 0 + 0 2 ^ C H 3 C H O + H 0 2 6.03E+10 0.0 1648.0 [44] 204 C 2 H 5 0 2 + C H 2 0 ^ C 2 H 5 0 2 H + HCO 2.00E+12 0.0 11663.0 [44] 205 C 2 H 5 0 2 + H 0 2 ^ C 2 H 5 0 2 H + 0 2 1.75E+11 0.0 -1696.0 [44] 206 C 2 H 5 0 2 + C H 4 ^ C H 3 + C 2 H 5 0 2 H 8.10E+11 0.0 18580.0 [44],a 207 C 2 H 5 0 2 + C 2 H 4 ^ C 2 H 4 0 + C 2 H 5 0 2.82E+12 0.0 17105.0 [44] 208 C 2 H 5 0 2 + C 2 H 4 ^ C 2 H 3 + C 2 H 5 0 2 H 1.12E+12 0.0 30411.0 [44] 209 C 2 H 5 0 2 + C 2 H 6 ^ C 2 H 5 + C 2 H 5 0 2 H 1.70E+13 0.0 20450.0 [44] 210 C 2 H 5 0 2 + C 3 H 8 ^ n C 3 H 7 + C 2 H 5 0 2 H 1.70E+13 0.0 20450.0 [44] 211 C 2 H s 0 2 + C 3 H 8 ^ i C 3 H 7 + C 2 H s 0 2 H 1.70E+13 0.0 20450.0 [44] 212 C 2 H 5 0 2 H ^ C 2 H 5 0 + OH 4.00E+15 0.0 42920.0 [143] Continued on next page Chapter 3. Ignition of Methane with Higher Alkane Additives A cj cj fl 3 fl T—H 3 fl (S T—I ^ 1 ^ 4 ^ 4 - ^ 4 ^ 4 T—I Tt< T—1 T — I T — I T—< T — I T — I T — I T — I T — I I — I O L O o o o C O o O N o o o C O o o V O C O ^ 00 ^ _ _ I T ) U\\ 00 0\\ ON N N in ^ ON ^ ' O N C N T—I T—I L O T—I O O O O ^ § ° 00 C O C O o m O N O N L ^ T—I O N T—I T—I o o o o o o o o o o o o o o o o LO V O V O V O CM C O C O O O O N O N ^ O O O O O V O O O O J ^ - O O O O O O O T - H T - J V O ' O O O T - H T - H O O V O O O O O T - H T — I V O L O C O C N C O C N C M V O L \\ C O C O O \\ C M V O L O C M C M T—I T—I O + + o + + + + o o + + + + + + + + + + + V O V O 00 T—I o v o C N CO VO 00 K o a + o o T f a a CN CM o o L O O O O O C O O O v D C O O O C O C N O L O O O p p o N p T t j p T - j p G N v o v o L q ' ^ v o o q N ^ n ^ o \\ N m o N N H ^ H O \\ i n d o o CO a o + a CN o o o + X o + CO CN a a I N I N o o O M I N CM I N a o o o + a o H, 10 o a a CN oo O Q ffi I + W m o c5 U U I N CN o o a a I N I N o o co ^ T—I T—I CM a CO o 11 g o o CN CO a a u o + + CN CN o o a a CO CO o o fl fl 0) 00 fl PH •*-> X 0 fl o 0) .s 00 O N O T - I S CM CN CO CO r , CM CM CM CM U 0) 0) o a fl c« +J c 0) O) o u C 0) a 0) fl o> >. O 4p fl H- t c Ol o> o C J Si c/) cn Ol ! H OH CU 75 Chapter 3. Ignition of Methane with Higher Alkane Additives C N C N C N C N C N C N C N C N C N C N C N C N C N C N C N C N C N 1^ C N C N T—1 H i-H i—l l-H i-H i—l 1—1 l—l i—l i—l i—l i—l i—l i-H i-H C N C N o o O O O o o o o o o o o O O O O O C O C O o o L O L O o o o o o o o o o o o o o o C N C N o o C O C O V C vC 00 00 C O C O vC vC V C o V C o IN . IN . vc vC C N C N C O C O L O L O o o o o C N C N C O C O C O C O o o O N CT\\ 00 00 V C o O o I N . o I N . O O T—l l-H i i C O C O T-H i-H T-H 1—1 C N C N C N C N C N i—l C N i-H i—l i-H o o o o o o o o o o O o o O o o o o o o o o o o o o o o o o o o o o o o o o o o C N C N o o C O C O C N C N i-H i-H C O C O C O C O C O C N C O C N o o 1—1 i-H 1—1 I-H 1—1 i-H 1—1 1—1 1—1 i-H 1—1 1—1 i-H i-H 1—1 T-H l-H l-H + + + + + +. + + + + + + + +' + + + + + + w pq PL) PL) P-l p-l P-l pq p-l pq pq pq pq pq pq pq pq pq pq o o L O L O C O C O o o o o I-H r H o o o o o o l-H T-H v© 00 I N . IN . 1—1 C O C O i—i 1—1 o O L N I N . L N o L N o o o L O < N 1—5 i — ! 1—5 i-5 vc vC 00 00 C O C O r-5 T-5 i—i C N i—5 C N L O L O o CO a o o o a + a CN O u TL o o o CN CO a a + + CN CM o o b- b-a a CO CO o o C N C O C O C O C N C N a + ^ a 0 « I v O 1 I I b-« a % o u u 2 o a K + + 3 a *tf L O C O C O C N C N a cs K t- o a *~ co td o i=l + CO a cs o u cs o b-a + CO a CN o 1L CN O b-a a « a O cs b- O a *-CO ffi CJ « a cj + a a o CJ PI + + a a cs cs CJ cj o a cs a o « o a « a CJ CO a cs o 5 g U 1L o a CO CJ a CO CJ el + + a a o o CO CO a a o o vC I N . 00 O N C O C O C O C O C N C N C N C N + + H 1L ° 9 . + + a a o o O i—H C N C N y a cj el a o b-a CO o ei + cs a C N C N a CN o b-a CO CJ + a a CN cj + a CN o a CO o CD a o cj a + CO CS 9 o ^ b-+ a cs r^ a y a ^ M a a a cs cs cs o o o b- t- t-a a a CO CO CO vj CJ CJ el el ei a a a CS CO CO o o o + + + cs cs cs o o o b- b- b-a a a CO CO CO o o o ei el el a £ + + a a CS CN o o b- b-a a CO CO o o 'u TL 00 00 a a CO CO CJ CJ + + cs cs o o b- b-a a CO CO o o cj o + + a m ^ i , co a \"3 CO if) vO N C N C N C N C N C N OO O N C N C N CU 00 U 1L & E H H X, H H CD (M CS w o o 5 + + g 00 oo 1-^ 3 a a C2H5+CH302H C2H6+H02<=>H202+C2H5 C2H40<=>CH3+HCO CH302+CH4<=>CH302H+CH3 _CH3O2+CH3<=>CH3O+CH30 CH3+02<=>CH302 2CH3(+M)<=>C2H6(+M) 2H02<=>02+H202(R116) 2H02<=>02+H202(R115) ^ H02+CH4<=>CH3+H202 =i CH3+02<=>0+CH30 H02+CH3<=>OH+CH30 OH+CH20<=>HCO+H20z^ OH+CH4<=>CH3+H2C 20H(+M)<=>H202(+M) H 975K I I 1250K H+CH4<=>CH3+H2= H+02<=>0+0H -0.3 -0.2 -0.1 0 0.1 0.2 0.3 dz/dk Figure 3.6: Sensitivity of ignition delay with respect to rates of elemen-tary reactions for the 3.7% C 2 H 6 / 9 6 . 3 % C H 4 mixture 87 Chapter 3. Ignition of Methane with Higher Alkane Additives nation of two methyl radicals 2CH 3 (+M) =± C 2 H 6 ( + M ) (R158) is a main chain termination step at high temperatures [149]. The pres-ence of ethane addition at the beginning of reaction shifts the equilib-r ium of R158 so that less methyl is consumed in this reaction. More im-portantly, as shown by Westbrook et al. [150], the H-atom abstraction from ethane and the subsequent decomposition of ethyl radical is very efficient in producing H atoms, which lead to a rapid chain branching through H + 0 2 ^ O H + O. (R38) This mechanism is supported by Fig. 3.6, which shows that the sensi-tivity to R38 increases drastically with increasing temperature. While the ignition-promoting effect of ethane at high temperatures is rela-tively well understood, the above mechanism does not account for the enhanced effect of ethane at lower temperatures, since the sensitivi-ties to both R158 and R38 decrease significantly at lower temperatures. Previous work [28] shows that at temperatures below HOOK, the rate of methylperoxy formation C H 3 + 0 2 =± C H 3 0 2 (R179) rises rapidly. Reaction 179 starts to compete with R158 for methyl rad-icals. Since the reaction channel formed by R179 and C H 3 0 2 + CH3 =± 2 C H 3 0 (R182) is very efficient in accelerating the oxidation of methyl, R158 becomes less rate-limiting at low temperatures. In this regime, the overall reac-tion rate is mainly limited by the depletion of OH radicals. This trend 88 Chapter 3. Ignition of Methane with Higher Alkane Additives can be seen from Fig. 3.6 which shows a higher (negative) sensitivity with respect to a main OH formation step 20H(+M) ^ H 2 0 2 ( + M ) • (R85) at 975K compared with the sensitivity at 1250K. The increasing promoting effect of ethane appears to coincide with the increasing formation rate of methylperoxy with reducing tempera-ture. Correspondingly, the ignition delay shows an increased sensitiv-ity to C 2 H 6 + C H 3 0 2 ^ C 2 H 5 + C H 3 0 2 H (R193) at the lower temperature. As reported by Tsang et al. [141], the rate of R193 is significantly higher than the rate of reaction C H 4 + C H 3 0 2 ^ C H 3 + C H 3 0 2 H , (R183) thus more methylhydroperoxide ( C H 3 0 2 H ) is produced through re-action with ethane. Subsequent reaction following R183 and R193 in-volves the decomposition of C H 3 0 2 C H 3 0 2 H ^ O H + C H 3 0 , (R181) which forms an OH radical. To see whether this mechanism is sufficient to account for the enhanced effect of ethane at reduced temperature, a series of reaction flow analyses was conducted to examine the contri-bution from elementary reactions to OH generation at various stages of the induction period. As shown in Fig. 3.7, although the overall con-tribution from R181 is significantly less than R85 and R119 C H 3 + H 0 2 ^ C H 3 0 + OH, (R119) it does dominate the OH formation in the initial stage of the ignition process. The species involved in the R179-R193-R181 branch are all 89 Chapter 3. Ignition of Methane with Higher Alkane Additives 0.9 t [ms] Figure 3.7: Cumulative contribution to O H formation. The analysis was carried out for the fuel mixture of 3.7% C 2 H 6 and 96.3% C H 4 (mixture No . 1) at 40 bar and 975K. 90 Chapter 3. Ignition of Methane with Higher Alkane Additives present in the upper portion of the main reaction path, which explains the higher efficiency of this branch at the beginning of the reaction. In the later stage, the concentration of hydroperoxy (HO2) radicals in-creases, which leads to a faster accumulation of hydrogen peroxide (H 2 0 2 ) through 2 H 0 2 =± 0 2 + H 2 0 2 . (R115,R116) Subsequently, R85 catches up and replaces R181 as a main source for O H radicals. Another way to examine the contribution of R179-R193-R181 path is to reduce the reaction rate of R193 to be identical to that of R183. It is found that the reduction in ignition delay disappears at tem-peratures below HOOK as this accelerating mechanism is turned off. It should be noted that the formation and decomposition of ethylhy-droperoxide ( C 2 H 5 0 2 H ) has a similar effect of generating O H radical as that of C H 3 0 2 H ; however, the overall contribution from C 2 H 5 0 2 H is less significant than C H 3 0 2 H due to the relatively low concentration of ethane present in this work. The above analysis shows that the sensitization mechanism for me-thane ignition with minor ethane addition changes with temperature. A t higher temperatures (T HOOK), the H-atom abstraction of ethane leads to a quicker chain branching that accelerates ignition; at lower temperatures, the presence of ethane leads to a more efficient produc-tion of C H 3 0 2 H and C 2 H 5 0 2 H than reactions with methane. The sub-sequent decomposition of these alkylhydroperoxide species enhances the production of O H radicals at the early stage of the induction period, which significantly promotes the overall reaction rate. The transition from the high-temperature to the low-temperature mechanism occurs with the rise of C H 3 0 2 formation at around HOOK. This could account 91 Chapter 3. Ignition of Methane with Higher Alkane Additives for the change of efficiency with changing temperature that has been observed in this study. 3.4.3 Kinetic Analysis—Methane/Propane Mixtures A comprehensive kinetics study of methane ignition with propane ad-dition at high temperatures was reported by Frenklach et al. [22]. The ignition-promoting effect of propane is mainly attributed to its rapid decomposition: C H 3 + C 2 H 4 + (M) + H ^ C 3 H 8 ( + M ) , (R221 + R74) which provides active H radicals leading to the chain branching. The efficiency of propane, however, is lower than that of ethane because H -atom abstraction of propane forms propyl radicals which can decom-pose to produce methyl radicals [150] as that shown in Eq. R221+R74. As stated above, the recombination of two methyl radicals (Eq. R158) leads to the chain termination immediately. Similar to ethane, it is found that the reaction between propane and methylperoxy needs to be taken into account to address the enhanced promoting effect at lower tem-peratures: C H 3 0 2 + C 3 H 8 ^ C H 3 0 2 H + n C 3 H 7 (R224) and C H 3 0 2 + C 3 H 8 ^ C H 3 0 2 H + i C 3 H 7 . (R225) The rates of R224 and R225 are comparable to that of R193 so that pro-pane shows a similar promoting effect to that of ethane at lower tem-peratures. 92 Chapter 3. Ignition of Methane with Higher Alkane Additives 3.4.4 Analyt ical M o d e l Figure 3.8 presents the main reaction path during the ignition in a methane/ethane/propane mixture, which is identified based on sen-sitivity and integrated reaction flow analyses. Figure 3.8 reveals two facts regarding the ignition mechanism of methane with minor ethane-/propane additives: 1. The presence of small fractions of C2H6-C3H8 does not modify the main reaction path of methane (this is also true at higher temper-atures [22]). 2. The ignition of methane fuel in this study is dominated by a l im-ited number of key elementary reactions. It is thus possible to describe the current reaction system with a rel-atively simple analytical model. One advantage of analytical models comes from the explicit description of the relation between an inte-grated system and individual factors. It allows a deeper insight into the general nature of the system that are otherwise hard to identify with numerical approaches. In the past, such analytical models based on skeletal mechanisms for methane ignition have been reported by several researchers [149,151,152]; however, none of the previous mod-els has taken into account the alkylperpxy chemistry, which is essential under the conditions of this study. In this work, this new mechanistic feature was included to address the transitional behavior of natural gas at intermediate temperature. The derivation of an analytical model requires further simplifica-tion from the main reaction path identified in Fig. 3.8. Following the description of Fig. 3.8 and the previous analysis, the key reactions in the methane oxidation path can be sorted into three main steps: 93 Chapter 3. Ignition of Methane with Higher Alkane Additives OH C H 3 0 2 ^ C H 3 HO CH, CH 4 C 2 H 6 C 3 H 8 C H 3 0 2 H CH.O HCCO Figure 3.8: M a i n reaction path for methane oxidation during the induc-tion period with the presence of minor ethane and propane additives. The \" R 0 2 , R 0 2 H path\" represents reactions related to the formation and decomposition of C 2 H 5 0 2 , C 2 H 5 0 2 H , C 3 H 7 O 2 and C 3 H 7 O 2 H radi-cals. 94 Chapter 3. Ignition of Methane with Higher Alkane Additives 1. Chain initiation C H 4 + 0 2 ^ C H 3 + H 0 2 (R118) 2. Chain propagation (a) Ma in propagation path: C H 4 + O H ^ C H 3 + H 2 0 (R98) C H 3 + H 0 2 ^ C H 3 0 + OH (R119) C H 3 0 + 0 2 ^ C H 2 0 + H 0 2 (R170) C H 2 0 + O H ^ CHO + H 2 0 (R101) C H 2 0 + H 0 2 ^ CHO + H 2 0 2 (R121) CHO + 0 2 ^ CO + H 0 2 (R168) The oxidation of C H 3 0 radicals occurs at a significantly faster rate than their formation; R119 and R170 can be essentially written as C H 3 + 0 2 ^ C H 2 0 + OH. (S3) In the lower part of the reaction chain, hydrogen peroxide from R121 rapidly decomposes to OH radicals via R85. The conversion between the H 0 2 and OH radicals occurs thro-ugh the path of R115, R116 and R85. These radical-radical conversions have much higher rates than the main chain propagation steps, so that R101, R121 and R168 can be read-ily approximated with a one-step reaction C H 2 0 + | o 2 ^ C O + i H 2 0 + OH. (S4) 95 Chapter 3. Ignition of Methane with Higher Alkane Additives (b) Methylperoxy chemistry: C H 3 + 0 2 =± C H 3 O 2 (R179) CH3O2 + C H 3 ^ 2 C H 3 0 . (R182) Interestingly, a one-step approximation of R179, R182 (with the consideration of R119 and R170) takes the same form as that of S3. 3. Chain termination 2 C H 3 =i C 2 H 6 . (R158) In summary, the pre-ignition reactions of methane can be simplified into five steps (with some reactions renumbered): C H 4 + 0 2 ->• C H 3 + H 0 2 (Si) C H 4 + OH —> C H 3 + H 2 0 ' (S2) C H 3 + 0 2 -* C H 2 0 + OH (S3) C H 2 0 + ^ 0 2 CO + ^ H 2 0 + OH (S4) 2CH3 —» C 2H6- (S5) For simplicity, all the reactions are assumed to be non-reversible. The rate laws for intermediate species concentrations in the reduced mechanism are given by = -k2[CHA][OH] + k3[CH3][02] + kA[CH20}[02f4 (3.2) d[CH3] I S I T T 1 ;„ r p u i r n 1 o7„ \\nv 12 dt d[CH2Q] dt k2[CH4}[02] - k3[CH3}[02] - 2h[CH3}2 (3.3) = k3[CH,}[02] - kA[CH20][02fl\\ (3.4) 96 Chapter 3. Ignition of Methane with Higher Alkane Additives where k\\ to k5 are the rate constants of reactions S i to S5 respectively. As stated above, the rate k3 represents a combined effect of the main propagation steps and the branch through methylperoxy. Some \"conventional\" simplifications for the ignition analysis can be applied to the current system [151,152]. During the induction period, the temperature and pressure in the reaction system do not change sig-nificantly so that constant reaction rates can be used. Meanwhile, the concentrations of methane and oxygen remain nearly unchanged prior to ignition [153]; one can approximate the concentration of methane and oxygen with their initial values. Because OH is an active radical, its rate of consumption is signif-icantly faster than major intermediates, while its net concentration is much lower. One can readily apply the steady state assumption [62] to OH radicals so that their concentration can be calculated from l n m „ k3[CH3}[02} + k4[CH20]{02}^ [ 0 H ] k^cm • (3-5) Strictly speaking, the steady assumption can only be applied to CH 2 0 when its consumption rate is significantly higher than that of its generation (this implies a large coefficient, kA, for the stabilization term.) This is not the case in the early phase of ignition delay because the conversion from methoxy to formaldehyde occurs at a relatively high rate (since R170 does not rely on radicals). In the latter phase, with the increase of OH and H 0 2 radicals, the concentration of CH 2 0 does level off rapidly. Since the initiation phase represents a relatively small portion of the induction period, one can roughly assume partial-equilibrium between the generation and consumption of C H 2 0 for most of the induction period. With the above assumptions and simplifica-tions, the rate law of methyl radicals is given by 97 Chapter 3. Ignition of Methane with Higher Alkane Additives « h[CH4}[02} + k3[CH3}[02] - 2k5[CH3}2. (3.6) Since rate change of C H 3 concentration best reflects the rate of con-sumption of methane, it can be used to characterize the system progress in the induction period. A close look at Eq. 3.6 shows that it well rep-resents the first and second phases of ignition [22] prior to thermal ex-plosion. The terms on the right hand side of Eq. 3.6 designate the initia-tion, propagation and termination steps respectively. It is interesting to note that Eq. 3.6 resembles the characteristic equation reported by Za-mansky for non-branching reactions system with no promoters [151]; the main difference is the presence of the second term in the current equation, which represents the self-promoting (feedback) mechanism in the chain propagation steps. A similar equation is also reported by Frenklach [22] for methyl radicals at higher temperatures and lower pressures, although the derivation is slightly different from this work. Mathematically, Eq. 3.6 is a first-order Ricatti's differential equation with constant coefficients (see Appendix C for details); its closed form solution is given by [ C H 3 ] = 4 + Q + l ) \" \" ' ( 3 ' 7 ) where o- = KUO2})2 + SkMCH4}[02}]1/2; (3.8) b = -2k5; (3.9) k3[Q2] - [(k3[Q2})2 + 8kMCH,}[02)}1/2 c = L- *—. (3.10) -2k5 The asymptotic (or equilibrium) concentration, which is a steady state solution of Eq. 3.7, can be found from 1/2 l r n ] fc3[02]+[(fc3[02])2 + 8fc1/c5[Cg4][02]J [CH3\\equi = — , (3 .11) 98 Chapter 3. Ignition of Methane with Higher Alkane Additives x 10 K— Ignition Delay ->! t [ms] Figure 3.9: Definition of ignition delay in the analytical model. The system approaches this concentration at the end of the second phase of the induction period, which represents a quasi-steady state where the kinetic growth of methyl is balanced by the second-order re-combination step. A breakthrough from this state with the presence of thermal feedback leads to the third phase of ignition - the chain explo-sion. The third phase, although important, takes only a small portion of the induction period [22]. It thus seems reasonable to approximate the ignition delay time with a characteristic time scale that represents the kinetic delay up to the end of the second phase. In this work, this time scale is defined by extrapolating the maximum growth rate of CFf 3 to the asymptotic concentration as shown in Fig. 3.9. With this definition, the ignition delay can be calculated from the formula 99 Chapter 3. Ignition of Methane with Higher Alkane Additives r = aU + ( 3 _ 1 2 ) —a a Table 3.4 lists the rate coefficients of reactions involved in the above skeletal mechanism. 100 Table 3.4: The skeletal mechanism used in the analyti-cal model N o Reaction A B E(cal) 4.04E13 0.0 56807 2.00E06 1.6 3120 1.08E12 0.0 22000 5.00E12[O 2] a 2 5 0 40000 4.24E14 -1.0 620 6.64E42 -8.24 93553 5.50E16 0.0 84416 1.20E14[M]9 0.0 14931 1.10E13/[O2]°-5 0.0 42996 3.54E05 2.1 870 3.16E06 1.8 994 1 2 3 4 5 6(C 2 H 6 ) 6(C 3 H 8 ) 7 8 9(C 2 H 6 ) 9(C 3 H 8 ) C H 4 + 0 2 -C H 4 + OH C H 3 + 0 2 -C H 3 + H 0 2 » C H 3 + H 2 0 C H 2 0 + OH C H 2 0 + | 0 2 -* CO + | H 2 0 + OH 2 C H 3 C 2 H 6 C 3 H 8 C 2 H 6 2 C H 3 . C 2 H 5 + C H 3 • C 2 H 5 + C H 3 0 2 H + ± H 2 0 C H 3 0 2 H + ± 0 2 -> C H 2 0 + 20H C 2 H 6 + OH - . H 2 0 + C 2 H 5 C 3 H 8 + OH - f H 2 0 + n C 3 H 7 C 2 H 6 + C H 4 + | 0 2 9total mixture concentration Chapter 3. Ignition of Methane with Higher Alkane Additives They are mostly taken from the rate-limiting reactions in the de-tailed mechanism. The pre-exponential coefficients of these reactions are adjusted to account for factors such as changing reaction order or effects from reverse-reactions. The activation energy and the tempera-ture exponents are unchanged to maintain the characteristic tempera-ture responses of the reaction rates. A comparison of the ignition delay calculated from Eq. 3.12 and that from predicted by the detailed mec-hanism is shown in Fig. 3.10. Fig. 3.11 shows a comparison of species concentrations prior to ig-nition based on the analytical model with the results from the detailed chemistry at 1000 and 1300K: It can be seen that the analytical solutions predict both the ignition delay and the methyl concentration very well. As expected, the agreement for the concentration of C H 2 0 radicals is very good at the end of the second phase where the partial-equilibrium assumption becomes more valid, but less good in the earlier phase. When ethane is present in the system, four skeletal steps need to be added to address the two promoting mechanisms aforementioned. The first is the decomposition of ethane to form methyl radicals: C 2 H 6 -> 2 C H 3 (SE6) To model the reactions between methylperoxy and ethane, the forma-tion of methylperoxy is approximated with a one-step reaction by com-bining the effects of R118, R85, R98 and R179: C H 4 + ^ 0 2 - C H 3 0 2 + ^ H 2 0 ; further combining this reaction with R193 froms C H 3 0 2 H C 2 H 6 + C H 4 + -02 -+ C 2 H 5 + C H 3 0 2 H + ^ H 2 0 . (SE7) 102 Chapter 3. Ignition of Methane with Higher Alkane Additives 1000K7T Figure 3.10: Comparison of ignition delay calculated from the model with that from the detailed chemistry at 16 bar (upper lines) and 40 bar (lower lines). 103 Chapter 3. Ignition of Methane with Higher Alkane Additives 0.5 Detailed Mechansim Analytical 1.5 2 t [ms] 2.5 3.5 (a) 40 bar, 1000K 10 -4 10 - 6 cj 1 10\" c o B 10 CO C J •10 o CJ 10\"\" 10\" CH 2 0 US O H -i\\ Detailed Mechansim Analytical 0.05 0.1 0.15 t [ms] 0.2 0.25 0.3 (b) 40 bar, 1300K Figure 3.11: Species concentrations predicted by the analytical model compared with the calculation from the detailed kinetic model for the stoichiometric methane/air mixture. 104 Chapter 3. Ignition of Methane with Higher Alkane Additives This is followed by a series of decomposition and oxidation reactions (R181, R170 and R85), which can be approximated by C H 3 0 2 H •+ -> C H 2 0 + 20H. (S8) Finally, the oxidation of ethane mainly goes through the reaction with OH radicals: C 2 H 6 + OH -> C 2 H 5 + H 2 0 . (SE9) With these extended reactions, the updated rate laws for the inter-mediate species are = -k2[CH4][OH] + k3[CH3][02] + h{CH20]{02f4 (3.13) at + 2k8[CH302H][02] 1/ 2 - k9[C2H6][OH] = h[CHA][02] + k2[02)[OH] - k3[CH3][02] (3.14) -2k5[CH3}2 + 2ke[C2He] d[Cft2°] = h[CH3}[02] - k,[CH3}[02fA + 2k&[CH302H][02] 1' 2 (3.15) d[CHf*H] = k7[CH4][C2H6][02}^ - U02Y'2. (3.16) It should be noted that, unlike methane, the concentration of ethane changes significantly during the induction period [153]. The detailed kinetic model shows that more than 70% of added ethane has decom-posed or been oxidized compared to around 10% of methane prior to ignition at 1300K. It is thus not appropriate to simplify the calculation wi th a constant ethane concentration. Instead, the concentration of et-hane must be calculated from its rate expression ^ 1 ^ 1 = -(k7[CH,][02f4 + k6 + k9[OH])[C2H6] + k5[CH3}2. (3.17) The contribution from the second term of Eq. 3.17 is negligible com-pared with that of the first term for the concentration of ethane present 105 Chapter 3. Ignition of Methane with Higher Alkane Additives in this study. An approximate solution of Eq. 3.17 can be obtained thr-ough asymptotic analysis. The change of ethane concentration is given by [C2He] = [C2H6]0exp j - (k7[CHA] [02fA + k6 + 2^^A[CH3]equ^ i j (3.18) where [C 2H 6] 0 is the initial concentration of C 2 H 6 . Subsequently, a dif-ferential equation for C H 3 0 2 H can be solved to get its concentration at time t [CH302H] = e-nt[e{n-k)t _ 1 ] ; ( 3 1 9 ) n — k where k = k7[CH4][02f\" + k6 + 2-^}^[CH3}equi] K2 [L>n4\\ m = k6[CH4][02fi-n = k&{02]1'2. The solution to this equation indicates that the concentration of C H 3 0 2 H should experience an initial growth followed by an exponential decay with the depletion of ethane. With the C 2 H 6 and C H 3 0 2 H concentrations known, one can pro-cess Eq. 3.13 using a similar steady-state assumption technique as that introduced earlier. This leads to a new rate law for the methyl concen-tration: +H-^5U) [ c g i ]-w-Comparing Eq. 3.20 with Eq. 3.6, one can find that the initiation term (the terms in the parentheses on the right hand side) of Eq. 3.6 is en-106 Chapter 3. Ignition of Methane with Higher Alkane Additives hanced through the decomposition of ethane and methylhydroperox-ide. The large rate constant k$ at low temperatures accounts for the higher promoting effect in that regime. However, ethane shows a neg-ative effect on the propagation term. It can be seen that if the concen-tration of ethane is significantly lower than that of methane, the second part of the propagation term becomes negligible, i.e. M f r g e ] < < x ( 3 2 1 ) (k2[CH4] + k9[C2H6})k3[02} ^ x' ^ l ) and the promoting effect from the initiation term dominates. With an increase in ethane concentration, the ignition promoting efficiency of ethane starts to decrease due to the increasing negative effect of its concentration on the propagation term. The physical meaning of this reduction is related to the competition between ethane and methane for OH radicals. A maximum reduction is achieved when the positive and negative effects of ethane are balanced. This trend in the change of ethane efficiency found by simply analyzing Eq. 3.20 agrees with the experimental results from this work as well as those reported in the literature. Also, since the concentrations of C2H6 and C H 3 0 2 H are two exponential functions which vanishes as they approach infinity, the asymptotic concentration of methyl does not change with ethane addition. The above characteristics revealed by the analytical equation are consistent with the numerical analyses presented earlier regarding the promoting mechanism of ethane. First, the kinetic effect of ethane pri-marily comes from its contribution to the initiation phase of the induc-tion period. Second, while additional ethane does not change the main reaction path of methane ignition (ethane does not lower the critical concentration of methyl required for the transition to the third phase), 107 Chapter 3. Ignition of Methane with Higher Alkane Additives it does change the rate of progress of different phases in the induction period. The skeletal mechanism for propane addition is very similar to that for ethane addition. For propane/methane mixtures, reaction SE7 is replaced by C 3 H 8 + CFf4 + ^ 0 2 n C 3 H 7 + C H 3 0 2 H + ^ H 2 0 (SP7) The decomposition of propane replaces that of ethane in SE6 C 3 H 8 - + C 2 H 5 + C H 3 . (SP6) The oxidation of propane goes through C 3 H 8 + O H n C 3 H 7 + H 2 0 . (SP9) The governing equation for methyl concentration with additional pro-pane takes a form analogous to that with ethane. Note that the decom-position reaction of propane contributes only one methyl radical com-pared with two of ethane; this could partially account for the weaker promoting effect of propane at relatively high temperatures. Fig. 3.12 presents the ignition delay calculated from the analytical model for two methane/ethane and methane/propane mixtures at 40 bar. Fig-ures 3.13 and 3.14 compare the predicted species concentrations of ethane/propane and mefhylhydroperoxide by the analytical model w-ith those calculated from the detailed chemistry. It can be seen that, de-spite the simplicity of the analytical model, the results agree reasonably wel l with the detailed chemistry and the experiments. The key features of the ignition system with the presence of ethane/propane are well captured by the analytical model. In particular, the model reproduces the transition of the promoting mechanism in the correct temperature range. 108 Chapter 3. Ignition of Methane with Higher Alkane Additives 5 10\" t-> 10 0.8 0.85 n 100% C H 4 1 1 1 Q 95% CH 4 , 5% C 3 H g 0.9 0.95 1 1000K/T 1.05 1.1 (b) 5% C 3 H 8 / 95% C H 4 / 40 bar Figure 3.12: Ignition delay with ethane/propane additives: comparison of analytical model and experimental results. 109 Chapter 3. Ignition of Methane with Higher Alkane Additives (b) 40 bar, BOOK Figure 3.13: Species concentrations predicted by the analytical model compared with the calculation from the detailed kinetic model for mix-ture No.3 (90% C H 4 , 1 0 % C 2 H 6 ) . 110 Chapter 3. Ignition of Methane with Higher Alkane Additives I i i i U 0 0.5 1 1.5 2 t [ms] (a) 40 bar, 1000K t i i i J 0 0.05 0.1 0.15 0.2 t [ms] (b) 40 bar, 1300K Figure 3.14: Species concentrations predicted by the analytical model compared with the calculation from the detailed kinetic model for mix-ture No.6 (95% CFf 4 , 5% C 3 H 8 ) . I l l Chapter 3. Ignition of Methane with Higher Alkane Additives 3.5 C o n c l u s i o n s Experimental and analytical studies have been conducted on the igni-tion of homogeneous methane/ethane/air and methane/propane/air mixtures at conditions relevant to practical combustion devices, par-ticularly IC engines and gas turbines. A complex behavior has been observed for the promoting effect of ethane/propane, which varies sig-nificantly with temperature and pressure. A detailed chemical kinetic model, which includes elementary reactions between ethane/propane and methylperoxy, have been used to study the new reaction system. The predicted ignition delay time based on this mechanism agrees well with the experimental results. The kinetic study reveals a switch of the dominant ignition-promoting mechanism for ethane/propane with changing temperature. The enhanced promoting effect of ethane and propane at low temperatures (T<1100K) is attributed to the increasing significance of methylperoxy and methylhydroperoxide chemistry. An analytical solution was obtained for the ignition delay of met-hane enriched with ethane and/or propane by systematically reduc-ing the skeletal mechanism identified from the reaction-flow analysis using the steady state assumption. The analytical model agrees well with the detailed-kinetic model for both ignition delay and the concen-trations of main intermediate species. Both the analytical model and kinetic analysis of the detailed mechanism show that the addition of ethane/propane does not change the main reaction path of pure me-thane during ignition. The promotion of ignition is realized through accelerating the initiation phase in the induction period. 112 Chapter 4 Ignition of Methane with Hydroqen Addition10 y u r u g e n 4.1 Introduction Hydrogen enriched methane fuels are being increasingly studied due to their economic and environmental benefits. Experimental and com-putational results reported recently [45-50] show that by extending the lean-burn limit or increasing the fraction of exhaust gas recircu-lation (EGR) with hydrogen addition, a significant reduction of N O x and hydrocarbon emissions from natural-gas-fueled internal combus-tion engines can be achieved with unaffected or slightly higher eng-ine efficiency. The gain in combustion efficiency mainly comes from the enhanced reaction rate and the higher flame velocity. Potential ap-plications of premixed methane/hydrogen fuels are also found in the aerospace [51] and gas turbine [52] industries. A n understanding of the autoignition behavior of these mixtures is clearly needed. Shock tube studies of high-temperature ignition in C H 4 - H 2 - 0 2 mixtures have been reported by Lifshitz et al. [10] as well as by Cheng and Oppenheim [10]. In both cases, the reactants were diluted with io A version of this chapter has been accepted for publication. Huang. J., Hi l l , P. G., Bushe, W. K. and Munshi, S. R. (2005) Experimental and kinetic study of shock initiated ignition in homogeneous methane-hydrogen-air mix-tures at engine-relevant conditions, Int. J. Chem. Kine., In Press 113 Chapter 4. Ignition of Methane with Hydrogen Addition 90 percent argon. The data of Lifshitz et al. [10] measured at a fixed pressure of 185 torr and covered temperatures from 1600 to 1800K. A thermal-based-promotion theory was proposed to account for the effects of hydrogen addition. Cheng and Oppenheim [12] conducted experiments for temperatures from 800 to 2000 K and pressures from 1 to 3 atm. They correlated the ignition delay of pure methane, pure hydrogen and their mixtures with the formula r = TCH4{1-e)W (4.1) where e is the mole fraction of hydrogen in the total fuel and T C H 4 and TH2 are the ignition delay times of pure methane and pure hydrogen under the same conditions. Fotache et al. [53] investigated the ignition delay of hydrogen-enriched methane by heated air using a counter-flow reactor. They identified three ignition regimes depending on the mole fractions of hydrogen. Methane ignition was found to benefit from hydrogen addition mainly due to the kinetic interactions between the two fuels. The modeling study shows that the promoting effect is enhanced by the spatial separation of the branching and termination steps resulting from the high diffusivity of atomic and molecular hy-drogen. Ju and Kiioka [51] conducted a numerical study of hydrogen-/methane ignition in a supersonic mixing layer. They found a reduc-tion of ignition delay proportional to the fraction of hydrogen addition. They suggested that the extra H and 0 radicals from H 2 oxidation was the main reason for the increased reaction rate. In reviewing the literature, it is clear that most of the prior work on the ignition of methane/hydrogen/air or methane/hydrogen/oxygen mixtures was at relatively low pressures (< 8 atm) [10,12,53] and was focused on high temperatures [10,12]; systematic studies of methane-114 Chapter 4. Ignition of Methane with Hydrogen Addition hydrogen-air ignition and the underlining chemical kinetics under con-ditions pertinent to internal combustion engines have rarely been re-ported. In this work, experimental and computational investigations have been conducted to address this specific issue. This study focuses on ignition in homogeneous, quiescent methane/hydrogen/air mix-tures behind a reflected shock [133] so that fluid dynamic effects can be isolated from chemical kinetics to simplify the analysis. 4.2 Experiments The experiments were conducted in a shock tube facility introduced in detail previously [28]. The stainless-steel shock tube has a circular cross section with an inner diameter of 59 mm. The length of the driver section is 3.18 m; the driven section is 4.25 m long. Ignition is initi-ated behind the reflected shock. A PCB pressure transducer mounted at the end plate and two photomultipliers (one equipped with a 427 nm interference filter) were used to detect the ignition based on pressure and optical emission signals. The experimental temperature was calcu-lated using the incident shock velocity, which is measured by 5 flush-mounted PCB pressure transducers with a minimum response time of 2 ji s. The uncertainty in the calculated temperatures is around 15K. To extend the experimental time for ignition to lower temperatures, the tailored-interface technique [27] was adopted; this technique allows a maximum experimental time of up to 5 ms to be achieved with the current facility. The compositions of gas mixtures in this work are presented in Ta-ble 4.1. 115 Chapter 4. Ignition of Methane with Hydrogen Addition Table 4.1: Molar fractions of main constituents in the test mixtures (4> = 1 for all mixtures, balanced by air (79%N2, 21%02)) Mixture # X H 2 X c H 4 1 15% 85% 2 35% 65% The focus of this work was on stoichiometric fuel-air mixtures so that the equivalence ratio was fixed at unity. The test gases had the fol-lowing purities: methane (99.95%), hydrogen (99.995%) and air (99.9%; 21%02, 79%N2). Ignition was defined by extrapolating the maximum value of dP/dt (measured at the end plate of the driven sections) back to the pressure right after shock reflection; this method has been widely adopted in shock tube ignition studies [9,30]. 4.3 Kinetic Model In the range of intermediate temperature (1000-1300K) and high pres-sure (16-230bar), previous studies show that the alkylperoxy chemis-try, especially reactions related to the formation and consumption of CH 30 2, are important to better describe the ignition characteristics of methane [28,29]. In this work, a mechanism from a previous work [28] was adopted; it is a modified version of a mechanism proposed by Petersen et al. [29]. The main part of the mechanism is identical to that of GRI-Mech 1.2 [39] with major CH 30 2 and C 2H 50 2 chemistry included. The numerical models based on this mechanism are able to reproduce experimental ignition delay results of methane/oxygen-/argon and methane/oxygen/nitrogen mixtures over a wide range of 116 Chapter 4. Ignition of Methane with Hydrogen Addition conditions [28]. It should be noted that the hydrogen subset of this mechanism is identical to that in GRI-Mech 1.2, which has been validated substan-tially using low-pressure, high-temperature experimental data, but not optimized for conditions similar to those of this study. A recent inves-tigation conducted by Akbar et al. [154] shows a significant discrep-ancy between the numerical results calculated using various hydro-gen mechanisms and experimental data for hydrogen/air ignition at pressures above 1 atm. To improve the performance of the new mech-anism at higher pressures, rate coefficients of several reactions key to hydrogen ignition were updated using recent values reported in the literature. These reactions were identified based on the results of pre-liminary sensitivity and reaction flow analyses. Among the changes, the rates of reactions related to H + 0 2 + M = i H 0 2 + M (reactions R36-R38) were replaced using values proposed by Baulch et al. [70]. The reaction with methane as a third body (R37) was added to the mechanism. For the rate of reaction C H 4 + H = i C H 3 + H 2 , (R54) the experimentally obtained values by Sutherland [155] were used. The rate constants of reaction OH + H 2 =± O + H 2 0 (R85) were adopted from the work of Isaacson [156]; his rate coefficients were supported by latest kinetic experiments of Smith and Cr im [157]. Due to the significance of C H 3 0 2 chemistry in the low-intermediate tem-perature range, reactions between hydrogen and methylperoxy were 117 Chapter 4. Ignition of Methane with Hydrogen Addition added: CH3O2 + H ^ C H 3 O + H (R194) and CH3O2 + H 2 ^ C H 3 O 2 H + H . (R195) The rate coefficients of these reactions were taken from Tsang et al. [141]. The final mechanism includes 40 species and 195 reversible reactions. A list of reactions extended or modified, including their rate coeffi-cients and source, is given in Table 4.2. 118 Table 4.2: Modified and extended reactions in the new model (based on GRI-Mech 1.2) Rx# Reaction A B E(cal) Ref. R1-R33 R1-R33, GRI-Mech 1.2 - - - [39] R34 H + 2 0 2 - H 0 2 + 0 2 1.41E+18 -0.8 0.0 [70] R35 H + 0 2 + H 2 0 ^ H 0 2 + H 2 0 9.38E+18 -0.8 0.0 [70] R36 H + 0 2 + N 2 ^ H 0 2 + N 2 2.20E+18 -0.8 0.0 [70] R37 H + 0 2 + C H 4 ^ H 0 2 + C H 4 4.40E+18 -0.8 0.0 u R38 H + 0 2 + A R ^ H 0 2 + A R 7.00E+17 -0.8 0.0 [70] R39-R53 R38-R52, GRI-Mech 1.2 - - - [39] R54 H + C H 4 ^ C H 3 + H 2 3.99E+03 3.2 8760.0 [155] R55-R84 R54-R83, GRI-Mech 1.2 - - - [39] R85 OH + H 2 ^ H + H 2 0 1.81E+09 1.2 4707.0 [156] R86-R157 R85-R156, GRI-Mech 1.2 - - - [39] R158 H 0 2 + C H 4 ^ C H 3 + H 2 0 2 4.48E+13 0.0 24629.0 [38] R159-R178 R158-R177, GRI-Mech 1.2 - - - [39] R179 C H 3 + OH ^ C H 2 0 + H 2 8.00E+12 0.0 0.0 [29] R180 C H 3 + 0 2 ^ C H 3 0 2 2.13E+58 -15.0 17018.0 [28] R181 C 2 H 5 0 ^ C H 2 0 + C H 3 1.00E+15 0.0 21606.0 [29] R182 C H 3 0 + H 0 2 ^ C H 2 0 + H 2 0 2 1.20E+13 0.0 0.0 [29] Continued on next page Rate estimated in this work Table 4.2 - continued from previous page Rx# Reaction A B E(cal) Ref. R183 C H 3 0 + C H 3 =i C H 2 0 + C H 4 2.41E+13 0.0 0.0 [29] R184 C H 3 0 2 H =± CH3O + OH 4.00E+15 0.0 42996.0 [143] R185 C H 3 0 2 + C H 3 =± CH3O + CH3O 3.00E+13 0.0 -1200.0 [29] R186 C H 3 0 2 + H 2 0 2 =i C H 3 0 2 H + H 0 2 2.40E+12 0.0 9942.0 [29] R187 C H 3 0 2 + C H 2 0 =i C H 3 0 2 H + HCO 2.00E+12 0.0 11663.0 [29] R188 C H 3 0 2 + C H 4 C H 3 0 2 H + C H 3 1.80E+11 0.0 18475.0 [29] R189 C 2 H 5 + 0 2 =i C 2 H 5 0 2 1.00E+12 0.0 0.0 [29] R190 C 2 H 5 0 2 + C H 2 0 C 2 H 5 0 2 H + HCO 2.00E+12 0.0 11663.0 [29] R191 C 2 H 5 0 2 H =± C 2 H 5 0 + OH 1.00E+16 0.0 42977.0 [29] R192 C H 3 0 2 + H 0 2 =i C H 3 0 2 H + 0 2 4.60E+10 0.0 -2600.0 [29] R193 C H 3 0 2 + C H 3 0 2 =± 0 2 + C H 3 0 + CH3O 3.70E+11 0.0 2200.0 [29] R194 C H 3 0 2 + H =± CH3O + OH 9.60E+13 0.0 0.0 [141] R195 C H 3 0 2 + H 2 =± C H 3 0 2 H + H 1.50E+14 0.0 26030.0 [141] The forward rate constant k f is calculated using the formula k f = AT^exp(E/RT); A , E are in cgs units. The reverse rate constant k r is calculated from the forward rate constant and the equilibrium constant K e q using k r = k f / K e q . A l l the reactions in the mechanism are reversible. Chapter 4. Ignition of Methane with Hydrogen Addition The numerical simulation of ignition behind the reflected shock was performed using a zero-dimensional model with constant-density and adiabatic boundary conditions. A system of ordinary differential equa-tions (ODEs) containing mass and energy conservation equations was solved using V O D E [148]. To be consistent with the base mechanism, the thermodynamic properties of reactants were mostly taken from the GRI thermo database; species properties not available in the GRI database were taken from those of Burcat [146] and Curran et al. [158]. The gas-phase reaction rates and mixture properties were calculated using Chemkin II [147]. To validate the modification to the base mechanism, computations were performed for hydrogen/oxygen/argon and hydrogen/air igni-tions at various pressures. Figure 4.1 presents the results of calculated ignition delay time using the modified mechanism together with ex-perimental data from the literature [30,159-161]. The results from the base mechanism (GRI-Mech 1.2) are also presented for comparison. It can be seen that the modified mechanism basically maintains or slightly improves the accuracy of GRI-Mech 1.2 in predicting ignitions in hydrogen/oxygen/argon mixtures, while it significantly improves the agreement between numerical and experimental results for igni-tions in hydrogen/air mixtures. A sensitivity analysis shows that the improvement can be largely attributed to the slightly higher rate of H + 0 2 + N 2 ^ H 0 2 + N 2 (R36) in the new mechanism at intermediate temperatures compared with the original rate. As we w i l l show later, R36 is a main H-scavenging step for hydrogen ignition; the increase in the rate of R36 typically leads to a reduction in the overall reaction rate. To examine if the modi-121 Chapter 4. Ignition of Methane with Hydrogen Addition 0.9 1 1000K/T (a) Ignition Delay in H2/02/Ar Mixtures 0.75 0.8 0.85 0.9 0.95 lOOOKyT 1.05 1.1 (b) Ignition Delay in H 2 / A i r Mixtures Figure 4.1: Calculated ignition delay of hydrogen/oxygen/argon mix-tures with experimental data from the literature. Symbols designate experimental data. In plot a: + - 2%H 2 + 1%02 + 97%Ar [159]; o -8%H 2 + 2%0 2 + 90%Ar [160]; x - 2%H 2 + 1%02 + 97%Ar [30];D -2%H 2 .+ 1%02 + 97%Ar [30]. In plot b: 15%H 2 + 85%Air [161] for all mixtures. Lines designate calculated delay. 122 Chapter 4. Ignition of Methane with Hydrogen Addition fication has caused any significant deterioration of the base mechanism in the regime where it has been extensively optimized, extra computa-tions have been performed for the ignition of methane under a variety of conditions and compared results with the experimental data from the literature as well as the predictions using GRI-Mechl.2 and GRI-Mech 3.0 [41]. It was found that above modifications did not signifi-cantly affect the calculated methane ignition results as demonstrated in Fig. 4.2, given as an example. 4.4 Results and Discussion Table 4.3 lists the experimental conditions and measured ignition delay time; the results are also presented graphically in Figure 4.3. 123 Chapter 4. Ignition of Methane with Hydrogen Addition 1 1 1 1 1 • • Experimental (Seery and Bowman) Current Mechanism - - - GRI-Mech v 1.2 GRI-Mech v3.0 I I I I I I I 'I 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 1000K/T (a) Low-pressure C H 4 Ignition Delay 1 1 40 atm s ' s ' ' 130 atm • Experimental (Huang et al.) O Experimental (Petersen et al.) Current Mechanism - - - GRI-Mechv 1.2 - - GRI-Mech v3.0 y\\ 1 1 • ! 0.7 0.75 0.8 0.85 0.9 0.95 1 1000K/T (b) High-pressure C H 4 Ignition Delay Figure 4.2: Predicted methane ignition delay using the new mechan-ism. In plot a, the experimental data were taken from Seery and Bow-man [9] for a C H 4 - 0 2 - Ar (9.1%-18.2%-72.7%) mixture with a con-centration of 14.6 mol /m3 (1.8-2.1 atm). In plot b, the experimental data were taken from Huang et al. [28] for a CH 4 -a i r mixture at 40 atm (0 = 1) and Petersen et al. [29] for a C H 4 - 0 2 - Ar mixture at 130 atm (0 = 6) 124 Chapter 4. Ignition of Methane with Hydrogen Addition Table 4.3: Experimental conditions and ignition delay results for methane/hydrogen mixtures mol H 2 /molf u e i T ( K ) P (bar) mol H 2 /mo l f u e i T (K) P (bar) T(/IS) 0.15 1297 36.4 332 0.15 1307 17.0 462 0.15 1259 40.2 405 0.15 1284 17.2 552 0.15 1212 40.3 588 0.15 1239 16.0 738 0.15 1183 42.7 672 0.15 1213 16.4 840 0.15 1116 39.2 752 0.15 1178 17.2 1170 0.15 1083 39.4 1164 0.15 1155 14.8 1014 0.15 1056 33.4 1356 0.15 1146 17.5 1152 0.15 1028 37.5 1614 0.15 1096 18.8 1676 0.35 1318 37.4 158 0.35 1318 17.0 260 0.35 1276 39.8 264 0.35 1276 15.9 372 0.35 1235 40.2 390 0.35 1201 16.9 876 0.35 1201 41.1 726 0.35 1134 15.6 1098 0.35 1156 41.1 759 0.35 1134 15.9 864 0.35 1103 37.6 1056 0.35 1103 15.9 1044 0.35 1072 36.0 882 0.35 1079 16.6 1736 0.35 1045 36.6 1500 0.35 1051 16.8 2189 0.35 1043 38.3 1986 0.35 1036 16.2 2016 0.35 1035 39.2 2166 0.35 1032 16.8 1822 0.35 1019 16.0 3250 For the mixture with a high hydrogen fraction (mixture #2), the promoting effect of hydrogen addition is more evident at high tem-peratures; however, a rapid reduction in the difference between pure-methane and methane/hydrogen mixtures is observed as the tempera-ture decreases. Specifically, at 1300K and 40 bar, the hydrogen addition in mixture #2 reduces the ignition delay time of pure methane by a fac-tor of 1.5, but the difference is hardly distinguishable for temperatures below 1200K. A similar trend is observed at 16 bar but the overall dif-ference between the two mixtures is more evident. For a low hydrogen fraction (in mixture #1), the effect of hydrogen appears to be largely 125 Chapter 4. Ignition of Methane with Hydrogen Addition 10 i 10\" io; / i ^ K o <\\ ^ r < [ O 0 -o < V ^ O > ^ / 0 o • 6 5 % C H 4 , 3 5 % H 2 8 5 % C H 4 , 1 5 % H 2 1 0 0 % C H , 4 0.7 0 .75 0.8 0 .85 1 0 0 0 K / T 0.9 0 .95 (a) 16 bar 10 O 6 5 % C H 4 , 3 5 % H 2 8 5 % C H 4 , 1 5 % H 2 1 0 0 % C H . 0 .85 0.9 1 0 0 0 K / T 0.95 (b) 40 bar Figure 4.3: Measured and calculated ignition delay of test mixtures. Symbols designate experimental results: • - Pure methane ($=1); o -Mixture #1 (85%CH 4 ,15%H 2); 0 - Mixture #2 (65%CH 4 , 35%H2). Lines represent calculated ignition delay from the new model. A horizontal error bar represents the typical uncertainty in the experimental tem-perature. 126 Chapter 4. Ignition of Methane with Hydrogen Addition inhibited. The reduction in ignition delay time with hydrogen in mix-ture #1 is small at 16 bar and can be barely seen at 40 bar. The agreement between the numerical results and experimental data is less than ideal; the new mechanism under-predicts the ignition delay of methane/hydrogen mixtures at relatively high temperatures, which highlights the need to further refine the elementary reaction rates. De-spite the discrepancy, the new model shows a reduction in the promot-ing effect with decreasing temperature, which is in agreement with the experimental observation described above. The model also predicted the slightly stronger effect of hydrogen addition at 16 bar compared to that at 40 bar. A n empirical correlation with the form r = A[CH4}a[02}b[H2}cexp(E/RT) (4.2) was fit to the experimental results of this work (see Fig.4.4). Table 4.4 lists the coefficients in Eq. 4.2 obtained by a linear regression method. Table 4.4: Correlation coefficients for Eq. 4.2 A(s . (cc /mol ) a + b + c ) a b c E (kcal/mol) 2.10E-8 2.31 -2.83 -0.0055 19.06 The activation energy E of the current correlation is 19.06 kcal /mol , which is consistent with value of 19.0 kca l /mol reported by Petersen et al. [29] for pure methane ignition under similar temperature condi-tions. The overall effect of hydrogen on methane ignition observed in this work is weak as indicated by a small, negative value of the coeffi-cient c in the above correlation (note that c < a). 127 Chapter 4. Ignition of Methane with Hydrogen Addition Correlation 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1000K/T Figure 4.4: Comparison of measured ignition delay with formula 4.2 Fotache et al. [53] attribute the significantly stronger effect of hydr-ogen on methane/air ignition in a diffusive system to the separation of hydrogen from other species due to its high diffusivity. In a homoge-neous system, the scavenging effect of methane on hydrogen radicals can be much stronger, leading to the lesser effect of hydrogen addition on the reaction system. For a methane/air mixture with a high hydro-gen mole fraction, Fotache et al. [53] report that the ignition is effected through radical rather than thermal explosion. To confirm this point, we removed reaction heat of all elementary reactions from the energy equation. The results show an increase of ignition delay for both meth-ane and methane/hydrogen mixtures; however, the relative difference between mixtures does not change significantly. This suggests that the thermal feedback in the current mixtures, being important to the over-all reaction rate, has similar contributions to the oxidation of methane 128 Chapter 4. Ignition of Methane with Hydrogen Addition and hydrogen. In other words, the global activation energy of pure methane and methane/hydrogen systems may be similar under the current conditions. The promoting effect from hydrogen appears to be mainly kinetic based. A \"brute force\" sensitivity analysis was conducted to investigate the significance of elementary reactions to the overall reaction rate. The normalized sensitivity is calculated using the formula _ T(2kj) - r(0.5fcj) h = r (2fc < ) - r (0 .5fc i ) k i ~ 2ki - 0.5ki r(ki) 1.5r(fci) ^ ' ' The results are presented in Fig. 4.5 for pure methane and 65%CH 4 — 35%H 2 mixtures respectively. Comparing the results between the two mixtures, one immediately observe that the sensitivities to reactions 2 C H 3 ( + M ) =i C 2 H 6 ( + M ) (R159) and H 0 2 + C H 3 =i OH + C H 3 0 (R120) are reduced significantly with the presence of hydrogen at 1300K, while, the sensitivities to reactions O H + H 2 =± H + H 2 0 (R85) and H + 0 2 =± O + OH (R39) increase drastically. The recombination of two methyl radicals in R159 is a well-known chain termination reaction at high temperatures. Re-action R120 is a main competitor of R159 for methyl radicals; it also generates an active hydroxyl radical which is a crucial species in the 129 Chapter 4. Ignition of Methane with Hydrogen Addition ! CH302+CH4<=>CH302H+CH3 ! CH3+02<=>CH302 2CH3(+M)<=>C2H6(+M)S ! H02+CH4<=>CH3+H202 ! CH3+02<=>0+CH30 I H02+CH20<=>HCO+H202 ! H02+CH3<=>OH+CH30 2H02<=>02+H202(R116) U H M H B = I 2H02<=>02+H202(R115) mm ^ = ^ H + C H 4 < = > C H 3 + H 2 0 OH+H202<=>H02+H20 OH+H02<=>02+H20 £=• ^ ^ ^ ^ ^ ^ 20H(+M)<=>H202(+M) • I 1000K I I 1300K H+CH4<=>CH3+H2f ! 02+CH20<=>H02+HCO -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.; Normalized Sensitivity (a) 100%CH4 I CH302+H2<=>CH302H+H CH302+CH3<=>CH30+CH30L — 2CH3(+M)<=>C2H6(+M) •mm =* H02+CH4<=>CH3+H202 H02+CH20<=>HCO+H202 H02+CH3<=>OH+CH30 2H02<=>02+H202(R116) S 2H02<=>02+H202(R115) I O H + C H 4 < = > C H 3 + H 2 0 ^ ™ [ OH+H02<=>02+H20i ^ 20H(+M)<=>H202(+M) ^.OH+H2<=>H+H20 =^H+CH4<=>CH3+H2 H+02<=>0+OH H+02+N2<=>H02+N2 S mm IOOOK I I 1300K -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 Normalized Sensitivity 0.3 0.4 0.5 (b) 65%CH 4,35%H 2 Figure 4.5: Normalized sensitivity of ignition delay with respect to in-dividual reaction rate for methane and methane/hydrogen mixtures at 40 bar. 130 Chapter 4. Ignition of Methane with Hydrogen Addition main oxidation path of methane. With the presence of hydrogen, the production rate of hydrogen radical through OH + H 2 =t H + H 2 0 (R85) increases. Hydrogen radicals are more active than OH radicals, they branch through H + 0 2 =i O + OH, (R39) which rapidly increases radical concentrations and hence promotes the ignition. A t 1000K, the rate limiting mechanism changes. As shown in our previous study [28], with the reduction of temperature, the total reac-tion rate is increasingly limited by the depletion of hydroxyl. This trend is demonstrated by the increased sensitivity to a main OH generating path through H 0 2 + C H 4 =± C H 3 + H 2 0 2 (R158) and 2 0 H + (M) =i H 2 0 2 ( + M ) . (R86) A t 1000K, the rate of H radical production through R85 is low, as is the branching efficiency of R39, which leads to a reduced effect of hydro-gen on methane ignition. However, the sensitivity to the reaction rate of C H 3 0 2 + H 2 =± C H 3 0 2 H + H (R195) is higher at 1000K due to the increased concentration of C H 3 0 2 . Re-action R195 produces an H radical directly; also the subsequent de-composition of methylhydroperoxide ( C H 3 0 2 H ) species generates an O H radical which is key to intermediate-temperature methane ignition. The rate of R195 is higher than the rate of C H 3 0 2 + C H 4 =± C H 3 0 2 H + C H 3 (R188) 131 Chapter 4. Ignition of Methane with Hydrogen Addition so that the overall production rate of C H 3 O 2 H is higher with hydrogen replacing some of methane. Compared with the branching-dominated promoting mechanism at higher temperature, the low-temperature pro-moting mechanism is much weaker. This accounts for the reduced dif-ference between pure methane and hydrogen-enriched mixtures at low temperatures. It is interesting to note the sign change of sensitivity to reaction H + C H 4 =± C H 3 + H 2 (R54 ) in the two different mixtures. Reaction R54 is usually considered a ma-jor H-scavenging step in methane/hydrogen systems [51,53]. However, in this work, R54 shows a promoting effect on ignition at 1300K when a high fraction of hydrogen presents. This is because the high H 2 con-centration in the system changes the equilibrium of R54 and leads to a reverse of it primary direction. This reaction is beneficial for increas-ing the overall reaction rate in two ways. First, it produces an active H radical that leads to chain branching; second, it reduces the termina-tion efficiency of R159 (combination of two methyl) by competing for C H 3 radicals. It should be noted that the shift in direction of R54 can only occur when the concentration of hydrogen molecule is sufficiently high. This appears to be consistent with the much weaker effect when a smaller fraction of hydrogen is present. By comparing the sensitivity analyses between pure methane and methane/hydrogen mixtures, one can identify several reactions whose uncertainties in their rate coefficients could be responsible for the dis-crepancy between the new model and experimental results. We shift our attention away from those reactions with rates to which pure met-hane or pure hydrogen ignition (for which the delay time is relatively 132 Chapter 4. Ignition of Methane with Hydrogen Addition well predicted) are sensitive and focus instead on the rates that to which only the ignition in methane/hydrogen mixtures is sensitive. The new numerical model over-predicts the reduction of ignition de-lay in methane/hydrogen mixtures at high temperature. It is clear that at 1300K, the ignition delay with hydrogen addition is sensitive to the rate of R54 ( C H 4 + H C H 3 + H 2 ) , but the ignition of pure methane is less sensitive to this reaction. Since the reaction involves both methane and hydrogen, it cannot affect the pure-hydrogen submechanism ei-ther. A review of the literature related to the rate of this reaction shows that the data (experimental or theoretical) reported by different authors vary by a factor of 3, which could account for a change of 20% to 30% in the calculated delay time. A t lower temperatures, the uncertainty in the rate of R195 ( C H 3 0 2 + H 2 =± C H 3 0 2 H + H) could be responsible for the under-predicted effect of hydrogen. There are no experimental measurements of the rate of this reaction in the temperature range of this study, and the data reported in the literature shows an uncertainty of a factor of ten. The sensitivity analysis reveals that the effect of hydrogen on me-thane ignition is closely tied to the production and consumption of hydrogen radicals. In Fig. 4.6, the concentrations of hydrogen radicals during the induction period in two different mixtures are presented at 1000K and 1300K respectively. It can be seen that, for both tempera-tures, the concentration of H radicals increases significantly with hydr-ogen replacing some methane in the fuel; as expected, the difference is more prominent at higher temperatures, where a stronger promoting effect is observed. Further insight into the function of hydrogen can be obtained by tracking the contribution of elementary reactions to the generation and consumption of hydrogen radicals during the induc-133 Chapter 4. Ignition of Methane with Hydrogen Addition 10 ' i i i • i i 0 0.2 0.4 0.6 0.8 Normalized time, tlx Figure 4.6: H concentration during the induction period in different mixtures 134 Chapter 4. Ignition of Methane with Hydrogen Addition 158. H+CH20(+M)<=>CH30(+M 185. OH+H2<=>H+H20 154. H+CH4<=>CH3+H2 170. |H+CH30H<=>CH30+H2 0.5 1 1.5 65%CH4,35%H2,1000K 58. H+CH20(+M)<=>tt;H30(+M) 1168. I H C O + M < = > H + C O + M 0 0.5 1 1.5 100%CH ,1000K 4 158. H+CH20(+M)<=>CH30(+M) 185. OH+H2<=>H+H20 154. H+CH4<=>CH3+H2 1168. HCO+M<=>H+CO+M 0.5 1 1.5 65%CH4,35%H2,1300K 58. H+CH20(+M)<=>CHf30(+M) HCO+M<=>H+CO+M 0 0.5 1 1.5 2 100%CH ,1300K Figure 4.7: Integrated contribution of major H-generation reactions for methane and methane/hydrogen mixture at 1000K and 1300K respec-tively. The analysis was performed at P=40 bar. tion period. A s shown in Figs.4.7 and 4.8, for the temperatures investigated, the third-order reactions H + 0 2 + M =± Ff0 2 + M (R33-R38) are mainly responsible for the consumption of H radicals. But in the methane/hydrogen mixture at 1300K, a relatively higher fraction of hydrogen radicals is consumed in the branching reaction H + 0 2 ^ O + OH, (R39) which indicates that this is the main promoting mechanism. It is also noticeable that, in agreement with the sensitivity analysis, R54 (hydr-135 Chapter 4. Ignition of Methane with Hydrogen Addition 136. |H+02+N2<=>H02+N2 137. |H+02+CH4<=>H02+CH4 134. |H+202<=>H02+02 133. |H+02+M<=>H02+M 0.5 1 65%CH4,35%H2,1000K 1.5 |H+02+N2<=>H02+N2 |H+02+CH4<=>H02+CH4 34. |H+202<=>H02+02 33. |H+02+M<=>H02+M 0.5 1 100%CH ,1000K 4 1.5 136. lH+02+N2<=>H02+N2 139. lH+02<=>0+OH 137. lH+02+CH4<=>H02+CH4 134. lH+202<=>H02+02 133. lH+02+M<=>H02+M 0.5 1 65%CH ,35%H ,1300K 4' 2 1.5 36. lH+02+N2<=>H02+N2 IH+CH4<=>CH3+H2 lH+02+CH4<=>H02+CH4 139. H+02<=>0+OH H+202<=>H02+02 0.5 1 100%CH ,1300K 4 1.5 Figure 4.8: Integrated contribution of major H-consumption reactions for methane and methane/hydrogen mixture at 1000K and 1300K re-spectively. The analysis was performed at P=40 bar. 136 Chapter 4. Ignition of Methane with Hydrogen Addition 10 o Recombination Branching 10 -10 0 0.2 0.4 0.6 Normalized time, t/x 0.8 1 Figure 4.9: Comparison of rates of progress of recombination reac-tions (sum of R33-R38) and the chain branching reaction (R39) in the 0.35H 2 /0.65CH 4 mixture at 40 bar during the induction period. ogen reacting with methane) is a secondary reaction for the consump-tion of hydrogen radicals in pure methane mixtures at 1300K, but ap-pears as a major H generating reaction with the presence of a large ini -tial H 2 concentration. Figures 4.7 and 4.8 suggest that the competition between chain termination (R33-R38) and chain branching (R39) reac-tions is key to the effect of hydrogen in methane ignition. As shown in Fig.4.9, at 1000K, the recombination reactions dominate. As a re-sult, branching through R39 is largely suppressed. A t 1300K, the dif-ference between the rate of chain branching and that of recombination decreases significantly, which leads to a stronger promoting effect. Figure 4.10 presents the main oxidation path during the induction period for the hydrogen/methane mixture based on integral reaction 137 Chapter 4. Ignition of Methane with Hydrogen Addition flow analyses [62]. It shows that the oxidation of hydrogen is mainly through the path H 2 -» H -> H 0 2 H 2 0 2 —> OH —> H 2 0 . Two of the most important intermediate species for methane oxida-tion - OH and H 0 2 - are among the main intermediates in the reaction pathway of hydrogen; however, the presence of hydrogen addition in the system does not change the fundamental reaction path of meth-ane discussed in previous work [28,29], except for an additional C H 3 0 oxidation channel through reaction H + C H 3 O H =± C H 3 O + H 2 . (R66) This new branch, C H 3 O -»• C H 3 O H -> C H 2 O H C H 2 0 , is in parallel with the main oxidation path through H + C H 2 0 ( + M ) =± C H 3 0 ( + M ) (R58) and C H 3 O + 0 2 =± H 0 2 + C H 2 0 . (R171) The significance of this extra path to the overall reaction rate is only secondary, because the oxidation of C H 3 0 is not a main rate-limiting step for methane ignition [28]. Since the third order reactions R33-R38 are more sensitive to the change of concentration than the second-order reaction R39, it is rea-sonable to believe that to reduce the pressure in the reaction system should favor the chain branching and thus more rapid ignition. This is consistent with the experimental observation. A slightly stronger pro-moting effect was observed at 16 bar compared with that at 40 bar. To 138 Chapter 4. Ignition of Methane with Hydrogen Addition OH CH 3 Oy H CH, ( X o2 C H 3 t * — H O . HO, H 2 0 2 H 2 0 2 | O H H 2 i \" H 2 0 OH C H 3 O 2 ^ O : C H 3 H , C H , HO, CH, H, - C 2 H 6 C R O ^ O L O H , OH 3 M OH C H 3 0 2 H C R O A C E O H OH , C H 3 I H T C H O C O Figure 4.10: Integral reaction flow analysis showing main reaction pathway for mixture #2 at 1300k and 40 bar. Smaller-font letters desig-nate main intermediate species participating the reaction steps. 139 Chapter 4. Ignition of Methane with Hydrogen Addition Figure 4.11: Calculated ignition limit of hydrogen for a given ignition delay in this work. The thin dash lines indicate the pressure range stud-ied in this work. clearly illustrate the pressure dependence of hydrogen, a traditional \"ignition l imit\" calculation has been performed for mixture #2 with methane treated as an inert gas (rates of reactions related to methane were set to zero). Figure 4.11 shows the pressure and temperature con-ditions required by the hydrogen-only systems to obtain ignition at a given delay time. The well-defined reversed S curves are observed, which show the (three) classical ignition limits of hydrogen. Figure4.11 shows that, over the experimental range of this work, particularly for ignitions with a delay time below one millisecond, we are primarily at the second ignition limit where the ignition temperature shows a negative pressure dependence. Namely, the third-order recombination reactions (R33-R38) are the rate-limiting steps in this region. This leads to a more discernible difference in ignition delay between the pure me-140 Chapter 4. Ignition of Methane with Hydrogen Addition thane and methane/hydrogen mixtures at the lower pressure. 4.5 Conclusions Shock tube experiments on the ignition of two stoichiometric methane hydrogen air mixtures under high pressure and moderate temperature conditions have been conducted. It has been observed that the pro-moting effect of hydrogen decreases with decreasing temperature. The difference between pure methane and methane/hydrogen mixtures is more prominent at 16 bar than that at 40 bar. A low fraction of hydro-gen addition shows only weak effects on the ignition delay of methane under the conditions explored. A numerical study of methane/hydrogen/air ignition under the ex-perimental conditions has been conducted using a detailed chemical kinetic mechanism. The mechanism was modified to obtain an im-proved agreement between the model and experimental results. The effect of hydrogen on methane ignition was primarily related to the generation and consumption of H radicals. A t high temperatures, the rapid oxidation of hydrogen molecules through R85 (OH + H 2 =± O + H 2 0 ) and the fast branching reaction R39 ( H + 0 2 =± O+OH) are mainly responsible for the stronger ignition promoting effect. The rates of both R85 and R39 decrease rapidly with decreasing temperature. A t lower temperatures, reactions between H 2 and C H 3 0 2 account for a weak ef-fect of hydrogen on methane ignition due to the production of extra H radicals. The effect of hydrogen in the current system exhibits a nega-tive pressure dependence, which implies that the third-order recombi-nation reactions (R33-R38) are rate-limiting steps under the conditions of this work. Further experimental and kinetic studies are necessary to improve 141 Chapter 4. Ignition of Methane with Hydrogen Addition our understanding of the ignition characteristics of methane/hydrogen mixtures under high-pressures and intermediate-temperatures. This should include studies on key elementary reactions such as R54 and R195 as well as ignition studies at other equivalence ratios and with a wider range of hydrogen concentrations. 142 Chapter 5 Evaluation of Detailed Nitrogen Oxides Mechanisms in Methane Flames12 5.1 Introduction Computational studies of the formation of oxides of nitrogen (NOx) in natural gas fueled internal combustion engines using detailed kinetic mechanisms have been reported recently in the literature [162-164]. The agreement between the numerical results and experimental data is sometimes compromised by the under-performance of the underly-ing detailed kinetic mechanism for combustion, which has often not been sufficiently validated under engine-relevant conditions [162]. A l -though the results under certain operating conditions can possibly be improved by empirically adjusting the key reaction rate constants, this approach is clearly questionable given the lack of justification on sci-entific grounds, and could lead to reduced credibility of the detailed mechanism. In our earlier work, we have developed a detailed reaction mecha-nism for the ignition of methane including minor higher hydrocarbons 1 2 A version of this chapter has been submitted for publication. Huang, J. and Bushe, W. K. Evaluation of detailed nitrogen oxides mechanisms in met-hane/ air flames, Submitted to Combust. Flame 143 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames (C2 and C3). The mechanism has been validated against the shock tube ignition data for methane and various synthetic natural gas fuels over a wide range of pressure and temperature conditions [28,165], par-ticularly under conditions relevant to practical combustion devices. Unfortunately, the earlier versions of the mechanism do not contain sub-models for NOx, which are considered important pollutants from combustion processes; given that, in most practical applications, natu-ral gas is burned in air, it is important to include a sub-model for NOx. The objective of the work described in this chapter is to find an appropriate sub-model for NOx to supplement our existing mechan-ism, which would make the resulting mechanism capable of predicting ignition, combustion and pollutant formation for all major gas-phase species considered to be pollutants under engine-relevant conditions. In the work described in this paper, we evaluate detailed NOx mecha-nisms from the literature used in conjunction with our natural gas re-action mechanism in predicting experimental data for NOx formation in methane/air combustion systems. 5.2 Experimental Data Two sets of experimental data on NOx formation have been selected in this study: the measurement of NOx and N 2 0 profiles in a stirred reactor by Steele et al. [66] and the measurement of the NO profile in a methane/air two-stage counterflow flame by L i and Williams [36]. Three reasons are behind the selection of the above two data sets: 1. They represent typical premixed and non-premixed reaction sys-tems. 2. The premixed experiment [66] focused on fuel-lean conditions 144 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames {(f) = 0.41 0.67) with relatively low temperature (T = 1415 -—1845^), conditions which are closely related to those found in H C C I engines. The counter-flow diffusion flame experiment [36] is more relevant to combustion in direct-injection natural gas en-gines. 3. The N O x control techniques studied in both experiments - lean combustion, staged flame and H 2 0 / C 0 2 addition - are among the most important and widely used in practical combustion de-vices. 5.3 NOx M e c h a n i s m s A set of five N O x mechanisms chosen from the literature has been tested in this work. Information regarding these mechanism is pro-vided in Table 5.1 Table 5.1: Selected N O x mechanisms from the literature Source Nitro Componds N O x Related Reactions Glarborgl [74] 17 78 Glarborg2 [84] 24 207 L i and Williams [36] 16 52 GRI-Mech 2.11 [40] 17 102 GRI-Mech 3.0 [41] 19 105 It should be emphasized that the above mechanisms are subsets of the detailed hydrocarbon reaction mechanisms from which they were taken. The intent of this work is to examine their performance when used in conjunction with our mechanism for natural gas combustion, 145 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames not to judge the accuracy of the individual N O x mechanisms, which cannot be justified without referring to their parent mechanisms. 5.4 Results and Discussion A perfectly stirred reactor (PSR) model was used to simulate the exper-iment of Steele et al. [66]. The experimentally measured reactor temper-ature was used as an input to the PSR model. Figures 5.1 and 5.2 show the comparison between the experimental data and model results us-ing the detailed chemistry with the different N O x mechanisms intro-duced above. The mass fraction of N O x is the sum of the mass fraction of N O and that of N 0 2 . It can be seen that the model results using the Glarborg2 mechanism show the best agreement with the experimental data for both N O x and N 2 0 mass fractions. The L i and Williams me-chanism and GRI-Mech 2.11 give reasonable prediction for N O x but underpredict the mass fraction of N 2 0 , while the predictions from the Glarborgl mechanism and GRI-Mech 3.0 deviate significantly from the experimental results for both N O x and N 2 0 . The counter flow diffusion flame was modeled using FlameMas-ter [166], which solves the transformed conservation equations of a two-dimensional potential flow in similarity coordinates. The fuel and oxidizer inlet temperatures are constant at 300K. The fuel stream con-tains methane/air mixtures at different equivalence ratios. Figures 5.3 and 5.4 present the numerical results incorporating different N O x sub-mechanisms compared with the experimental data. The predicted temperature profiles agree fairly well with the ex-perimental measurements. The fact that all the temperature profiles calculated with the different N O x sub-mechanisms collapse into one line shows that the N O x reactions have little influence on the original 146 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames 5.5 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 x 10 • Experimental • +Glarborg] O +Glarborg2 > +Li&Williams 0 +GRI2.11 + +GRI3.0 Hr + -H-1600 1650 1700 1750 T[K] 1800 1850 (a) f = 3.42ms 4 3.5 3 c •2 2.5 x 10 O z • Experimental • +Glarborg' < +Glarborg2 > +Li&Williams • 0 +GRI 2.11 + +GRI 3.0 1.5 0.5 > mm > > . + + 1600 1620 1640 1660 1680 1700 1720 1740 1760 T[K] (b) r = 6.49ms Figure 5.1: NOx prediction in simulated PSR in comparison with exper-imental data [66]. In subplot (a), the average residence time is 3.42 ms; in subplot (b), the average residence time is 6.49 ms. $ 0.517 — 0.673. 147 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames 2.5 x 10 O Z • Experimental • +GIarborg1 < +Glarborg2 • 0 +Li&Williams 0 +GRI2.11 + +GRI 3.0 1.5 0.5 i - -t-32 00 > > ++ CP 00 ++ 1600 1650 1700 1750 T[K] 1800 1850 (a) r = 3.42ms 9 x l 0 \" 6 2|, • Experimental 1.8 - • 4-Glarborg1 -1.6 < +Glarborg2 - > +Li&Williams 0 +GRI2.11 -c 1.4 - + +GRI3.0 -5 Fractio 1.2 + + + + + ,0 Mas; 1 0.8 0.6 0.4 • Z . 0 a • 0 0 < • • 0 f a 0 0.2 t> t> i i i 1> : i i i i i i i i 1600 1620 1640 1660 1680 1700 1720 1740 1760 T [ K ] (b) f = 6.49ms Figure 5.2: N 2 0 prediction in simulated PSR in comparison with exper-imental data [66]. In subplot (a), the average residence time is 3.42 ms; in subplot (b), the average residence time is 6.49 ms. $ « 0.517 — 0.673. 148 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames .x 10 (a) NO (b) Temperature Figure 5.3: NO prediction in a laminar counterflow diffusion flame in comparison with experimental data [36]. In subplot (a), the predicted NO profile; in subplot (b), the predicted temperature profile. The equiv-alence ratio in the fuel stream is 1.5. 149 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames 2 ^ • O Experimental Model(+Glaborg') Model(+Glaborg2) - 9 - Model(+Li&Williams) x [mm] (a) NO (b) Temperature Figure 5.4: NO prediction in a laminar counterflow diffusion flame in comparison with experimental data [36]. In subplot (a), the predicted NO profile; in subplot (b), the predicted temperature profile. The equiv-alence ratio in the fuel stream is 2.5. 150 Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames methane/air reaction system. As for the predicted N O profile, the re-sults from GRI-Mech 2.11 and Glagborg 2 show the best agreement with the experimental data. The other three mechansms tend to overpredict the NO mole fraction in the flame. Since the main reactions in the thermal and prompt. NOx routes have been weir established in the past decades [63,72], it is intuitive to examine the rates of these reactions for the possible cause of the dif-ferent performance among the tested NOx mechanisms. In Table 5.2, a comparison is made for the forward reaction rates of several key NOx formation steps at 1600 and 2000K respectively. A n immediate obser-vation is that the rates of key reactions from the two Glarborg mecha-nisms are very similar; the only significant difference is that the main NO reburn reaction [84] NO + C H 3 = i H C N + H 2 0 ( R R 1 ) in methane/air flames is not involved in the Glarborgl mechanism. Under the experimental conditions of Steele et al. [66], a high concen-tration of methyl radicals exists in the reactor, thus the absencence of key reburn reactions could lead to a significant over-prediction of NO concentration. In the counterflow diffusion flame, the peak flame tem-perature is around 2000K. The reaction rate of the key prompt NO for-mation step C H + N 2 =± H C N + N, ( R T I ) calculated from L i and Williams mechanism and GRI-Mech 3.0 is sig-nificantly higher than the rates from Glarborg2 mechanism and GRI-Mech 2.11, which could contribute to the over-prediction of NO mole fraction by the former mechanisms. 151 Table 5.2: Comparison of key reaction rates at 1600 and 2000K T[K] Reaction 1 3 Glagborgl Glarborg2 Li&Williams GRI 2.11 GRI 3.0 1600 R T I 3.02E+14 3.02E+14 6.99E+0314 3.15E+13 2.41E+13 R T 2 1.42E+12 1.42E+12 1.42E+12 3.54E+11 1.86E+12 R T 3 3.80E+13 3.80E+13 3.80E+13 5.15E+13 2.98E+13 R P I 2.63E+09 1.93E+09 4.36E+09 1.57E+09 3.65E+09 Rp2 8.40E+11 8.40E+11 1.09E+10 6.64E+11 1.22E+12 R R I N / A 8.53E+09 5.27E+09 1.11E+10 1.11E+10 2000 R T I 3.23E+14 3.23E+14 8.52E+05 3.22E+13 2.47E+13 RT2 2.63E+12 2.63E+12 2.63E+12 5.29E+11 3.50E+12 RT3 3.80E+13 3.80E+13 3.80E+13 5.53E+13 3.05E+13 R P I 6.19E+09 9.78E+09 1.74E+10 7.23E+09 1.58E+10 Rp2 2.07E+12 2.07E+12 2.57E+10 1.64E+12 3.00E+12 R R I N / A 2.40E+10 1.45E+10 6.82E+10 6.82E+10 \"Reactions RT1 R T 3 designate thermal NOx route. R T 1 N + NO =^ O + N 2 / R T 2 0 2 + N ;= NO + O, R T 3 N + OH ;= NO + H. Reactions R P i - R P 2 designate prompt NOx route. R P i CH + N 2 =^ HCN + N, R P 2 HCN + O =^ HCO + H. Reaction R R I is a main reburn reaction NO + C H 3 = i HCN + H 2 0 . 14Reverse reaction rate Chapter 5. Evaluation of Detailed NOx Mechanisms in Methane Flames 5.5 Conclusions NOx formation mechanisms selected from the literature in conjunc-tion with a natural gas reaction mechanism from our previous work have been examined to predict experimental NOx data from a stirred reactor [66] and a counterflow diffusion flame [36]. The numerical re-sults using the mechanism proposed by Glarborg et al. [84] and GRI-Mech 2.11 show the best agreement with the experimental data. The reburn of NO by methyl radicals is found to be important in reaction systems at relatively low temperatures such as those in this work. A consistent trend has been identified between higher predicted NO con-centration and higher rate of a key reaction for prompt NO formation, which suggests that the prompt NO mechanism is important under the conditions of this study. 153 Chapter 6 Trajectory Generated Low Dimensional Manifold in C H 4 - Air Combustion Systems15 6.1 Introduction Combustion simulations incorporating detailed chemical kinetic mech-anisms are being increasingly used in studying reactive flow problems; however, for many practical combustion systems, the increase of C P U load as a result of using detailed chemistry could be prohibitively high. A detailed chemical kinetic mechanism for a combustion process typi-cally involves tens or even hundreds of reactive scalars with hundreds or thousands of chemical reactions, each with their own different time scales, which give rise to stiffness in the governing ordinary differen-tial equations (ODEs). To solve such a stiff system of ODEs is very time-consuming since the smallest time scale must be resolved for the nu-merical solution to be stable [57,62]. Therefore, there is clearly a need to reduce the dimensionality and the stiffness in the detailed chemistry to reduced the computatonal time for combustion simulations. The intrinsically low-dimensional manifold (ILDM) method pro-1 5 A version of this chapter has been submitted for publication. Huang, J. and Bushe, W. K. Validation of trajectory generated low-dimensional manif-old in methane/air combustion systems, Submitted to the 31th International Symposium on Combustion 154 Chapter 6. Validation of TGLDMin CR4-Air Combustion Systems posed by Maas and Pope [90] reduces the detailed chemistry from an n-dimensional composition space !R n to a ns-dimensional (ns < n) manif-old governed by the slow processes based on the concept of local time-scale separation. Following Warnatz et al. [62], the governing ODEs of a conserved chemical reaction system (i.g. a constant-pressure, isoen-thalpic reactor) can be written in general vector form f = F ( Y ) , (6.1) where Y denotes an n-element vector of species mass fractions. The local time scales of the reaction system can be obtained by performing a Schur decomposition of the Jacobian matrix as in J = g=VAV, (6.2) where V is the right Schur vector matrix of the Jacobian J ; A is an upper triangular matrix whose diagonal entries are the eigenvalues of J , which represents the inverse of the characteristic time scales; V de-notes inverse matrix of V . The eigenvalues and corresponding Schur vectors in A , V and V are sorted in descending order with the nega-tive eigenvalues of greatest magnitude representing the fastest process at the bottom of A . If there exists a distinct difference between time scales associated with the n s th and ns + 1th process, i.e. || A(ns, ns)\\\\ « \\\\A(na + l,ns + 1)||, the fast processes corresponding to n — ns large negative eigenvalues quickly relax to quasi-steady state, leading to the convergence of the reaction system onto a n s-dimensional manifold. Mathematically, the low-dimensional manifold is expressed as a ns-dimensional subspace of W which satisfies V / F ( Y ) = 0, (6.3) 1 5 5 Chapter 6. Validation of TGLDM in CR^-Air Combustion Systems where V / is an (n — ns) x n matrix taken from the last n — ns rows of V . The physical meaning of Eq. 6.3 is that the projection of the reaction vector along the slow manifold onto the fast subspace vanishes. Eq. 6.3 can be solved by adding ns parametric equations representing system constrains (i.e. element mass conservation, energy conservation, etc.) and a small number of progress variables. The final form of the I L D M equations can be written as where denotes ns parametric equations, and r denotes the values of parameters. A particularly attractive feature of the I L D M is that the species mass fractions and the reaction rates on the low-dimensional manifold can be tabulated prior to the reactive flow calculation; they can be quickly retrieved from the table based on the values of progress variables from the flow calculation to provide closure for the chemical source term. In this way, an orders-of-magnitude reduction in compu-tational time can be achieved, which makes the I L D M method an excel-lent choice for complex turbulent combustion models. Another advan-tage of the I L D M method over other mechanism reduction techniques, particularly the quasi-steady-state-assumption (QSSA) method, is that the reduction can be performed automatically; the modeller is not re-quired to have a priori knowledge of the reaction mechanism. Consid-erable success has been achieved in implementing the I L D M method in direct numerical simulation (DNS) and probability-density-function (PDF) modeling of premixed and non-premixed turbulent flames [167-169]. A variety of extensions of the original I L D M model to address the dynamic effects from the flow field on the construction of manifold has also been reported in the literature [170-172]. (6.4) 156 Chapter 6. Validation of TGLDM in C H V A i r Combustion Systems Nafe and Maas [92] discussed several issues associated with the original I L D M model. Besides the numerical difficulties associated with obtaining a converged solution using Eq. 6.4, the I L D M calculated us-ing Maas and Pope method represents an inertial manifold only when the time scale separation approaches infinity. Hence, the reaction vec-tors along the I L D M are mostly not in the tangent plane of the ma-nifold. Singh et al. [173] used the analysis of Davis and Skodje [174] to show that the I L D M of a reaction system is only an approxima-tion of the more fundamental slow invariant manifold (SIM). A n im-provement of the I L D M using a modified Fraser algorithm proposed by Davis and Skodje [174] has been studied by Nafe and Maas [92]. In that case, the I L D M calculated using the method of Maas and Pope [90] was used as an initial guess for the calculation of the SIM for the same system. A n alternative approach to obtain an inertial manifold is to use the Trajectory Generated Low-Dimensional Manifold (TGLDM) method developed by Pope and Maas [94]. The T G L D M method generates the manifold by simply integrating Eq.6.1 from selected initial states in the composition space, therefore a converged solution is guaranteed. Since the manifold is generated along the trajectory, the reaction vector is al-ways in the tangent plane of the manifold, which eliminates the need for back projection as required by the original I L D M method [94]. The-oretically, if an attractive manifold does exist in a dynamic system, the trajectory whose initial state is on the manifold stays on the manifold. Thus, such a trajectory w i l l be equivalent to a SIM. Therefore, for the T G L D M method, the global optimization problem is reduced to one of local optimization, i.e. to locate the optimal initial state for each trajec-tory. 157 Chapter 6. Validation of TGLDM in CH 4 -Air Combustion Systems For a conserved reaction system, the initial state exists in a subspace of 5Rn defined by the linear equations KID™=' = = = i - n ' - (6-5) where M(i,j) denotes the mass fraction of the ith element in the jth species; P denotes a coefficient matrix for the parametric equations; Yc(i) denotes the mass fraction of ith element, Yp(k) denotes the as-signed value of the fcth parameter; ne denotes the total number of el-ements in the system; ns denotes the total number of species; np de-notes the total number of parameters, ie. the dimension of the T G L D M . The attractive point is naturally the equilibrium state in this subspace. Such a constrained equilibrium state can be calculated readily using the method of Lagrange multipliers [175,176] or Gibbs function con-tinuation [177]. Thus, the constrained-equilibrium state is proposed in this work as the initial point for the trajectories in building a T G -L D M . Figure 6.1 shows a comparison between the resulting T G L D M and an I L D M calculated using the Maas and Pope [90] method for a C O / H 2 / a i r reaction system. In this case, the agreement between the two manifolds is very good. For trajectories starting from a region where the slowest chemical time scale in the fast subspace is comparable to or longer than the physical time scale - ie. during the induction period of a fuel jet in IC engines - constrained equilibrium or the original I L D M method can lead to significant transient state error [173]. The error can be reduced by using a higher dimensional manifold, but the task of storing a high dimensional table or switching I L D M dimensions during the flow cal-culation is very complex. A compromise is to use T G L D M with ini-tial conditions that resemble the possible states in the flow field, for 158 Chapter 6. Validation of TGLDM in CH 4-Air Combustion Systems Manifold with MP method Manifold with current method Reaction trajectories 4 c|) c o^mol/kg] 6 8 (a) H 2 0 - C 0 2 4 4>co5[mol/kgf (b) H - C 0 2 Figure 6.1: Comparison of a one-dimensional T G L D M with optimal initial states and a conventional I L D M . The manifold is projected onto H 2 0 - C 0 2 and H - C 0 2 planes; YU,H2O and Yco2 > Yu,co2> where YUtH2o and YUjCo2 designates the mass fraction of H 2 0 and C 0 2 in the unreacted mixture. This guarantees that the trajec-tories starting from the initial mixing states are included in the manif-old. The reaction rates were calculated using Chemkin II [147], and the governing O D E systems were solved using a stiff O D E solver [148]. A s pointed out by Pope and Maas [94], the procedure of construct-ing a T G L D M introduces a natural parameterization, i.e. the trajecto-ries can be defined using the normalized length rt and angle of the initial point, 9t, wi th respect to the equilibrium state. The domain of the T G L D M in parameter space thus becomes a circle with unit ra-dius. Each trajectory is transformed into a straight line starting from the edge of this circle and evolving towards the center, which is the equilibrium state. The implementation of the manifold parameterized 160 Chapter 6. Validation of TGLDM in CEU-Air Combustion Systems in this way is still difficult, particularly when it comes to transform-ing a perturbation in the physical coordinates defined by the progress variables to the rt — 9t coordinates. This is because, in the physical coor-dinates, the trajectories bunch as they approach the equilibrium state, but in the rt — 6t coordinates, they remain separated. This leads to a singularity in the transformation matrix - the transformation from a low dimensional space to a high dimensional space is not uniquely de-fined. In this work, two progress variables, Yco2 a n d YH2O, have been used as the table entry directly rather than parameterizing the T G L D M in transformed coordinates. The perturbation and reaction vectors are projected onto the same plane. Since there is no parameterization of the T G L D M , the table lookup is realized through a subprogram which performs a Delaunay triangulation, an interior point search and a sur-face interpolation on a two-dimensional unstructured grid [178]. Fig-ure 6.2 presents an example of the triangulated T G L D M projected into the YCo2 — YH2O plane. One advantage of the T G L D M is that the den-sity of mesh varies naturally with the change of stiffness in the reaction system since the stiff O D E solver used to generate the trajectories can respond to such a change by providing more solution points locally. Thus the mesh has an adaptive resolution which provides a higher ac-curacy in terms of rate interpolation. 6.3 Results and Discussion The performance of the T G L D M in three different methane/air reac-tion systems has been investigated: a unstrained premixed laminar flame, a non-premixed flamelet, and a perfectly stirred reactor. De-pending on the reaction system, the savings in C P U time by using the T G L D M method range from two to three orders of magnitude when 161 Chapter 6. Validation of TGLDM in CH-^-Air Combustion Systems 0.14 C02 Figure 6.2: Delaunay-triangulated T G L D M in Yco2 ~ YH2O plane. The unburned C H 4 / a i r mixture temperature (at the lower-left corner) is 400K, P = latm, and $ = 1. The mesh density is reduced for the pur-pose of demonstration. 162 Chapter 6. Validation of TGLDM in CH 4-Air Combustion Systems compared with the method of direct integration using the detailed ch-emistry. 6.3.1 Unstrained Premixed Laminar Flame The governing equation for the species mass conservation of a one-dimensional, premixed, free laminar flame is given by dYi dYi d pLW+pu^i = - d x - { p Y i V i ) + ^ ( 6 - 6 ) where u designates the bulk convection velocity, V designates the dif-fusion velocity, and Co designates the rate change of mass fraction due to chemical reactions. For the mixture averaged diffusion velocity Vit which incorporates both mass diffusion and thermal diffusion veloci-ties, one can assume Fickian diffusion [57]: V; = - i A v X , - | | v T . (6.7) where D is the mixture averaged diffusion coefficient, X is the mole fraction, and DT is the thermal diffusivity at temperature T. The energy conservation equation is given by + P U C ^ = Yx oo. (6.10) ox ox The subscript u designates the unburned mixture. The above govern-ing equations are solved using a Strang operator splitting [179] scheme 163 Chapter 6. Validation of TGLDM in CH 4-Air Combustion Systems that separates the solution of chemical source terms from the convecti-on/diffusion term. Each numerical time step is split into a convection-/diffusion time step and a chemical time step. In the physical time step, the convection/diffusion terms in the governing PDEs are first solved with the reaction source term suppressed. In the chemical time step, the flow field obtained from the physical time step is frozen and each computational cell is treated as being spatially homogeneous; the O D E system that governs the chemical reactions is then solved. Since the T G L D M is an inertial manifold, the solution in the chemical time step stays on the manifold so that no projection is required. The physical perturbation caused by the convection/diffusion terms, on the other hand, is not confined by the low-dimensional manifold. In this work, orthogonal projection [94] is used to return the solution to the manifold at the end of the physical time step. Figure 6.3 compares the steady-state solutions of the flame structure for a premixed methane/air laminar flame obtained using detailed ch-emistry and the T G L D M method respectively. For the profiles of tem-perature and major species, the agreement between the T G L D M and detailed chemistry is very good. For minor species, the peak hydro-gen radical mass fraction in the flame front is over-predicted by about 15 percent, while the peak O H radical mass fraction is over-predicted by around 65 percent. One factor that could contribute to the over-prediction of these minor radicals is their high diffusivity, which leads to a higher mass-transport rate from the flame. In the current T G L D M method, differential diffusion effects are cannot be accounted for, as was mentioned earlier. Away from the flame front, the agreement be-tween T G L D M and detailed chemistry for minor species improves qu-ickly as the burned mixture approaches the equilibrium state. 164 Chapter 6. Validation of TGLDM in CH4 -A i r Combustion Systems • T G L D M C0 2 • T G L D M H z O T G L D M C H 4 • T G L D M 0 2 • T G L D M T — Detailed Chemistry 0.25 ^ 0.2 0.15 0.1 0.05 0 0 W t M M M 3 x [mm] (a) Temperature and Major Species * TGLDM H • TGLDM OH • TGLDM CH — Detailed Chemistry (b) Minor Species Figure 6.3: Comparison of calculated flame structure of a one-dimensional premixed laminar methane/air flame. The flame prop-agates from right to left into an unburned C H 4 / a i r mixture at 400K. P = latm, $ = 1. 165 Chapter 6. Validation of TGLDM in CH 4-Air Combustion Systems 150 - - T G L D M — Full Chemistry A Experimental 100 E o • 50 • ) i . . , > , . _ .00 350 400 450 500 550 600 T[K] Figure 6.4: Laminar flame velocity of a stoichiometric methane/air mixture calculated using T G L D M method. The experimental data shown are those of Warnatz et al. [62]. Figure 6.3.1 shows the laminar burning velocities calculated using the T G L D M for the same mixture as a function of unburned mixture temperature. It can be seen that the agreement between the detailed chemistry and T G L D M is satisfactory for the conditions investigated and both give reasonable predictions of the experimentally measured values. 6.3.2 Laminar Flamelet The basic flamelet equations [96] in mixture fraction space are d Y xd2Yi (6-11) 166 Chapter 6. Validation of TGLDM in CR^-Air Combustion Systems where Z denotes the mixture fraction, and x denotes the instantaneous scalar dissipation rate, which is defined by X - m ^ f , (6.12) where Xj denotes the direction normal to the flame interface. The scalar dissipation describes the mixing time scale in a diffusion flame. A larger value of x indicates more rapid turbulent mixing. Law and Chung [102] and Peters [96] proposed using the scalar dissipation at the stoichio-metric surface, Xst = x(Zst), as the representative parameter for diffu-sion flames. From the analytical solution of a steady-state counter flow configuration, it can be shown that X(Z) = -^p{-2[er fc - 1 (2Z)] 2 } , (6.13) where erfc - 1 designates the inverse of the error function. The func-tional dependence of x(Z) on Xst can be derived from Eq. 6.13 x(Z)=Xstf(Z)/f(Zst), (6.14) where f(Z) designates the exponential term in Eq. 6.13. Unlike the premixed flame, which corresponds to a single manifold at one equivalence ratio, the solution of the flamelet equations requires a complete spectrum of low-dimensional manifolds that spans the full range of mixture fraction. In this work, a discrete set of T G L D M s were generated at fifty two points in mixture fraction from zero to unity. The distribution of node points was made dense in the lean and stoichio-metric regions where the gradients in reactive scalars are high. Figs. 6.5 and 6.6 show a comparison of steady-state flamelet solu-tions obtained using the detailed chemistry and the T G L D M at Xst = 100s _ 1 and Xst = 300s - 1 respectively. The oxidizer temperature at the 167 Chapter 6. Validation of TGLDM in CH4 -A i r Combustion Systems left boundary (Z = 0) is 1400K, and the fuel temperature at the right boundary (Z = 1) is 300K. The pressure is 30 bar throughout the do-main. The two values of scalar dissipation selected represent moder-ate perturbation and high (near-quenching) perturbation cases for this system. In general, for Xst = 100s - 1 , the predictions of temperature, major and minor species profiles using the T G L D M method agree well with the results from the detailed chemistry, but some deviations of temper-ature and methane mass fraction on the rich side of stoichiometric can be observed. For mixtures with high equivalence ratios, the chemical time scales associated with the consumption of methane are long due to the relatively low mixture temperature and the high methyl concen-tration. The reduction in the separation of chemical and physical time scales could be a major factor that leads to the under-prediction of me-thane by the T G L D M . The improved agreement of minor species in the flame zone compared to that in the premixed flame is partially at-tributed to our having neglected differential diffusion: the same scalar dissipation used for all the species in the flamelet equations. For the high scalar dissipation case at 300s - 1 , an increase in the de-viation between the T G L D M results and those from detailed chemistry is observed, although the overall agreement remains reasonably good. In particular, the peak flame temperature and peak C 0 2 mass fraction is slightly under-predicted. It appears that the flame profile obtained from the detailed chemistry is more \"resistant\" to a high level of per-turbation than that obtained from the low-dimensional manifold. 168 Chapter 6. Validation of TGLDM in C H 4 - A i r Combustion Systems 0.25 0.15 0.05' TGLDM C 0 2 TGLDM H 2 0 TGLDM 0 2 TGLDM CH, 4 TGLDM CO TGLDM T Detailed Chemistry 3000 2500 2000 & 1500 1000 0.2 (a) Temperature and Major Species 0.01 0.009 0.008 0.007 0.006 JH\" 0.005 0.004 0.003 0.002 0.001 o TGLDM H - > TGLDM OH < TGLDM CH 2 0 - 0 TGLDM H0 2 ^ > / 0 H Detailed Chemistry - J r --$*% \\> HO x 150 z \\o \\ / f \\o \\ / AL \\ 0 \\ / • -M V / \\ CH Ox20 / mm $• Hxio w v y< < / \\ A 7 / 4> \\ , 9 ^ f N Y < \" — — • — —S * i 0.05 0.1 z 0.15 0.2 (b) Minor Species Figure 6.5: Comparison of calculated temperature and species profile of a steady state methane/air flamelet at 30 bar in the mixture fraction space. Xst = 100s - 1 169 Chapter 6. Validation of TGLDM in CH 4-Air Combustion Systems 0.25 0.15 0.05' TGLDM C 0 2 TGLDM H 2 0 TGLDM 0 2 TGLDM CH, 4 TGLDM CO TGLDM T Detailed Chemistry 3000 2500 2000 £ 1500 0.2 1000 (a) Temperature and Major Species o TGLDM H > TGLDM OH 4 TGLDM CH 20 TGLDM HO, Detailed Chemistry (b) Minor Species Figure 6.6: Comparison of calculated temperature and species profile of a steady state methane/air flamelet at 30 bar in the mixture fraction space. Xst — 300s\"1 170 Chapter 6. Validation of TGLDM in CE.4-Air Combustion Systems 6.3.3 Perfectly Stirred Reactor Using a T G L D M to model transient process such as ignition and flame quenching is more challenging than the previous two cases tested be-cause more chemical time scales are in a range comparable to the phys-ical time scales under these conditions. The perfectly stirred reactor (PSR) model was chosen to test performance under transient condi-tions because, in the PSR, the physical time scale can be specified ex-plicitly. The governing equation for a perfectly stirred reactor is given by dY _ = F ( Y ) - W ( Y - Y , ) , (6.15) where Y j designates the mass fraction vector of the flow entering the PSR, and cu is the inverse of residence time which represents the physi-cal time scale. In this work, the PSR model for a stoichiometric methane-/air system at 30 bar constant pressure was investigated. Figure 6.7 presents the change in temperature and species mass fraction during the induction period of the reactants in the PSR. Two inflow temperatures, 1229K and 1024K, were examined. They are equiv-alent to the temperatures when stoichiometric mixtures are formed with methane at 300K and air at 1400K and 1150K, respectively. The traces marked cu = 0 represents an un-perturbed reaction system. The slight phase difference between the results with the T G L D M and those from the detailed chemistry in this case is mostly due to interpolation errors in getting the instantaneous reaction rate. When an inflow of un-reacted mixture is introduced as a perturbation, the ignition delay time is retarded. For a moderate perturbation, the phase shift predicted by the T G L D M method is in reasonable agreement with detailed chemi-stry. For instance, at an inflow temperature of 1229K, the error in the 171 Chapter 6. Validation of TGLDM in C H 4 - A i r Combustion Systems predicted phase shift is about 10 percent for a retard in ignition delay of nearly 100 percent. The error is larger for the case with lower in-flow temperature. Reasonable predictions of phase retard are only ob-tained for relatively small perturbations, as shown in Fig.6.7(a). Note that the value of u in Fig.6.7(a) is significantly smaller than that of that in Fig.6.7(b). This is clearly due to the larger number of long chemical time scales associated with lower temperatures which leads to a sig-nificant deviation of reaction trajectories subject to large perturbations from the low-dimensional manifold. A combustion system can be quenched when the physical pertur-. bation is such that the chemical reaction rate cannot be sustained due to an excessive heat and radical loss. Figure 6.3.3 shows the prediction of temperature and species mass fraction history during a quenching process using the T G L D M method. In simulating the quenching pro-cess, an unperturbed reaction system is first evolved to the equilibrium state; the perturbation term is then introduced with a large value of to. It can be seen that the calculated results using the T G L D M method essentially match those obtained with detailed chemistry. However, it has been noticed in this work that the performance of T G L D M in the near quenching region is not ideal. In particular, the current T G L D M tends to under-predict the quenching limit: the model employing the T G L D M w i l l predict quenching at a lower level of perturbation than is observed with detailed chemistry. Future work on improving the per-formance of manifold methods in the region of competing chemical and physical time scales is clearly warranted. 172 Chapter 6. Validation of TGLDM in CH4 -A i r Combustion Systems 0.05 r x 10 (a) T in=1024K 0.1 0.05 1 1 1 1 co=0 s\"1 \\ / \\ V T co=3xl04 s\"1 \\ I \\ 1 1 ~7~~f OHxlO -Detailed Chemistry - - T G L D M 3000 2500 2000 Hi 500 1000 £ 500 0 0.5 1.5 t[s] 2.5 x 10 (b) T in=1229K Figure 6.7: Calculated temperature and species mass fraction history during the induction period for a stoichiometric methane/air mixture in a PSR. The solid line designates results from detailed chemistry; the dash line designates results from the T G L D M . 173 Chapter 6. Validation of TGLDM in CH 4 -Air Combustion Systems Figure 6.8: Calculated temperature and species mass fraction history during a quenching process of a stoichiometric methane/air mixture in a PSR. The inverse of the characteristic residence time, cu, is set at 6 x 10 6 s _ 1 ; the inflow mixture temperature is 1024K. 174 Chapter 6. Validation of TGLDM in CH 4 -Air Combustion Systems 6.4 Conclusions A n improved T G L D M method [94] that uses optimized initial states for the trajectories has been discussed. The difference between the result-ing T G L D M and the I L D M generated using the method of Maas and Pope [90] has been investigated. Several merits of the T G L D M method have been identified: the construction of the manifold is less mathe-matically involved than the conventional I L D M method; the reaction vector of a T G L D M is always in the tangent plane of the manifold. A n alysis of a hypothetical reaction system showed that if the initial state of the reaction system is known, the T G L D M with optimized initial conditions gives a better approximation of detailed chemistry than the I L D M . The performance of the T G L D M in predicting reactive scalar pro-files for a methane/air mixture has been investigated in modeling lam-inar premixed flames, laminar flamelets, and the transient processes in a perfectly stirred reactor. In general, reasonable agreement between the results with the T G L D M and those from detailed chemistry has been found. When the physical time scales approaches the chemical time scales in the fast subspace, relatively large errors occur using the T G L D M . Further study on improving the performance of manifold methods in these regions has been suggested. 175 Chapter 7 Simulation of Turbulent Reactive Methane Jets Under Engine Relevant Conditions 7.1 I n t r o d u c t i o n Direct-injection natural gas (DING) engines offer a combined bene-fit of low pollutant emissions and high thermal efficiency [4-8]. The combustion process in a D I N G engine, particularly that of a transient turbulent natural gas jet, is a subject under intense study. Thanks to progress in the study of the chemical kinetics that are important to natural gas combustion at elevated pressures and reduced tempera-tures, the combustion chemistry for natural gas under engine relevant conditions is now better understood [28,29,38]. However, the imple-mentation of detailed chemistry in a multi-dimensional simulation of a turbulent reactive natural gas jet remains a challenging task. In par-ticular, two fundamental problems have to be addressed in order to simulate turbulent combustion using detailed chemistry. First, for any computational fluid dynamic (CFD) models that do not resolve all of the turbulence scales, a turbulent combustion model is required to ac-count for the effects of turbulent fluctuations on the chemical reaction rates. Second, the combustion system described by detailed chemistry 176 Chapter 7. Simulation of Turbulent Reactive Methane Jets is usually very stiff: the chemical time scales associated with different reaction scalars vary drastically. To solve such a stiff system directly in most practical turbulent reactive flow simulations is still beyond the reach with the existing computational power. Thus methods to reduce the CPU time required for computing the reaction rates from detailed chemistry must be used. In this work, a combustion model is proposed that provides solu-tions to the above two problems. The model incorporates the Condi-tional Source-term Estimation (CSE) [127] and Trajectory Generated Low-Dimensional Manifold (TGLDM) [94] methods, which can pro-vide a closure for the chemical source term at the level of the first mo-ment with a relatively low computational cost. The model is applied to simulate the ignition and combustion of experimental methane jets under engine relevant conditions. The simulation results are validated using experimental data obtained from a high-pressure shock tube fa-cility as well as data reported in literature. 7.2 C S E - T G L D M Combustion Model 7.2.1 Conditional Source-term Estimation The Conditional Source-term Estimation (CSE) method [127] seeks clo-sure of the chemical source term using the conditional average of the reactive scalars in a manner similar to the first-order Conditional Mo-ment Closure (CMC) method. The solution for the conditional aver-ages, however, is not obtained by solving transport equations as is done in CMC. Based on the assumption that the conditional averages of scalars have very low gradients in space, Bushe and Steiner pro-posed a method to obtain the conditional averages of reactive scalars by inverting an integral equation for an ensemble of discrete points on 177 Chapter 7. Simulation of Turbulent Reactive Methane Jets a computational grid [127]. Mathematically, the unconditional average of a random event can be calculated from its conditional average if the local probability density function (PDF) is known. {Y(x,t))= f' P(x,t;r))(Y(x,t)\\V)dV. (7.1) Jo For a selected spatial ensemble of N points, it is assumed that (Yix^fM^iYlr))^ (7.2) where the superscript n denotes the nth point in the ensemble; the sub-script A indicates that the average is taken over the ensemble of A dis-crete grid points in space. Thus a discrete set of TV integrals {Y(x^,t)}= f P(x,t;V)(Y\\r1)AttdV. (7.3) Jo By approximating the integral using a numerical quadrature, one ob-tains a linear system of equations which can be inverted for the con-ditional average of the scalar of interest using a method such as linear regularization. The CSE method was initially implemented in a Large Eddy Simu-lation of a piloted methane/air diffusion flame with encouraging suc-cess in predicting experimental measurements [128]. Later, the concept of CSE was tested in conjunction with the unsteady laminar flamelet model, in which the conditional averages of reaction scalars were cal-culated using a linear combination of flamelet solutions [129]. In this incarnation of CSE, one decomposes the conditional average into a linear combination of laminar flamelets; one then inverts an integral equation to determine the optimal linear combination of flamelets for each ensemble. The method has been used in the context of RANS mo-del to study turbulent methane jet ignition under engine-relevant con-ditions with some success [130,131]. In order to address the issue raised 178 Chapter 7. Simulation of Turbulent Reactive Methane Jets by ill-posedness in Eq. 7.3 as well as to provide temporal continuity in the solution, Grout [131] proposed using a general Tikhonov regular-ization method for the inverting process, given by mm{\\\\n(Y\\VY - {Yf\\\\ + A||(Y|T7)* - W ^ ' H } , (7.4) where Q is the original coefficient matrix for the discrete integral equa-tion; the superscripts t and t — At are the times when the scalars are evaluated; A is a weighting coefficient specified by the modeler. The regularization term A|| (Y\\r))t—(Y\\r])t~At\\\\ was added to limit the change of conditional average between two consecutive time steps as well as to stabilize the solution. In Grout's work, A was set to the ratio of the trace of Cl to the trace of an identity matrix with the same size. In this work, efforts have been made to further improve the reg-ularization method by including spatial continuity for the conditional scalar field. The incremental limiter, || ( Y ^ ) * — ( Y|77) t _ A i || in Eq. 7.4 is re-placed by ||(Y|?7)f — (Y*|T7)*|| where (Y*|r7)* is calculated from the trans-port equation {Y*\\vy - (Y\\r,y-*t ^ dpjuiMWv)*-\" ~~At ^ • ( 7 - 5 ) It has been reported from experimental studies and D N S simulations that for steady, axi-symmetric jet flames, the cross-stream variations of conditional means are not significant [127,180,181]. Thus the convec-tion of the conditional means described in Eq. 7.5 is only considered in the axial direction of the jet. The conditional mean velocity at the ax-ial location x, {ux\\rf), is obtained by cross-stream averaging along the isopleth , , v _ JQR(u(x,r))P(mx,r)dr {Ux\\V) JR— (7.6) J 0 P(r};x,r)dr where R denotes the radius of the jet. 179 Chapter 7. Simulation of Turbulent Reactive Methane Jets The main advantage of the CSE method is that the computational cost is substantially lower than that of C M C . In addition, it does not involve constraining assumptions such as those employed in laminar flamelet models, and is thus applicable to a wide range of turbulent flames. 7.2.2 Trajectory Generated Low-Dimensional Manifold Manifold methods for reducing detailed chemistry are based on the separation of chemical time scales associated with different reaction scalars. If the time scale separation is large enough, fast processes with short time scales approach a quasi-steady state rapidly; these can be decoupled from slow processes to reduce the total dimensionality of the reaction system. The remaining low-dimensional manifold can be used to approximate the detailed chemistry with a high degree of accu-racy. For a two-dimensional manifold, for example, the instantaneous rates of reaction scalars Y can be obtained from the manifold using the formula dY dY(u,v)du dY(u,v)dv ~dt= Wu dt + dv dt ( ' ' where u and v are progress variables used to parameterize the manif-old. Maas and Pope [90] proposed a mathematical model for comput-ing the Intrinsic Low-Dimensional Manifold (ILDM) by minimizing the reaction vector projected into the fast subspace, which is defined by eigenvectors associated with large negative eigenvalues of the Ja-cobian matrix. The manifold generated by this method is somewhat optimal globally; however, the implementation of the method is very involved. In the T G L D M method proposed by Pope and Maas [94], the man-ifold is generated along reaction trajectories so that the construction of 180 Chapter 7. Simulation of Turbulent Reactive Methane Jets the manifold is significantly simpler than that of the I L D M . The bound-ary formed by the initial states of the trajectories, which is called the manifold generator, can be obtained using the extreme-value-of-major-species method [94] to achieve a maximum overlap between the TG-L D M and I L D M . It should be note that the initial states of T G L D M s can be tailored to match the initial mixture compositions in the flow field to provide a better approximation of detailed chemistry. The parameterization of the T G L D M can be realized using the nor-malized trajectory length and the initial locations of the trajectory with respect to some reference. However, in locations where the reaction trajectories bunch, the projecting matrix which maps the perturbation from physical space to the manifold space becomes nearly singular. To avoid this problem in this study, two reaction scalars, YCo2 a n d YK2o, are used as progress variables for the manifold without parameteriza-tion. The projected T G L D M in the YCo2 — ^ H 2 O plane is triangulated using the Delaunay method to form an unstructured mesh. A subpro-gram is then used to perform the interior point search and interpola-tion on the manifold surface based on the instantaneous value of YCQ2 and Y H 2 O -7.3 Experiments A n experimental study of the ignition and combustion of a transient turbulent methane jet has been conducted in a shock tube facility; de-tailed experimental results have been reported by Sullivan et al. [182]. A schematic of the experimental setup used in this study is shown in Fig. 7.1. The shock tube has a circular cross section with an inner diame-ter of 59 mm. The lengths of the driver and driven section are 3.11 m 181 Chapter 7. Simulation of Turbulent Reactive Methane Jets Data Acquisition System Driver J u u Double Diaphragms >40cm He and Air Injector i Contact Surface ' I (During Experiment) 5.9 cm Figure 7.1: Schematic of the shock tube facility 182 Chapter 7. Simulation of Turbulent Reactive Methane Jets and 4.79 m respectively. The incident shock velocity is measured using five flush mounted PCB pressure transducers along the driven section. The minimum response time of these transducers is 2 ps. The signal from the pressure transducers were acquired using a Wavebook 512 data acquisition system sampling at a rate of 140 k H z per channel. The temperature behind the reflected shock is calculated based on normal shock relations using the measured incident shock velocity [27]. The uncertainty in the experimental temperature, calculated based on the uncertainty in the measured shock velocity, is around 14K. A stainless steel section equipped with three 199 x 15 mm fused quartz windows was attached to the end of the driven section to pro-vide optical access to the experimental area. Methane (99.97% purity) was injected into the shock tube along its axis through a 1.1 m m diame-ter single nozzle natural gas injector equipped with a magnetostrictive actuator. The response time of the injector is around 250 ps with a stan-dard deviation of 20 ps. The injection was initiated between 100 and 800 ps after the shock reflection by a customized controller unit trig-gered using the pressure signal sensed at the end plate of the driven section. The total injection duration was 1.5 ms. A Vision Research Phantom high speed C M O S digital camera operated at a rate of 31000 fps was used to capture broad band luminosity during the ignition and combustion process of the methane jet behind the reflected shock wave. The location of initial ignition kernel is defined as the emergence of a non-contiguous flame region, which is able to develop into a fully fledged jet flame. The ignition delay is thus defined as the time delay from the start of gas injection to the emergence of the ignition kernel. Tailored interface conditions [27] were used to obtain a maximum experimental time of around 5 ms. The calculated contact surface po-183 Chapter 7. Simulation of Turbulent Reactive Methane Jets sition when stopped by the reflected shock wave was more than 40 cm away from the end plate of the driven section, which is considerably longer than the maximum penetration length of the methane jet dur-ing the experimental period. This ensured that the temperature, pres-sure and gas composition during the experiment were not affected by the driver gas. A t the end of the experimental time, the reaction in the experimental region was quenched due to the arrival of the rarefaction wave that rapidly reduces the temperature and pressure in the test sec-tion. After each experiment, the complete contents of the shock tube were vented through an impacter-type filter into a 400 L carbon-impregnated polyolefin sampling bag. The sample was then analyzed using an A P I 200E Chemiluminscent NOx analyzer to obtain the concentration of NO and N 0 2 in the combustion products. Table 7.1 summarizes the experimental conditions and main param-eters in these experiments. Table 7.1: Experimental conditions of this work No. Runs dinjector [mm] Piie [bar] P c 1 7[bar] Ti 1 8[ms] % 1 9 r o 2 0 K 28 1.1 75 30 1.5 300 1150-1400 16Injection pressure 1 7Back pressure (pressure behind reflected shock) 18injection duration 1 9Fuel temperature 20Temperature behind reflected shock 184 Chapter 7. Simulation of Turbulent Reactive Methane Jets 7A Model Formulation and Validation 7.4.1 C F D Model Formulation For the experimental conditions just described, the jet is choked at the nozzle exit. The flow field in the region close to the nozzle exit has a high Mach number and thus need to be treated as being fully com-pressible. The Reynolds averaged transport equations for mass, mo-mentum and energy in the cylindrical coordinates with an axisymmet-ric configuration are Continuity Symbol ~ denotes Favre average and p is the mean of density. Momentum / dur _ dur _ dur \\ dp f i d , . drrz \\ _ _ _(duz duz duz\\ dp (Id drzz\\ P \\ - d J + U ^ + U ^ ) = - ^ [ r d - z ^ + -dz-) ' ( 7 - 1 0 ) where r designates the stress tensor whose components are given by r r r = ^ [ 2 ^ - ? ( V - v ) ] (7.11) r , , = * [ 2 ^ - f ( V . v ) ] q- — ~- — l?duz dur] rzr — rrz — pt z— h TT~ L dr dz 185 Chapter 7. Simulation of Turbulent Reactive Methane Jets Energy „ {df „ df „ df\\ {dp _ dp „dp\\ ^ U + ^ * + ^ ) \" U + + ^ J (712) i d {dr\\ d 2n Lr/<%.\\2 t u r \\ 2 /<%:x2' = A r dr \\ dr J dz2 + Pt<2 dr J V r j V dz {dur duz\\2 2 / 1 <9 . . duz\\2\\ + {-Y + -dV) - 3 { r d - r i r U r ) + ^z-)l+UJ Body force and radiation effects have been neglected in these equa-tions. A l l the chemical species are assumed to diffuse at the same rate and the Lewis number is unity. The turbulent viscosity pt is calculated using the standard k — e model [183] k2 Pt = pCu —, (7.13) e where k is the turbulent kinetic energy; e is the dissipation rate of tur-bulent kinetic energy; the value of coefficient Cu is 0.09. 7.4.2 Combustion Model Formulation In the CSE method, a probability density function for mixture frac-tion is constructed from its local mean and variance. The closure for transport equations for the mean and variance of mixture fraction is achieved by employing a gradient transport hypothesis: • Mean mixture fraction J,(Z\"2) PDF co) H,o) q CSE (yco)n),(YHJn) 1 (1200K) / the H -atpm abstraction from ethane and the subsequent decomposition of ethyl radical is very efficient in producing H atoms, which enhance the chain branching and accelerate methane ignition. The promoting effect of ethane and/or propane at low temperature (T<1100K) is attributed to the increasing significance of methylperoxy and methylhydroperoxide chemistry. In particular, additional ethane and propane was found to increases rate of production of methylhy-droperoxide (CH3O2H), which subsequently decomposes to generate an O H radical. Since the depletion of O H radical is the main rate-limiting mechanism for methane ignition at temperature below HOOK, the contribution from the above mechanism is more prominent at low temperatures. To further examine the kinetic interaction between methane and higher alkanes in this study, an analytical solution was obtained for the ignition delay of the test mixtures by systematically reducing the skeletal mechanism identified from the reaction-flow analysis using the steady state assumption. The analytical model agrees well with the detailed-kinetic model for both ignition delay and the concentrations of main intermediate species. Both the analytical model and kinetic analysis of the detailed mechanism show that the addition of ethane-/propane does not change the main reaction path of the methane sys-tem. The promotion of ignition is realized through accelerating the ini-216 Chapter 8. Conclusions and Future Work nation phase in the induction period. Another fuel of interest, whose ignition mechanism was studied in this work, is hydrogen-enriched methane. This is motivated by the po-tential of achieving higher combustion efficiency and lower emission level through hydrogen enrichment. Shock tube experiments on the ignition of two stoichiometric CH4-H2-air mixtures under high pres-sure and moderate temperature conditions have been conducted. It has been observed that the promoting effect of hydrogen is, in general, less than that of ethane and propane; it also decreases with decreas-ing temperature. The difference between pure methane and methane-/hydrogen mixtures is more prominent at 16 atm than that at 40 atm. A low fraction of hydrogen addition shows only weak effects on the ignition delay of methane under the conditions explored. A numerical study of methane/hydrogen/air ignition under the ex-perimental conditions of this work has been conducted using the de-tailed chemical kinetic mechanism developed from the work described in Chapter 3. The mechanism was further optimized based on the lat-est kinetic data from the literature to obtain an improved agreement between the model and experimental results. The effect of hydrogen on methane ignition was mainly related to the generation and con-sumption of H radicals. A t high temperature, the rapid oxidation of hydrogen molecules through R85 {OH + H2 ^ O + H20) and the fast branching reaction R39 {H + 02 ^ O + OH) are mainly responsible for the stronger ignition promoting effect. The rates of both R85 and R39 decrease rapidly with decreasing temperature. A t lower temper-ature, reactions between H2 and CH302 account for a weak effect of hydrogen on methane ignition due to the production of extra H radi-cals. The effect of hydrogen in the reaction system exhibits a negative 217 Chapter 8. Conclusions and Future Work pressure dependence, which implies that the third-order recombina-tion reactions (R33-R38) are rate-limiting steps under the conditions of this work. The new mechanism was completed by adding a NOx submecha-nism selected from the literature. A series of five candidate NOx mech-anisms from the literature was added to the new natural gas reaction mechanism developed in this work; they were used to predict exper-imental NOx data from a stirred reactor [66] and a counterflow diffu-sion flame [36]. The numerical results using the mechanism proposed by Glarborg et al. [84] and GRI-Mech 2.11 show the best agreement with the experimental data. The reburn of NO by methyl radicals is found to be important in reaction systems at relatively low tempera-tures (such as those in the experiment of Steele et al. [66]. A consistent trend has been identified between higher predicted N O concentration and higher rate of a key reaction for prompt N O formation, which sug-gests that the prompt N O mechanism is important under the condi-tions of this study. The manifold methods were used to address the second question at the beginning of this chapter, i.e. to find an effective method to use de-tailed chemistry with the maximum reduction of computational time and reasonable accuracy. While the Intrinsically Low Dimensional M a -nifold method proposed by Maas and Pope [94] is plausible for its rig-orous mathematical derivation, its solution is difficult to achieve due to the highly non-linear implicit construction equation. Furthermore, the fact that I L D M is a non-inertial manifold increases the difficulty of applying the method. The trajectory generated low-dimensional man-ifold method [94], on the other hand, offers an solution to the above two problems of I L D M . However, the method has never been system-218 Chapter 8. Conclusions and Future Work atically validated or put in direct comparison with I L D M to illustrate the difference. The above issues of T G L D M were addressed in this the-sis work. For a one-dimensional manifold, the T G L D M method with con-strained equilibrium states as the onset for the trajectories yields an es-sentially identical manifold as does the I L D M method. This result sup-ports the concept proposed in this thesis work. That is the trajectory method can be used to reduce the problem of solving for an inertial manifold from a global optimization problem to a local optimization problem by optimizing the onset of trajectories. The difference between the general T G L D M and the I L D M gener-ated using the method of Maas and Pope [90] has been investigated. Several merits of the T G L D M method have been identified: the con-struction of the manifold is less mathematically involved than the con-ventional I L D M method; the reaction vector of a T G L D M is always in the tangent plane of the manifold. A n analysis of a hypothetical re-action system showed that if the initial state of the reaction system is known, the T G L D M with optimized initial conditions gives a better approximation to detailed chemistry than the I L D M . The performance of the T G L D M in predicting reactive scalar pro-files for a methane/air mixture has been investigated in modeling lam-inar premixed flames, laminar flamelets, and the transient processes in a perfectly stirred reactor. In general, reasonable agreement between the results using the T G L D M and those from the detailed chemistry has been reached. When the physical time scales approaches the chem-ical time scales in the fast subspace, relatively large errors occur using the T G L D M . Combustion systems modeled using the detailed chemi-cal kinetic mechanism appears to be more resilient to large perturba-219 Chapter 8. Conclusions and Future Work tions than those modeled with the T G L D M method. Finally, to address the last question as one of the objectives of this thesis work, a C S E - T G L D M combustion model has been proposed. The model was used in simulating transient turbulent natural gas com-bustion under engine-relevant conditions to examine the performance of T G L D M derived from the new reaction mechanism in conjunction with a turbulence closure method . The model is able to predict cor-rectly the experimental ignition delay and ignition kernel locations measured in this study as wel l as ignition data reported in the liter-ature. The formation of N O , which is a slow process, was decoupled from the main combustion process. The fact that the concentration of NO approaches equilibrium at a rate significantly slower than most of the species related to NO formation suggests that the formation rate of NO converges significantly faster than its mass fraction on a manif-old. The C S E - T G L D M model is then adjusted to calculate the uncon-ditional rate of N O rather than the unconditional mass fraction as was used for other species. The total NOx calculated from the simulation is less than the experimental measurement. The discrepancy could be explained by the finite cooling time of the rarefaction wave which was not reflected in the numerical model. For the development of the jet flame after ignition, two different modes have been identified based on temperature. A t high tempera-ture, the ignition starts prior to the end of injection. A thin diffusion flame confined to the outer skin of the jet was observed. During the in-jection, the strain close to the exit of the jet is high; the diffusion flame is not able to survive in that area. The end of injection leads to a rapid entrainment of air at the tail of the jet, and the spread of the diffu-220 Chapter 8. Conclusions and Future Work sion flame locally. A t lower temperature, a partially premixed mode combustion due to a significant mixing of fuel and oxidizer prior to ignition was found. The flame is able to propagate rapidly through the core of the jet that leads the jet to expand rapidly. For both combustion modes, the favorable conditions for the propagation or spread of flame are those with the right stoichiometry and low strain-rate. Through the course of this thesis work, the three fundamental ques-tions regarding the natures of natural gas combustion under engine-relevant conditions and the methods of exploring these natures have been addressed consecutively. The study of chemical kinetics for me-thane - additive ignition and NOx formation increased our knowl-edge regarding the governing factors which can be used to control the combustion of natural gas. In particular, a better understanding have been obtained regarding the effects of the environmental factors such as temperature and pressure as well as the effects from the fuel con-stituents on the ignition characteristics of natural gas. The develop-ment of T G L D M method and C S E - T G L D M model provides a powerful tool for applying the knowledge of combustion chemistry in detailed modeling of a turbulent combustion process. This study has shown that the T G L D M method, although being relatively simple in its for-mulation, can provide a reasonable approximation to detailed chem-istry in modeling many reactive flow problems with significantly re-duced computational cost. Such modeling work in conjunction with experimental investigations offers a promising approach for obtaining critical information regarding the combustion process in natural gas engines, and can eventually be used to direct engine designs and im-prove combustion control strategies. 221 Chapter 8. Conclusions and Future Work 8.2 Future Work Further experimental and kinetic studies are necessary to improve our understanding of the ignition characteristics of methane+additive mix-tures under high-pressure and intermediate-temperature conditions. There are still a significant number of elementary reactions in the de-tailed reaction mechanism whose rate constants are less well studied. Some of these rate constants are essential for determining the overall reaction rate of the system. There is a high priority to obtain reliable ex-perimental or theoretical data on the rate constants of these reactions. Given the need to achieve accurate NOx predictions in modeling practical combustion devices, detailed reaction flow and sensitivity analyses are warranted to understand the major factors contributing to the NOx formation under various operating conditions. In particu-lar, there are strong interests in understanding the reason behind the high N 0 2 measurement from a directly-injected natural gas engine un-der high EGR conditions. It is important to further our study on improving the performance of manifold methods in the regions where the physical time scale are comparable with the time scale in the fast subspace. There is also an urgent need to study the generation of manifolds that can be used in the combustion simulation for internal combustion engines. The cur-rent manifolds are generated subject to a constrain of constant pressure and adiabatic boundary conditions, while in an IC engine, the pressure varies constantly and the heat transfer to the wal l of the combustion chamber is significant. A higher dimensional manifold is thus neces-sary to account for these factors. With the increase of dimensions of manifold, the storage and retrieving methods soon becomes a critical 222 Chapter 8. Conclusions and Future Work issue. One possible solution to this problem is to use in situ tabulation of the manifold that eliminates the requirement of storing manifolds or sections of manifold that are never used by a specific problem. While the method in its concept is very attractive, its realization, limitation, efficiency and accuracy needs to be carefully studied before it can be applied to turbulent combustion modeling for IC engines. 223 Bibliography [1] S. Soylu, Examination of combustion characteristics and phasing strategies of a natural gas hcci engine, Energy Conservation and Management 46 (2005) 101-119. [2] A . Ishida, A . Nishimura, M . Uranishi, R. Kihara, A . Nakamura, P. Newman, The development of the ecos-ddf natural gas eng-ine for medium-duty trucks: Exhaust emission reduction against base diesel engine, JSAE Review 22 (2001) 237-243. [3] P. Ouellette, High-pressure direct injection of natural gas in diesel engines, proceeding of ngv 2000: Ngv- transportation for the new century, Seventh International Conference and Exhibi-tion on Natural Gas Vehicles, International Association For Nat-ural Gas Vehicles, Inc., Yokohama, Japan (2000) 235-241. [4] G . McTaggart-Cowan, S. Rogak, P. H i l l , W. Bushe, S. Munshi, Effects of operating condition on particulate matter and nitro-gen oxides emissions from a heavy-duty direct injection natural gas engine using cooled exhaust gas recirculation, International Journal of Engine Research 5 (2004) 499-511. [5] G . McTaggart-Cowan, W. Bushe, P. H i l l , S. Munshi , The applica-tion of exhaust gas recirculation to a heavy duty direct injection of natural gas engine, International Journal of Engine Research 5 (2004) 175-191. [6] G. McTaggart-Cowan, W. Bushe, S. Rogak, P. H i l l , S. Munshi , The effects of varying egr test conditions on a direct injection of nat-ural gas heavy-duty engine with high egr levels, S A E Techni-cal Paper 2004-01-2955. SAE Transactions, Journal of Engines 113 (2004)1500-1509. [7] G . McTaggart-Cowan, W. Bushe, P. H i l l , S. A . Munshi , Super-charged single cylinder heavy duty engine for high pressure di-rect injection of natural gas, International Journal of Engine Re-search 4 (2004)315-329. [8] G . McTaggart-Cowan, W. Bushe, S. Rogak, P. H i l l , S. Munshi , In-jection parameter effects on a direct injected, pilot ignited, heavy duty natural gas engine with egr, S A E Technical Paper 2003-01-3089. S A E Transactions, Journal of Fuels and Lubricants 112 (2003)2103-2109. 224 Bibliography [9] D. J. Seery, C. T. Bowman, A n experimental and analytical study of methane oxidation behind shock wave, Combustion and Flame 14 (1970) 37-47. [10] A . Lifshitz, k. Scheller, A . Burcat, G. B. Skinner, Shock-tube in-vestigation of ignition in methane-oxygen-argon mixtures, Com-bustion and Flame 16 (1971) 331. [11] T. Tsuboi, H . G. Wagner, Homogeneous thermal oxidation of me-thane in reflected shock waves, Proceedings of the Combustion Institute 15 (1974) 883. [12] A . K. Cheng, R. K. and Oppenheim, Autoignition in methane-hydrogen mixtures, Combustion and Flame 58 (1984) 125-139. [13] A . Gril lo, M . W. Slack, Shock tube study of ignition delay times in methane-oxygen-nitrogen-argon mixtures, Combustion and Flame. 27 (1976) 377-381. [14] D . W. Walker, L. H . Diehl, W. A . Strauss, R. Edse, Investigation of ignition properties of flowing combustible gas mixtures, Report No. AFAPL-TR-69-82, U S A F Report. [15] A . Turbiez, A . ElBakali, J. F. Pauwels, A . Rida, P. Meunier, Exper-imental study of a low pressure stoichiometric premixed meth-ane, methane/ethane, methane/ethane/propane and synthetic natural gas flames, Fuel 83 (2004) 933-941. [16] M . Mbarawa, E. W. Mureithi, Modell ing the effects of natural gas composition on dual-fuel combustion under constant-volume conditions, Journal of the Institute of Energy 76 (2003) 2-9. [17] J. Hiltner, R. Agama, F. Mauss, B. Johansson, M . Christensen, Ho-mogeneous charge compression ignition operation with natural gas: Fuel composition implications, Journal of Engineering for Gas Turbine and Power-Transaction of A S M E 125 (2003) 837-844. [18] E. B. Khal i l , G . A . Karim, A kinetic investigation of the role of changes in the composition of natural gas in engine applications, Journal of Engineering for Gas Turbine and Power-Transaction of A S M E 124 (2002) 404-411. [19] A . Burcat, K. Scheller, A . Lifshitz, Shock-tube investigation of comparative ignition delay times for cl-c5 alkanes, Combustion and Flame 16 (1971) 29. 225 Bibliography [20] R. W. Crossley, E. A . Dorko, K . Scheller, A . Burcat, The effect of higher alkanes on the ignition of methane-oxygen-argon mix-tures in shock waves, Combustion and Flame 19 (1972) 373-378. [21] C. K. Westbrook, W. J. Pitz, Effects of propane on ignition of methane-ethane-air mixtures, Combustion Science and Technol-ogy 33 (1983)315-319. [22] M . Frenklach, D. E. Bornside, Shock-ignition in methane-propane mixture, Combustion and Flame 56 (1984) 1-27. [23] L. J. Spadaccini, M . B. Colket III, Ignition delay characteristics of methane fuels, Progress in Energy and Combustion Science 20 (1994) 431^60. [24] A . E l Bakali, P. Dagaut, L. Pillier, P. Desgroux, J. F. Pauwels, A . Rida, P. Meunier, Experimental and modeling study of the oxidation of natural gas in a premixed flame, shock tube, and jet-stirred reactor, Combustion and Flame 137 (2004) 109-128. [25] N . Lamoureux, C. E. Paillard, Natural gas ignition delay times behind reflected shock waves: Application to modeling and safety, Shock Waves 13 (2003) 57-68. [26] Annual energy outlook 2006 with projections to 2030, Energy In-formation Administration Report DOE/EIA-0383(2006), http : //www.eia.doe.gov/oiaf/aeo/. [27] A . G. Gaydon, I. R. Hurle, The shock tube in high temperature chemical physics, Chapman and Ha l l Ltd. , 1963. [28] J. Huang, P. G . H i l l , W. K . Bushe, S. R. Munshi , Shock-tube study of methane ignition under engine-relevant conditions: experi-ments and modeling, Combustion and Flame 136 (2004) 25-42. [29] E. L. Petersen, D. F. Davidson,\" R. K . Hanson, Kinetics model-ing of shock-induced ignition in low-dilution ch4/ o2 mixtures at high pressures and intermediate temperatures, Combustion and Flame 117 (1999) 272-290. [30] E. Petersen, D. Davidson, R. Rohrig, M . and Hanson, High-pressure shock-tube measurements of ignition times in stoichio-metric h2/o2/ar mixtures, In 20th Symp. Int. Shock Waves (1996) 941-946. [31] D. C. Bull , J. E. Elsworth, G . Hooper, Susceptibility of methane-ethane mixtures to gaseous detonation in air, Combustion and Flame. 34 (1979) 327. 226 Bibliography [32] C. S. Eubank, M . J. Rabinowitz, W. C. J. Gardiner, R. E. Zellner, Shock initiated ignition of natural gas-air mixtures, Eighteenth Symposium (International) on Combustion. The Combustion In-stitute, Pittsburgh (1981) 1767. [33] R. M . R. Higgin, A . Williams, A shock-tube investigation of the ignition of lean methane and n-butane mixtures with oxygen, Sixteenth Symposium (International) on Combustion. The Com-bustion Institute, Pittsburgh (1968) 579. [34] R. Zellner, K . J. Niemitz, J. Warnatz, W. C. J. Gardiner, C. S. Eu-bank, J. M . Simmif, Prog, in Astro and Aero 88 (1981) 252. [35] J. R. Griffiths, D. Coopersthwaite, C. H . Philips, C. K . West-brook, Auto-ignition temperature of binary mixtures of alkane in a closed vessel: Comparisons between experimental mea-surements and numerical predictions, 23rd Symposium (inter-national) on Combustion (1990) 1745. [36] S. C. L i , F. A . Williams, Nox formation in two-stage methane-air flames, Combustion and Flame 118 (1999) 399^14. [37] K . Hughes, T. Turnyi, A . Clague, P. M . J., Development and test-ing of a comprehensive chemical mechanism for the oxidation of methane, International Journal of Chemical Kinetics 33 (2001) 513-538. [38] T. B. Hunter, H . Wang, L. T. A . , F. M . , The oxidation of methane at elevated pressure: experiments and modeling, Combustion and Flame 87 (1991) 365-370. [39] M . Frenklach, H . Wang, M . Goldenberg, G . P. Smith, D. M . Golden, B. C. T , R. K . Hanson, W. C. Gardiner, V. Lissianski, Gri-mech-an optimized detailed chemical reaction mechanism for methane combustion, Report No. GRI-95/0058. [40] C. T. Bowman, R. K . Hanson, D. F. Davidson, W. C. J. Gardiner, V. Lissianski, G . P. Smith, D. M . Golden, F. M . , M . Goldenberg, Gri-mech 2.11, http : //www.me.berkeley.edu/grimech/. [41] G . P. Smith, D. M . Golden, M . Frenklach, N . W. Moriarty, B. Eite-neer, M . Goldenberg, C. T. Bowman, R. K. Hanson, S. S. Song, W. C. Gardiner, V. V. Lissianski, Z . Qin, Gri-mech3.0, http : j/www. me.berkeley.edu/'grimech/'.. [42] E. Ranzi, A . Sogaro, P. Gaffuri, C. Pennati, C. K . Westbrook, W J. Pitz, A new comprehensive reaction mechanism for combustion of hydrocarbon fuels, Combustion and Flame 99 (1994) 201-211. 227 Bibliography [43] C. K . Westbrook, A n analytical study of the shock tube ignition of methane and ethane, Combustion Science and Technology 20 (1979) 5. [44] T. B. Hunter, T. A . Litzinger, H . Wang, M . Frenklach, Ethane ox-idation at elevated pressures in the intermediate temperature regime: Experiments and modeling, Combustion and Flame 104 (1996)505-523. [45] I. R. Roger, Natural gas/hydrogen mixtures for low noxious emissions, Journal of Scientific and Industrial Research 62 (2003) 64-70. [46] A . Raju, B. Ramesh, A . and Nagalingam, Effect of hydrogen in-duction on the performance of a natural-gas-fueled lean-burn si engine, Journal of the Institute of Energy 73 (2000) 143-148. [47] C. G . a. F. T. W. Bauer, Effect of hydrogen addition on the perfor-mance of methane-fueled vehicles, part i , effect on si engine per-formance, International Journal of Hydrogen Energy 26 (2001) 55-70. [48] R. Sierens, E. Rosseel, Variable composition hydrogen/natural gas mixtures for increased engine efficiency and decreased emis-sions, Journal of Engineering for Gas Turbines and Power-Transactions of the A S M E 122 (2000) 135-140. [49] S. Allenby, A . Chang W. C. and Megaritis, M . L. Wyszynski, H y -drogen enrichment: a way to maintain combustion stability in a natural gas fueled engine with exhaust gas recirculation, the potential of fuel reforming, proceedings of the institution of me-chanical engineers, part d, Journal of Automobile Engineering 215 (2001) 405^18. [50] J. L. Gauducheau, B. Denet, G . Searby, A numerical study of lean ch4/h2/air premixed flames at high pressure, Combustion Sci-ence and Technology. 137 (1998) 81-99. [51] Y. Ju, T. Ni ioka , Ignition simulation of methane/hydrogen mix-tures in a supersonic mixing layer, Combustion and Flame 102 (1995) 462-470. [52] J. Y. Ren, T. T. Egolfopoulos, F. N . and Tsotsis, Nox emission control of lean methane-air combustion with addition of met-hane reforming products, Combustion Science and Technology 174 (2002)181-205. 228 Bibliography [53] C. G . Fotache, G . Dreutz, T, C. K . Law, Ignition of hydrogen en-riched methane by heated air, Combustion and Flame 110 (1997) 429^40. [54] D. Davidson, M . DiRosa, A . Chang, R. Hanson, C. Bowman, A shock tube study of methane decomposition using laser absorp-tion by ch3, Twenty-fourth Symposium (International) on Com-bustion, The Combustion Institute (1992) 589-596. [55] C. L. Yu, C. Wang, M . Frenklach, Chemical kinetics of methyl oxidation by molecular oxygen, Journal of Physical Chemistry 99 (1995) 14377 - 14387. [56] R. A . Marcus, Unimolecular dissociation and free radical recom-bination reactions, Journal of Chemical Physics 20 (1952) 359-364. [57] K . K . Kuo, Principles of combustion, John Wiley and Son, New York, 2005. [58] E. F.N. , P. Cho, C. Law, Laminar flame speeds of methane-air mixtures under reduced and elevated pressures, Combustion and Flame 76 (1989) 375. [59] A . Garforth, C. Rallis, Laminar burning velocity of stoichiometric methane-air: pressure and temperature dependence, Combus-tion and Flame 31 (1978) 53. [60] G . Andrews, D. Bradley, Burning velocity of methane/air mix-tures, Combustion and Flame 19 (1972) 275. [61] T. Iijima, T. Takeno, Effects of temperature and pressure on burn-ing velocity, Combustion and Flame 65 (1986) 35. [62] J. Warnatz, U . Maas, W. Dibble, R, Combustion, 2nd Edition, Springer-Verlag, New York, 1999. [63] Y. B. Zeldovich, The oxidation of nitrogen in combustion and ex-plosives, Acta Physicochim USSR 21 (1946) 577. [64] J. S. Bernstein, A . Fein, J. B. Choi , C. Terrill, A . Sausa, S. L. Howard, R. J. Locke, A . W. Miziolek, Laser-based flame species profile measurements: A comparison with flame model predic-tions, Combustion and Flame 92 (1993) 85-105. [65] D. Heard, J. Jeffries, G . Smith, D. Crosley, Lif measurements in methane/air flames of radicals important in prompt-no forma-tion, Combustion and Flame 88 (1992) 137. 229 Bibliography [66] R. Steele, R Malte, D. Nichol , J. Kramlich, Nox and n2o in lean-premixed jet-stirred flames, Combustion and Flame 100 (1995) 440. [67] W. Bartok, V. Engleman, R. Goldstein, E. delValle, Basic kinetic studies and modeling of nitrogen oxide formation in combustion processes, AIChE Symposium Series 68 (1972) 30-38. [68] B. Williams, J. Fleming, Comparative species concentrations in ch4/ o2/ar flames doped with n2o, no and no2, Combustion and Flame 98 (1994) 93. [69] T. Etzkorn, S. Muris , J. Wolfrum, C. Dembny, H . Bockhorn, P. Ne l -son, A . Attia-Shahin, J. Warnatz, Destruction and formation of no in low pressure stoichiometric ch4/o2 flames, Twenty-fourth Symposium (International) on Combustion, The Combustion In-stitute (1992) 925-932. [70] D. L. Baulch, C. J. Cobos, R. A . Cox, P. Frank, G . Hayman, T. Just, J. A . Kerr, T. Murrells, M . J. Pil l ing, J. Troe, J. Walker, R. W. and Warnatz, Evaluated kinetic data for combustion mod-eling, supplement i , ref. data, Journal of Physical Chemistry 23 (1994)847-1033. [71] J. B. Heywood, Internal Combustion Engine Fundmentals, McGraw-Hi l l , 1988. [72] C. P. Fenimore, Studies of fuel-nitrogen species in rich flame gases, 17th Symposium (intl) on combustion, the combustion in-stitute, Pittsburgh (1979) 661-670. [73] C. P. Fenimore, Formation of nitro oxide in premixed hydrocar-bon flames, 13th Symposium (intl) on combustion, the combus-tion institute, Pittsburgh (1970) 373-380. [74] P. Glarborg, J. A . Miller, R. J. Kee, Kinetic modeling and sensitiv-ity analysis of nitrogen oxide formation in well-stirred reactors, Combustion and Flame 65 (1986) 177-202. [75] Y. Matsui, A . Yuuki , Radical concentrations and prompt no for-mation in hydrocarbon-air premixed flames, Jpn. J. A p p l . Phys. 24. [76] D. Lindackers, M . Burmeister, P. Roth, Perturbation studies of high temperature c and ch reactions with n2 and no, Symp. Int. Combust. Proc. 23 (1991) 251 - 257. 230 Bibliography [77] J. Miller, S. Walch, Prompt no: theoretical prediction of the high-temperature rate coefficient for ch + n2 —>• hen + n, International Journal of Chemical Kinetics 29 (1997) 253 - 259. [78] J. A . Miller, C. T. Bowman, Mechanism and modeling of nitrogen chemistry in combustion, Progress in Energy and Combustion Science 15 (1989) 287-338. [79] A . H . Lefebvre, The role of fuel preparation in low-emission combustion, Journal of Engineering for Gas Turbine and Power-Transaction of A S M E 117 (1995) 617-654. [80] S. C. L i , I. N , F. A . Williams, Reduction of nox formation by wa-ter sprays in strained two-stage flames, Journal of Engineering for Gas Turbines and Power-Transactions of the A S M E 119 (1997) 836-843. [81] M . U . Alzueta, P. Glarborg, K. Dam-Johansen, Low temperature interactions between hydrocarbons and nitric oxide: A n experi-mental study, Combustion and Flame 109 (1997) 25-36. [82] P. G . Kristensen, P. Glarborg, K . Dam-Johansen, Nitrogen chemi-stry during burnout in fuel-staged combustion, Combustion and Flame 107 (1996) 211-222. [83] P. Glarborg, P. G . Kristensen, S. H . Jensen, A flow reactor study of hnco oxidation chemistry, Combustion and Flame 98 (1994) 241-258. [84] P. Glarborg, M . U . Alzueta, K . Dam-Johansen, A flow reactor study of hnco oxidation chemistry, Combustion and Flame 115 (1998) 1-27. [85] M . Frenklach, Complex chemical reaction systems, mathematical modeling and simulation, Springer series in chemical physics 47 (1987) 2-16. [86] H . Wang, M . Frenklach, Detailed reduction of reaction mechan-ism for flame modeling, Combustion and Flame 87 (1991) 365-370. [87] S. Lam, D. A . Goussis, Research on computational singular per-turbation, Progress Report #1, Report #1749(a)-MAE, Depart-ment of Mechanical and Aerospace Engineering, Princeton U n i -versity, 1987. 231 Bibliography [88] S. Lam, D. A . Goussis, The csp method for simplifying kinetics, International Journal of Chemical Kinetics 26 (1994) 461^486. [89] T. L u , Y. Ju, C. K. Law, Complex csp for chemistry reduction and analysis, Combustion and Flame 126 (2001) 1445-1455. [90] U . Maas, S. B. Pope, Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space> Combustion and Flame 88 (1992) 239-264. [91] W. Press, B. Flannery, S. Teukolsky, W. Vetterling, Numerical recipes: The art of scientific computing, Cambridge University Press, New York, NY, 1986. [92] J. Nafe, U . Maas, Modeling of no formation based on i ldm re-duced chemistry, Proceeding of the Combustion Institute 29 (2002) 1379-1385. [93] S. J. Fraser, The steady and equilibrium approximation: a geo-metrical picture, Journal of Chemical Physics 88 (1988) 4732. [94] S. B. Pope, U . Maas, Simplifying chemical kinetics: trajectory-generated low-dimensional manifolds, F D A 93-11. [95] J. P. Ignizio, T. M . Cavalier, Linear programming, Prentice Ha l l International Series in Industrial and Systems Engineering, Pren-tice Ha l l , Englewood Cliffs, NJ , 1994. [96] P. N . , Laminar diffusion flamelet models in non-premixed turbu-lent combustion, Progress in Energy and Combustion Science 10 (1982)369-399. [97] A . Y. Klimenko, R. W. Bilger, Conditional moment closure for tur-bulent combustion, Progress in energy and combustion science 25 (1999)595-687. [98] S. B. Pope, Pdf methods for turbulent reactive flow, Progress in Energy and Combustion Science 11 (1985) 119-192. [99] A . R. Masri , S. B. Pope, Pdf calculations of piloted turbulent non-premixed flames of methane, Combustion and Flame 81 (1990) 13-29. [100] A . T. Norris, S. B. Pope, Modeling of extinction in turbulent-diffusion flames by the velocity-dissipation-composition pdf method, Combustion and Flame 100 (1996) 211-220. 232 Bibliography [101] V. Saxena, S. B. Pope, Pdf calculation of major and minor species in a turbulent piloted jet flame, Twenty-Seventh Symposium (International) on Combustion, the Combustion Institute (1998) 1081-1086. [102] C. K. Law, Steady-state diffusion flame structure with lewis number variation, Combustion Science and Technology 29 (1982) 129-145. [103] W. P. Jones, W. J. H , Calculation methods for reacting turbulent flows: A review, Combustion and Flame 48 (1982) 1-26. [104] T. Saitoh, O. Y., Unsteady behavior of diffusion flames and pre-mixed flames for counter flow geometry, Combustion Science and Technology 12 (1976) 135. [105] F. MauB, D. Keller, N . Peters, A lagrangian simulation of flamelet extinction and re-ignition in turbulent jet diffusion flames, Twenty-third symposium (Int.) on combustion, the combustion institute (1990) 693-698. [106] N . Peters, Turbulent combustion, Cambridge University Press, Cambridge, 2000. [107] P. J. Coelho, N . Peters, Unsteady modeling of a piloted meth-ane/air jet flame based on the eulerian particle flamelet model, Combustion and Flame 124 (2001) AM-A65. [108] H . Barths, C. Antoni , N . Peters, Three-dimensional simulation of pollutant formation in a di diesel engine using multiple interac-tive flamelets, S A E technical paper, No. 982459. [109] C. A . Hergart, H . Barths, N . Peters, Modeling the combustion in a small-bore diesel engine using a method based on representative interactive flamelets, S A E technical paper, No. 1999-01-3550. [110] C. Hasse, H . Barths, P. N . , Modeling the effect of split injections in diesel engines using representative interactive flamelets, S A E technical paper, No. 1999-01-3547. [ I l l ] H . Barths, C. Hasse, P. N . , Computational fluid dynamics mod-elling of non-premixed combustion in direct injection diesel en-gines, International Journal of Engine Research 1 (2000) 249-267. [112] S. Rao, C. J. Rutland, A flamelet time scale model for non-premixed combustion including chemical kinetic effects, Com-bustion and Flame 133 (2003) 189-191. 233 Bibliography [113] N . Swaminathan, Flamelet regime in non-premixed combustion, Combustion and Flame 129 (2002) 217-219. [114] A . Y. Klimenko, Multicomponent diffusion of various admix-tures in turbulent flow, Fluid Dynamics 25 (1990) 327-334. [115] R. W. Bilger, Conditional moment methods for turbulent reacting flow using crocco variable conditions, Report T N F-99, Depart-ment of Mechanical Engineering, University of Sydney. [116] R. Lee, J. H . Whitelaw, T. S. Wung (Eds.), Aerothermodynamics in combustion, Springer, Berlin, 1991. [117] N . Smith, R. Bilger, C. Carter, A comparison of cmc and pdf mod-eling predictions with experimental nitric oxide lif/raman mea-surements in a turbulent h-2 jet flame, Combustion Science and Technology 105 (1995) 357-375. [118] M . Roomina, R. Bilger, Conditional moment closure modeling of turbulent methanol jet flames, Combustion Theory and Model-ing 3 (1999) 689-708. [119] M . Roomina, R. W. Bilger, Conditional moment closure (cmc) prediction of a turbulent methane-air jet flames, Combustion and Flame 125 (2001) 1176-1195. [120] M . Fairweather, R. M . Woolley, First-order conditional moment closure modeling of turbulent, non-premixed methane flames, Combustion and Flame 138 (2004) 3-19. [121] C. Devaud, K . N . C. Bray, Assessment of the applicability of con-ditional moment closure to a lifted turbulent flame: first order model, Combustion and Flame 132 (2003) 102-114. [122] S. K i m , K . Y. Huh , L. Tao, Application of the elliptic conditional moment closure model to a two-dimensional nonpremixed methanol bluff-body flame, Combustion and Flame 120 (2000) 75-90. [123] S. K i m , K. Y. Huh , Use of the conditional moment closure model to predict no formation in a turbulent ch4/h-2 flame over a bluff-body, Combustion and Flame 130 (2002) 94-111. [124] A . Kronenburg, R. Bilger, J. Kent, Modeling soot formation in tur-bulent methane-air jet diffusion flames, Combustion and Flame 121(2000)24-40. 234 Bibliography [125] D. Bradley, D. Emerson, P. Gaskell, Mathematical modeling of turbulent non-premixed piloted-jet flames with local extinction, Proceedings of the combustion institute 29 (2002) 2155-2162. [126] E. Mastorakos, R. Bilger, Second-order conditional moment clo-sure for the autoignition of turbulent flows, Physics of Fluids 10 (1998)1246-1248. [127] W. Bushe, H . Steiner, Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows, Physics of Fluids 11 (1999) 1896-1906. [128] H . Steiner, W. K. Bushe, Large eddy simulation of a turbulent reacting jet with conditional source-term estimation, Physics of Fluids 13 (2001) 754-769. [129] W. Bushe, H . Steiner, Laminar flamelet decomposition for condi-tional source-term estimation, Physics of Fluids 15 (2003) 1564-1575. [130] C. Blair, Implementation of conditional source term estimation for prediction of methane ignition, Master Thesis, Department of Mechanical Engineering, University of British Columbia, 2003. [131] R. Grout, Combustion modeling with conditional source term estimation and laminar flamelet decomposition, Master Thesis, Department of Mechanical Engineering, University of British Columbia, 2004. [132] D. Flowers, S. Aceves, C. K . Westbrook, J. R. Smith, D. R., Detailed chemical kinetic simulation of natural gas hcci com-bustion: Gas composition effects and investigation of control strategies, Journal of Engineering for Gas Turbine and Power-Transaction of A S M E 123 (2001) 433-439. [133] G. Ben-Dor, O. Igra, T. Elperin, In Handbook of Shock Waves, Vol. 1, Academic Press, New York, 2001. [134] J. E. G . Lapworth, K. C. and Townsend, Temperature and pres-sure studies in the reservoir of a reflected-shock hypersonic tun-nel, A . R . C 26 (1964) 100. [135] H . Mark, The interaction of a reflected shock wave with the boundary layer in a shock tube, N S C A T M (1958) 1418. [136] L . Davies, The interaction of the reflected shock with the bound-ary layer in a shock tube and its influence on the duration of hot flow in the reflected-shock tunnel, A . R . C (1967) 880-881. 235 Bibliography [137] G. H . Markstein, Flow disturbances induced near a slightly wavy contact surface or flame front traversed by a shock wave, J. Aero-nautic Sci. 24 (1957) 238-239. [138] G. I. Taylor, The instability of l iquid surface when accelerated in a direction perpendicular to their planes 1, Proc. Roy. Soc. (1950) A 201. [139] J. P. Holman, Experimental methods for engineers, seventh Edi -tion, McGraw-Hi l l , New York, 2001. [140] G. B. Skinner, Limitation of the reflected shock technique for studying fast chemical reactions, Journal of Chemical Physics 31 (1959)268. [141] W. Tsang, R. Hampson, Chemical kinetic data base for combus-tion chemistry, part i . methane and related compounds, Journal of Chemical Physics: Reference Data 15 (1986) 1087. [142] H . J. Curran, P. Gaffuri, W. J. Pitz, C. K. Westbrook, A comprehen-sive modeling study of iso-octane oxidation, Combustion and Flame 129 (2002) 253-280. [143] D. L. Baulch, C. J. Cobos, R. A . Cox, C. Esser, P. Frank, T. Just, J. A . Kerr, M . J. Pil l ing, J. Troe, R. W. Walker, J. Warnatz, Evalu-ated kinetic data for combustion modeling, Journal of Physical Chemistry 21 (1992) 411-429. [144] N . Cohen, Are reaction rate coefficients additive? revised transi-tion state theory calculations for oh + alkane reactions, Journal of Chemical Kinetics 23 (1991) 239^17. [145] A . Burcat, B. McBride, 1994 ideal gas thermodynamic data for combustion and air-pollution use, Technion Aerospace Engi-neering Report No. T A E 697. [146] A . Burcat, Third millennium ideal gas and condensed phase ther-mochemical database for combustion, technion aerospace engi-neering report no. tae 867. [147] R. J. Kee, F. M . Rupley, J. A . Miller, Chemkin-II: A fortran chemi-cal kinetics package for the analysis of gas phase chemical kinet-ics, SAND89-8009, UC401. [148] P. N . Brown, G. D. Byrne, H . A . C , Vode: a variable coefficient ode solver, S I A M J. Sci. Stat. Comput. 10 (1988) 1038-1051. 236 Bibliography [149] W. C. Gardiner, V. Lissianski, V, M . Zamanski, V, Reduced chem-ical reaction mechanism of shock-initiated ignition of methane and ethane mixtures with oxygen, 19th International Sympo-sium on Shock Waves, Marseille (1993) 154-159. [150] C. K . Westbrook, W. J. Pitz, H . C. Curran, J. Boercker, Chemical kinetic modelling study of shock tube ignition of heptane iso-mers, International Journal of Chemical Kinetics 33 (2001) 868-877. [151] V. M . Zamansky, A . A . Borisov, Promotion of high-temperature self-ignition, Energy Combust. Sci 18 (1992) 297-235. [152] S. C. L i , F. A . Williams, Reaction mechanisms for methane igni-tion, Journal of Engineering for Gas Turbine and Power - Trans-action of A S M E 124 (2002) 471-480. [153] D. F. Cooke, A . Wiliams, Shock-tube study of the ignition delay and combustion of ethane and slightly rich methane mixtures wi th oxygen, 13th Symposium (Int) on Combustion, The Com-bustion Institute (1970) 757-765. [154] R. Akbar, M . Kaneshige, E. Schultz, J. E. Shepherd, Detonations in h2—n2 o—chA—nh3—o2—n2 mixtures, Technical Report FM97-3, GALCIT . [155] J. W. Sutherland, M . C. Su, J. V. Michael, Rate constants for h + ch4,ch3 + h2, and c/ i 4 dissociation at high temperature, In-ternational Journal of Chemical Kinetics 33 (2001) 669-684. [156] A . D., Isaacson, Harmonic and anharmonic rate constants and transmission coefficients obtained from ab initio data, Journal of Chemical Physics 107 (1997) 3832-3839. [157] S. I. W. M . and C r i m F. E , The chemical kinetics and dynamics of the prototypical reaction: oh + h2 —> h2o + h, Journal of Physical Chemistry Chemical Physics 4 (2002) 3543-3551. [158] H . J. Curran, W J. Pitz, N . M . Marinov, C. K . Westbrook, P. Da-gaut, M . Boettner, J-C. and Cathonnet, A wide range modeling study of dimethyl ether oxidation, Journal of Chemical Kinetics 30(1998)229-241. [159] R. Craig, A shock tube study of the ignition delay of hydrogen-air mixtures near the second, explosion limit, Technical Re-port, AFAPL-TR-66-74, A i r Force Aero-Propulsion Lab, Wright-Patterson. 237 Bibliography [160] G. Skinner, G. and Ringrose, Ignition delays of a hydrogen-oxygen-argon mixture at relatively low temperature, Journal of Chemical Physics 42 (1966) 2190-2192. [161] R. Blumenthal, K . Fieweger, K. Komp, G. Adomeit, B. Gelfand, Self-ignition of h2-air mixtures at high pressure and low temper-ature, In 20th Symp. Int. Shock Waves (1996) 935-940. [162] A . B. Bendtsen, P. Glarborg, K . Dam-Johansen, Chemometric analysis of a detailed chemical reaction mechanism for meth-ane oxidation, Chemometric and intelligent laboratory systems 44 (1998)353-361. [163] C. Mansour, A . Bounif, A . Aris , F. Gaillard, Gas-diesel (dual-fuel) modeling in diesel engine environment, International Journal of Thermal Sciences 40 (2001) 409-424. [164] Q. P. Zheng, H . M . Zhang, D. F. Zhang, A computational study of combustion in compression ignition natural gas engine with separated chamber, Fuel 84 (2005) 12-13. [165] J. Huang, W. K . Bushe, Experimental and kinetic study of au-toignition in methane/ethane/air and methane/propane/air mixtures under engine-relevant conditions, Combustion and Flame. [166] H . Pitsch, Flamemaster v3.1: a C++ computer program for. 0d combustion and Id laminar flame calculations, available from http://www.stanford.edu/ hpitsch. [167] O. Gicquel, D. Thevenin, M . Hi lka , Direct numerical simulation of turbulent premixed flames using intrinsic low-dimensional manifolds, Combustion Theory and Modell ing 3 (1999) 479-502. [168] D. Thevenin, Three-dimensional direct simulations and struc-ture of expanding turbulent methane flames, Proceedings of the Combustion Institute 30 (2005) 629-637. [169] P. Nooren, H . Wouters, T. W. J. Peeters, Monte carlo pdf modeling of a turbulent natural-gas diffusion flame, Combustion Theory and Modeling 1 (1997) 79-96. [170] H . Bongers, J. A . Van Oijen, L. De Goey, Intrinsic low-dimensional manifold method extended with diffusion, Pro-ceedings of the Combustion Institute 29 (2002) 1371-1378. 238 Bibliography [171] J. A. VanOijen, L. De Goey, Modelling of premixed laminar flames using flamelet-generated manifolds, Combustion Science and Technology 161 (2000) 113-137. [172] B. Fiorina, O. Gicquel, L. Vervisch, Premixed turbulent combus-tion modeling using tabulated detailed chemistry and pdf, Pro-ceeding of the Combustion Institute 30 (2005) 867-874. [173] S. Singh, J. Powers, S. Paolucci, On slow manifolds of chemically reactive systems, Journal of Chemical Physics 117 (2002) 1482-1496. [174] M. Davis, R. Skodje, Geometric investigation of low-dimensional manifolds in systems approaching equilibrium, Journal of Chemical Physics 111 (1999) 859-874. [175] S. Gordon, B. J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311. [176] W. C. Reynolds, The element potential method for chemical equilibrium analysis: implementation in the interactive program stanjan, Stanford University Report ME 270 HO No.7. [177] S. B. Pope, Gibbs function continuation for the stable computa-tion of chemical equilibrium, Combustion and Flame 139 (2004) 222-226. [178] R. Renka, Algorithm 751: Tripack: A constrained two-dimensional delaunay triangulation package, A C M Transactions on Mathematical Software 22 (1996) 1-8. [179] G. Strang, On construction and comparison of differential schemes, SIAM Journal of Numerical Analysis 5 (1968) 506. [180] S. S. H., R. W. Bilger, K. M. Lyons, J. H. Frank, M. B. Long, Con-served scalar measurements in turbulent-diffusion flames by a raman and rayleigh ribbon imaging method, Combustion and Flame 99 (1994) 347-354. [181] Y. Chen, M. S. Mansour, Measurements of scalar dissipation in turbulent hydrogen diffusion flames and some implications on combustion modeling, Combustion Science and Technology 126 (1997)291-313. [182] G. D. Sullivan, J. Huang, T. X. Wang, W. K. Bushe, S. N. Rogak, Emissions variability in gas fuel direct injection compression ig-nition combustion, SAE Technical Paper No. 2005-01-0917. 239 Bibliography [183] B. E. Launder, D. B. Spalding, Mathematical model of turbulence, Academic Press, 1972. [184] B. J. P., D. L. Book, Flux-corrected transport i. shasta, a fluid transport algorithm that works, Journal of Computational Physics 11 (1973) 38-69. [185] P. G. Hill, P, Ouellette, Transient turbulent gaseous fuel jets for diesel engine, Journal of Fluid Engineering 121 (1999) 93-101. [186] J. S. Turner, The 'starting plum' in neutral surrounding, Journal of Fluid Mechanics 13 (1962) 356-368. [187] F. P. Ricou, D. B. Spalding, Measurement of entrainment by asymmetrical turbulent jets, Journal of Fluid Mechanics 11 (1961) 21-32. [188] J. D. Naber, D. L. Siebers, J. A. Caton, C. K. Westbrook, S. S. Di Julio, Natural gas autoignition under diesel conditions: exper-iments and chemical kinetic modeling, SAE Paper No. 942034. [189] J. Iaconis, An investigation of methane autoignition behav-ior under diesel engine-relevant conditions, Master Thesis, Department of Mechanical Engineering, University of British Columbia. [190] D. C. Lay, Linear algebra and its applications, Addison-Wesley, 2003. 240 Appendix A Derivations of Key Equations in CSP For a reaction system involving N species and K elementary reactions, the rates in which the mass fractions of species change are given by Gi = ^ = S<,;F;> (A- 1 ) where y is the array of mass fraction of N species; £ = l...JV;SisaNxK matrix of stoichiometric coefficient; j = 1...K, and F is the rate array of K reactions. In the Computational Singular Perturbation (CSP) method [87,88], the composition space, $tN, is split into a fast subspace, G™ast, and a slow subspaces, G™low < G%st U Gnslow >= Span(G7ast U G ^ J = MN (A.2) < G / L t n Gliow >=o. The superscripts m and n denote the dimensions of the fast subspace and slow subspace respectively. Let be a set of basis vectors span-ning GJast; let bl be the dual vectors of a ; that satisfy the orthonormal relation b^- = 5), (A.3) in which 8 is the Kronecker delta. The basis-vector matrix, A, whose columns are the basis vectors of Gfast, defines a linear transformation 241 Appendix A. Deriva tions of Key Equa tions in CSP (projection) onto the fast subspace. G%st = Af, (A.4) where f is a vector whose elements represent the magnitude of fast reactions. Mul t ip ly both sides of Eq. A.4 by a dual-vector matrix, B , whose rows are bj. One get VGJast = 6i^ (A.5) B A f = f. By definition, the nul l space (or kernel) of transformation B is the slow subspace B G l r a = 0. (A.6) The CSP method reduces the computational time for solving stiff ODEs by examining the contribution of G™ a s t, which is represented by the magnitude of f. When the contribution of Gfast becomes trivial com-pared with that of G™low, the processes in the fast subspace is in quasi-steady state. Large time steps is then used to speed up the numerical integration of Eq. A . l because the system is now governed only by the slow processes. The equation governing the change of f can be ob-tained by differentiating Eq. A.5 df _ dGfast dB - d t - B - d T + G ^ H ( A - 7 ) _udGdy dB ~alfydt +(*fast~dt _ -p dG ( r i m r „ , _,m dB — a~f^\\K:'fast ^slow) ^fast ^ = ( B J A + f A ) f , in which J = dG/dy is the Jacobian matrix. Note that Eq. A.6 was used to derive the second step in the above equation. 242 A p p e n d i x B Eigen Decomposition and Schur Decomposition For a n x n square matrix 2 2 A, if A has non-degenerative eigenvalues, there exists a decomposition of A so that A = P D P 1 , (B.l) where P is a matrix of right eigenvectors, D is a diagonal matrix whose diagonal entries are the eigenvalues of A. The matrix decomposition of Eq. B.l is called eigen decomposition. By definition, a general eigenvalue problem regarding a linear sys-tem A can be written as Ax = Ax, (B.2) where x is the right eigenvector corresponding to eigenvalue A. Equa-tion B.2 can be rearranged into (A - l)x = 0, (B.3) were I is an n x n identity matrix. Analytically, the eigenvalues of A can be computed from det(A - I) = 0, (B.4) 2 2 A matrix whose vertical and horizontal dimensions are the same 243 Appendix B. Eigen Decomposition and Schur Decomposition which requires solving a polynomial of nth order. Numerically, to solve Eq. B.4 is highly inefficient especially for high-dimensional problems. In practice, A is usually transformed into a Hessenberg matrix through a similarity transformation [190] (e.g. householder transformation) S A S ^ S z = XSx (B.5) By = Ay. Matrix B is a Hessenberg matrix, which is subsequently decomposed to B = V A V T , (B.6) where V is a unitary matrix,A is an upper triangular matrix whose di-agonal entries are the eigenvalues of B (which shares the same eigen-values as A ) . The decomposition given in Eq. B.6 is called Schur de-composition and the resulting matrix V contains a set of orthonor-mal basis vectors of A , which are called Schur vectors [190]. Given the eigenvalues, the eigenvectors of A can be computed using the iterative QR factorization method. In many occasions, the Jacobian matrix of a reaction system does not have a full rank as a result of nearly liner dependence among the rows in the Jacobian matrix. Consequently, the eigenvalue decomposi-tion may not be computed, but a Schur decomposition always exists. For this reason, Maas and Pope [90] proposed using Schur vectors in-stead of eigenvectors to form the basis of the Jacobian matrix. 244 Appendix C Solution of Ricatti's Differential Equation The Ricatti's equation with constant coefficients is given by du = m + ny + py . (C.l) dt The initial condition for the current problem is yt=o = 0. Namely, no methyl radical appears in the initial mixture before reaction. To solve Eq. C . l , we rewrite the equation as d(y + s) dt = a(y + s) + b(y + sY, (C.2) where m = as + bs2 (C3) n = a + 2bs p — b. Solving the equation system C.3, we get n±^/n2 - 4pm S = 2p ( C 4 ) b = p. 245 Appendix C. Solution ofRicatti's Differential Equation Let y + s = u, Eq. C.2 can be rewritten as du , =au + bu2. (C.5) at Rearranging Eq. C.5, we get d 1 1 -£ = -a--b. (C.6) Let z = we rewrite Eq C.6 as — + az + b = 0. (C.7) This is a first-order homogeneous equation, which can be readily solved to get 2 = - - + c e - a * . (C.8) a where c is an integral constant, which need to be determined from the initial condition of the problem. Given yt=0 = 0, we get ^ = 0 = i = - - + C ) (C.9) s a or c = - + - . (CIO) s a The final form of z can be obtained by substituting Eq. C.10 into Eq C.8 to get z = - b - + ( \\ + b - ) e - < * . ( C H ) a s a Thus the solution for y is given by V=-l + p + l y * - ' - a v s a> The current analytical model determines the ignition delay time based on extrapolation of the slope at the inflection point on the y — t curve. 246 Appendix C. Solution ofRicatti's Differential Equation The slope is determined by the first derivative of y with respect to t, given by y> = ^ (C.13) [ - ! + [ ( i + ! ) e - ] 2 The second derivative of y with respect to t is The inflection point can be calculated by setting y\" = 0, which yields = a V s (C.15) a The equilibrium value of y can be obtained at t —> oo from Eq. C.12 = - S- (C.16) The time from the inflection point to the intersection of extrapolated slope and equilibrium value of y can be calculated from tex = V e q u i ~ V t = t i . (C.17) Vt=ti Substituting Eqs. C15 ,C13 and C.12 into Eq. C.17 and solve for tex, we get tex = -. (C.18) a Finally, the ignition delay is calculated by combining tex with ^ T = ill < (C.19) 247 Appendix D Flux Corrected Transport Flux-Corrected Transport (FCT) [184] is a conservative, monotone tech-nique for integrating generalized continuity and hydromagnetic equa-tions. It is especially useful for solving compressible-flow problems, particularly those involving shock and rarefaction waves and contact discontinuities. FCT accomplishes this objective by combining integra-tion schemes with low and high orders of spatial accuracy. The low-order scheme provides a monotone solution, usually by the introduc-tion of diffusive numerical fluxes, while the high-order scheme pro-vides high accuracy in regions of smooth flow and shallow gradients. The high-order solution is obtained by \"antidiffusing\" the low-order, monotone solution, but only to such an extent that no new extrema are created and no existing extrema are accentuated. This is done by limit-ing, or \"correcting,\" the antidiffusive fluxes of the high-order scheme, hence the name Flux-Corrected Transport. A key feature of the FCT scheme is to use flux limiters to remove non-physical extremes. We use the Euler equations to illustrate the application of such flux limiters. dpib dFx, <9FX 9 <9FX, „ , . „ S 248 Appendix D. Flux Corrected Transport Using a second-order central time advance scheme, Eq. D . l can be dis-cretized temporally into PJ^ = pJ^_^dF\\ 2At 2At ^ dxi K ' ' i=i 1 For a Courant-Friedrichs-Lewy (CFL) number less than unity, the lower limit of updated scalar tp must satisfy £±%iP J cv where VWn is the local minimum of the transported scalar. A sufficient condition that satisfies Eq. D.3 is 2At Axip ^t + M - > / Font + (D.4) J cv where f Fout is the total flux leaving the control volume. The flux l im-iter ratio is thus defined as ^t+At _ ^t+At Ax~p lev F ° u t + Ip* A t oout r y min /r-> c\\ >J ~ 2At r T 7 i . / + _ A + \\ V J - D ) If f3°ut is greater than 1, then no flux correction is necessary; otherwise, the flux limiter should be applied to the out flux ^ = m m ( l , p™\\)max(0, F{) + m m ( l , p°ut)min(0, Fi) (D.6) where i is the node at which /3°ut is evaluated. In the current study, the FCT scheme is used to resolve the fine structures of the flow field as wel l as to enhance the numerical stability - the existence of a large pressure drop at the nozzle exit and non-linear nature of the chemical reaction may lead to large local gradients of reaction scalars, which may cause non-physical oscillations due to the numerical dispersion with ordinary second-order schemes. 249 Appendix E Numerical Schlieren The Schlieren image represents the density gradient integrated along the transverse direction of the test section. In a cylindrical coordinates, the transverse-direction density gradient can be obtained from dp(r, 9) ^dpdr | dp 86 • dx dr dx d6 dx The differential term with respect to 9 drops due to the axisymmetry of the flow field. We integrate Eq. E . l along the transverse direction of the test section Jr=R ordx Jr=R dr where / is proportional to the light intensity of the Schlieren image. Note that we only integrate dp in first quadrant to take advantage of the symmetry of the round jet. 250 "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2006-05"@en ; edm:isShownAt "10.14288/1.0080745"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mechanical Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Natural gas combustion under engine-relevant conditions"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/18215"@en .