CONDITIONAL SOURCE-TERM ESTIMATION METHODS FOR TURBULENT REACTING FLOWS by BEI JIN BE Tsinghua University, 2000 ME Tsinghua University, 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA December 2007 © Bei Jin, 2007 Abstract Conditional Source-term Estimation (CSE) methods are used to obtain chemical closure in turbulent combustion simulation. A Laminar Flamelet Decomposition (LFD) and then a Trajectory Generated Low-Dimensional Manifold (TGLDM) method are combin- ed with CSE in Reynolds-Averaged Navier Stokes (RANS) simulation of non-premixed autoigniting jets. Despite the scatter observed in the experimental data, the predictions of ignition delay from both meth- ods agree reasonably well with the measurements. The discrepancy between predictions of these two methods can be attributed to differ- ent ways of generating libraries that contain information of detailed chemical mechanism. The CSE-TGLDM method is recommended for its seemingly better performance and its ability to transition from au- toignition to combustion. The effects of fuel composition and injection parameters on ignition delay are studied using the CSE-TGLDM met- hod. The CSE-TGLDMmethod is then applied in Large Eddy Simulation of a non-premixed, piloted jet flame, Sandia Flame D. The adiabatic CSE-TGLDM method is extended to include radiation by introducing a variable enthalpy defect to parameterize TGLDM manifolds. The re- sults are compared to the adiabatic computation and the experimental data. The prediction of NO formation is improved, though the predic- tions of temperature andmajor products show no significant difference from the adiabatic computation due to the weak radiation of the flame. The scalar fields are then extracted and used to predict the mean spec- tral radiation intensities of the flame. Finally, the application of CSE in turbulent premixed combustion is explored. A product-based progress variable is chosen for condi- tioning. Presumed Probability Density Function (PDF) models for the progress variable are studied. A modified version of a laminar flame- based PDF model is proposed, which best captures the distribution of the conditional variable among all PDFs under study. A priori tests are performed with the CSE and presumed PDF models. Reaction rates of turbulent premixed flames are closed and compared to the DNS data. The results are promising, suggesting that chemical closure can be achieved in premixed combustion using the CSE method. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Chapter 1. Introduction and Thesis Outline . . . . . . . . . . . . . . . . . . . . 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Chapter 2. Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 CFD Tools for Turbulent Combustion . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Reynolds Averaged Navier Stokes . . . . . . . . . . . . . . . . . . 11 2.2.3 Large Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Non-premixed Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Infinitely Fast Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3.2 Laminar Flamelet Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.3 Conditional Moment Closure . . . . . . . . . . . . . . . . . . . . . . 24 2.3.4 Conditional Source-term Estimation . . . . . . . . . . . . . . . . 28 2.3.5 Including Chemistry in CSE . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Premixed Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.4.1 Eddy-Break-Up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4.2 The Bray-Moss-Libby Model . . . . . . . . . . . . . . . . . . . . . . . 40 2.4.3 Models Based on Surface Area Estimation . . . . . . . . . . 44 2.4.4 Models Based on the G-equation . . . . . . . . . . . . . . . . . . . 45 2.5 Reduction of Detailed Reaction Mechanism . . . . . . . . . . . . . . 46 2.5.1 Intrinsic Low-Dimensional Manifold . . . . . . . . . . . . . . . 47 2.5.2 Trajectory Generated Low-Dimensional Manifold . . 50 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 iii Chapter 3. IgnitionDelay ofMethane Jets with Conditional Source- term Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.2 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.2.1 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.3.1 Simulation of a Non-reactive Jet . . . . . . . . . . . . . . . . . . . . 57 3.3.2 Simulations with CSE-TGLDM and CSE-LFD . . . . . . 61 3.3.3 Effects of Different Fuel Composition . . . . . . . . . . . . . . 68 3.3.4 Effects of Different Injection Parameters . . . . . . . . . . . . 76 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Chapter 4. Modelling Radiative Effects of a TurbulentNon-premixed Flame using CSE-TGLDM . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2 LES of Sandia Flame D using CSE-TGLDM . . . . . . . . . . . . . . . 81 4.3 CSE-TGLDM with Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.1 Radiation Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4 Radiation Intensities Analysis using CSE . . . . . . . . . . . . . . . . 101 4.4.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Chapter 5. A Progress Variable and Presumed Probability Den- sity Function for Premixed Combustion . . . . . . . . . . . . 108 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.2 The PDF Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.2.1 β-PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.2.2 A Laminar Flame-based PDF . . . . . . . . . . . . . . . . . . . . . . 112 5.2.3 Modified Laminar Flame-based PDF . . . . . . . . . . . . . . 114 5.3 Comparison with DNS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3.1 DNS Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3.2 Comparisons of PDF Models . . . . . . . . . . . . . . . . . . . . . . 118 5.3.3 Effects of Chemical Kinetics . . . . . . . . . . . . . . . . . . . . . . . 132 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Chapter 6. Closure of Reaction Rates in Turbulent PremixedCom- bustion with CSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.2 The CSE Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 iv 6.3 Results and Discussion : RANS tests . . . . . . . . . . . . . . . . . . . . . 142 6.3.1 Conditional Averages: Mechanism I . . . . . . . . . . . . . . . 142 6.3.2 Reaction Rate Predictions: Mechanism I . . . . . . . . . . . 146 6.3.3 Comparison with DNS: Mechanism II . . . . . . . . . . . . . 147 6.4 Results and Discussion : LES tests . . . . . . . . . . . . . . . . . . . . . . . 151 6.4.1 Comparison with DNS: Mechanism I . . . . . . . . . . . . . . 151 6.4.2 Comparison with DNS: Mechanism II . . . . . . . . . . . . . 161 6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 Chapter 7. Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . 166 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 v List of Tables 3.1 Autoigniting jets with various fuel composition . . . . . . . . . . . . . . 61 4.1 Coefficients for Planck mean absorption coefficients. . . . . . . . . . 86 vi List of Figures 2.1 Regimes in non-premixed turbulent combustion . . . . . . . . . . . . 19 2.2 Experiment measurement of scattered and conditional average temperature of Sandia D-flame at x/d = 30 . . . . . . . . . . . . . . . . . . 25 2.3 Borghi diagram for premixed combustion . . . . . . . . . . . . . . . . . . . 35 2.4 Regimes in premixed turbulent combustion . . . . . . . . . . . . . . . . . 37 2.5 Bray-Moss-Libby model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6 Comparing one dimensional manifold generated by ILDM and TGLDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.1 Penetration length of a cold jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2 TGLDM for a methane/ethane jet . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 LFD library of a methane/ethane jet . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4 Maximum increase of 〈T |ζ〉 for methane/ethane jets . . . . . . . . 65 3.5 Comparison of CSE-TGLDM and CSE-LFD predictions . . . . . 67 3.6 Comparison of TGLDM and laminar flamelet library . . . . . . . . 68 3.7 Ignition Delay τ vs T0 of methane/nitrogen jets . . . . . . . . . . . . . 69 3.8 Ignition Delay τ vs T0 of methane and methane with additives jets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.9 Conditional averaged temperature evaluated using TGLDMs of different fuel composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.10 Comparison of Zst contours in different jets . . . . . . . . . . . . . . . . . 73 3.11 Temperature fields with Z contours at t = 0.9ms . . . . . . . . . . . . 74 3.12 Temperature fields with Z contours at t = 1.1ms . . . . . . . . . . . . 75 3.13 Ignition delay τ and pressure ratio . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.14 Ignition delay τ and injection duration ti . . . . . . . . . . . . . . . . . . . . 78 4.1 Centerline profiles of temperature and species . . . . . . . . . . . . . . 92 4.2 Temperature radial profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 Radial profiles of YCO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.4 Radial profiles of YH2O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.5 Radial profiles of YNO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.6 Conditional averaged CO2 at different locations . . . . . . . . . . . . . 97 4.7 Conditional averaged H2O at different locations . . . . . . . . . . . . . 98 4.8 Conditional averaged NO at different locations . . . . . . . . . . . . . . 99 4.9 Spectral radiation intensities in Flame D . . . . . . . . . . . . . . . . . . . 104 4.9 Spectral radiation intensities in Flame D (cont.) . . . . . . . . . . . . 105 4.10 Spectral radiation intensities in Flame D based on adiabatic and non-adiabatic computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 vii 4.10 Spectral radiation intensities in Flame D based on adiabatic and non-adiabatic computations (cont.) . . . . . . . . . . . . . . . . . . . . . . . . 107 5.1 Probability Density Function in a turbulent premixed flame 111 5.2 Progress variable in an unstrained laminar flame . . . . . . . . . . . 114 5.3 Probability Density Function (PDF) and Cumulative Distribution Function (CDF) in a premixed turbulent flame . . . . . . . . . . . . . 120 5.4 PDF on planes in a premixed turbulent flame . . . . . . . . . . . . . . 121 5.5 CDF on planes in a premixed turbulent flame . . . . . . . . . . . . . . 122 5.6 PDF in different cubes in a premixed turbulent flame . . . . . . 124 5.7 CDF in different cubes in a premixed turbulent flame . . . . . . 125 5.8 PDF and CDF in the whole domain . . . . . . . . . . . . . . . . . . . . . . . . 127 5.9 Test II: PDF on different planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.10 Test II: CDF on different planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.11 Test II: PDF in different cubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.12 Test II: CDF in different cubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.13 PDF and CDF in the whole domain using the right and wrong chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.14 PDF on planes using the right and wrong chemistry . . . . . . . 134 5.15 CDF on planes using the right and wrong chemistry . . . . . . . 135 5.16 PDF in different cubes using the right and wrong chemistry 136 5.17 CDF in different cubes using the right and wrong chemistry 137 5.18 Comparison of laminar flames with different chemistry . . . . 138 6.1 Conditional averages calculated in RANS paradigm . . . . . . . 144 6.2 Conditional averages calculated using β-PDF . . . . . . . . . . . . . . 145 6.3 Conditionally averaged chemical reaction rates calculated in RANS paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 6.4 Conditionally averaged chemical reaction rates calculated in RANS paradigm using β-PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 6.5 Reaction rates on planes in a premixed turbulent flame . . . . 150 6.6 Test II: Conditional averages calculated in RANS paradigm 152 6.7 Test II: Conditionally averaged chemical reaction rates calculated in RANS paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.8 Test II: Chemical reaction rates on planes in a premixed turbulent flame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 6.9 Conditional averages calculated in LES paradigm . . . . . . . . . . 155 6.10 Conditionally averaged chemical reaction rates calculated in LES paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 6.11 Chemical reaction rates calculated in LES paradigm . . . . . . . . 158 6.12 Chemical reaction rates calculated in LES paradigm, with denser grid for conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 viii 6.13 Chemical reaction rates calculated in LES paradigm, conditional means extracted from DNS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 6.14 Test II: Conditional averages calculated in LES paradigm . . 162 6.15 Test II: Conditionally averaged chemical reaction rates calculated in LES paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 6.16 Test II: Chemical reaction rates calculated in LES paradigm 164 ix Nomenclature Symbol Meaning Units c Progress variable deq Equivalent diameter m D Diffusivity m2· s−1 Da Damköhler number m2· s−1 E Activation Energy cal·mol−1 F Body force h Specific Enthalpy J·kg−1 I Identity Matrix J Jacobian Matrix J h Enthalpy diffusivity J/(m2·s) J k Species molecular diffusivity kg/(m2·s) Ka Karlovitz Number lF Premixed flame length scale m lref Non-premixed reference length scale m p Pressure bar P Probability Density Function Pr Prandtl number Re Reynolds number sL Laminar flame velocity m/s Sc Schmidt number τc Chemical time scale µs,ms or s tF Premixed flame time scale µs,ms or s tt Turbulent time scale µs,ms or s tk, tη Kolmogorov time scale µs,ms or s uref Non-premixed reference velocity m/s vη Kolmogorov velocity scale m/s v velocity vector m/s V Eigen(Schur) Basis Matrix T Temperature K x Cartesian coordinates Y Mass Fraction x Symbol Meaning Units Z Mixture Fraction 〈Y |c = c∗〉 Conditional Average of Ywith Condi- tion c = c∗ 〈Y |Z = ζ〉 Conditional Average of Ywith Condi- tion Z = ζ Z ′′2 Variance of Mixture Fraction BML Bray-Moss-Libby Model CDF Cumulative Distribution Function CMC Conditional Moment Closure CSE Conditional Source-term Estimation DNS Direct Numerical Simulation ILDM Intrinsically Low Dimensional Manif- old LES Large Eddy Simulation LFD Laminar Flamelet Decomposition PDF Probability Density Function QSSA Quasi-steady State Approximation RANS Reynolds Averaged Navier-Stokes RTE Radiative Transfer Equations SFL Steady Flamelet Library SPM Stochastic Particle Model TGLDM Trajectory-Generated Low Dimensio- nal Manifold TRI Turbulent-Radiation Interaction UFL Unsteady Flamelet Library η Kolmogorov length scale m Λ Diagonal Matrix of Eigenvalues s−1 ρ Density kg·m−3 τ Ignition Delay or Time Duration µs,ms or s ω̇ Chemical Source Term s−1 χ Scalar Dissipation Rate s−1 ζ Sample space for mixture fraction xi Acknowledgments The thesis work described here would not be possible without the sup- port and help of many people. I would like to thank all my teachers, old and new friends in China and Canada who have seen me through this long journey. I would like to thank Drs. M. H. Davy, J. J. Feng and S. N. Ro- gak, who kindly agree to be my research committee despite their busy schedules. I would like to thank my colleagues and friends in my research group. Thanks to Drs. Jian Huang and Mei Wang for their help on my project. Thanks to Drs. Gord McTaggart-Cowan and Andrea Frisque and Heather Jones for being supportive and encouraging all along. I would like to thank my supervisor, Dr. W. Kendal Bushe, for his advice, patience and continuous support in every way throughout my project. My PhD study has been a great experience thanks to the pre- cious opportunities and guidance he gave me. Special thanks to my family - my parents, sister, brother-in-law and nephew, who always believe in me and give me their unconditional love. Lastly, I would like to dedicate this work to my father, who inspired and motivated me every step of the way. I wish my efforts would have made him proud and I will always cherish the memory of him. xii Chapter 1 Introduction and Thesis Outline 1.1 Introduction The last century has witnessed soaring gas prices, deteriorating air quality and alarming global climate changes [1]. In recent years, in- creasing concerns have been raised with respect to the environmental impacts of energy consumption. For centuries, the combustion of fos- sil fuels has been a major form of energy consumption that has been widely used in stationery power generation and transportation indus- tries. Transportation in particular has been identified to be one of the largest sources of energy consumption, greenhouse gases and air pol- lutants emissions. As a result, governments now set more and more stringent standards [2, 3] for vehicles powered by internal combustion engines. Great efforts have been made to improve fuel efficiency and reduce harmful emissions of traditional IC engines. One promising so- lution is to substitute conventional fuels such as gasoline and diesel with clean-burning fuels. Natural gas is a viable alternative fuel to replace diesel in compres- sion ignition engines because of its relative abundance, lower emis- sions and lower cost [4]. In terms of availability, Canada is one of the world’s largest producers of natural gas [4]. As for emission, natural 1 gas as a transportation fuel has the potential of significantly reducing harmful emissions. For example, it has been demonstrated [5, 6] that medium- and heavy-duty natural gas engines have over 90% reduc- tions of carbon monoxide (CO) and particulate matter and more than 50% reduction in nitrogen oxides (NOx) relative to commercial diesel engines. Natural gas has high antiknock performance, eliminating the need to add antiknock ingredients to the fuel. It is safe to use if han- dled properly. Conventional diesel engines can be converted to run on natural gas with relatively minor changes while maintaining similar maximum horsepower. Despite these benefits, using natural gas does not automatically guarantee lower emissions. The performance of the engine depends on the type of engine technology deployed and the implementation. Currently, the most common use of natural gas as a transportation fuel is in the form of compressed gas (CNG). Spark-ignited CNG en- gines involve mixing low-pressure natural gas with air before the mix- ture enters the combustion chamber. Alternatively, pilot diesel may be injected for ignition. Recently, High Pressure Direct Injection (HPDI) technology enables direct injection of CNG into the cylinder, allowing diesel engines to operate on clean-burning CNG. To implement these technologies properly, it is crucial to understand the ignition and com- bustion of natural gas. This has motivated recent studies of the com- bustion of methane and methane mixed with additives. Among different research approaches such as experimental and the- oretical studies, numerical simulation provides an efficient way of ana- lyzing combustion phenomena of natural gas-like mixtures. Turbulent 2 combustionmodelling is known to be challenging. It involves a variety of closure problems including the traditional problems of turbulence modelling and unclosed terms caused by the coupling between differ- ent mechanisms such as turbulence/chemistry and turbulence/radia- tion interactions. It also involves computational challenges such as the numerical solution of stiff ODE systems and the incorporation of com- plicated chemical kinetics in CFD. This thesis focuses on the numerical study of turbulent combustion. The objective is to apply the Conditional Source-term Estimation (CSE) method [7] to study properties of non-premixed and premixed com- bustion of natural gas-like mixtures. Although the cases under study all include methane as a major component in fuel, the modelling tech- nique can be applied to other gaseous fuels with proper adaption. 1.2 Outline The thesis starts with a comprehensive review in Chapter 2 to discuss combustion simulation subjects relevant to the thesis work. Basic tools of CFD are described first. Different combustion models are then in- troduced and their applications in non-premixed and premixed com- bustion are discussed. Some commonly used techniques for including detailed chemistry in CFD are introduced and the chapter ends with a list of questions to be studied in the thesis work. In Chapter 3, Reynolds Averaged Navier-Stokes (RANS) simula- tions of transient jet flames are performed to understand the autoigni- tion of methane under engine-relevant conditions. Two different meth- ods are combined with CSE to incorporate full chemical mechanisms 3 and their performances are compared head-to-head. The CSEwith Tra- jectory Generated LowDimensionalManifold (TGLDM)method is rec- ommended and used thereafter to study non-premixed jet flames. The predictions of ignition delay are compared to experimental data; the effects of additives and injection parameters on the autoignition of me- thane are analyzed. In Chapter 4, Large Eddy Simulation (LES) of a non-premixed steady jet flame is performed with the recommended CSE-TGLDM method. An extended version of the CSE-TGLDM method is proposed to in- clude a radiation model in the simulation. The effects of radiation heat loss on the temperature field and species formation are analyzed. The mean spectral radiation characteristics are then analyzed in a post- processing fashion where the effectiveness of using the CSE method to model turbulence/radiation interaction is further proved. Chapters 5 and 6 discuss the application of the CSE method in pre- mixed flames. In Chapter 5, a conditioning variable is introduced and the presumed probability density function (PDF) of this variable is studied. PDF models from previous studies are reviewed and an ex- tended laminar flame based PDF is proposed. The distribution pre- scribed by different PDFs are compared to Direct Numerical Simula- tion (DNS) results. The precision of PDF models is analyzed and the importance of chemical kinetics in the PDF computation is discussed. In Chapter 6, a priori tests of CSE with various presumed PDF mod- els are performed for turbulent premixed flames. Reactive scalars such as conditional averages of species are evaluated and chemical reaction rates are closed using the CSE method; comparison to the DNS data 4 shows the potential of the CSE method as a tool to get closure in tur- bulent premixed combustion. 5 Chapter 2 Literature Review 2.1 Introduction Turbulent reacting flows are widely observed in nature and in indus- trial devices. Numerical simulation of turbulent combustion is attrac- tive because it is affordable and informative. Great progress has been made thanks to the development of modern computer technology that enables the simulation of complex turbulent combustion phenomena. Despite this progress, it is still challenging to simulate turbulent re- acting flows. A wide variety of problems are coupled and need to be solved: • It is computationally expensive to directly solve for all the scales in a turbulent flow, which introduces the traditional closure prob- lem of turbulence modelling. • It is important to use detailed chemical kinetics to describe com- bustion processes, especially in the study of ignition, extinction and pollutant formation. Such mechanisms usually include large number of species and reactions over a wide range of time scales, which often leads to a large system with extreme stiffness. • In some cases, two or three phases of mixture might exist in the 6 system, which requires modelling of complicated interactions be- tween different phases. For example, to study the liquid fuel in- jection and ignition in diesel engines, the breakdown and vapor- ization of the liquid droplets, turbulent mixing and combustion need to be modelled. • Radiation as an inherent heat transfer mechanism influences the flame structure, soot and NOX formation. Radiation modelling is known to be complicated; yet accurate calculation of radiation effects is necessary especially in luminous flames. • The problems listed above are coupled in turbulent combustion simulation due to the interactions between turbulence and other mechanisms. For example, turbulence/chemistry and turbulence /radiation interactions lead to typical closure problems in turbu- lent combustion modelling. Clearly, turbulent combustion modeling is a complicated and broad subject. This chapter will address a few aspects of this subject that are relevant to the thesis work. The review starts with basic Compu- tational Fluid Dynamics (CFD) methods, then introduces some com- bustion models for non-premixed and premixed turbulent combustion and ends with a discussion of different methods used to incorporate detailed chemical kinetics in CFD. 2.2 CFD Tools for Turbulent Combustion Numerical study of turbulent combustion involves the solution of ba- sic balance equations such as theNavier-Stokes, continuity, species and 7 energy equations. Given the chemical mechanism and transport mod- els, such equations may be solved directly without modelling turbu- lence or interactions between turbulence and other mechanisms. Al- ternatively, the equations may be averaged or filtered, resulting in un- closed terms that need to be modelled. Major CFD tools for turbulent combustion simulation include Direct Numerical Simulation (DNS), Reynolds-Averaged Navier-Stokes (RANS) methods and Large Eddy Simulation (LES). 2.2.1 Direct Numerical Simulation As indicated by the name, DNS method directly solves the governing equations of turbulent reactive flows [8]: • Mass conservation (continuity): ∂ρ ∂t + ∂ρuj ∂xj = 0, (2.1) where ρ and uj are density and velocity respectively. • Momentum: ∂ρui ∂t + ∂ρujui ∂xj = − ∂p ∂xi + ∂τij ∂xj + Fi, (2.2) where τij is the viscous stress tensor and Fi is a body force. • Species transport equation: ∂ρYk ∂t + ∂ρujYk ∂xj = −∂Jj k ∂xj + ω̇k, (2.3) where Yk is the mass fraction of the species k, Jj k is the species diffusive flux and ω̇k is the species chemical reaction rate. 8 • Energy equation: ∂ρht ∂t + ∂ρujht ∂xj = ∂p ∂t + ∂ ∂xj (Jj h + uiτij) + ujFj, (2.4) where ht is the total specific enthalpy defined by ht = h+ uiui/2, Jj h is the enthalpy diffusion, uiτij and uiFj denote the effects of viscous and body forces respectively. For Newtonian fluids, τij is defined by: τij = µl( ∂ui ∂xj + ∂uj ∂xi )− 2 3 µlδij( ∂uk ∂xk ), (2.5) where µl is the molecular viscosity. Species molecular diffusivity Jj k can be approximated by Fick’s law: Jj k = − µl Sck ∂Yk ∂xj . (2.6) Here Sck is the Schmidt number of species k defined by Sck = µl ρDk ; it characterizes the simultaneous momentum and mass diffusion con- vection processes in turbulent flows. Dk is the molecular diffusivity of species k relative to the major species. The enthalpy diffusion Jj h can be approximated by Fourier’s law: Jj h = − µl Pr [ ∂h ∂xj + N∑ k=1 ( Pr Sck − 1)hk ∂Yk ∂xj ] , (2.7) where N is the number of species. The Prandtl number, Pr, compares the diffusive transport of momentum and temperature; it is defined using the thermal diffusivity λ and constant pressure specific heat Cp as Pr = µCp λ . The Lewis number of species k compares thermal and mass diffu- sivities and is defined by: Lek = Sck Pr = λ ρCpDk . (2.8) 9 Eq. 2.7 can be simplified under the assumption of unity Lewis num- ber; Eq. 2.3 and 2.4 are formally identical under the assumption of low Mach number. These assumptions may be applied to simplify turbu- lent combustion modelling under proper circumstances. Besides those numbers already discussed, listed below are some other dimensionless numbers that are commonly used to characterize turbulent reactive flows. The turbulent Reynolds number is defined as: Re = u′lt ν , (2.9) where u′ is the velocity root mean square related to turbulent kinetic energy k by u′ = √ 2/3k, lt is the turbulence integral length scale and ν is the kinematic viscosity. The Reynolds number, Re, compares tur- bulent transport with viscous forces. It indicates the ratio of the large scales of turbulent motions to the Kolmogorov microscales. The turbulent Damköhler number is defined by: Da = τt τc , (2.10) where τt and τc are the turbulent and chemical time scale respectively. The Damköhler number measures how important the turbulence/ch- emistry interaction is in comparison to other processes. The Karlovitz number is used to characterize interactions between turbulent motions and flame structure; it is defined as: Ka = τc τk , (2.11) where τk ≡ ( ν ǫ )1/2 is the Kolmogorov time scale, ǫ is the average rate of energy dissipation per unit mass. 10 In DNS, the balance equations Eq. 2.1- 2.7 are solved directly after discretization. Meshes used in DNS need to be fine enough to resolve the smallest eddies in the turbulent flow. DNS of turbulent combustion is often prohibitively expensive. Firstly and mostly, the computational cost is high to resolve all time and length scales which might differ by orders of magnitude. To simulate a real flow with large Re, the num- ber of cells needed can be overwhelming [9]. Moreover, the number of scalars that need to be solved for in turbulent combustion is often very large. A detailed chemical kinetic system might include tens or hun- dreds of species and reactions. Lastly, the boundary and initial condi- tions of a practical system are often difficult to define in DNS. Con- sequently, DNS is mostly used to analyze turbulent flames in simple configurations [10–12] such as homogeneous isotropic turbulent flows and mixing layers. 2.2.2 Reynolds Averaged Navier Stokes Reynolds-averagedNavier-Stokes (RANS) approaches obtain themean flow field by solving the averaged governing equations. Each quantity Q is decomposed into a mean value Q and a deviation Q′ away from the mean: Q = Q+Q′, (2.12) where Q′ = 0. The instantaneous balance equations may be ensemble or time averaged to derive the transport equations forQ. This Reynolds averaging procedure leads to unclosed terms such as u′Q′ that need to be modelled. This method has been widely used in simulations of non- reacting turbulent flows. In turbulent combustion, because of the fluctuations in density caused 11 by heat release, Reynolds averaging generates unclosed terms includ- ing fluctuation correlations involving density. To avoid explicitly mod- elling such correlations, a Favre (mass weighted) averaging of a quan- tity Q is introduced: Q = Q̃+Q′′, (2.13) where Q̃ = ρQ ρ and Q̃′′ = ρ(Q− eQ) ρ = 0. The Favre averaged equations are: • Mass conservation (continuity): ∂ρ ∂t + ∂ρũj ∂xj = 0 (2.14) • Momentum: ∂ρũi ∂t + ∂ρũjũi ∂xj = −∂ρũ ′′ i u ′′ j ∂xj − ∂p ∂xi + ∂τ̃ij ∂xj + F̃i (2.15) • Species transport equation: ∂ρỸk ∂t + ∂ρũjỸk ∂xj = −∂ρũ ′′ jY ′′ k ∂xj − ∂Jj k ∂xj + ω̇k (2.16) • Energy equation: ∂ρh̃t ∂t + ∂ρũjh̃t ∂xj = −∂ρũ ′′ jh ′′ t ∂xj + ∂p ∂t + ∂ ∂xj (Jj h+uiτij)+ujFj. (2.17) Unclosed terms are evaluated by modelling or deriving balance equations and solving them. The Reynolds stresses ũ′′i u ′′ j are oftenmod- elled using turbulence models developed for non-reacting flows, such as the k − ǫ model [13], without explicitly including heat release ef- fects. Species (ũ′′jY ′′ k ) and temperature (ũ ′′ jT ) turbulent fluxes are usu- ally modelled under a gradient transport hypothesis: ρũ′′jY ′′ k = − µt Sckt ∂Ỹk ∂xj , (2.18) 12 where µt is the turbulent viscosity estimated with turbulence models and Sckt is a turbulent Schmidt number of species k. At high Reynolds numbers, laminar diffusive fluxes such as Jkj and J h j are usually small compared to turbulent transport. Mean species chemical reaction rates ω̇k need to be closed with turbulent combustion models, which will be discussed in greater detail later in this chapter. RANS simulation solves for the averaged quantities, providing the statistical means of scalars associatedwith turbulent flames. Because of the strong unsteady mixing effects observed in turbulent combustion, such mean properties are not always sufficient to describe the flames. As an alternative, Large Eddy Simulation may be used. 2.2.3 Large Eddy Simulation Large Eddy Simulation [14] is a compromise between extremely ex- pensive DNS and (perhaps overly) simple RANS methods. It explic- itly computes the large scales that are larger than the mesh size, while modelling the small scales. Although its application in turbulent com- bustion [15] is still at an early stage, it is a promising method because of the following features: • Compared to DNS, LES only solves for large structures, which is computationally more affordable. • The unsteady large scale mixing is simulated, instead of being averaged as in RANS. • LES might capture the instabilities caused by the interactions be- tween heat release, hydrodynamic flow and acoustic waves. 13 In LES, a quantity Q is filtered in the spectral space where compo- nents larger than a given cut-off frequency are suppressed, or in the physical space where weighted averaging of the quantity is performed in a given volume. The filtering operation is defined by: Q(x) = ∫ Q(x∗)F (x− x∗)dx∗ (2.19) where F is the LES filter, normalized such that:∫ +∞ −∞ ∫ +∞ −∞ ∫ +∞ −∞ F (x1, x2, x3)dx1dx2dx3 = 1. (2.20) Typical examples include a cut-off filter in the spectral space, a box fil- ter in the physical space and a Gaussian filter in the physical space [14]. In turbulent combustion, a mass weighted, Favre filtering is usually used: ρQ̃(x) = ∫ ρQ(x∗)F (x− x∗)dx∗ (2.21) Any quantity Q may be decomposed into a filtered component Q and a fluctuating component Q′ such that Q = Q+Q′. The filtered balance equations are: • Mass conservation (continuity): ∂ρ ∂t + ∂ρũj ∂xj = 0 (2.22) • Momentum: ∂ρũi ∂t + ∂ρũjũi ∂xj = − ∂ ∂xj [ρ (ũiuj − ũiũj)]− ∂p ∂xi + ∂τij ∂xj + Fi (2.23) • Species transport equation: ∂ρỸk ∂t + ∂ρũjỸk ∂xj = − ∂ ∂xj [ ρ ( ũjYk − ũj ) Ỹk ] + ω̇k (2.24) 14 • Energy equation: ∂ρh̃t ∂t + ∂ρũjh̃t ∂xj = − ∂ ∂xj [ ρ ( ũjht − ũjh̃t )] + ∂p ∂t + ∂ ∂xj ( Jj h + uiτij ) +ujFj. (2.25) The unclosed terms include Reynolds stresses, the species and en- thalpy fluxes, filtered laminar diffusion fluxes and filtered chemical re- action rates. Reynolds stresses are usually evaluated using a subgrid scale (SGS) turbulencemodel. Common SGSmodels include Smagorin- sky models [16] and dynamic Smagorinsky models [17–19]. In a stan- dard dynamic procedure, the coefficients of the SGS model are solved by assuming that they are the same at the grid and test-filter level. The dynamic subgrid models usually have better performance. LES has been successfully used to study a wide variety of turbu- lent reacting flows. Since neither RANS nor LES resolves the smallest scales, at which the rate-controlling molecular mixing processes and chemical reactions occur, the formulations of these two methods both require modelling of these scales. Therefore, most of the RANS com- bustionmodelsmay bemodified and applied in LES. A comprehensive review of LES for turbulent combustion is presented by Pitsch [20]. 2.3 Non-premixed Combustion Non-premixed combustion is observed inmany industrial applications, where fuel and oxidizer are initially separated; both mixing and burn- ing processes occur in the combustion chamber. Examples of non-pre- mixed combustion include diesel engines and furnaces where fuel is supplied separately. In non-premixed combustion, the time for turbu- lent mixing is usually much longer than the chemical time scale [21]. 15 As a result, the mixing process is often the rate-controlling process. It is observed that the reaction rates of non-premixed combustion are closely related to the state of mixing [8]. A conserved scalar, the mixture fraction Z, is often used to describe the state of mixing in non-premixed flames. The mixture fraction can be defined in various ways. Different definitions of Z [21] essentially indicate the fraction of themixture that originated from the fuel stream. It is normalized such that Z = 0 in the oxidizer stream and Z = 1 in the fuel stream. One definition is given by Bilger [22]: Z = Yi − Yi,O Yi,F − Yi,O , (2.26) where Yi is the mass fraction of element i in the local mixture, and Yi,F and Yi,O are the the mass fractions of element i in the fuel and oxidizer streams respectively. Assuming all Lewis numbers to be unity, Z is independent of the selection of element i. The mixture fraction Z is a conserved scalar that is not affected by chemical reactions. The balance equation of Z is: ∂ρZ ∂t + ∂ρuiZ ∂xi = ∂ ∂xi ( ρDZ ∂Z ∂xi ) . (2.27) The diffusion coefficientDZ can be approximated with the thermal dif- fusivity. The mixture fraction Z is attractive because there is no chemi- cal source term in the balance equation. It is found that in non-premixed flames, chemical reaction rates are crucially related to mixture fraction. Many non-premixed models are developed based on Z. The means and variance of mixture fraction describe the mixing in non-premixed reacting flows. These two variables Z̃ and Z̃ ′′2 might be modelled or solved for using their transport equations. For example, 16 in RANS simulation, the following transport equations can be derived: ∂ρZ̃ ∂t + ∂ρũjZ̃ ∂xj = ∂ ∂xj ( ρD ∂Z ∂xj − ρũ′′iZ ′′ ) , (2.28) ∂ρZ̃ ′′2 ∂t + ∂ρũjZ̃ ′′2 ∂xj = ∂ ∂xj ( ρũ′′iZ ′′2 ) − 2ρũ′′iZ ′′ − ρχ̃. (2.29) Here the flux ũ′′iZ ′′ might be modelled as: ũ′′iZ ′′ = −Dt ∂Z̃ ∂xi ; (2.30) The mean scalar dissipation rate defined by: χ̃ = ˜2D| ▽ Z|2 (2.31) needs to be modelled. Scalar dissipation rate is important in the study of non-premixed combustion. It also appears in the laminar flamelet equations that will be discussed in detail later. In laminar flows, the instantaneous dissi- pation rate χ is used: χ = 2D| ▽ Z|2 (2.32) A diffusion time scale τχ may be defined using the inverse of the dis- sipation rate χ and this is often used to characterize extinction of non- premixed flames. To study non-premixed turbulent combustion phenomena, differ- ent regime diagrams have been developed using dimensionless num- bers. Some diagrams useDa andRe [23] to characterize different regimes. At largeDa numbers the chemistry is sufficiently fast and the flame has a structure similar to the laminar flame; at smallDa numbers where the chemical time scale is relatively large, extinction might occur. 17 Peters proposed a regime diagram to include the effects of turbu- lent mixing [21]. The diagram is defined using the ratio χq/χ̃st and Z ′st/(△Z)F . Here, χq is the extinction scalar dissipation rate, χ̃st is the conditional Favre mean scalar dissipation rate at the stoichiometric mixture fraction, Z ′st is the fluctuation around the mean value at the stoichiometric mixture fraction, (△Z)F is the diffusion thickness in the Z space defined as: (△Z)F = | ▽ Z|stlD. (2.33) The diffusion length lD in physical space is defined by: lD = ( Dst α )1/2 , (2.34) where Dst is the diffusion coefficient and α is the strain rate. (△Z)R is a reaction zone thickness in mixture fraction space defined as: (△Z)R = ε(△Z)F , (2.35) where ε is a scaling factor. It is found [21] to be proportional to χ 1/4 st in a methane-air diffusion flame with a four-step mechanism and χ 1/3 st in a diffusion flame with a one-step mechanism. As shown in Fig.2.1, extinction occurs at high scalar dissipation rates where χq/χ̃st < 1. Beyond this, three regimes are defined accord- ing to the fluctuations in mixture fraction. In the separated flamelets regime, the fluctuations extend to sufficiently lean and rich mixture and the diffusion layers around the reaction zone are separated. In the connected flame zones regime, the fluctuation is relatively small due to intense mixing or partial premixing of the fuel stream. In the connected reaction zones regime, the mixture fraction fluctuations are 18 Figure 2.1: Regimes in non-premixed turbulent combustion. smaller than the thickness of the reaction zones and the mixture frac- tion field is almost homogeneous. The contour of the mean stoichio- metric mixture fraction in a jet diffusion flame is imposed on the plot, which shows the relevance between different zones and the local con- ditions in the jet. 2.3.1 Infinitely Fast Chemistry In non-premixed turbulent combustion, the time needed for convec- tion and diffusion is usually much longer than that of chemical reac- tions. Therefore, it is appropriate to assume infinitely fast chemistry in the study of global properties [24]. With this assumption, the reactive scalars in the flame can be modeled as a function of Z. For example, non-premixed flames can be modelled using a pre- sumed Probability Density Function (PDF) and the infinitely fast che- 19 mistry model (IFCM). Assuming an infinitely fast single step chemical reaction [25]: F +O → P, (2.36) the reactive scalars are related to the mixture fraction by: YF = Y IFCM F (Z); YO = Y IFCM O (Z); T = T IFCM(Z). (2.37) Given a presumed PDF of mixture fraction, the averages of reactive scalars can be evaluated by integration: Ỹ IFCMF = ∫ 1 0 Y IFCMF (ζ)P̃ (ζ;x, t)dζ (2.38) Ỹ IFCMO = ∫ 1 0 Y IFCMO (ζ)P̃ (ζ;x, t)dζ (2.39) T̃ IFCM = ∫ 1 0 T IFCM(ζ)P̃ (ζ;x, t)dζ, (2.40) where ζ is the sample space of mixture fraction, P̃ (ζ;x, t) is the pre- sumed PDF at (x, t). A β-PDF prescribed by Z̃ and Z̃ ′′2 is usually used to model the air/fuel mixing in non-premixed combustion [26]. The means and variances of mixture fraction can be obtained by solving their transport equations. The infinitely fast chemistry with a presumed PDF is attractive be- cause it calculates the averages of reactive scalars without solving their transport equations. A chemical equilibrium condition can be adopted in case of multiple-step chemistry [27]. The IFCM method provides valuable information about the global flame structure and the maxi- mum possible heat release. However, in turbulent flows where the lo- cal diffusion time scales are not so large, the fast chemistry assump- tion is no longer appropriate and non-equilibrium effects need to be included. Moreover, relatively slow chemistry is poorly predicted in 20 the IFCM method; therefore, it should not be used to study formation of particulate matter, CO or NOx. 2.3.2 Laminar Flamelet Model Flamelet models [28] introduce the scalar dissipation rate χ as a second parameter besides Z to account for the departure from IFCM equilib- rium states. It is assumed that the local balance between diffusion and reaction in a turbulent reactive flow is similar to that of a laminar flame with the state of mixing characterized by the same Z and χ. The solu- tions of such a prototype laminar flame can be obtained by solving the flamelet equations. The flamelet equations for non-premixed laminar flames are de- rived by applying the following transformation: ∂ ∂t = ∂ ∂τ + ∂Z ∂t ∂ ∂Z , (2.41) ∂ ∂x1 = ∂Z ∂x1 ∂ ∂Z , (2.42) ∂ ∂xk = ∂ ∂Zk + ∂Z ∂xk ∂ ∂Z (k = 2, 3), (2.43) where (x1, x2, x3, t) is the original coordinate system in the physical space; (Z1, Z2, Z3, τ) is the new coordinate system, defined as Z1 = ▽Z |▽Z| being locally normal to the stoichiometric flame surface, Z2 and Z3 be- ing tangential to the stoichiometric isosurface, and τ = t. Assuming unity Lewis numbers for all species, the balance equations of the lami- nar flame are rewritten in the mixture fraction space as below: ρ ∂Yi ∂τ = ρχ 2 ∂2Yi ∂Z2 + ω̇i −R(Yi), (2.44) ρ ∂T ∂τ = ρχ 2 ∂2T ∂Z2 − n∑ i=1 hi cp ω̇i + 1 cp {∂p ∂t + qR} −R(T ), (2.45) 21 where qR is the radiation heat loss. The flamelet equations are sim- plified by neglecting the radiation and pressure fluctuation terms in Eq. 2.45, as well as R(Yi) and R(T ), which are of a lower order com- pared to the other terms. The simplified flamelet equations are: ρ ∂Yi ∂τ = ρχ 2 ∂2Yi ∂Z2 + ω̇i, (2.46) ρ ∂T ∂τ = ρχ 2 ∂2T ∂Z2 − n∑ i=1 hi cp ω̇i. (2.47) The coordinate system of equations is reduced to (Z, τ) after the trans- formation. The scalar dissipation rate χ is used to relate the laminar flame solutions to the turbulent flow. A typical model for χ takes the form χ = χ0χ ∗(Z), where χ∗(Z) is some presumed function describing the relationship between Z and χ. Peters [28] developed an analytical solution of χ∗(Z) using the configuration of an infinite one dimensional mixing layer: χ∗(Z) = exp{−2[erfc−1(2Z)]2} exp{−2[erfc−1(2Z0)]2} . (2.48) Pitsch and Peters [29] obtained an alternate functional form for χ∗ us- ing a symmetric one-dimensional mixing layer configuration: χ∗(Z) = Z2lnZ Z20 lnZ0 . (2.49) Solutions of the flamelet equations are usually calculated in a pre- processing fashion and stored as flamelet libraries. Using such libraries, it is feasible to incorporate detailed chemical mechanisms in CFD com- putations. In flamelet libraries, reactive scalars are often tabulated for given values of the two control parameters as Yi(Z, χ). The solutions of flamelet equations without the time dependent terms are used in the Steady Laminar Flamelet Model (SLFM). The mean properties of a tur- bulent flame can be obtained by using a SLFM library combined with 22 a joint PDF model of Z and χ: Ỹi = ∫ Z∗ ∫ χ∗ Y SLFMi (Z ∗, χ∗)P̃ (Z∗, χ∗;x, t)dχ∗dZ∗, (2.50) where Y SLFMi (Z ∗, χ∗) is the SLFM solution and P̃ (Z∗, χ∗;x, t) is the pre- sumed joint PDF. SLFM models have been successfully used to study finite-rate che- mistry problems. However, neglecting time-dependence on the dissi- pation rate makes it insufficient in the study of transient phenomena such as extinction and reignition. The effects of unsteadiness are in- cluded in Unsteady Laminar Flamelet Models (ULFM) [30–32]. Two different approaches, i.e. the Largrangian Flamelet Model (LFM) and the Eulerian Particle FlameletModel (EPFM) have been usedwithULF- M. For example, Mauss et al. [33] studied extinction and re-ignition in a turbulent jet flame using unsteady flamelets. In his work, a Lagrangian time measured along the stoichiometric interface was used to account for the effects of time-evolving χ. Barths et al. [34, 35] used unsteady flamelet solutions by tracing imaginary marker particles. In his me- thod, flamelets were attached to these particles, which were used to trace different histories of scalar dissipation rates in the flow; a mod- eled Eulerian convective-diffusive equation was solved to evaluate the probability of finding these particles at specific locations. The flamelet models considerably relax the hypothesis required in the infinite fast chemistry model and have been widely used in non- premixed combustion modelling. Efforts have been made to improve flamelet models, such as including differential diffusion effects [36,37] and including diffusion in the direction tangential to the isosurface of Zst [38], etc. 23 2.3.3 Conditional Moment Closure Conditional Moment Closure (CMC) methods were first proposed by Bilger [39] and Klimenko [40] independently. CMC is based on the observation that in non-premixed systems, the reaction rates are non- linearly and crucially dependent on the stoichiometry; the fluctuations around mean properties are significantly reduced after conditioning with a stoichiometry variable. Therefore, accurate closure might be achieved in non-premixed combustion by using conditional moments. Mixture fraction Z is used as the conditioning variable in non-pre- mixed combustion. In the absence of extinction and reignition, condi- tioningwithZ is usually sufficient to reduce the impact of nonlinearity. For illustration, temperature at the station x/d = 30 in Sandia flame D [41] is plotted against mixture fraction Z in Fig. 2.2. The scattered data points all lie within close vicinity of the conditional average of temperature. Bilger [39, 40, 42] developed CMC equations with a decomposition approach. In his work, the conditional average of a scalar Y (x, t) is defined as: Q(ζ, x, t) ≡ 〈Y (x, t)|Z(x, t) = ζ〉 ≡ 〈Y |ζ〉, (2.51) where the angular brackets is an ensemble average over different real- izations of the flow and the vertical bar means “given the condition”. The basic decomposition of Y (x, t) can be expressed as: Y (x, t) = Q(Z(x, t), x, t) + Y ′′(x, t), (2.52) where Q is the conditional average of Y given Z(x, t) = ζ and Y ′′(x, t) is the fluctuation around the conditional means. The derivatives of Y 24 0 0.2 0.4 0.6 0.8 1 0 500 1000 1500 2000 2500 Mixture fraction Te m pe ra tu re (k) Figure 2.2: Experiment measurement of scattered and conditional av- erage temperature of Sandia D-flame at x/d = 30 (grey dots: scattered data; solid black line: conditional average temperature). 25 can be rewritten as: ∂Y ∂t = ∂Q ∂t + ∂Q ∂ζ ∂ζ ∂t + ∂Y ′′ ∂t , (2.53) ∇Y = ∇Q+ ∂Q ∂ζ ∇ζ +∇Y ′′. (2.54) The molecular diffusion term is given by: ∇ · (ρD∇Y ) = ∇ · (ρD∇Q) + ∂Q ∂ζ · (ρD∇Z) (2.55) + ρD(∇Z)2∂ 2Q ∂ζ2 + ρD∇Z · ∇∂Q ∂ζ +∇ · (ρD∇Y ′′). The CMC equation can be derived by substituting the above deriva- tives into the balance equations of Y (x, t); leading to: 〈ρ|ζ〉∂Q ∂t + 〈ρ|ζ〉〈uj|ζ〉 ∂Q ∂xj −〈ρ|ζ〉〈χ|ζ〉 2 ∂2Q ∂ζ2 (2.56) = 〈ρ|ζ〉〈ω̇|ζ〉+ eQ + eY , where 〈χ|ζ〉 = 2D〈∇Z · ∇Z|ζ〉 (2.57) is the conditional scalar dissipation rate that can be modeled [43–46]; other unclosed terms are: eQ ≡ 〈∇(ρD∇Q) + ρD∇Z · ∇∂Q ∂ζ |ζ〉, (2.58) and eY ≡ −〈ρ∂Y ′′ ∂t + ρu · ∇Y′′ −∇ · (Dρ∇Y′′)|ζ〉. (2.59) The term eQ represents contributions from molecular diffusion down the mean gradients of Q and Z; it is negligible at large Re numbers. 26 The other term eY needs to be modelled. Klimenko [42] proposed to close eY with: eY ≈ ∇ · (〈ρ|ζ〉〈u ′′ Y ′′|ζ〉P(ζ)) P (ζ)〈ρ|ζ〉 ≈ ∇ · (〈ρ|ζ〉(−Dt∇Q)P (ζ)) P (ζ)〈ρ|ζ〉 . (2.60) Similar to Eq. 2.56, a CMC equation can be derived for enthalpy. In a homogeneous field, the CMC equations reduce to: ∂Q ∂t − 〈χ|ζ〉 2 ∂2Q ∂ζ2 = 〈ω̇|ζ〉, (2.61) which resembles the flamelet equations. This shows the close relation between CMC and the flameletmodel where the same one-dimensional mixing assumption is often adopted for species diffusion. The chemical source terms are often closedwith the first order CMC hypothesis – that is, the conditional averages of the chemical reaction rate is approximated by evaluating the rate expressions with the con- ditional means of the species mass fraction vector, temperature and density: ω̇|ζ ≈ ω̇ ( T |ζ, Yi|ζ, ρ|ζ ) . (2.62) The first order CMC hypothesis is justifiable as long as the fluctua- tions around the conditional averages are small. In CMC, conditional means of reactive scalars are obtained by solving the above CMC trans- port equations. It has been successfully applied to study a wide variety of problems [43, 45, 47–50]. In cases with extinction and reignition, the first order CMC hypoth- esis might be imprecise and refinements to the closure hypothesis need to be applied. One solution is to use high order closure [51, 52]. Kim et al. [46] used the second order conditional moment closure of CMC to 27 study turbulent non-premixed hydrocarbon flames with significant lo- cal extinction and reignition. Although the predictions of the interme- diates and fuel reaction rates are improved with the second-order clo- sure method, the requirement of modelling additional unclosed terms make it less desirable. Another approach is to introduce a second con- ditioning variable. Cha et al. [53] suggested using the scalar dissipation rate as the second conditioning variable. The method predicts the ex- tinction reasonably well, but predicts the onset of reignition too early. Kronenbug [54] proposed to use a normalized enthalpy parameter to- gether with mixture fraction for double conditioning. The local extinc- tion and reignition can both be captured successfully. The traditional CMC methods are computationally expensive not only because of the unclosed terms in CMC transport equations, but also because of the need to directly solve CMC equations – each con- ditioning variable introduces a new dimension to the system. Some researchers have experimented with using coarse grids in the spatial dimensions to fix this problem. 2.3.4 Conditional Source-term Estimation The Conditional Source-term Estimation (CSE) method was first pro- posed by Bushe and Steiner [7]. It adopts the CMC first order closure hypothesis, but reduces the computational cost by eliminating the ad- ditional dimension introduced in CMC. Instead of directly solving the conditional averaged transport equations, the conditional means are obtained by inversion of an integral equation. 28 The integral equations for inversion are built based on the identity: f̃(x, t) = ∫ 1 0 f |ζ P̃ (ζ;x, t) dζ, (2.63) where f represents a reactive scalar such as density, temperature and species concentration, P̃ (ζ;x, t) is the density-weighted or Favre prob- ability density function (FPDF) defined (for LES implementation) by: P̃ (ζ;x, t) = 1 ρ(x, t) ∫ D ρ(x′, t)δ[ζ − Z(x′, t)]G(x− x′)dx′, (2.64) where ζ is the sample space of mixture fraction, δ is the delta function and Z(x′, t) is the instantaneous mixture fraction in the flow field. The scalar f̃(x, t) is obtained by solving a transport equation; the FPDF can be obtained by solving a FPDF-transport equation, or by using a pre- sumed function form such as the β-PDF. It is observed in experiments and DNS that the spatial gradient of f |ζ is small [7]. Therefore, spatial homogeneity in f |ζ might be as- sumed in many flows. For example, in an axisymmetric jet, variation in f |ζ on planes normal to the jet axis is usually very small. Gener- ally, it is assumed that some ensemble A may be selected in the flow field with statistical homogeneity in the conditional averages, i.e., the conditional means at an individual point j can be approximated by an ensemble conditional average: (f |ζ)j ≈ 〈 f |ζ 〉 A j ∈ A, (2.65) where 〈〉 denotes the average of ensemble A. Eq. 2.63 is then rewritten as: f̃ ≈ ∫ 1 0 〈 f |ζ 〉 A P̃ (ζ;xj, t) dζ, j ∈ A. (2.66) 29 Such an equationmay bewritten for each point in the ensemble. Us- ing a numerical quadrature of N points in Z space, one can discretize the integral above and construct a system of linear equations. The con- ditional means can then be obtained by inverting Eq. 2.66, which is a Fredholm equation of the first kind [55]. Inverting for the conditional means greatly reduces computational cost; however, Eq. 2.66 is usually ill-posed. To address this issue, Steiner and Bushe [56] proposed to simultaneously minimize the residual of Eq. 2.66 and the derivative of the conditional means with respect to mixture fraction. Later Grout et al. [57] proposed to use the Tikhonov regularization method that introduces a priori knowledge of the solu- tion; the inversion problem is reformulated as a minimization prob- lem. In this approach, the target is to minimize the residual of Eq. 2.66 as well as the deviation from the a priori information: min{||M〈Y |ζ〉 − Y ||+ λ||〈Y |ζ〉 − 〈Y |ζ〉0||}, (2.67) where M is the coefficient matrix of the discretized integral equations represented by Eq. 2.66, λ is a coefficient evaluated based on the char- acteristics of the inversion problem and 〈Y |ζ〉0 is the a priori informa- tion introduced. In Grout’s formulation [58], λ is evaluated with: λ = tr(M∗M)/n, (2.68) where tr(M∗M) is the trace of the n by n matrixM∗M . Using Eq. 2.68 for λ leads to a well-behaved solution 〈Y |ζ〉. Grout’s study also shows that the CSE inversion method is not sensitive to λ; increasing or de- creasing λ by a factor of 10makes little difference in the solution of con- ditional means. 〈Y |ζ〉0 is chosen to be the solution from a previous time 30 step or iteration, which is proved to stabilize the inversion process. Truncation and re-scaling is then applied to eliminate non-physical so- lutions. The Tikhonov regularization method has been applied in CSE with success [57–59]. It has also been found that optimizing the meshes in conditioning variable sample space and normalizing the inversion equations can improve the solutions. Compared to CMC, CSE is more computationally efficient. It does not require the constraining assumptions of the fast chemistry and laminar flamelet models. Therefore, it can be used to study various non-premixed flames and possibly premixed flames as long as an ap- propriate conditioning variable and an effective inversion algorithm can be selected. 2.3.5 Including Chemistry in CSE The closure method of CSE is desirable considering the reduced com- putational cost and precision of closure; however, to directly include detailed chemistry in CSE is still too expensive. The system of inver- sion equations will be too large and the species transport equations might be too stiff to be handled by available ODE solvers. To address this issue, different methods can be combined with CSE model. The CSE method was first implemented in LES to study a piloted methane/air diffusion flame with a simplified two-step chemistry [7]. The predictions of scalars such as temperature, mixture fraction and major species were consistent with the experimental data. Huang and Bushe [60] combined CSE with a chemistry reduction technique, the Trajectory Generated Low-Dimensional Manifold (TG- LDM) method using a purpose-built high-order finite-volume solver. 31 Details of the TGLDMmethod will be discussed later in this chapter. The Laminar Flamelet Decomposition(LFD) method [61] is an al- ternative to using reduced chemical kinetics, which incorporates de- tailed chemical mechanisms in combustion simulation using flamelet libraries. In the decomposition algorithm, conditional means are ap- proximated as a linear combination of basis functions: 〈f |ζ〉(t;A) ≈ Nf∑ i=1 aiΘi(ζ), (2.69) where Θi(ζ) is the basis function. The basis functions should capture the relationship between the conditional means and the conditioning variable. For non-premixed flames, one option is to use the laminar flame solutions, where a library of basis functions (a flamelet library) can be generated by solving the unsteady laminar flamelet equation: ρ ∂ψi ∂t = ρ Lei χ 2 ∂2ψi ∂Z2 + ω̇i, (2.70) with Z being the mixture fraction and χ(Z, t) being the scalar dis- sipation rate. Solutions at different flamelet times can be stored and used as the basis functions for decomposition. Alternatively, one can use the solutions of steady laminar flamelet equations or a combina- tion of both [62, 63]. The choices of scalar dissipation rate χ include using constant values or using a presumed functional form such as χ(Z, t) = χ0(t)χ(Z). According to Peters [28], the functional depen- dence on mixture fraction χ(Z) can be defined using the solution of a one-dimensional mixing layer prescribed by: χ(Z) = e−2(erfc −1(2Z)) 2 . (2.71) 32 The temporal function χ0(t) can be defined arbitrarily as: χ0(t) = c1e −t/c2 + c3, (2.72) where coefficients c1 and c2 are evaluated to approximate the behavior of mean scalar dissipation rate on the stoichiometric interface in the re- active flow under study; c3 is given different constant values to incor- porate flamelets of high dissipation rates, which correspond to tempo- rary fluctuations in scalar dissipation observed in turbulent combus- tion. Now, Eq. 2.66 can be rewritten as: f̃(x, t) ≈ ∫ 1 0 Nf∑ i=1 aiΘi(ζ) P̃ (x, t; ζ)dζ. (2.73) With Eq. 2.73, one solves for the coefficient vector ai rather than solv- ing for the conditional averages directly. In Grout’s implementation, the equations of inversion were constructed using selected scalars that were found to best differentiate between flamelet solutions at various stages of the ignition process in flamelet time, which included tem- perature and the mass fractions of CO and CH3OH [57]. Other scalars such as the reaction rates can be evaluated with Eq. 2.73, using solu- tions of the coefficient vector and the basis functions. This makes CSE with laminar flamelet decomposition more similar to flamelet meth- ods, rather than CMC methods. The CSE-LFD method was first studied in an a priori test by Bushe and Steiner [61]. Grout [58] implemented CSE-LFD in Fluent to study different flames. He used a steady flamelet library in a RANS simu- lation of a piloted steady jet flame, Sandia flame D and an unsteady flamelet library for autoigniting jets under engine-relevant conditions. 33 It was argued that the unsteady flamelet library better captured the transient characteristics of autoignition processes.Wang and Bushe [62, 63] implemented CSE-LFD in LES to study Sandia flame D. It was found that a mixed flamelet library consisting of both steady and un- steady flamelet solutions led to better predictions of the flame. Both studies showed that by optimizing the library of basis functions and the combination of scalars for inversion, the precision of closure with CSE-LFD could be improved. However, Wang and Bushe [63] found that the LES results were sensitive to the composition of the flamelet library, which suggests that one must be very careful in constructing the flamelet library if one is to use the CSE-LFD method. 2.4 Premixed Combustion For premixed combustion, the fuel and oxidizer are completely mixed before the combustion starts. Examples of premixed combustion in- clude lean-burn gas turbines and homogeneous charge spark ignition engines. Different from non-premixed combustion, premixed combus- tion phenomena have a characteristic velocity scale- the laminar burn- ing velocity- and a characteristic length scale- the flame thickness. These new parameters are used to define different regimes in premixed com- bustion. Regime diagrams have been proposed by numerous research- ers [21]. The diagrams are usually based on the comparison between the turbulent and chemical length scales and time scales- as was the case in non-premixed combustion. Regimes are defined based on dif- ferent non-dimensional parameters. For example, the Borghi diagram [64] uses a velocity ratio u′/uref 34 D a = 1 1 -1 log( ll lref ) Well Stirred Reactor, lo g ( u ′ u r e f ) Island Formation and Wrinkled Laminar Flame, Flamelet Regime T or n F la m e Fr on ts , P er tu rb ed F la m el et sThickened Flame R e T = 1 Laminar Flame Fronts Ka = 1 1 Figure 2.3: Borghi diagram for premixed combustion and a length ratio ll/lref . In a classic Borghi diagram, u ′ is the turbu- lent intensity, uref = sL is the laminar flame velocity, ll is the integral scale and lref = δL is the laminar flame thickness. As shown in Fig. 2.3, at Ka < 1, where island formation and wrinkled laminar flames are observed, the flamelet regime is defined. The laminar flamelet model (to be discussed later) works well for flames within this regime. At Da < 1, the thickened flame regime is defined. Here small scale eddies penetrate the reaction zones and the flamelet concept is inapplicable. In between is the torn flame front regime. Within this regime, the tur- bulence interacts with the flame without completely overwhelming it. Peters [30] proposed a regime diagram for premixed flames where, for scaling purposes, it is assumed that the diffusivities of all scalars are equal and the Schmidt number is unity. The length and time scales 35 of the flame are defined by: lF = D sL , (2.74) tF = D sL2 . (2.75) The turbulent Damköhler and Karlovitz number are defined as: Da = sLl υ′lF , (2.76) Ka = tF tη = l2F η2 = υ2η s2L , (2.77) where tη, η and vη are the Kolmogorov time, length and velocity scales respectively; υ′ is the turbulent intensity. A second Karlovitz number is defined by: Kaδ = l2δ η2 = δ2Ka, (2.78) where lδ is the thickness of the inner layer between the chemically in- ert preheat zone and the oxidation layer. The inner layer is where the reaction process of a premixed flame is kept alive. The thickness lδ is evaluated using lδ = δlF , where the value of δ is typically around one tenth. Different regimes can be definedwith the above characteristic scales as shown in Fig. 2.4. The ratios l/lF and υ ′/sL parameterize the regime diagram, which are related as: υ′ sL = Re ( l lF )−1 = Ka2/3 ( l lF )1/3 . (2.79) The line Re = 1 separates turbulent flames from laminar flames. The lines υ′ = sL, Ka = 1 and Kaδ = 1 divide premixed turbulent 36 10−1 100 101 102 103 104 10−1 100 101 102 103 l/lF υ ′/s L Re = 1 η = lδ Kaδ = 1 η = lF Ka = 1 broken reaction zones thin reaction zones corrugated flamelets wrinkled flamelets laminar flamelets Figure 2.4: Regimes in premixed turbulent combustion. flames intowrinkled and corrugated flamelet regimes, the thin reaction zones and the broken reaction zones regimes. In the wrinkled flamelets regime, even large eddies are not large enough to compete with the laminar flame propagation, which con- sequently dominates the combustion process. In the broken reaction zones, the inner layer is disturbed by the smaller eddies. Enhanced heat loss may cause temperature decrease and loss of radicals; as a re- sult, extinction might occur. For practical reasons, these two regimes are of less interest in industrial applications. The corrugated flamelets regime features Ka < 1 and υ′ ≥ sL ≥ υη. In this regime, small eddies interact with the advancing laminar flame. The thin reaction zones regime is characterized byRe > 1,Kaη < 1 and Ka > 1. The upper limit of this regime is defined by Kaδ = 1 and in 37 most cases Ka = 100 can be used as an approximation. In this regime, small eddies may increase scalar mixing, but cannot penetrate the in- ner layer. The corrugated flamelet and thin reaction zones regimes are separated by the line Ka = 1. Under such circumstances, the thick- ness of the flame is equal to the Kolmogorov length scale; tF = tη and υη = sL can be derived from Eq. 2.77. Many premixed combustion models are developed based on the characteristics of different regimes. 2.4.1 Eddy-Break-Up The Eddy-Break-Up (EBU)model is based on a phenomenological anal- ysis and is primarily used in premixed combustion with high Re and Da numbers. It is argued [65] that in cases where mixing determines the rate of chemical reactions, the cascade process from the integral to molecular scales of turbulent mixing also controls the reaction process. In case of excessive oxidizer, the mean reaction rate of fuel can be eval- uated using the turbulent mixing time τt: ω̇F = −ρCEBU (Ỹ ′′ F 2 )1/2 τt , (2.80) where Ỹ ′′ F 2 is the variance of the fuel mass fraction; CEBU is the Eddy- Break-Up constant (on the order of unity); τt is the turbulent mixing time estimated from the turbulent kinetic energy k and its dissipation rate ε by using τt = k/ε. A progress variable c for the chemical reaction is defined such that c = 0 in unburnt gas and c = 1 in fully burnt mixture. Examples of c include a normalized temperature or a normalized mass fraction: c = T − Tu Tb − Tu , (2.81) 38 c = YP YP,b , (2.82) where T , Tu and Tb are the temperatures of the local mixture, unburnt gas and fully burnt mixture, YP and YP,b are the mass fractions of prod- ucts in the local and fully burnt mixture. The reaction rate may be cal- culated using a progress variable c as: ω̇ = −ρCEBU √ c̃′′2 τt . (2.83) Eq. 2.83 leads to inconsistencies at c̃ = 0 and c̃ = 1 when the flame front is infinitely thin (i.e. c = 0 or c = 1). Under such conditions, c2 = c is valid and c̃′′2 is given by: ρc̃′′2 = ρ(c− c̃)2 = ρ(c̃2 − c̃2) = ρc̃(1− c̃). (2.84) Substituting Eq. 2.84 in the reaction rate equation, one can find dω̇/dc̃ is infinite at c̃ = 0 and c̃ = 1. A corrected version was proposed by Borghi [8]: ω̇ = CEBUρ c̃(1− c̃) τt . (2.85) The EBU model was modified by Magnussen and Hjertager [66], leading to the Eddy DissipationModel (EDM). In EDM, themeanmass fraction of the deficient reactant is used and the minimum of the reac- tion rates for fuel, oxidizer and products is used to evaluate the mean chemical source terms. The EBU model replaces the chemical time scale of an assumed one-step chemical reaction with the turbulent time scale τt and calcu- lates the reaction rates simply as a function of some known quantities. Because of its simplicity, it has been widely used in commercial CFD packages. However, because the chemical kinetics are eliminated from 39 0 1 α(x, t) β(x, t) Figure 2.5: Presumed probability density function in premixed turbu- lent combustion, prescribed by the Bray-Moss-Libby model. the modelled reaction rates, the model represents the fast chemistry limit only. The model constant CEBU sometimes needs to be tuned ac- cording to specific problems. It is also observed that the model tends to over-estimate the reaction rates. 2.4.2 The Bray-Moss-Libby Model The Bray-Moss-Libby (BML) [67]model implements the classical flamelet concept for premixed turbulent combustion. Named after its authors, themodel has been studied and improved by numerous researchers [68]. The BML model introduces a progress variable c to use the flamelet concept. The presumed probability density function of the progress variable c at a given time and location (x, t) is modeled with: P (c∗, x, t) = α(x, t)δ(c∗) + β(x, t)δ(1− c∗) + γ(x, t)f(c∗, x, t). (2.86) The presumed PDF prescribed by Eq. 2.86 is presented in Fig. 2.5. It describes the mixture at (x, t) as a sum of unburnt, fully burnt and burning gases, whose contributions are represented by the three terms 40 on the r.h.s. of Eq. 2.86 respectively. The coefficients α, β and γ denote the probabilities of having fresh, completely burnt and reacting mix- tures at (x, t). The Dirac delta functions δ(c∗) and δ(1− c∗) correspond to fresh mixture at c = 0 and complete combustion at c = 1 respec- tively. The functional form f(c∗, x, t) for the probability of burning gas can be defined such that: ∫ 1 0 f(c∗, x, t)dc∗ = 1 (2.87) with f(0) = f(1) = 0. The coefficients are subject to the normalization of the PDF:∫ 1 0 P (c∗, x, t)dc∗ = 1, (2.88) from which one can derive: α+ β + γ = 1. (2.89) The model is further simplified in cases where Re >> Da >> 1 . Under such circumstances, the premixed turbulent flame is extremely thin and the combustion is mainly controlled by the turbulent trans- port. Consequently, γ << 1, γ << α and γ << β are valid. The PDF is reduced to a quasi-bimodal form: P (c∗, x, t) = α(x, t)δ(c∗) + β(x, t)δ(1− c∗). (2.90) The coefficients α and β can be easily evaluated using the normaliza- tion of PDF prescribed by Eq. 2.88 and averages of the progress vari- able c̃: ρc̃ = ∫ 1 0 ρc∗P (c∗, x, t)dc∗, (2.91) 41 from which one can derive: β = ρc̃ ρb , (2.92) α = 1− β. (2.93) Here ρb is the density of fully burnt mixture. The reaction heat release factor τ = Tb/Tu − 1 is often introduced in the BML analysis. At constant pressure P , the following equations can be derived by assuming perfect gas behavior: τ = Tb Tu − 1 = ρu ρb − 1, (2.94) ρu = (1 + τ)ρb = ρ(1 + τ c̃). (2.95) Substituting these in Eq. 2.93 and Eq. 2.92 and one obtains: α = 1− c̃ 1 + τ c̃ , (2.96) β = (1 + τ)c̃ 1 + τ c̃ . (2.97) The heat release factor τ is determined by the specific chemical mec- hanism and is therefore fixed. The PDF P (c) depends solely on c̃ and τ . The governing equation of the progress variable is: ∂ρc ∂t +▽ · (ρuc) = ▽ · (ρD▽ c) + ω̇ (2.98) Under the assumptions of the BMLmodel, where the progress variable c is equal to zero or unity, an equation can be derived for reaction rates: 2ρD▽ c · ▽c = 2cω̇ − ω̇ (2.99) After averaging, Eq. 2.99 can be rewritten as: 2ρD▽ c · ▽c = (2cm − 1)ω̇, (2.100) 42 where cm = R 1 0 cω̇f(c)dcR 1 0 ω̇f(c)dc is a progress variable characterizing the chemical kinetics. The mean reaction rate is: ω̇ = 2 ρχ 2cm − 1 , (2.101) where ρχ is the scalar dissipation rate of the progress variable c: ρχ = ρχ̃ = ρD▽ c · ▽c = ρD ∂c ∂xi ∂c ∂xi . (2.102) The scalar dissipation rate χ̃ can be evaluated by solving a transport equation [69] or by modelling: ρχ ≈ ρc ′′2 τt , (2.103) where τt is the turbulent time scale. Assuming the flame front is in- finitely thin (i.e. c = 0 or c = 1), c2 = c is valid and c̃′′2 is given as: ρc̃′′2 = ρ(c− c̃)2 = ρ(c̃2 − c̃2) = ρc̃(1− c̃). (2.104) Eq. 2.101 is rewritten as: ω̇ = 2 2cm − 1 ρc̃(1− c̃) τt , (2.105) which is consistent with the EBU model represented by Eq. 2.85. The link between the BML and EBU models can be also established using flamelet analysis [8]. Using the progress variable, the average of any quantity may be evaluated by conditional averaging. For example, in the bimodal form of BML, the Favre average of a scalar Q can be calculated as: ρQ̃ = ρQ = ∫ 1 0 ρQP (c)dc = αρuQu + βρbQb, (2.106) whereQu andQb are the conditional means in the fresh gases and fully burnt mixture respectively. 43 2.4.3 Models Based on Surface Area Estimation Another main category of turbulent premixed combustion models is based on the concept of flame surface density Σ. By definition, Σ is the flame surface area per unit volume. The mean reaction rate of species i is modeled with: ω̇i = Ω̇iΣ, (2.107) where Ω̇i is the local reaction rate per unit of flame area integrated in the direction normal to the flame surface. It is related to the properties of the local flame structure and can be estimated using a prototype laminar flame in a pre-processing fashion. This method decouples the chemical kinetics described by Ω̇i and the turbulence/flame interaction described by Σ. Bray et al. suggested that the mean chemical source term must be proportional to the flamelet crossing frequency and proposed a mo- del [68]: ωc = ρus 0 LI0Σ, (2.108) where s0L is the laminar burning velocity of the unstretched flame and I0 is the stretch factor. The flame surface density Σ can be calculated using an algebraic expression such as [68]: Σ ≈ c(1− c) L̂y , (2.109) where L̂y is the crossing length scale; alternatively, it can be obtained by solving a transport equation. For example, an equation was pro- posed by Trouvé and Poinsot [70]: ∂Σ ∂t +▽ · (υ̃Σ) = ▽ · (Dt ▽ Σ) + C1 ε k Σ− C2sL Σ 2 1− c. (2.110) 44 The terms on the left hand side of the equation represent the local rate of change and convection, and the terms on the right hand side rep- resent turbulent diffusion, production by flame stretch and flame sur- face annihilation respectively. Other forms of the balance equation for Σ have been developed in various premixed combustion models [8] where the source and sink terms are modelled differently. 2.4.4 Models Based on the G-equation The premixed combustion models discussed previously involve the use of a progress variable. An alternative approach is the G-equation method, where a non-reacting scalar G is introduced. An iso-surface at a critical value G0 divides the flow field into regions of burnt gases and unburnt gases. The value ofG0 is fixed for a particular combustion event. A governing equation for G was proposed by Williams [71] as: ∂G ∂t + v · ▽G = sL| ▽G| (2.111) This G-equation is applicable to thin flames that propagate at a well- defined burning velocity and is therefore suited for describing pre- mixed turbulent combustion in the corrugated flamelets regime. Pe- ters [21] formulated a G-equation in the thin reaction zones regime: ∂G ∂t + v · ▽G = [▽ · (ρD▽ T ) + ωT ρ| ▽ T | ] | ▽G|. (2.112) A common level set equation was formulated for both regimes: ∂G ∂t∗ + v∗ · ▽∗G = sL,s νη | ▽∗ G| − D ν κ∗| ▽∗ G|, (2.113) where the independent variables and curvature are normalized by the Kolmogorov length, time and velocity scales: t∗ = t/tη, x ∗ = s/η, v∗ = v/vη, κ ∗ = ηκ, ▽∗ = η▽ . (2.114) 45 Analogous to the role of the conserved scalar Z in non-premixed flames, G and/or its variance are often used in premixed combustion models to evaluate various properties of premixed flames. For exam- ple [72], a flamelet library may be generated where the species, tem- perature and density of the prototype flames are stored as functions of distance from the G = 0 surface and strain rate; it is then incorporated in the combustion simulation using a presumed probability density function. 2.5 Reduction of Detailed Reaction Mechanism A detailed chemical reaction system generally includes many species and elementary reactions. To simulate turbulent reactive flows at af- fordable computational cost, it is usually necessary to reduce the full mechanisms. The most commonly used methods include the Quasi- steady State Approximation [73] (QSSA) and the partial equilibrium methods [74]. In QSSA, analytical approximations of reaction rates are derived by assuming equal consumption and formation rates for cer- tain species; in partial equilibrium methods, they are derived by as- suming small net reaction rates of fast reactions. These methods have been applied to study premixed and non-premixed flames with suc- cess [75–77]. Recently, new techniques have been developed and applied in tur- bulent combustion simulation, which typically involve computation with full chemistry in a pre-processing step and tabulation. Examples include the Intrinsic Low-Dimensional Manifold (ILDM) [78] method, Invariant Constrained-equilibrium Edge Pre-Image Curve (ICE-PIC) 46 method [79] and methods based on flamelet concepts [80, 81]. 2.5.1 Intrinsic Low-Dimensional Manifold The Intrinsic Low-DimensionalManifold (ILDM)methodwas first pro- posed by Maas and Pope [78]. The method is based on the fact that a wide range of time scales are involved in a chemical reaction system (from 10−10 s for radical reactions to more than 1 s for NO formation). Fast processes attenuate quickly and reaction systems converge onto manifolds of fewer dimensions on longer time-scales such as those rel- evant to mixing in a turbulent non-premixed flame. The ILDMmethod is used to find such a manifold, which is parameterized by progress variables and stored as a look-up table in combustion simulation. The first step is to define the low-dimensional manifold bymeans of a time scale analysis. The detailed chemical kinetic system is described with: ∂ψ ∂t = F(ψ) + Ξ(ψ,∇ψ,∆ψ), (2.115) where ψ = (h, p, ψ1, ψ2, ..., ψns) T is the scalar vector consisting of en- thalpy, pressure and all ns species in the detailed chemistry, ψk = Yk Wk is the specific mole number of species k ([mol/kg]), F(Ψ) is the chemical source-term and Ξ denotes physical processes such as convection, dif- fusion, etc. For a simple homogeneous, adiabatic and isobaric system, Eq. 2.115 reduces to: ∂ψ ∂t = F(ψ). (2.116) The time scales of the system can be separated through a Schur decom- position analysis: J = VΛṼ, (2.117) 47 where J is the Jacobian matrix of the chemical source term F(ψ) with Jij = ∂Fi ∂ψj , V is the n × n right Schur vector matrix of J, Λ is an upper- triangular matrix with the eigenvalues of J as its diagonal entries. The eigenvalues λk inΛ are sorted in descending order such that the eigen- value representing the fastest process is at the bottom; V and Ṽ are sorted accordingly. The characteristic chemical time scales of different processes are evaluated with: τk = λ −1 k . (2.118) Different time scales are grouped into nf fast time scales and nc slow time scales, where the fast processes quickly attenuate to quasi-steady states and the reaction system converges onto an nc-dimensional man- ifold defined by: ṼfF(ψ) = 0, (2.119) where Ṽf consists of the last n− nc rows taken from Ṽ, corresponding to the fast processes. Such a subspace can be found using Eq.2.119, to- gether with additional equations representing system constraints and those of selected progress variables. The final form of the ILDM equa- tion is: G(ψ, τ) = ṼfF (ψ) P (ψ, τ) = 0, (2.120) where P contains the parametric equations such as the adiabaticity, constant pressure and element conservation. Once the low dimensional manifold is defined, the detailed chemi- cal kinetic system can be parametrized by a small number of variables θ = (θ1, θ2, ..., θN ) (N < n). The n-dimensional system can be projected 48 onto the new manifold with the governing equations: ∂θ ∂t = S(θ) + Γ(ψ(θ),∇ψ(θ),∆ψ(θ)). (2.121) To apply ILDM in combustion simulations, the low-dimensional manifolds are usually calculated beforehand and stored as look-up ta- bles. These tables can be retrieved during the combustion simulation and the themokinetic properties are evaluated using the progress vari- ables solved in the flow computation. Since the calculation of ILDM is only based on the analysis of the reaction system, they can be used in various flame configurations as long as the element composition, pres- sures and enthalpies are a subset of the conditions defined in ILDM computation. Another attractive feature of ILDM is that the mathe- matical formulation of ILDM provides an automatic error control of the simplification scheme. The ILDM method has been used to study laminar and turbulent flames of various configurations [82, 83]. However, the ILDM method is based on the assumption that the fast chemical time scales can be well separated from the slow ones, which is not necessarily true for any chemical system. To solve this problem, Nafe andMaas [84] introduced the Slow Invariant Manifolds (SIM) methods, where the original ILDM computation was used as an initial guess and a post-processing step was applied such that only movements perpendicular to the manifold were allowed. The resultant manifold was found to be more accurate than the original ILDM. Another issue is that at relatively low temper- atures, the diffusion and chemical time scales are of the same order. Therefore, the diffusion processes need to be included in the computa- tion of manifolds. 49 Some researchers have combined ILDM with different ways of ma- nifold generation to solve this problem. For example, a related ap- proach called Flame Prolongation of ILDM (FPI) was studied by Gic- quel et al. [85]. In this method, the low temperature regime in ILDM was constructed by using the solutions of laminar premixed flames. Oi- jen et al. [81] introduced the Flamelet-Generated Manifold (FGM) me- thod based on the similarity of flame structure between a propagating three-dimensional flame and a one-dimensional flame. A combination of ILDM and FGM methods called Phase-space ILDM (PS-ILDM) was studied by Bongers et al. [86]. In the PS-ILDMmethod, the ILDM strat- egy was used to decouple the fast and slow processes in the solution of flamelet equations. Besides the efforts to improve the quality of manifolds, the tabula- tion of manifolds has been studied to reduce the requirements of stor- age. For instance, a method called in situ adaptive tabulation (ISAT) was combined with ILDM methods by Pope [87], which built the ma- nifold as needed during simulation. 2.5.2 Trajectory Generated Low-Dimensional Manifold An alternative technique of generating low dimensional manifolds is the Trajectory Generated Low-Dimensional Manifold (TGLDM) met- hod, originally proposed by Pope and Maas [88]. The TGLDMmethod involves solving for physically realizable initial states of the system and integrating Eq. 2.115 to evolve the reaction trajectories. Each tra- jectory consists of solutions in the composition space starting from a given initial state of the system and evolving toward the chemical equi- librium state. Theoretically, if a viable low dimension manifold does 50 exist for a chemical system, the trajectories with initial states on the manifold will stay on the manifold. In this sense, such trajectories are equivalent to a SIM. An important advantage over ILDM is that the al- gorithm used in TGLDM guarantees a converged solution. The global optimization problem in ILDM is reduced to a local optimization prob- lem of defining initial states for trajectories. As in ILDM, TGLDMman- ifolds can be calculated and tabulated in pre-processing and retrieved later during a combustion simulation. To select the initial states, Pope and Maas [88] adopted the extreme- value-of-major-species method, in which linear equations were solved based on the conservation of elements. For example, major species in the reaction systems may be selected and the mass conservation of el- ements can be written for these species: ns∑ i=1 Yc;j,iYi = Ye;j j = 1, ..., ne, (2.122) where ns is the number of selected species, ne is the number of ele- ments, Yc;j,i is the mass fraction of element j in species i, Yi is the mass fraction of species i and Ye;j is the mass fraction of element j in the fresh mixture. By solving such equations, initial states can be found with extreme values of mass fractions at the boundary of the composi- tion space. Huang and Bushe [60] used a slightly different approach. To con- struct an np dimensional manifold, constraints on progress variables are applied in addition to the element conservation equations: Pk,iYi = Yp;k; i = 1...ns, k = 1...np. (2.123) where np is the number of progress variables that parameterize the low 51 Figure 2.6: Comparing one dimensional manifold generated by ILDM and TGLDM. dimensional manifold, Pk,i is the coefficient for the parametric equa- tions, Yp;k is the assigned value of progress variable k. A natural choice of the constraints is to use the equilibrium state in the subspace, which can be calculated with the method of Lagrange multipliers [89, 90] or Gibbs function continuation [91]. Huang [83] used this method to gen- erate a one dimensional manifold which was compared to an ILDM and good agreement was observed. His analysis also showed that the TGLDM createdwith constrained equilibrium better approximated de- tailed chemistry than ILDM methods under moderate perturbation. The TGLDM generated was applied in three different methane/air re- action systems with success including an unstrained, premixed lami- nar flame, a non-premixed laminar flame and a perfectly stirred reac- tor. 52 Later Huang and Bushe [60] combined TGLDM with CSE to study a methane jet flame where fuel was injected into a shock-tube under engine-relevant conditions. In this study, a two dimensional manifold parameterized by YCO2 and YH2O was constructed. The manifolds were stored using Delaunay triangulation [92] and retrieved by surface in- terpolation over the two-dimensional unstructured grid [92]. Five ma- jor species (O2, CH4, CO, CO2 and H2O) were used to construct the equations for the initial states, including element conservations of C, H andO, and constrained equilibrium equations of the two progress vari- ables. The predictions of autoignition characteristics of the jets were found to be promising. In another work on autoigniting jets, Wang et al. [62, 93] discussed that at low pre-combustion temperatures, the ignition delay time be- comes exceedingly long such that it is extremely difficult for the ODE solver to obtain a solution for the trajectory. In previous studies, this problem was circumvented by mixing the unburned mixture with a small amount of mixture at equilibrium. The initial state point was therefore shifted slightly into a state of a higher equivalent temperature and the stiffness of the systemwas reduced. The resulting trajectory ch- emistry does not describe autoignition correctly. To solve this problem, the Stochastic Particle Model (SPM) [94–96] was used to generate the autoignition trajectory at low temperatures. In Wang’s study, the SPM was used to calculate several realizations of trajectories, which were filtered, averaged, and combined with TGLDM trajectories generated using a traditional continuum approach. The SPM-extended TGLDM was then used (in the same high-order code as Huang and Bushe [60]) 53 to study a methane jet and the fluctuations in ignition delays due to randomness of chemical reaction paths were analyzed. Wang and Bushe [59] also used CSE-TGLDM to study a piloted, non-premixed turbulent jet flame, Sandia flame D. The same approach as discussed in Huang and Bushe’s work [60] was adopted to close the chemical source terms in LES. The predictions of major species, temperature field as well as the conditional averages agreed well with the experiment. 2.6 Summary Turbulent combustion simulation is challenging because of the cou- pling of mechanisms such as fluid dynamics, chemical reactions, radi- ation and phase changes, etc. Depending on the characteristics of the combustion phenomena, different methods can be used to approach this problem. Among these methods, Conditional Source-term Estima- tion provides a closure method with a good balance between computa- tional cost and precision. In previous research, CSE-TGLDM and CSE- LFD have been successfully implemented in different studies. More work needs to be done to explore the applications of CSE: • Although CSE-TGLDM and CSE-LFD have been implemented in individual studies, these two methods have never been com- pared head-to-head. Such a study would help to clarify the un- derlying differences between these two methods and compare the performance of them on a fair basis. • So far, the CSE applications have never accounted for the effects of radiation heat loss, which is important especially in the pre- 54 diction of NO formation and might not be neglected in luminous flames. • In previous work, the application of CSE is limited to modelling non-premixed flames. The possibility of using CSE in premixed flames is worth exploring. These research problems are studied in this thesis work and will be discussed in the following chapters. 55 Chapter 3 Ignition Delay of Methane Jets with Conditional Source-term Estimation 3.1 Introduction Conditional source-term estimation (CSE) has been used to study the autoignition process in non-premixed jet flames of methane and me- thane mixed with additives, where conditional means of temperature and species are used to close chemical source-terms. CSE is combined with a Trajectory Generated Low-DimensionalManifold (TGLDM)me- thod and then with a Laminar Flamelet Decomposition (LFD) method to reduce the complexity of detailed chemistry. The predictions of ig- nition delay time are compared with measurements of shock-tube ex- periments. 3.2 Model Formulation 3.2.1 Numerical Scheme Reynolds-Averaged Navier-Stokes (RANS) simulation of autoigniting jets has been performed in an axisymmetric computational domain us- ing the standard k − ǫ [14] turbulence model. Besides continuity and momentum equations, Reynolds averaged transport equations [8] of 56 energy, selected species, the mean and variance of mixture fraction are solved in a cylindrical coordinate. In the current study, body force is neglected in themomentum equa- tions and radiation heat loss is neglected in the energy equation. It is assumed that the Schmidt number Sc is 0.7; the Lewis number Le is unity and the diffusivities of individual gaseous species are not signif- icantly different. The chemical source terms in the Reynolds averaged energy and species equations are closed with CSE. For details of the CSE-LFD and CSE-TGLDM methods, please refer to Chapter 2, Sec- tion 2.3.5 and 2.5.2. 3.3 Results and Discussion To validate and compare the CSE-TGLDM and CSE-LFD methods de- scribed above, these methods have been implemented in OpenFOAM [97] to simulate a variety of shock-tube experiments where the shock is used to compress air and fuel is injected along the centerline of the shock-tube with an injector mounted on the end-plate. The resulting jet ignites and the ignition delay time is measured using a high-speed camera, imaging the entire jet in the visual spectrum. Ignition is de- tected by identifying the first visible flame kernel that subsequently grows into a full jet flame. The experimental techniques used are de- scribed by Sullivan et al. [98]. 3.3.1 Simulation of a Non-reactive Jet A cold jet was simulated to validate solvers developed using Open- FOAM. The jet under discussion was studied experimentally and nu- merically by Ouellette [99] . In his experiments, cold methane jets were 57 injected at high Reynolds number through a nozzle and mixed with cold air in the shock-tube. Because of the high pressure ratio upstream of the chamber and the high exit pressure at the nozzle, underexpanded jets were observed. The penetration length was measured and aver- aged over three independent tests. A scalingmodel was proposed byOuellette [99] based on the vortex- quasi-steady-state assumption and the mass/momentum entrainment measurement from Ricou and Spalding [100]. According to the scaling model, the penetration of a self-similar transient jet can be modeled as:( Zt dn √ ρi ρ0 ) = Γ √ π 4 ( U0t dn √ ρi ρ0 )0.5 , (3.1) whereZt is the penetration length of the jet, dn is the nozzle diameter, ρi and ρ0 are the densities of the gas at the nozzle exit and in the ambient air respectively, U0 is the jet exit velocity, and Γ is an empirical penetra- tion number with a value around 3. An equivalent diameter defined by deq = dn √ ρi ρ0 is introduced to non-dimensionalize the penetration length Zt in post-processing. To test the CFD solvers developed with OpenFOAM, the predicted penetration length of the jet is compared to experimental data as well as predictions of Fluent [101] and the scaling models. The cold jet un- der study is injected from a 0.5 mm nozzle at the injection pressure 2.285 MPa and temperature 300 K into the shock-tube filled with cold air at the pressure of 1.494MPa. Simulations of the non-reactive jet have been performed in an axi- symmetric computational domain over a structured mesh. The mesh is refined near the exit of the nozzle to resolve the large local gradient. A test run with a doubly refined mesh shows no significant changes in 58 the solution. The inflow boundary conditions are defined by assuming a polytropic expansion of the fuel jet at the injector tip. The penetration length is evaluated using the method proposed by Ouellette [102]. The cold jet is also simulated using the commercial software Fluent. For comparison, the configurations of the OpenFOAM and Fluent compu- tations are set as close to each other as possible. Fig. 3.1 shows the non-dimensionalized penetration length of the methane jet as a function of time. The predictions of the OpenFOAM and Fluent simulations are compared with the experimental data and the prediction of the scaling model. The experimental data contain er- ror because of the injection delay between issuing the command to in- ject and the actual onset of injection. Therefore, it is more important to compare the simulation results to the scaling model prediction, which is calculated in two steps. Firstly, the scaling model is rewritten to in- clude an injection delay time and Γ is estimated by curve fitting the new function form to the experimental data; secondly, Eq. 3.1 is used to calculate the penetration length as a function of time (without the injection delay). The figure shows that the simulation results are generally consis- tent with the measurements. The profiles of the OpenFOAM and Flu- ent predictions are close to each other, while the OpenFOAM and the scaling model calculation are barely distinguishable. The comparison shows that for the two dimensional RANS simulation of this cold jet, the performance of the solver developed with OpenFOAM is at least as good as that of Fluent. On the other hand, the OpenFOAM solver completes the computation of an injection duration of around 3 ms in 59 0 0.5 1 1.5 2 2.5 3 x 10−3 0 10 20 30 40 50 60 70 80 90 t(s) Z t / d n Figure 3.1: Penetration length of a cold jet. : experimental data; — -: simulation with OpenFOAM, −−: simulation with Fluent; · · · · ·: prediction of the scaling model. 60 Table 3.1: Autoigniting jets with various fuel composition Parametera Pi (bar) ti (ms) Fuel T0 (K) Jet A 120 1.0 100%CH4 1200-1400 Jet B 120 1.0 90%CH4 + 10%C2H6 1200-1400 Jet C 120 1.0 80%CH4 + 20%N2 1200-1400 Jet D 60-120 1.0 100%CH4 1300 Jet E 120 1.0-2.5 100%CH4 1300 a Pi: injection pressure; ti: injection duration; T0: pre-combustion temperature; fuel composition is given by volume. about half the time that Fluent takes. The simulations of autoigniting jets hereafter are performed using OpenFOAM. 3.3.2 Simulations with CSE-TGLDM and CSE-LFD The autoignition of jets of pure methane and of methane mixed with additives has been simulated with OpenFOAM. For all the cases under study, the nozzle diameter is 0.28 mm and cold fuel jets at around 300 K are injected into the shock-tube with a pre-combustion pressure of 30 bar. Other important characteristics of the investigated jets are listed in Table 3.1 The experimental results are provided by Wu et al. [103, 104]. To compare the CSE-TGLDM and CSE-LFD methods, the two dif- ferent CSE approaches are implemented using the same code, numer- ical scheme, grid as well as the same initial and boundary conditions. A detailed chemical mechanism including 71 species and 379 reac- tions [105] is used in generating the TGLDMand LFDflamelet libraries. This mechanism is tuned to premixed shock-tube experiments. In the simulation with the CSE-TGLDM model, the mass fractions 61 of YCO2 and YH2O are selected as progress variables because of their rel- atively long formation time scales. As an example, Fig. 3.2 shows the TGLDMmanifold and its Delaunay triangulation for amethane/ethane jet in YCO2 and YH2O plane, generated at the stoichiometric mixture frac- tion, with a pre-combustion temperature of 1200 K and at a constant pressure of 30 bar. In the simulation with the CSE-LFD model, unsteady flamelet so- lutions are used as basis functions for decomposition. For example, Fig. 3.3 shows flamelet solutions of a methane/ethane jet extracted at various flamelet times (corresponding to different flamelet numbers), with a pre-combustion temperature of 1200 K and at a constant pres- sure of 30 bar. The flamelets shown are all taken prior to ignition in the laminar flamelet time. It is observed that before ignition, the temper- ature change is not sufficient to differentiate between flamelets; while the combination of T , YCO and YCH3OH might be used to detect differ- ent stages of the autoigniting process. Therefore, temperature and the mass fractions of CO and CH3OH are selected as scalars to construct the inversion equations of CSE in this study. It has been found in earlier work [57] that the simulated ignition de- lay time is largely insensitive to the definition used: the state changes so abruptlywhen ignition occurs thatmany different variables all change dramatically and almost simultaneously, such that it is possible to de- tect ignition by tracking any one of several different variables. For con- venience, the ignition delay is defined based on a rise in the conditional mean temperature. Fig. 3.4 shows the maximum increase in the condi- tional means of temperature in the reaction field plotted against time 62 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 YCO2 Y H 2 O (a) TGLDMmanifold 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 YCO2 Y H 2 O (b) Delaunay triangulation Figure 3.2: TGLDM for methane/ethane jet at stoichiometric mixture fraction, T0 = 1200 K, P0 = 30 bar. 63 0 0.2 0.4 0.6 0.8 1 200 400 600 800 1000 1200 1400 Z T (K ) 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 x 10 11 Z ω T (e r g / g ) 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 −5 Z Y C O 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 −6 Z Y C H 3 O H Figure 3.3: LFD library ofmethane/ethane jet, T0 = 1200K, P0 = 30 bar. —-: flamelet 5; · · · · ·: flamelet 10; − · − line: flamelet 20; −−: flamelet 25. 64 0 0.2 0.4 0.6 0.8 1 1.2 1.4 x 10−3 0 5 10 15 20 25 30 35 40 t(s) M ax . In cr ea se in 〈T |ζ〉 Figure 3.4: Maximum increase of 〈T |ζ〉 predicted with CSE-TGLDM for several different initial temperatures with Jet B. —-: T0 = 1400 K; −−: T0 = 1350 K; · · · · ·: T0 = 1300 K; − · −: T0 = 1250 K; ◦: T0 = 1200. at different pre-combustion temperatures for the methane/ethane jet calculated by CSE-TGLDM model. It can be seen that this maximum value is almost zero early in the simulation and later increases dramat- ically. Similar observations are found with CSE-LFD calculation and in the simulations of pure methane jets. In this study, the ignition delay time is calculated as the time when the dramatic change in conditional mean temperature begins. The autoignition of pure methane and methane/ethane jets is sim- ulated using both CSE-TGLDM and CSE-LFD methods. Predictions of ignition delay are compared with experiments, as shown in Fig. 3.5. Predictions with both methods appear to agree with the experimental results despite (if not because of) the scatter in the experimental data. 65 The scatter observed in the measurements can be attributed to fluctua- tions in turbulent mixing as well as to different realizations of chemical reaction paths - even in homogeneous, quiescent ignition, it has been found that the ignition delay time is a random variable [93, 96]. Al- though the present models capture the trend of ignition delay, further work is needed to account for these fluctuations. The ignition delay predicted with CSE-TGLDM is somewhat better than that predicted with with CSE-LFD – the predictions of CSE-LFD being universally shorter than those with CSE-TGLDM. This discrep- ancy is due to the different techniques of incorporating detailed ch- emistry in TGLDM and LFD. For comparison, the TGLDM and lam- inar flamelet library generated at the stoichiometric mixture fraction for the 1200 K methane/ethane jet are plotted in Fig. 3.6. For a given temperature, the rate of change of temperature (Ṫ ) in the unsteady flamelet library is generally higher than that of the TGLDM, especially in the temperature range of 1200–1500 K within which autoignition occurs. The TGLDM also includes more data points with lower change rates that are missing from the flamelet library. This difference between the predictions obtained with the two libraries appears to be a conse- quence of the modelling technique. For the unsteady flamelet library, the scalar dissipation rate is modelled according to the approximate behavior of the mean scalar dissipation along the stoichiometric inter- face in the jets. That the modelled scalar dissipation rate cannot exactly reproduce all possible realizations of the scalar dissipation rate in the real flame could be a source of the discrepancy. Another issue is that the solutions in the flamelet library are ex- 66 0.7 0.75 0.8 0.85 0.5 1.0 1.5 1000/T (1/K ) τ (m s) (a) Jet A 0.7 0.75 0.8 0.85 0.5 1.0 1.5 1000/T (1/K ) τ (m s ) (b) Jet B Figure 3.5: Comparison of CSE-TGLDM and CSE-LFD predictions to experimental results. +: experimental data; : calculation with CSE- TGLDM; △: calculation with CSE-LFD. 67 1000 1500 2000 2500 3000 100 102 104 106 108 1010 T (k) Ṫ (K / s ) Figure 3.6: Comparison of TGLDM and laminar flamelet library. +: laminar flamelet library; •: TGLDM library. clusively igniting flamelets: the flamelet method is not able to account for the transition from ignition to burning. The TGLDM library cov- ers a much wider range of conditions and is able to transition from igniting to burning conditions seamlessly. This, combined with the (ar- guably) better predictions of ignition delay, leads to the conclusion that the CSE-TGLDMmethod is the better of the two approaches available. 3.3.3 Effects of Different Fuel Composition The CSE-TGLDM method is adopted to simulate jets of different fuel composition. To study the effects of dilution on the ignition delay, met- hane was diluted with nitrogen and autoignition was studied in shock- tube experiments. The characteristics of the diluted jet are listed as Jet C in Table 3.1. The experimental results are discussed by McTaggart- 68 0.7 0.75 0.8 0.85 0.5 1.0 1.5 1000/T (1/K ) τ (m s ) Figure 3.7: Ignition Delay τ vs T0. +: experimental data; : predictions of CSE-TGLDM. Cowan et al. [104]. This series of experiments is simulated and the pre- dictions are compared to the experiments in Fig. 3.7. The predictions of ignition delay time capture the trend observed in the experimental data. To analyze the effects of different additives on the autoignition pro- cess, the measured and predicted ignition delay of three series of jets (from jet A to jet C) are compared. The effects of additives are not ob- vious in the experimental data because of the scatter in the measure- ments. There is as yet insufficient data from the methane/ethane ex- periments to be able to draw statistically significant conclusions about the effects of ethane addition on the ignition these jets. However, as is shown in Fig. 3.8, the CSE-TGLDM model predicts a modest reduc- tion in ignition delay time of the methane/ethane jet – in fact, the re- 69 0.7 0.75 0.8 0.85 0.5 1.0 1.5 1000/T (1/K ) τ (m s ) Figure 3.8: Ignition Delay τ vs T0, predictions of CSE-TGLDM. ◦: CH4; : 90%CH4 + 10%C2H6; △: 80%CH4 + 20%N2. duction is considerably smaller than the scatter in either the methane or methane/ethane data. At high pre-combustion temperatures, even this modest reduction becomes negligible. Further work is necessary to predict pollutant formation so as to study the trade-off between the re- duced ignition delay and increased emissions that result from adding ethane to the fuel. Similarly, it is hard to draw conclusions on the effect of dilution from the experiment because of the scatter in the measurements and the amount of available experimental data points. The comparison of ignition delay predictions in Fig. 3.8 shows that the addition of nitro- gen delays methane autoignition. The effect is most significant at the low pre-combustion temperature of 1200 K and becomes negligible at higher temperatures. Further work is necessary to predict the change 70 in NOx emission caused by the dilution. Mixing methane with additives can change the chemical kinetics of the flame. Not only does it change the stoichiometry of the flame, but it also affects the initial combustion rate by changing the energy density of the fuel. These effects are reflected in the TGLDM libraries. For example, in Fig. 3.9, a comparison of conditional mean temper- atures at the stoichiometric mixture fraction is presented. Given the same conditional means of progress variables, temperatures interpo- lated from the TGLDM vary with the fuel composition. The tempera- tures for methane/ethane jets are the highest in most cases, followed by pure methane jets and then methane/nitrogen jets. The most ob- vious difference is observed at small values of YCO2 and YH2O, corre- sponding to conditions at the start of ignition. For large values of the progress variables, the difference is small; 〈T |ζ〉 for pure methane and methane/nitrogen jets are almost the same, which corresponds to the conditions approaching equilibrium. As a result, it is expected that pre- dictions of ignition delay change with fuel composition. In this study, the CSE-TGLDM model predicts ignition delay being the shortest for methane/ethane jets and the longest for jets diluted with nitrogen. Additives also affect the mixing of fuel and oxidizer. Given the same injection parameters, more fuel is injected when C2H6 or N2 is added because of the larger molecular weight. The gaseous jet momen- tum is increased and mixing is enhanced. As a demonstration, the con- tours of stoichiometry Zst in the pure methane andmethane with addi- tives jets are studied at injection pressure 120 bar and pre-combustion temperature 1200 K. Contours of Zst shown in Fig. 3.10 are recorded 71 0 0.01 0.02 1300 1350 1400 YCO2 T (K ) 0 0.02 0.04 1550 1650 1750 YCO2 T (K ) 0 0.02 0.04 0.06 0.08 1800 1900 2000 2100 YCO2 T (K ) 0 0.05 0.1 2200 2400 YCO2 T (K ) Figure 3.9: Conditional averaged temperature evaluated using TGLDMs of different fuel composition, at stoichiometric mixture frac- tion and pre-combustion temperature 1200 K.−·−: 90%CH4+10%C2H6, —-: CH4, −−: 80%CH4 + 20%N2. Top left: YH2O = 0.02, top right: YH2O = 0.045, bottom left: YH2O = 0.07 , bottom right: YH2O = 0.095. 72 0 10 20 30 40 50 60 70 80 0 5 10 x/deq r / d e q (a) t=0.9ms 0 10 20 30 40 50 60 70 80 0 5 10 x/deq r / d e q (b) t=1.1ms Figure 3.10: Comparison of Zst contours in different jets. —-: pure me- thane; −−:80%CH4 + 20%N2; − · −: 90%CH4 + 10%C2H6. at t = 0.9ms (during injection), t = 1.1ms (after injection stops). The figures show that the methane/additive jets penetrate further than the pure methane jet. The Zst of methane/nitrogen jets is larger than those of methane and methane/ethane jets. When the contours at the same mixture fraction are compared for the three jets, the methane/nitrogen jet always penetrates furthest and expands widest, followed by meth- ane/ethane and then methane jets. The temperature fields of the different jets are compared during in- jection and after injection stops. Contours of mixture fraction are im- posed on temperature fields to compare the mixing of fuel and air. A similar observation as discussed previously is found in Fig. 3.11- 3.12, where enhanced mixing is exhibited in jets of methane mixed with ad- ditives. This comparison is of interest because of the possible applica- tion in compression ignition engines, where pilot fuel injection can be 73 10 20 30 40 50 60 0 5 10 x/deq r / d e q 400 600 800 1000 1200 (a) CH4 10 20 30 40 50 60 0 5 10 x/deq r / d e q 400 600 800 1000 1200 (b) 80%CH4 + 20%N2 10 20 30 40 50 60 0 5 10 x/deq r / d e q 400 600 800 1000 1200 (c) 90%CH4 + 10%C2H6 Figure 3.11: Temperature fields with Z contours at t = 0.9 ms. Black: Z = 0.1, 0.2, 0.3, ..., 0.9; white: Z = Zst. 74 10 20 30 40 50 60 70 80 0 5 10 x/deq r / d e q 600 800 1000 1200 (a) CH4 10 20 30 40 50 60 70 80 0 5 10 x/deq r / d e q 600 800 1000 1200 (b) 80%CH4 + 20%N2 10 20 30 40 50 60 70 80 0 5 10 x/deq r / d e q 600 800 1000 1200 (c) 90%CH4 + 10%C2H6 Figure 3.12: Temperature fields with Z contours at t = 1.1 ms. Black: Z = 0.1, 0.2, 0.3, ..., 0.9; white: Z = Zst. 75 used to ignite the gas near the tip of the main jet and the shape of the jet thus affects the engine performance. Altogether, the effects of additives are a result of the competition between chemical kinetics and mixing. The changes in ignition delay are related to the percentage and properties of the additives, as well as the pre-combustion conditions. At low pre-combustion tempera- tures, ignition delay seems to be mainly affected by changes in chem- ical kinetics, where ethane addition accelerates ignition and nitrogen retards ignition. At high temperatures, the formation of the ignition kernel seems to be mainly limited by the mixing of fuel and air, where the chemical kinetics of methane and methane with additives are fast enough and the changes in mixing caused by additives are not signifi- cant enough tomake any noticeable difference because of the relatively short time for mixing prior to ignition. 3.3.4 Effects of Different Injection Parameters The effects of different injection parameters are also studied with the CSE-TGLDM method. Simulations are performed for pure methane jets at various pressure ratios. In this series of simulations, injection pressure Pi is varied and predictions of ignition delay are compared to the measurements of Wu [103]. The characteristics of the jets are listed in Table 3.1 as Jet D. The predicted and measured ignition delays are plotted against the pressure ratio in Figure 3.13. The results show that the autoignition de- lay is reduced with increased pressure ratio. For this series of jets, the pre-combustion temperature of 1300 K is relatively high and turbulent mixing seems to be the rate-limiting process of autoignition. Increas- 76 1.5 2 2.5 3 3.5 4 4.5 5 5.5 0.5 1.0 τ (m s ) Figure 3.13: Ignition delay τ and pressure ratio. +: experimental data; : predictions of CSE-TGLDM. ing injection pressure helps to enhance mixing and therefore reduces ignition delay time. This decrease becomes negligible at relatively high pressure ratios, when the effect of enhanced mixing is no longer the dominant factor and the autoignition process seems to be limited by chemical kinetics. The effect of varying injection duration is studied at the pre-com- bustion temperature of 1300 K using the CSE-TGLDM method. The characteristics of the jet are listed in Table 3.1 as Jet E. The predicted ignition delays are compared with the experimen- tal data. Despite the scatter in the measurements, the numerical pre- dictions show that at the specified high pre-combustion temperature, ignition delay is insensitive to the injection duration. In fact, the pro- files of maximum 〈T |ζ〉 for all injection durations coincide with one 77 1 1.5 2.0 2.5 0.5 1.0 1.2 ti (ms) τ (m s ) Figure 3.14: Ignition delay τ and injection duration ti. Cross: experi- mental data; : predictions of CSE-TGLDM. other during the autoignition process and only slight differences are observed later in the burning stage. The irrelevance of injection duration to ignition delay, however, is inconclusive. For the specified temperature 1300 K, the ignition delay is relatively short and autoignition starts before the end of injection in all cases under study. At lower temperatures, ignition might start after injection stops and the ignition delay might be affected by injection duration because of the changes in mixing behavior. 3.4 Conclusions This chapter has presented simulations of autoigniting jets of meth- ane and methane mixed with additives. The model predictions of the 78 ignition delay are consistent with experimental results. However, the scatter observed in the experiments is not accounted for in the present models. The discrepancy between the predictions of the two methods appears to be due to the different modeling techniques used for gen- erating the libraries of full chemistry. Ultimately, it is found that the CSE-TGLDM approach is superior to the CSE-LFD approach. Adding ethane to the methane jet results in a modest reduction in the ignition delay time at lower initial temperatures. However, this reduction is sig- nificantly smaller than the scatter in the experimental results. Diluting methanewith nitrogen delays autoignition at low pre-combustion tem- peratures and the effects become less obvious at high temperatures. In- creasing pressure ratio reduces ignition delay and the effects become negligible at high pressure ratios. The injection duration seems irrel- evant when autoignition starts before the end of injection; its signifi- cance under other circumstances needs further study. 79 Chapter 4 Modelling Radiative Effects of a Turbulent Non-premixed Flame using CSE-TGLDM 4.1 Introduction Unlike RANS, Large Eddy Simulation of turbulent combustion solves large scale motions and might capture the transient behavior of tur- bulent flames. Conditional Source-term Estimation (CSE) has been ap- plied to study non-premixed turbulent flames in LES. As discussed in Chapter 2, CSE-LFD and later CSE-TGLDM [62] has been used in LES to model the Sandia Flame D – a piloted, non-premixed turbulent jet flame that has been widely studied numerically and experimentally. In spite of the general good agreement between the predictions of CSE and the experimental data, there are some unanswered questions: • The radiation properties of Sandia Flame D remain to be mod- eled. The models to be used should be a compromise among pre- cision, complexity and computational cost. • The possibility of using CSE to describe Turbulence/Radiation Interaction (TRI) is to be tested. Similar to the Turbulence/Che- mistry Interaction, TRI is complicated by the nonlinearity of the 80 radiation process. Theoretically, TRI might be modeled by CSE in a similar way to that with which it closes chemical source terms. • The effects of neglecting radiation heat loss on the predictions of temperature and species formation in Flame D need to be ana- lyzed. Previous studies with CSE are all based on the assumption of adiabatic combustion. The CSE methods are to be extended to model radiation. In this chapter, the CSE-TGLDM method is extended and applied in LES of Sandia Flame D to answer these questions. 4.2 LES of Sandia Flame D using CSE-TGLDM The Sandia D-flame is a piloted, non-premixed CH4/Air flame. The main jet is composed of 25% CH4 and 75% air. The nozzle diameter is 7.2 mm and the Reynolds number is 22400. The pilot flame burns a pre- mixture of C2H2, H2, CO2, N2 and air having nominally the same equi- librium composition and enthalpy as the main jet. The jet is enclosed by an air co-flow with the velocity of 0.9 m/s. A detailed description for the configuration of this flame and experiment documentation are provided by Barlow et al. [41]. Themodellingwork discussed is implemented in the LES paradigm, using the LES code developed by Steiner and Bushe [56]. The filtered balance equations of continuity, momentum, energy and species in LES have been discussed in Chapter 2. In this work, the unresolved subgrid- scale convective fluxes are modeled using the eddy viscosity and eddy 81 diffusivity models [56, 59]: τ̃ = −ρ̄(ũu− ũũ)) = ρ̄ν̃t((∇ṽ) + (∇ṽ)T − 2 3 (∇ · ṽ)I), (4.1) −ρ̄(ũYj − ũỸj) = ρ̄D̃t,j ▽ Ỹj, (4.2) −ρ̄(ũh− ũh̃) = k̃t ▽ T̃ , (4.3) where ν̃t is the eddy viscosity, D̃t,j is the eddy diffusivity of species j, k̃t is the thermal eddy conductivity and I is an identity matrix. The turbulent transport coefficients ν̃t, D̃t,j and k̃t are computed using the dynamic model for compressible flows by Moin et al. [17, 18]. LES is performed in spherical coordinates on a computational grid of 192 × 84 × 48 control volumes in the streamwise, cross-stream and azimuthal directions respectively [56,59]. A finite volumemethodwith second order accuracy in space is applied. A predictor-corrector-projec- tion scheme is used to integrate in time, with a second-order Adams- Bashford scheme in the predictor step and a semi-implicit Crank-Nichol- son scheme in the corrector step. A timestep splitting technique is used in the computation of conditional averaged source terms to overcome the stiffness of the chemical kineticmechanism. The inflow boundary is defined according to the experimental data, with velocity fluctuations superimposed on the measured mean velocity profile. A convective boundary condition [106] is used for the outflow boundary to stabilize the simulation. A traction-free condition [107,108] is used at the lateral boundary to allow for a flux of ambient fluid across the boundary of the spreading jet. The code is developed using theMessage-Passing In- terface (MPI) [109] for parallel computation on clusters. The chemical source-terms are closedwith the extended CSE-TGLDMwith radiation model. 82 4.3 CSE-TGLDMwith Radiation 4.3.1 Radiation Models Different models can be adopted to describe radiation in combustion, such as the optically thin approximation [110], discrete ordinates [111] and discrete transfer methods [112]. For turbulent flames, the radiation properties might be evaluated using the filtered values of species and temperature. However, similar to Turbulence/Chemistry interactions, the fluctuations around the means of relevant scalars affect the radia- tion quantities. This effect, or Turbulence/Radiation interaction (TRI), leads to closure problems similar to chemical closure. Coelho [113,114] studied Flame D using various models to evaluate radiation characteristics of gaseous species, combined with different methods to account for TRI. The pros and cons of a few radiation mod- els, gas radiation property models and TRI modelling methods were discussed. In terms of flame radiation, it was suggested that the opti- cally thin approximation (where flame absorption is neglected) might be adopted in the simulation of non-luminous flames; more accurate calculationsmight be done by solving radiative transfer equationswith the discrete ordinates method (DOM), where the absorption of radia- tion was described as a function of local conditions. As for gas radi- ation properties, the spectral line-based weighted-sum-of-gray-gases (SLW) model was suggested as an alternative to the weighted-sum- of-gray-gases (WSGG) model. Comparison was made among calcula- tions based on a variety of combinations: assuming an adiabatic flame or optically thin flame, or solving radiative transfer equations (RTE) combined with a SLW radiation property model or a model where ab- 83 sorption coefficients are approximated by curve fitting; combined with methods that fully or partially account for TRI effects. It was found that the optically thin approach tends to over-predict the radiation heat loss; a complete RTE method only slightly reduces the overprediction when the medium is treated as gray and Planck mean absorption coef- ficients are adopted. Accurate predictions are obtained with the com- bination of complete RTE and SLW models. Overall, the temperature difference between the adiabatic calculation and those with various ra- diation models was found to be less than 150 K in all calculations. 4.3.2 Formulation To include radiation effects in the combustion simulation using the CSE-TGLDMmodel, several problems need to be solved. To beginwith, the radiative properties of gaseous species need to be modelled with precision at an affordable cost. Secondly, a radiation model must be chosen. To use complicated models such as those solving the Radiative Transfer Equation (RTE) may significantly increase the computational cost of LES. A relatively simple radiation model with acceptable accu- racy is needed to evaluate the radiation heat loss term in the energy equation. Thirdly, Turbulence/Radiation Interaction (TRI) needs to be accounted for in the calculation. Similar to the Turbulence/Chemistry Interaction, TRI is caused by turbulent fluctuations. Lastly, the TGLDM method for reducing detailed chemical kinetics has so far been based on the assumption of adiabatic combustion. A TGLDMmethod includ- ing the effects of radiation needs to be developed. Considering the performance of different models and additional computational cost, an optically thin flame model is used in this work, 84 combined with a gas radiation property model derived from the RAD- CAL model by Grosshandler [115]. It is assumed that the flame is opti- cally thin and each radiation point source has an unimpeded isotropic view of the cold surroundings. The radiation heat loss is calculated with [116]: Q = Q(T, Yi) = 4σ ( nr∑ 1 pi · ap,i ) · (T 4 − T 4∞) , (4.4) where Q[W/m3] is the radiation heat loss rated per unit volume, σ = 5.669×10−8[W/m2K4] is the Steffan-Boltzmann constant, nr is the num- ber of the species included in the calculation, pi is the partial pressure of species i in atmospheres (equal to mole fraction times local pres- sure), ap,i[m −1atm−1] is the Plankmean absorption coefficient of species i and T∞ is the background temperature. In the context of the optically thin model, the inclusion of species such as CH4 and CO is not essential for methane flames and may be neglected. Moreover, previous study [117] shows thatCO2 andH2O are the major radiating species in Sandia Flame D. Only these two species are considered in the computation of radiation heat loss. The Planck mean absorption coefficients of the gaseous species are evaluated as a function of temperature by curve fitting within the temperature range of 300 K to 2500 K based on the RADCAL model: ap = c0+c1∗(1000 T )+c2∗(1000 T )2+c3∗(1000 T )3+c4∗(1000 T )4+c5∗(1000 T )5, (4.5) with the constant coefficients listed in Table 4.1 [116]. The CSE method is proposed to account for TRI where the nonlin- earity is partially accounted for with conditional averaging. The radia- 85 Table 4.1: Coefficients for Planck mean absorption coefficients. species c0 c1 c2 c3 c4 c5 H2O -0.2309 -1.1239 9.4153 -2.9988 0.51382 −1.8684× 10−5 CO2 18.741 -121.31 273.50 -194.05 56.310 -5.8169 tion heat loss in the turbulent flame is calculated with: 〈Q|ζ〉 = Q(〈T |ζ〉, 〈Yi|ζ〉), (4.6) Q̃(x, t) ≈ ∫ 1 0 〈Q|ζ〉 P̃ (ζ;x, t) dζ, (4.7) where the radiation properties of gaseous species are evaluated using conditional means. The radiative heat loss is presented in the energy equation as a sink term. The effects of radiation heat loss also need to be considered in the TGLDM calculation. If the radiation term is directly included in the en- ergy equation in the TGLDM computation, the trajectories will deviate from the desired equilibrium states and tail off to extinction. Moreover, radiation quantities at a particular location in a flame is related to local conditions as well as global conditions at distant locations. Therefore, it is physically unreasonable to develop trajectories with an arbitrary heat loss term or one evaluated using only the instantaneous solutions of species and temperature. A simple solution is to assume that radiation of the flame is not strong enough to significantly affect the trajectories. This simplifica- tion might be acceptable since the overall fraction of radiative heat loss in Flame D is only around 5%. The formulation of such a me- thod would be almost the same as the adiabatic CSE-TGLDM, where 86 an adiabatic TGLDM library will be used and the only difference is the sink term added to the energy equation in the combustion simu- lation. Alternatively, an additional parameter, the enthalpy defect can be introduced in the TGLDM calculation to describe the relationship between the manifolds and the local radiation conditions in the flame. The concept of enthalpy defect has been used in laminar flamelet models to account for non-adiabatic effects [118]. It is the difference between the actual and the adiabatic enthalpy. The specific enthalpy of a gaseous mixture can be calculated with: h = ns∑ 1 Yihi = ns∑ 1 ( h0i + ∫ T T 0 cp,idT ) , (4.8) where ns is the number of species in the mixture and h 0 i is the enthalpy of formation of species i at the standard reference temperature T 0. The adiabatic enthalpy is: had = ho + Z(hf − ho), (4.9) where ho and hf are the specific enthalpies of the oxidant and fuel in the fresh mixture respectively. The enthalpy defect△h is defined as: △h = h− had. (4.10) The calculation of TGLDM with enthalpy defect as an additional parameter involves the computation ofmultiple adiabatic TGLDMs for mixtures at various pre-combustion temperatures. Firstly, an adiabatic TGLDM can be solved at the given pre-combustion temperature T0 of fuel and oxidant for a sampled mixture fraction, where h0 = had(T0) and△h0 = 0. Another adiabatic TGLDM can be generated for the same mixture fraction at a lower initial temperature T1, where h1 = had(T1). 87 For this TGLDM, the enthalpy defect is: △h1 = h1 − had(T0) = had(T1)− had(T0). (4.11) Several TGLDMs with different enthalpy defect values can be gener- ated by changing the initial temperature of fuel and oxidant. Multi- ple TGLDMs can be generated at temperatures T0 > T1 > T2 > ... > Tj , where j is the number of TGLDM libraries and Tj should be low enough to account for themaximum radiation heat loss observed in the flame under study. These TGLDMs are then tabulated and retrieved later during combustion simulation. The modified CSE-TGLDMmethod can now evaluate reacting scal- ars using manifolds with the same enthalpy defects as those observed locally in the flame. The local enthalpy defect of the flame can be calcu- lated directly using the instantaneous solutions of species and temper- ature, as shown in Eq. 4.8 and Eq. 4.10. Once the enthalpy defect △h is calculated, interpolation needs to be performed based on progress variable (as in the adiabatic CSE-TGLDM) as well as△h. More specifi- cally, TGLDMs of two initial temperatures can be found such that: △hj1 ≤ △h ≤ △hj1+1, (4.12) where△hj1 is the enthalpy defect for the manifold of an initial temper- ature Tj1. Interpolation along the enthalpy defect dimension can be performed implicitly, assuming that reactive scalars can be retrieved from the TG- LDM using the progress variables. At the same time, 〈T |ζ〉 can be eval- uated using the CSE inversion method: T̃ ≈ ∫ 1 0 〈 T |ζ 〉 A P̃Z, j (ζ) dζ, j ∈ A. (4.13) 88 Supposing the local enthalpy defect in the flame corresponds to a man- ifold of an initial temperature Tr, ideally, the temperature interpolated using this manifold 〈T |ζ〉(Tr) should be equal to the inverted one at convergence: 〈T |ζ〉 ≈ 〈T |ζ〉(Tr). By comparing the interpolated and inverted 〈T |ζ〉, one can decide which manifold has the same enthalpy defect as that of the local mixture and should be used to evaluate chem- ical source terms. In practice, TGLDMs are generated and stored for selected initial temperatures T0, T1, ..., Tj , which might not include the desired tem- perature Tr. In this case, interpolation in the dimension of enthalpy defect is needed. Specifically, two TGLDM libraries are selected such that: 〈T |ζ〉(Tj1) ≤ 〈T |ζ〉 ≤ 〈T |ζ〉(Tj1+1), (4.14) Here, 〈T |ζ〉(Tj1) = g(〈YCO2|ζ〉, 〈YH2O|ζ〉), (4.15) is the conditional averaged temperature recovered by simple interpo- lation over the manifold generated at the initial temperature of Tj1. Eq. 4.14 indicates Tj1 ≤ Tr ≤ Tj1+1. The conditional mean of a reactive scalar 〈f |ζ〉 can then be calculated with: 〈f |ζ〉 = 〈f |ζ〉(Tj1) ∗ (1− ψ) + 〈f |ζ〉(Tj1+1) ∗ ψ, (4.16) where 〈f |ζ〉(Tj1) is the conditional mean interpolated using the TG- LDM manifold generated at the initial temperature Tj1. Interpolation in the direction of enthalpy defect is performed between manifolds of Tj1 and Tj1+1, where the interpolation factor is evaluated with: ψ = 〈T |ζst〉 − 〈T |ζst〉(Tj1) 〈T |ζst〉(Tj1+1)− 〈T |ζst〉(Tj1) , (4.17) 89 where ζst is the stoichiometric mixture fraction; or by solving the fol- lowing problem: min { nz∑ 1 ( 〈T |ζ〉(Tj1) ∗ (1− ψ) + 〈T |ζ〉(Tj1+1) ∗ ψ − 〈T |ζ〉 )2} , (4.18) where nz is the number of sample mixture fractions. 4.3.3 Results and Discussion Large Eddy Simulation of Sandia Flame D is performed using the ex- tended CSE-TGLDM method. GRI-Mech 2.11 [119] is used instead of GRI-Mech 3.0 [120] because it gave more accurate predictions of NO in previous study. A transport equation is solved for NO, with chemical source terms retrieved from the TGLDM library. This leads to a more accurate prediction compared to extracting NO from the TGLDM li- brary because of the relatively slow process of NO formation. The sim- ulation results are compared to experimental data as well as Wang et al.’s simulation [59] where adiabatic combustion was assumed. Figs. 4.1 to 4.5 show the centerline profiles and radial profiles of temperature and species. Generally speaking, the predictions of both adiabatic CSE-TGLDM [59] and CSE-TGLDM with radiation are in good agreement with the experimental data. The comparison of tem- perature andmajor products such as CO2 andH2O showsminor differ- ences between these twomethods, which indicates that the approxima- tion of adiabatic combustion is likely acceptable for the weakly radiat- ing Flame D. The predictions of CO2 and H2O are slightly improved in the far field, but the improvement is slight. On the other hand, the predictions of NO are quite different. As ex- pected, the adiabatic CSE-TGLDM over-predicts NO formation, espe- 90 cially in the far field where radiation is relatively strong. CSE-TGLDM with radiation gives better predictions of NO, which agree with the ex- perimental data well. Slight over-prediction is observed at x/D = 15 and under-prediction at other stations. The over-prediction at x/D = 15might be a result of accumulated error in the near field where the ef- fects of the pilot flame are not precisely described. The under-prediction of NO at other stations is expected because of the optically thin model used for radiation, which tends to over-predict heat loss in radiative CH4 flames due to the strong absorption by CO2 around the 4.3 µm wave length. This under-prediction might be improved if a more accu- rate radiation model including gaseous species absorption is used to evaluate the energy sink term in combustion simulation. Figs. 4.6 to 4.8 show the conditional means of progress variables CO2 and H2O as well as NO. As observed in the unconditional means, the predictions of 〈YCO2|ζ〉 and 〈YH2O|ζ〉 are slightly improved in the mid- and far-field using CSE-TGLDM with radiation; however, the difference is small. Fig. 4.8 shows the conditional means of NO. The under-prediction of NO caused by the optically thin model is more ob- vious in the conditional means. Overall, the adiabatic computation and that including radiation both capture the trends observed in the exper- iment. Compared to the adiabatic CSE-TGLDM method, the requirement of storage for CSE-TGLDM with radiation increases dramatically be- cause of the use of multiple TGLDMs. In terms of speed, the com- putation is slowed down because interpolation over each TGLDM is needed in order to evaluate the enthalpy defect level and to close the 91 0 20 40 60 0 2 4 6 8 x/D T/ T0 0 20 40 60 0 0.05 0.1 0.15 0.2 x/D Y C O 2 0 20 40 60 0 0.05 0.1 0.15 0.2 x/D Y H 2o 0 20 40 60 0 0.5 1 1.5 x 10−4 x/D Y N O Figure 4.1: Centerline profiles of Flame D. ∆: experimental data; — -: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE- TGLDM with radiation. 92 0 2 4 6 0 2 4 6 r/D T/ T0 x/D=7.5 0 2 4 6 8 0 2 4 6 r/D T/ T0 x/D=15 0 2 4 6 8 0 2 4 6 8 x/D=30 r/D T/ T0 0 5 10 0 2 4 6 8 x/D=45 r/D T/ T0 Figure 4.2: Temperature radial profiles at different stations down- stream. ∆: experimental data; —-: predictions of adiabatic CSE- TGLDM; − · −: predictions of CSE-TGLDM with radiation. 93 0 2 4 6 0 0.02 0.04 0.06 0.08 0.1 r/D Y C O 2 x/D=7.5 0 2 4 6 8 0 0.02 0.04 0.06 0.08 0.1 r/D Y C O 2 x/D=15 0 2 4 6 8 0 0.05 0.1 0.15 0.2 x/D=30 r/D Y C O 2 0 5 10 0 0.05 0.1 0.15 0.2 x/D=45 r/D Y C O 2 Figure 4.3: Mass fraction of CO2 radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE- TGLDM; − · −: predictions of CSE-TGLDM with radiation. 94 0 2 4 6 0 0.02 0.04 0.06 0.08 0.1 r/D Y H 2O x/D=7.5 0 2 4 6 8 0 0.02 0.04 0.06 0.08 0.1 r/D Y H 2O x/D=15 0 2 4 6 8 0 0.05 0.1 0.15 0.2 x/D=30 r/D Y H 2O 0 5 10 0 0.05 0.1 0.15 0.2 x/D=45 r/D Y H 2O Figure 4.4: Mass fraction of H2O radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE- TGLDM; − · −: predictions of CSE-TGLDM with radiation. 95 0 2 4 6 0 2 4 6 x 10−5 r/D Y N O x/D=7.5 0 2 4 6 8 0 2 4 6 8 x 10−5 r/D Y N O x/D=15 0 2 4 6 8 0 0.2 0.4 0.6 0.8 1 x 10−4 x/D=30 r/D Y N O 0 5 10 0 0.5 1 1.5 x 10−4 x/D=45 r/D Y N O Figure 4.5: Mass fraction of NO radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE- TGLDM; − · −: predictions of CSE-TGLDM with radiation. 96 0 0.5 1 0 0.05 0.1 z < Y C O 2|ξ > x/D=7.5 0 0.5 1 0 0.05 0.1 z < Y C O 2|ξ > x/D=15 0 0.5 1 0 0.05 0.1 z < Y C O 2|ξ > x/D=30 0 0.5 1 0 0.05 0.1 z < Y C O 2|ξ > x/D=45 Figure 4.6: Conditional averages of CO2 radial profiles at different sta- tions downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 97 0 0.5 1 0 0.05 0.1 z < Y H 2o |ξ> x/D=7.5 0 0.5 1 0 0.05 0.1 z < Y H 2o |ξ> x/D=15 0 0.5 1 0 0.05 0.1 z < Y H 2o |ξ> x/D=30 0 0.5 1 0 0.05 0.1 0.15 z < Y H 2o |ξ> x/D=45 Figure 4.7: Conditional averages of H2O radial profiles at different sta- tions downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 98 0 0.5 1 0 2 4 6 x 10−5 z < Y N O |ξ> x/D=7.5 0 0.5 1 0 2 4 6 8 x 10−5 z < Y N O |ξ> x/D=15 0 0.5 1 0 2 4 6 8 x 10−5 z < Y N O |ξ> x/D=30 0 0.5 1 0 0.5 1 x 10−4 z < Y N O |ξ> x/D=45 Figure 4.8: Conditional averages of NO radial profiles at different sta- tions downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 99 reaction rates accordingly. The additional dimension of enthalpy de- fect also leads to more unsmoothness in the interpolated reaction rates, which might cause the ODE solver to choose smaller sub-timesteps to converge solutions in the time-splitting. In this study of Flame D, the radiation is not strong and only three TGLDMs generated at differ- ent enthalpy defect levels are used; the increase of computational cost is still manageable. However, a more dense grid in enthalpy defect is necessary for highly radiating flames and CSE-TGLDMwith radiation needs to be optimized. One suggestion is to combine the parallel com- putation technique with the extended CSE-TGLDM method, where a segment of the flame simulated on each individual processor has a rel- atively small range of enthalpy defect and accordingly only a small number of TGLDMs (that might differ from processor to processor) need to be loaded and used on each processor. Altogether, the adiabatic combustion and optically thin flame as- sumptions are both acceptable in the study of Sandia Flame D. The adiabatic CSE-TGLDM method can be used to study non-premixed turbulent flames without strong radiation, such as Sandia Flame D; the method predicts temperature and major species (for which trans- port equations are solved) reasonably well. The radiation effects can be incorporated in CSE-TGLDM by introducing enthalpy defect as an ad- ditional variable to parameterize the manifold and interpolating based on the progress variables as well as the enthalpy defect. The extended CSE-TGLDM with radiation method improves the prediction of NO formation. Optimization of CSE-TGLDMwith radiation should be fur- ther studied to deal with the increase of storage requirement and com- 100 putational cost. 4.4 Radiation Intensities Analysis using CSE 4.4.1 Formulation The scalar fields predicted with CSE-TGLDM are used to evaluate the mean spectral radiation intensities in the Flame D. The radiation trans- fer equations are solved for various radiation paths with RADCAL, a program based on a narrow band radiation model [115]. To account for TRI in the analysis of radiation intensities, a stochastic time and space time series simulation (TASS) method might be adopted [121]. TASS analyzes statistics of radiation properties using many realizations of the flame under study. In this work, the CSE model is combined with the RADCAL model to simulate TRI. In the RADCALmodel, the radiative transfer equation (RTE) for an absorbing and emitting medium is solved by breaking the line-of-sight into uniform elements. The solution of the spectral intensity i′λ is [115]: i′λ(l) = i ′ λ,we −κλ(l) + ∫ κλ(l) 0 ib,λ(l ∗)exp[−(κλ(l)− κλ(l∗))]dκλ(l∗) (4.19) where ib,λ is the Planck blackbody function, κλ is the optical thickness defined as κλ = ∫ l 0 aλ(l ∗)dl∗, λ is the wavelength and the subscript w refers to the boundary condition of a wall. A rough estimate of themean spectral intensities can be obtained by using the means of relevant scalars in the RADCAL model. In such a calculation, due to the nonlinearity of Eq. 4.19, the effects of fluctuating scalar fields on the radiation intensities are neglected. To correct for this, the CSE model can be combined with the RADCAL model in a way similar to the closure for the reaction rates. The contributions of 101 uniform elements to the overall radiation intensity can be evaluated using conditional means: i′λ(l ∗, l) = ∫ 1 0 〈i′λ(l∗, l)|ζ〉P (ζ)dζ, (4.20) where 〈i′λ(l∗, l)|ζ〉 is evaluated with conditional averages of reactive scalars. The mean radiation intensity i′λ(l) is the summation of the con- tributions from all elements. 4.4.2 Results and Discussion The radiation characteristics of the pilotedmethane-air diffusion flame, Sandia Flame D are studied. Large Eddy Simulation of the Flame D has been done with the CSE-TGLDMmodel discussed previously. The code has been revised and re-run to extract the reactive scalars needed in the RADCAL model. An example of the conditional means at dif- ferent stations downstream of the nozzle has been shown in Figs. 4.6 and 4.7. The mean spectral radiation intensities along diametric paths at x/D = 30, x/D = 45 and x/D = 60 are calculated using the mean properties method and the CSE method as described in the previous section. The results are compared with the experimental data of Zheng et al. [122] in Fig. 4.9. Predictions of both methods appear to capture the trends in the measurements. The spectra show that the main sources of radiation in Flame D are water vapor and CO2. The strongest radiation intensity is observed at x/D = 45, which is close to the stoichiometric flame height of x/D = 47. The predictions of the CSE method are generally higher than the mean value calculation as a result of accounting for TRI. The radiation quantities predicted with the CSE method seem to be 102 closer to the experimental data only at x/D = 60. However, this does not necessarily mean that using the mean values is the better choice. In the work of Zheng et al. [122], the experimental data was used in the mean value method, and this lead to as much as a 50% underpredic- tion. The seemingly good performance of the mean properties method here might be as much attributable to a slight discrepancy between the predicted mean scalar fields and those of the real flame. Zheng et al. [122] also provided predictions of TRI using the com- plicated TASS method with the measured scalars as input data. The CSE predictions here are 15% to 22% higher than those of the mean properties method around peak values; this magnitude of discrepancy is similar to that observed with the TASS method. Ultimately, the CSE method seems promising for modelling TRI. For comparison, the spectral radiation intensities are also evalu- ated [117] with the CSE method using scalar fields extracted from the adiabatic calculation by Wang [59]. As shown in Fig. 4.10, the discrep- ancy is negligible, which might be explained by the small difference in temperature, CO2 and H2O predictions of the adiabatic CSE-TGLDM and CSE-TGLDM with radiation methods. It again proves that Flame D is not strongly radiating and an adiabatic assumption is acceptable for such a flame. This section has presented the analysis of spectral radiation inten- sities in the turbulent, piloted, Sandia Flame D. The CSE method has been used to model TRI, predicting higher radiation intensities than the mean value method. The model predictions are consistent with the experimental data. 103 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 7000 λ,µm I λm e a n , W /m 2 _ sr _µ m (a) x/D = 30 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 7000 8000 λ,µm I λm e a n , W /m 2 _ sr _µ m (b) x/D = 45 Figure 4.9: Spectral radiation intensities in Flame D. ◦: measured; —-: CSE method; −−: mean properties method. 104 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 λ,µm I λm e a n , W /m 2 _ sr _µ m (c) x/D = 60 Figure 4.9: Spectral radiation intensities in Flame D (cont.). ◦: mea- sured; —-: CSE method; −−: mean properties method. 4.5 Conclusion The CSE method has been used to analyze the radiation properties of the Sandia Flame D. The study shows that CSE can be used to model the Turbulence/Radiation interaction (TRI) effects. The CSE-TGLDM method has been extended with the enthalpy defect as an additional dimension in the manifolds. It has been combined with an optically thin radiation model and applied in the LES of Flame D. The predic- tions of NO are improved using the extended CSE-TGLDM method. The mean spectral radiation intensities in the flame have been ana- lyzed by solving radiative transfer equations and modelling TRI with CSE. The results agree well with experimental results. 105 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 7000 λ,µm I λm e a n , W /m 2 _ sr _µ m (a) x/D = 30 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 7000 8000 λ,µm I λm e a n , W /m 2 _ sr _µ m (b) x/D = 45 Figure 4.10: Spectral radiation intensities in Flame D, ◦: measured, —-: CSE-TGLDM with radiation, −−: adiabatic CSE-TGLDM. 106 1 1.5 2 2.5 3 3.5 4 4.5 5 0 1000 2000 3000 4000 5000 6000 λ,µm I λm e a n , W /m 2 _ sr _µ m (c) x/D = 60 Figure 4.10: Spectral radiation intensities in Flame D (cont.). ◦: mea- sured; —-: CSE-TGLDM with radiation; −−: adiabatic CSE-TGLDM. 107 Chapter 5 A Progress Variable and Presumed Probability Density Function for Premixed Combustion 5.1 Introduction As discussed in Chapter 2, partially premixed and premixed combus- tion also have wide applications in industry. Numerical simulation promises to assist in the design of premixed combustion devices [34, 123]. One major category in the models for turbulent premixed flames is statistical models that use a presumed probability density function (PDF) [124–128]. The PDF approaches are not limited to the thin flame- lets regime and can be used to account for finite rate chemistry. De- tailed chemistry can be included by using tabulated data. Because of these advantages, PDF models have been applied in RANS and LES of premixed flames [129]. Recently, efforts have been made to introduce the well-established Conditional Moment Closure (CMC) model [39, 40], which has been widely used to model non-premixed flames [43, 45, 46, 48, 50], to pre- mixed flame simulations. Klimenko and Bilger developed the theory of using CMC for premixed flames [42]. The method was applied in the study of premixed hydrogen-air flames by Swaminathan and Bil- 108 ger [130] and in modelling lean premixed gas turbines by Martin [131]. Although effective, CMC is computationally expensive because it in- troduces an additional dimension while solving the transport equa- tions for conditional means. In previous chapters, Conditional Source-term Estimation (CSE) has been shown to be a promising closuremethod for non-premixed flames [7, 57, 59, 60]. In CSE, conditional means are obtained by inverting in- tegral equations and chemical source terms are closed by invoking the first conditional moment hypothesis from CMC. It only requires the CMC closure hypothesis and the assumption of statistical homogene- ity, which do not exclude its application in premixed flames. To apply CSE to premixed flames, a suitable conditioning variable is needed. Since the fluctuations in temperature and mass fractions in premixed flames are mostly related to a combustion progress vari- able, this seems to be a logical choice for this purpose. Bray et al. [132] studied three different PDFs of a species-based progress variable c us- ing DNS with single step chemistry. The definition of c was such that c = 0 in reactants and c = 1 in products. A β-PDF, a twin delta func- tion PDF and a PDF based on unstrained laminar flame properties were adopted to predict the mean reaction rates in turbulent premixed flames with finite rate chemistry. Although the β-PDF and the lam- inar flamelet based PDF agreed with the DNS when the normalized variance g = c′2/[c(1 − c)] was around 1, the precision of modelled PDFs under other conditions was not tested. It was also suggested that the laminar flamelet based PDF be used only in the thin flamelet regime, where the Kolmogorov scales of the turbulence are larger than 109 the characteristic scales of the laminar flame. In this study, a product-based progress variable is introduced. Dif- ferent presumed PDF models are studied and a modification is sug- gested to improve the precision of the laminar flame-based PDFmodel of Bray et al. [132]. The DNS data from Grout’s work [133] is used to evaluate the precision of different PDF models in the RANS and LES paradigms. 5.2 The PDF Models A combustion progress variable can be defined with the mass fraction of products: c(x, t) = YP(x, t)/YP,eq, (5.1) where YP,eq is the mass fraction of products at equilibrium. By defi- nition, c = 0 in pure reactants and c = 1 in fully burnt mixture at equilibrium. The presumed PDF takes the form: P (c∗;x, t) ≈ p(c∗; c, c2), (5.2) where, presumably, a transport equation would be solved for c and possibly also for c2 in a real implementation in a CFD code. The re- action rates are zero at c = 0 and c = 1, which is similar to what is seen in non-premixed flames: in non-premixed flames, reaction rates are typically also zero at a mixture fraction Z = 0 and Z = 1. An example of P (c∗) taken from the DNS data of Grout [133] is il- lustrated in Fig. 5.1, showing the distribution of c at different locations in a turbulent premixed flame. Near the leading edge, where combus- tion has barely started, and the trailing edge, where the reactions ap- 110 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 35 c p(c ) x=10 x=30 x=60 x=90 x=120 Figure 5.1: Probability Density Function in a turbulent premixed flame. proach equilibrium, a bimodal distribution is observed and might be captured by a β-PDF. In the interior of the flame, because of the inter- mediate reactions, the formation of products might be incomplete even if the fuel is mostly consumed. A flame-based PDF might capture the distribution under such conditions if the premixed turbulent flame is extremely thin. 5.2.1 β-PDF Although the β-PDF approximates the real distribution only in bimodal distributions, it is tested here because of its wide use in both premixed and non-premixed flames. A β-PDF of c is defined as: P (c∗;x, t) = C (c∗)a−1 (1− c∗)b−1 , (5.3) where a, b and C are functions of local statistical moments of c: a = c ( c(1− c) (c2 − c2) − 1 ) (5.4) 111 b = a(1− c) c (5.5) C = 1 β(a, b) , (5.6) where β(a, b) = ∫ 1 0 (c∗)a−1 (1− c∗)b−1dc∗ is the usual β function. 5.2.2 A Laminar Flame-based PDF A laminar flame-based PDF model has been proved to work well in the thin flamelet regime [132]. The model is based on the observation that in the limit of a thin flame, the iso-surfaces of c are all parallel. The interior of the PDF is assumed to be inversely proportional to | ▽ c|. The PDF is defined as: P (c∗;x, t) = Aδ(c∗) +Bf(c∗) + Cδ(1− c∗), (5.7) where f(c∗) is calculated using the solution of an unstrained laminar flame: f(c∗) = 1 ▽c∗ . (5.8) The coefficients A, B and C are obtained based on the first three moment equations of c: 1 = ∫ 1 0 P (c∗)dc∗ = A+B ∫ 1 0 1 ▽c∗dc ∗ + C (5.9) c = ∫ 1 0 c∗P (c∗)dc∗ = B ∫ 1 0 c∗ ▽c∗dc ∗ + C (5.10) c2 = ∫ 1 0 c∗2P (c∗)dc∗ = B ∫ 1 0 c∗2 ▽c∗dc ∗ + C. (5.11) 112 From these equations the coefficients can be easily obtained: B = c− c2 I1 − I2 (5.12) C = c2I1 − cI2 I1 − I2 (5.13) A = 1−B ∗ I0 − C, (5.14) where I0, I1 and I2 are defined as: I0 = ∫ 1 0 1 ▽c∗dc ∗ (5.15) I1 = ∫ 1 0 c∗ ▽c∗dc ∗ (5.16) I2 = ∫ 1 0 (c∗)2 ▽c∗ dc ∗. (5.17) This PDF model has an unsolved problem. The PDF is constructed using the solution of an unstrained laminar flame. Ideally, the domain for the planar laminar flame would extend from negative infinity to positive infinity. If one were to have an “ideal” simulation of the lam- inar flame with a domain extending from negative infinity to positive infinity, the PDF one would obtain from this would just be two deltas at 0 and 1. In reality, the computational domain of the laminar flame is finite and this affects the precision of the PDF as c → 0 and c → 1. Moreover, by definition, the profile of this PDF is determined by a lam- inar flamewhich is independent of the local conditions in the turbulent flame. Consequently, the PDF always rises dramatically towards infin- ity as c→ 0 and c→ 1, even if c is not near 0 or 1 and c2 is very small. 113 −Inf xmin x1 x2 xmax +Inf 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 c Figure 5.2: Progress variable c in an unstrained laminar flame. 5.2.3 Modified Laminar Flame-based PDF To address the issues in the laminar flame-based PDF, some modifi- cations to the original model are proposed. The idea is to find a part of the laminar flame that approximates the conditions in the turbulent flame characterized by c and c2. Fig. 5.2 shows the profile of c in an unstrained laminar flame. Recognizing that the PDF will be stored at discrete points in c-space, a domain of interest for the laminar flame can be determined by the discretization scheme. A common approach in implementing presumed PDFs is to store the PDF in a look-up table; it is usually easier to store the integral (Pdc)i = ∫ ci+∆c/2 ci−∆c/2 Pdc (5.18) at discrete points i in c-space. It is possible to refine the spacing locally in c-space to better resolve the reaction zone. However, for this simple description (and in the later tests), it is assumed that there will be a 114 constant spacing ∆c in c-space. The first and last intervals in c-space must be treated differently, since the space is bounded by 0 and 1, such that (Pdc)0 = ∫ ∆c/2 0 Pdc (5.19) and (Pdc)n = ∫ 1 1−∆c/2 Pdc (5.20) where (Pdc)n is the last interval in c-space; this way Eq. 5.18 is used over the interval 1 < i < n− 1. This discretization can be used to provide bounds for using the laminar flame solution in assembling a PDF for a turbulent flame. A lower bound can be defined at cmin = ∆c/2 and an upper bound at cmax = 1 − ∆c/2. The positions xmin and xmax correspond to locations in the laminar flame where c = cmin and c = cmax. The modelled PDF would then comprise δ functions at 0 and 1 (at (Pdc)0 and (Pdc)n in the table) and the function f(c∗) over all of the other intervals, as was done in the original formulation of Bray et al. [132]. Further, however, it is possible to eliminate one or both of the δ functions from the PDF and truncate the function f(c∗) as required in the event that the local c and c2 in the turbulent flame indicate that this is needed. Given the local c and c2 in a turbulent flame, a truncated area [x1, x2] (with corresponding values of the progress variable [c1, c2]) can be lo- cated in the laminar flame such that the means and variance of c sam- pled within [x1, x2] match those in the turbulent flame. The functional form of the modified PDF changes depending on the locations of x1 and x2 relative to xmin and xmax. Four cases are possible: 115 • Case 1: xmin < x1 < x2 < xmax P (c∗;x, t) = 0 if c∗ < c1 or c∗ > c2B1f(c∗) if c1 ≤ c∗ ≤ c2 • Case 2: x1 < xmin < x2 < xmax P (c∗;x, t) = A2δ(c∗) +B2f(c∗) if c∗ ≤ c20 if c∗ > c2 • Case 3: xmin < x1 < xmax < x2 P (c∗;x, t) = 0 if c∗ < c1B3f(c∗) + C3δ(1− c∗) if c∗ ≥ c1 • Case 4: x1 < xmin < xmax < x2 P (c∗;x, t) = A4δ(c ∗) +B4f(c ∗) + C4δ(1− c∗). (5.21) Under all four circumstances, the first three moment equations of c can be used to locate the optimal area [x1, x2] within the laminar flame and evaluate the unknown coefficients. The unknowns are solved such that 1 = ∫ 1 0 P (c∗)dc∗ is satisfied and the deviation from the targeted c and c2 is minimized. The coefficients in the equations can be evaluated as: B1 = 1 I0 (5.22) B2 = c I1 (5.23) A2 = 1−B2I0 (5.24) 116 B3 = c− 1 I1 − I0 (5.25) C3 = 1−B3I0 (5.26) B4 = c− c2 I1 − I2 (5.27) C4 = c2I1 − cI2 I1 − I2 (5.28) A4 = 1−B4 ∗ I0 − C4. (5.29) I0, I1 and I2 are redefined as: I0 = ∫ xu xl 1 ▽c∗dc ∗ (5.30) I1 = ∫ xu xl c∗ ▽c∗dc ∗ (5.31) I2 = ∫ xu xl c∗2 ▽c∗dc ∗, (5.32) where xl = max(x1, xmin) and xu = min(x2, xmax). 5.3 Comparison with DNS 5.3.1 DNS Configuration The DNS [133] was performed in an inflow-outflow configuration us- ing the code SENGA [134], where the fully compressible Navier-Stokes equations were solved explicitly with the 3rd order Runge-Kutta time advances. The unburnt mixture consisted of 5%methane, 20% oxygen and 75% nitrogen at 900K. The flamewas solved in a 1283 domain, with 117 a physical edge domain of 6 mm, and the inflow velocity was forced to provide a decaying, isotropic, homogeneous turbulent field through which the flame propagates. The chemical mechanism used in the DNS of laminar and turbulent flames is a two-step mechanism: F +O → I + P (R1) I +O → 2P (R2) where F and O are model fuel and oxidizer which are analogous to CH4 and O2; I is a model intermediate which is a combination of CO andH corresponding to the full mechanism; P is a combination ofH2O and CO2. The rates of the two reactions, in [mol/cc/s] are calculated as: ω̇R1 = 7.8(10 14)e −7640 T [F ][R] (5.33) ω̇R2 = 8.9(10 16)T−1.8[O][R] (5.34) where [R] is the radical concentration given by: [R] = 4.5(105)e 2300 T [P ] [O]0.5[I]1.5e−λ √ 15 4 |F | |O| (5.35) The parameter λ is related to the ratio of the rate of the elementary fuel consumption reaction to the chain branching reaction. Reducing λ results in a disproportionate increase in the flame speed relative to its effect on the local gradient and giving a more easily resolvable flame front. 5.3.2 Comparisons of PDF Models The proposed PDFs are tested against the DNS of a premixed turbulent flame in the both the RANS and LES paradigms. For the RANS test, 118 the distributions of c are studied over the whole domain and on planes orthogonal to the direction of the mean flow – for this flow, the tur- bulence on these planes should be statistically homogeneous. The PDF profiles from the DNS are compared with the three presumed PDFs. Fig. 5.3 shows the distribution of c over the whole domain at 40 evenly spaced discrete points in c-space. The β-PDF captures the bi- modal shape of distribution, but the magnitude and location of the peaks are nearly symmetrical, which fails to capture the shape of the true PDF. The laminar flame-based PDF does a better job than the β- PDF, especially for large values of c. The original version of the laminar flame-based PDF requires that there be a δ at c∗ = 1, which does not match the behavior of the true PDF; the modified version of the lami- nar PDF more closely matches the true PDF, although it over-predicts the magnitude of the local maximum at c∗ ≈ 0.97. The difference be- tween PDFs from the original and the modified versions is more ob- vious in the Cumulative Distribution Function (CDF): the CDF calcu- lated with the original version departs from the DNS as c∗ increases, while the CDF calculated with the modified version almost perfectly reproduces the CDF from the DNS. Figs. 5.4 and 5.5 show the distributions of c on planes at different lo- cations. Overall, the β-PDF and the original laminar flame PDF predic- tions contain significant errors. On the other hand, the modified PDF agrees well with the DNS, although it generally fails to capture the exact location and magnitude of the peak on the products side of the distribution. For the LES tests of the PDFs, the distributions of c in cubes consist- 119 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 14 16 18 20 c∗ p( c∗ ) (a) Probability Density Function 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 c∗ ∫ p( c∗ )d c∗ (b) Cumulative Distribution Function Figure 5.3: Presumed Probability Density Function and Cumulative Distribution Function in a premixed turbulent flame. —-: DNS; −−: modified PDF model; − · −: original PDF by Bray et al.; · · · · ·: β-PDF. 120 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 c∗ p( c∗ ) (a) x/L = 0.234 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 c∗ p( c∗ ) (b) x/L = 0.391 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.99 0 5 10 15 20 25 c∗ p( c∗ ) (c) x/L = 0.625 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 10 20 30 40 c∗ p( c∗ ) (d) x/L = 1 Figure 5.4: Probability Density Function on planes in a premixed tur- bulent flame. —-: DNS; −−: modified PDF model; − · − line: original PDF by Bray et al.; · · · · ·: β-PDF. 121 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) x/L = 0.234 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) x/L = 0.391 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) x/L = 0.625 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) x/L = 1 Figure 5.5: Cumulative Distribution Function on planes in a premixed turbulent flame. —-: DNS; −−: modified PDF model; − · −: original PDF by Bray et al.; · · · · ·: β-PDF. 122 ing of 4 × 32 × 32 grid cells are analyzed. The results of the original laminar flame PDF are not presented since the improvements in the modified version are similar to those observed in the RANS tests. For this test, PDF and CDF tables were generated for 40 grid points in the space of the conditioning variable c. The tables are calculated using the modified laminar flame model for given values of c and c2. As a demonstration, the distributions in four cubes are shown in Figs. 5.6 and 5.7. The distribution calculated with the modified lami- nar flame PDF is obtained by look-up and linear interpolation in the pre-calculated tables. Observations from different cubes show that the β-PDF model predictions usually contain major errors. The interpo- lated values from tables of the modified PDF model agree well with the DNS, with the biggest discrepancy again being around the peak in the distribution towards the product side. Among the three presumed PDFs, the modified version of the lam- inar flame PDF model best captures the c distribution in the premixed turbulent flame. In the PDF plots, the maximum of the PDF from the DNS on the products side of the distribution tend to appear rounded off, likely due to turbulent mixing. This is typically the only feature that the modified PDF fails to capture. One way to improve this might be to add a filtering operator to the PDF: f̂(c∗) = ∫ +∞ −∞ f(c+)g(c+, c∗)dc+. (5.36) Further work would be necessary to define an appropriate filter kernel g. The discrepancy between themodified PDF and the DNS seems less significant in the CDF calculation. In CSE, the presumed PDFs are used 123 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 c∗ p( c∗ ) (a) c = 0.089, c2 = 0.048 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 c∗ p( c∗ ) (b) c = 0.406, c2 = 0.298 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 c∗ p( c∗ ) (c) c = 0.565, c2 = 0.451 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 c∗ p( c∗ ) (d) c = 0.800, c2 = 0.674 Figure 5.6: Probability Density Function in different cubes in a pre- mixed turbulent flame. —-: DNS; −−: modified PDF model; − · −: β- PDF. 124 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) c = 0.089, c2 = 0.048 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) c = 0.406, c2 = 0.298 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) c = 0.565, c2 = 0.451 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) c = 0.800, c2 = 0.674 Figure 5.7: Cumulative Distribution Function in different cubes in a premixed turbulent flame. —-: DNS; −−: modified PDF model; − · −: β-PDF. 125 to discretize integral equations: f ≈ J∑ i=1 〈f |c∗j〉P (c∗j; c, c2)dc∗j, (5.37) where J denotes the number of grid points in the conditioning vari- able space. In these equations, the precision of P (c∗j; c, c2)dc ∗ j is more of concern. It is more related to the CDF. Considering the good CDF predictions of the modified PDF, the discrepancy might be acceptable in the context of the CSE model. This argument will be explored in the CSE tests in Chapter 6. Another series of tests was performed to further analyze the preci- sion of the PDF models. The same RANS and LES tests are repeated for a second set of DNS data. This set of DNS data is of a turbulent premixed flame with lower turbulence intensity with a different two- step simplified chemical mechanism, where the coefficients in the Ar- rhenius expressions are tuned to accelerate the DNS computation. For clarity, the DNS data obtained using the original chemical kinetics will be referred to as solution ‘I’, and data obtained using the modified chemical kinetics as solution ‘II’ hereafter. The results of the a priori tests with data set II are shown in Fig. 5.8 to Fig. 5.12. As observed in the first series of tests, the modified laminar flame-based PDF model generally captures the distribution of c in the turbulent flame; overall, its predictions are better than those of the β- PDF and the original flamelet PDF model. The discrepancies between the modified PDF and the DNS from data set II are in similar locations to those found using set I. This suggests that this modified PDF may work well regardless of the chemical kinetics, although considerably more work would be required to state this conclusively. 126 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 c∗ p (c ∗ ) (a) Probability Density Function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p (c ∗ )d c∗ (b) Cumulative Distribution Function Figure 5.8: PDF and CDF in the whole domain. —-: DNS;−·−: original laminar flame-based model; −−: modified PDF model; · · · · ·: β-PDF. 127 0 0.2 0.4 0.6 0.8 1 0 5 10 15 c∗ p( c∗ ) (a) x/L = 0.234 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 c∗ p( c∗ ) (b) x/L = 0.391 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 c∗ p( c∗ ) (c) x/L = 0.625 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 70 c∗ p( c∗ ) (d) x/L = 1 Figure 5.9: Probability Density Function on different planes. —-: DNS; − · −: original laminar flame-based model; −−: modified PDF model; · · · · ·: β-PDF. 128 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) x/L = 0.234 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) x/L = 0.391 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) x/L = 0.625 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) x/L = 1 Figure 5.10: Cumulative Distribution Function on different planes. —-: DNS; − · −: original laminar flame-based model; −−: modified PDF model; · · · · ·: β-PDF. 129 0 0.2 0.4 0.6 0.8 1 0 20 40 60 80 100 c∗ p( c∗ ) (a) c̄ = 0.089,c̄2 = 0.046 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 c∗ p( c∗ ) (b) c̄ = 0.413,c̄2 = 0.282 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 16 18 c∗ p( c∗ ) (c) c̄ = 0.461,c̄2 = 0.311 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 c∗ p( c∗ ) (d) c̄ = 0.689,c̄2 = 0.505 Figure 5.11: Probability Density Function in different cubes. —-: DNS; −−: modified PDF model; − · −: β-PDF. 130 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) c̄ = 0.089,c̄2 = 0.046 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) c̄ = 0.413,c̄2 = 0.282 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) c̄ = 0.461,c̄2 = 0.311 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) c̄ = 0.689,c̄2 = 0.505 Figure 5.12: Cumulative Distribution Function in different cubes. —-: DNS; −−: modified PDF model; − · −: β-PDF. 131 5.3.3 Effects of Chemical Kinetics To study the sensitivity of the proposed PDF model to the chemical kinetics, the distribution of c in the second turbulent flame is evalu- ated with the modified PDF model; the predictions using the ‘right’ chemistry (based on laminar flame solution II) and the ‘wrong’ chemi- stry (based on laminar flame solution I) are compared to the DNS. The results are shown in Figs. 5.13 to 5.17. The results of both the RANS and LES tests show that the preci- sion of the modified laminar flamelet PDF is strongly affected by the chemical kinetics. When the laminar flame solution is based on a dif- ferent chemical mechanism, the performance of the PDF model seems to be far less satisfying, especially in the calculation of the CDF. The predictions of the β-PDF are not shown for better visibility; however, the modified PDF predictions are still closer to the DNS than the β- PDF even when the ‘wrong’ chemistry (i.e., laminar flame solution I) is used. The discrepancy between the PDF model predictions and the DNS is increased when different chemical kinetics are used for the solu- tions of the one dimensional laminar flame and the premixed turbulent flame. According to the modified PDF, the distribution of c between [0, 1] is proportional to f , the inverse of product gradient as defined in Eq. 5.8. In Fig. 5.18, comparison is made between normalized product gradient of the laminar flame solution I and II corresponding to two different chemical mechanisms. It is shown that with different chemi- stry, the profiles of normalized f differ significantly, with a shift in the location of peak and changes in the slopes of the curve. Such differ- 132 0 0.2 0.4 0.6 0.8 1 0 5 10 15 c∗ p (c ∗ ) (a) Probability Density Function 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p (c ∗ )d c∗ (b) Cumulative Distribution Function Figure 5.13: PDF and CDF in the whole domain. —-: DNS; − · −: mod- ified PDF model with laminar flame solution I;−−: modified PDF mo- del. 133 0 0.2 0.4 0.6 0.8 1 0 5 10 15 c∗ p( c∗ ) (a) x/L = 0.234 0 0.2 0.4 0.6 0.8 1 0 5 10 15 c∗ p( c∗ ) (b) x/L = 0.391 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 c∗ p( c∗ ) (c) x/L = 0.625 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 70 c∗ p( c∗ ) (d) x/L = 1 Figure 5.14: Probability Density Function on different planes.—-: DNS; −·−: modified PDFmodel with laminar flame solution I;−−: modified PDF model. 134 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) x/L = 0.234 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) x/L = 0.391 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) x/L = 0.625 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) x/L = 1 Figure 5.15: Cumulative Distribution Function on different planes. —-: DNS; − · −: modified PDF model with laminar flame solution I; −−: modified PDF model. 135 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 70 c∗ p( c∗ ) (a) c̄ = 0.089,c̄2 = 0.046 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 c∗ p( c∗ ) (b) c̄ = 0.413,c̄2 = 0.282 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14 16 18 c∗ p( c∗ ) (c) c̄ = 0.461,c̄2 = 0.311 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 c∗ p( c∗ ) (d) c̄ = 0.689,c̄2 = 0.505 Figure 5.16: Probability Density Function in different cubes. —-: DNS; −·−: modified PDFmodel with laminar flame solution I;−−: modified PDF model. 136 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (a) c̄ = 0.089,c̄2 = 0.046 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (b) c̄ = 0.413,c̄2 = 0.282 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (c) c̄ = 0.461,c̄2 = 0.311 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ ∫ p( c∗ )d c∗ (d) c̄ = 0.689,c̄2 = 0.505 Figure 5.17: Cumulative Distribution Function in different cubes. —-: DNS; − · −: modified PDF model with laminar flame solution I; −−: modified PDF model. 137 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c∗ m a x (f )/ f chem−1 chem−2 Figure 5.18: Comparison of laminar flames with different chemistry ences in f lead to the error in the PDF model when the ‘wrong’ chem- istry is used. It is also observed that the error introduced by different kinetics is less significant at relatively small values of the variance c′′2. Under such conditions, the modified PDF might predict a non-zero distribution over a small range of c∗ around the mean value c̄ to match c′′2 in the turbulent flame. Only a relatively small segment of the lam- inar flame solutions (corresponding to c∗ within this range) needs to be used. The error introduced by the laminar solution using a different chemical mechanism therefore might be less significant. 5.4 Conclusions Different forms for the presumed PDF of a product-based reaction prog- ress variable have been tested against DNS of turbulent premixed flames. On the whole, the β-PDF rarely captures the distribution of the vari- 138 able. A modification is suggested to extend the laminar flame based PDF proposed by Bray et al.. [132]. Comparison shows improvements in PDF and CDF predictions with the modified laminar flame PDF. The precision of the modified PDF model is affected by the chemical kinet- ics. The same mechanism used in the turbulent combustion simulation should be adopted to solve the prototype laminar flame that would be used to calculate the presumed PDF if this is to best match the subse- quently modelled turbulent flame. 139 Chapter 6 Closure of Reaction Rates in Turbulent Premixed Combustion with CSE 6.1 Introduction In previous studies, closure with CSE was obtained for non-premixed flames using the mixture fraction Z for conditioning. For premixed flames, chemical closure might be obtained by using the product based progress variable c for conditioning. In this chapter, the mean reaction rates in turbulent premixed flames are evaluated using the CSE model and the presumed PDFs proposed in Chapter 5. The DNS data from Grout’s work [133] is used to evaluate the precision of closure with CSE. 6.2 The CSE Formulation The formulation of CSE with c as the conditional variable is similar to the one with Z. The first conditional moment closure hypothesis of CMC is adopted, which evaluates chemical source terms using quanti- ties conditional on the progress variable c: ω̇|c∗ ≈ ω̇(T |c∗, Yi|c∗, ρ|c∗), (6.1) 140 where c∗ is the sample space for the progress variable. As discussed in previous chapters, instead of solving for the CMC transport equations for the conditional averages, CSE obtains condi- tional averages of different quantities by inverting integral equations, taking advantage of the statistical homogeneity of conditional aver- ages on pre-determined surfaces in the flow field. Similar to non-pre- mixed combustion simulation, ensembles of points with homogeneity may be selected in premixed flames, where the conditional means at different points are equal to those of the ensemble: f |c∗(x, t) ≈ 〈f |c∗〉(t;A). (6.2) Here, A is the ensemble of homogeneity and x is a point within that ensemble. The CSE inversion equations can be constructed in a way similar to the non-premixed cases if the probability density function of the new conditioning variable c is known. For a selected ensemble of statistical homogeneity, such integrals may be written for different points: f(x, t) ≈ ∫ 1 0 〈f |c∗〉(t;A)P(x, t; c∗)dc∗, (6.3) where P (x, t; c∗) is the local probability density function at point x. Assuming that the probability density function of c can be well- approximated by the PDF models described previously, Eq. 6.3 be- comes a relatively simple linear system that can be inverted so long as some kind of regularization is provided. Thus, the conditional av- erages of different scalars can be approximated. Chemical closure can then be achieved by invoking the CMC hypothesis, Eq. 6.1. The chem- ical source term in the transport equations of the unconditional means is closed by integrating over the sample space of c∗, using Eq. 6.3. 141 In the inversion process, a linear regularization method can be used to overcome ill-conditioning of the problem. In the simulation of real flows, it has been found that temporal continuity in the conditional av- erages can be assumed and conditional means obtained at the previous time step used for regularization [57]. In an a priori test, such informa- tion is unavailable. Alternately, twomethods are suggested for the CSE inversion computation: • Use the laminar flame solutions for regularization. The argument is that the premixed turbulent flame is extremely thin and simi- larity is expected between the mean flow and the laminar flame. • Solve the inversion equation iteratively from some arbitrary ini- tial condition and use the solutions of a previous iteration for reg- ularization. The computation stops when the solutions become converged orwhen somemaximumnumber of iterations are reach- ed. The idea is to mimic the time-advancement in a CSE-based simulation with this iteration process. 6.3 Results and Discussion : RANS tests 6.3.1 Conditional Averages: Mechanism I The mean reaction rates are closed using the presumed PDF and the CSE method. Tests are performed using the turbulent flame solutions I, with the two-step mechanism defined previously by Eqs. 5.33 to 5.35. In the RANS tests, the mean reaction rates on parallel planes or- thogonal to the direction of the main flow are extracted from the DNS data base. The reaction rates are also evaluatedwith scalars conditional 142 on c, using the presumed PDFs and the CSE formulation described pre- viously. An assumption of statistical homogeneity in the conditional aver- ages over the flow field is adopted to construct the linear equation sys- tem for inversion. As discussed previously, a suitable regularization method is important for evaluating conditional means in CSE. The ef- fectiveness of the two proposedmethods are compared, using the same presumed PDFs in the CSE inversion computation. The results presented in Fig. 6.1 are obtained using the modified laminar flame PDF. The method of regularizing with the laminar flame solution yields good estimates of the conditional means for all reactive scalars. The iteration method which started from an initial condition of zero for each scalar either fails to obtain converged solutions within 2000 iterations, or yields solutions less desirable than those of the first method. The most significant errors are observed in the iteration solu- tions of 〈YI |c∗〉. Since the second reaction in the two-step mechanism is highly sensitive to the change in intermediates, inaccuracy here causes even bigger errors in the reaction rate calculation. The inversion methods are also tested with the β-PDF to explore the sensitivity of CSE method to the presumed shape of the PDF. The results are presented in Fig. 6.2. The solutions with the iteration me- thod are extremely poor, and are therefore not presented. Regulariza- tion with the laminar flame solution yields an acceptable approxima- tion for the conditional means, but the deviation from the DNS is larger than that observed in the solutions with the modified PDF. The laminar flame solutions used for regularization are also plotted 143 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.01 0.02 0.03 0.04 0.05 c∗ < Y F |c∗ > (a) 〈YF |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c∗ < Y O |c∗ > (b) 〈YO|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 c∗ < Y I |c∗ > (c) 〈YI |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 500 1000 1500 2000 2500 3000 c∗ < T |c∗ > (d) 〈T |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 c∗ < ρ |c∗ > (e) 〈ρ|c∗〉 Figure 6.1: Conditional averages calculated in RANS paradigm, ob- tained by using modified laminar flame PDF and different regular- ization methods. —-: DNS; −−: regularization with laminar flame so- lution (CSE-1); − · −: regularization with iterations up to 2000 times (CSE-2); · · · · ·: laminar flame solution. 144 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.01 0.02 0.03 0.04 0.05 0.06 c∗ < Y F |c∗ > (a) 〈YF |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c∗ < Y O |c∗ > (b) 〈YO|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 c∗ < Y I |c∗ > (c) 〈YI |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 500 1000 1500 2000 2500 3000 c∗ < T |c∗ > (d) 〈T |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 c∗ < ρ |c∗ > (e) 〈ρ|c∗〉 Figure 6.2: Conditional averages calculated in RANS paradigm, ob- tained by using β-PDF and regularization with laminar flame solution. —-: DNS; −−: CSE. 145 in Fig. 6.1, which happen to coincides with the DNS conditional means in this test. This seems to cast questions on the effectiveness of the CSE inversion method, which will be clarified in the following discussion and later in another series of a priori tests. The problems with the iteration methods might be attributed to the lack of a stabilization mechanism that would be present in CSE simula- tions. In previous RANS and LES simulations using CSE, the changes in conditional means affect the reactive scalars through the chemical source terms, and changes in reactive scalars are directly reflected in the conditional means. The interactions between unconditional and conditional solutions appear to stabilize the CSE inversion process. Such a mechanism does not exist in the a priori tests, which may be leading to the deviation and fluctuations observed in the results from the iteration method. Regularization with laminar flame solutions works relatively well with both the β-PDF and the modified PDF. This observation indicates that a method similar to LFD [61] might work well in premixed flame simulations, in which a linear combination of laminar flame solutions at various strain rates might provide a better approximation to the flame. 6.3.2 Reaction Rate Predictions: Mechanism I The chemical reaction rates are modelled and compared with the DNS. Figs. 6.3 and 6.4 present the conditionally averaged reaction rates eval- uated using the modified laminar flame PDF and the β-PDF. The seem- ingly small errors in the solutions of conditionally averaged reactive scalars are amplified in the calculation of reaction rates. The CSE mo- 146 del with the modified PDF over-predicts the reaction rate 〈ω̇|c∗〉 by up to 20%, while the model with β-PDF predictions are obviously unac- ceptable. The chemical reactions rates on different planes are closed by in- tegration. Only predictions with the modified PDF are presented. As shown in Fig. 6.5, overall, the model predictions agree reasonably well with the DNS when the modified PDF is used. The peak reaction rates that occur in the middle of the flow field are captured well for both reactions. 6.3.3 Comparison with DNS: Mechanism II The a priori test described in the previous section is repeated for a sec- ond set of data, turbulent flame solution II. The conditional means of reactive scalars and mean reaction rates are evaluated using the CSE and presumed PDF method; the results of RANS tests are shown in Figs. 6.6 to 6.8. Similar observations are found in this series of tests. The results generally agreewith the DNS data, indicating the effectiveness of chem- ical closure with the CSE and the proposed PDF model in premixed flames. It is noteworthy that in this test, the conditional means from the DNS are different from the laminar flame solutions used for reg- ularization (method CSE-1); the conditional means obtained using the CSE inversion method are closer to the DNS than the a priori informa- tion from the laminar solutions. This clarifies the question raised in the first series of tests: good estimation of conditional means can be obtained using the CSE method, even when the a priori information in- troduced for regularization contains error. Altogether, the predictions 147 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 c∗ < ω̇ i |c∗ > (a) 〈ω̇i|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 < ω̇ ii |c∗ > c∗ (b) 〈ω̇ii|c∗〉 Figure 6.3: Conditionally averaged reaction rates 〈ω̇|c∗〉 calculated in RANS paradigm. Conditional means obtained by using the modified PDF and regularization with the laminar flame solution. —-: DNS;−−: CSE. 148 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 c∗ < ω̇ i |c∗ > (a) 〈ω̇i|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5 3 c∗ < ω̇ ii |c∗ > (b) 〈ω̇ii|c∗〉 Figure 6.4: Conditionally averaged reaction rates 〈ω̇|c∗〉 calculated in RANS paradigm. Conditional means obtained by using the β-PDF and regularization with the laminar flame solution. —-: DNS; −−: CSE. 149 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x/L ω̇ i (a) 〈ω̇i〉 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 ω̇ ii x/L (b) 〈ω̇ii〉 Figure 6.5: Reaction rates 〈ω̇〉 on planes in a premixed turbulent flame. Conditional means obtained by using the modified PDF and regular- ization with the laminar flame solution. —-: DNS; −−: CSE. 150 of the CSEmodel capture the trend observed in the DNS data, whereas the reaction rates tend to be over-predicted, as in the first test. 6.4 Results and Discussion : LES tests 6.4.1 Comparison with DNS: Mechanism I A priori tests are performed in the LES paradigm using the turbulent flame solution I. In the LES tests, the mean reaction rates in 4 × 32 × 32 cubes are estimated with CSE and the modified PDF. Again, the predictions are compared to the data extracted from the DNS. Again, the assumption of statistical homogeneity in the conditional averages throughout the flow field is adopted and the linear regularization with an unstrained laminar flame solution is used to solve for the conditional means. The solutions are presented in Fig. 6.9. As was seen in the RANS test, the predicted conditional means are close to the DNS; the laminar flame solution is barely distinguishable from the DNS. The chemical reaction rates are then modelled and compared with the DNS results. In Fig. 6.10, the conditionally averaged reaction rates are compared to the DNS. As in the RANS test, the reaction rates are over-predicted because of the errors in the conditional means of reac- tive scalars. The filtered chemical reactions rates in different cubes are evalu- ated by integration and plotted against the filtered DNS reaction rates. A constrained least squares analysis is performed to evaluate the pre- cision of the closure. As shown in Fig. 6.11, the modelled reaction rates are close to the DNS results. The slopes of the fitted lines for the two re- 151 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 c∗ < Y F |c∗ > (a) 〈YF |c∗〉 0 0.2 0.4 0.6 0.8 1 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 c∗ < Y O |c∗ > (b) 〈YO|c∗〉 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 c∗ < Y I |c∗ > (c) 〈YI |c∗〉 0 0.2 0.4 0.6 0.8 1 600 800 1000 1200 1400 1600 1800 2000 c∗ < T |c∗ > (d) 〈T |c∗〉 0 0.2 0.4 0.6 0.8 1 0.2 0.25 0.3 0.35 0.4 0.45 0.5 c∗ < ρ |c∗ > (e) 〈ρ|c∗〉 Figure 6.6: Conditional averages calculated in RANS paradigm, ob- tained by using modified laminar flame PDF and different regular- ization methods. —-: DNS; −−: regularization with laminar flame so- lution (CSE-1); − · −: regularization with iterations up to 2000 times (CSE-2); · · · · ·: laminar flame solution. 152 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 c∗ < ω̇ i |c∗ > (a) 〈ω̇i|c∗〉 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 c∗ < ω̇ ii |c∗ > (b) 〈ω̇ii|c∗〉 Figure 6.7: Conditionally averaged reaction rates 〈ω̇|c∗〉 calculated in RANS paradigm. Conditional means obtained by using the modified PDF and regularization with the laminar flame solution. —-: DNS;−−: CSE. 153 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/L ω̇ i (a) 〈ω̇i〉 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x/L ω̇ ii (b) 〈ω̇ii〉 Figure 6.8: Reaction rates 〈ω̇〉 on planes in a premixed turbulent flame. Conditional means obtained by using the modified PDF and regular- ization with the laminar flame solution. —-: DNS; −−: CSE. 154 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 c∗ < Y F |c∗ > (a) 〈YF |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c∗ < Y O |c∗ > (b) 〈YO|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 c∗ < Y I |c∗ > (c) 〈YI |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 500 1000 1500 2000 2500 3000 c∗ < T |c∗ > (d) 〈T |c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 c∗ < ρ |c∗ > (e) 〈ρ|c∗〉 Figure 6.9: Conditional averages calculated in LES paradigm, obtained by using the modified PDF and regularization with laminar flame so- lution. —-: DNS; −−: CSE; · · · · ·: laminar flame solution. 155 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 c∗ < ω̇ i |c∗ > (a) 〈ω̇i|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 c∗ < ω̇ ii |c∗ > (b) 〈ω̇ii|c∗〉 Figure 6.10: Conditionally averaged reaction rates 〈ω̇|c∗〉, calculated in LES paradigm, conditional means obtained by using the modified PDF and regularization with laminar flame solution. —-: DNS; −−: CSE. 156 actions are 1.222 and 1.649; the Pearson product moment coefficients of correlation are 0.9508 and 0.6924 respectively. The results seem promis- ing, although it is somewhat disappointing that there might be a posi- tive bias in the CSE prediction of reaction rates. Optimization of the PDF table might further improve the precision of closure. For illustration, a slightly denser grid in conditional vari- able space with 46 points was used. The same procedure of generating a PDF table, inverting for conditional means, invoking first conditional moment closure and integrating to obtain filtered reaction rates is re- peated. As shown in Fig. 6.12, the slopes of the fitted lines are reduced to 1.135 and 1.491, and the coefficients of correlation are increased to 0.9614 and 0.7167, indicating enhanced precision of closure. For further analysis, the mean reaction rates are calculated using the conditional means extracted from the DNS and the modified PDF. The results are plotted against the DNS reaction rates in Fig. 6.13. The slopes of the fitted lines of these data points are 0.9710 and 0.9684; the coefficients of correlation are 0.9602 and 0.9590 respectively. This analysis shows the error observed in Fig. 6.11 is mainly attributable to errors in the conditional means. As discussed previously, in real CSE modeling, the interaction between conditional means and the reactive scalars is likely to yield better solutions in the CSE inversion, which might improve the precision of closure as well. Therefore, the results shown in Fig. 6.13 give more confidence in the closure method with CSE and the modified PDF. 157 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 1.2 ω̇DNS ω̇ i (a) 〈ω̇i〉 −0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ω̇DNS ω̇ ii (b) 〈ω̇ii〉 Figure 6.11: Reaction rates 〈ω̇〉 calculated in LES paradigm. —-: ω̇ = ω̇DNS ; −−: ω̇i = 1.222 ∗ ω̇DNS,i and ω̇ii = 1.649 ∗ ω̇DNS,ii. 158 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 ω̇DNS ω̇ i (a) 〈ω̇i〉 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ω̇DNS ω̇ ii (b) 〈ω̇ii〉 Figure 6.12: Reaction rates 〈ω̇〉 calculated in LES paradigm, with denser grid in c-space. —-: ω̇ = ω̇DNS ; −−: ω̇i = 1.1345 ∗ ω̇DNS,i and ω̇ii = 1.4912 ∗ ω̇DNS,ii. 159 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ω̇DNS ω̇ i (a) 〈ω̇i〉 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ω̇DNS ω̇ ii (b) 〈ω̇ii〉 Figure 6.13: Reaction rates 〈ω̇〉 calculated in LES paradigm, evaluated by using the modified PDF and conditional mean reaction rates ex- tracted directly from the DNS. —-: ω̇ = ω̇DNS ; −−: ω̇i = 0.9710 ∗ ω̇DNS,i and ω̇ii = 0.9684 ∗ ω̇DNS,ii. 160 6.4.2 Comparison with DNS: Mechanism II The LES tests described in the previous section is repeated for the tur- bulent flame solution II. The conditional means of reactive scalars and unconditional means of reaction rates are compared to the DNS data as shown in Figs. 6.14 to 6.16. As observed in the RANS test II, the conditional means of the DNS are different from the laminar flame solutions; despite this imprecise a priori information, the conditional means obtained with the CSE inver- sion method are quite close to the DNS results. As observed in the LES test I, the predicted reaction rates are consistent with the DNS data. The slopes of the fitted lines in Fig. 6.16 are 1.0454 and 1.0799 respec- tively. The tendency of overprediction is still exhibited in this series of test; however, it is less significant than that of LES test I. 6.5 Conclusions Aproduct-based progress variable is introduced as a conditioning vari- able for the CSE model. The mean reaction rates in premixed turbulent flames are modelled using the CSE method and the proposed PDFmo- del. The closure method is tested in both RANS and LES paradigms. In the a priori tests, the conditional means of reactive scalars are ob- tained using the CSE inversion method. Good estimates of conditional means are obtained when the inversion equations are regularized with a laminar flame solution. The a priori tests further validate themodified version of the laminar flame based PDF model and the closure method of CSE. The predictions of the mean reaction rates agree reasonably well with the DNS in both RANS and LES tests. The application of the 161 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 c∗ < Y F |c∗ > (a) 〈YF |c∗〉 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 0.25 c∗ < Y O |c∗ > (b) 〈YO|c∗〉 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 c∗ < Y I |c∗ > (c) 〈YI |c∗〉 0 0.2 0.4 0.6 0.8 1 800 1000 1200 1400 1600 1800 2000 c∗ < T |c∗ > (d) 〈T |c∗〉 0 0.2 0.4 0.6 0.8 1 0.2 0.25 0.3 0.35 0.4 0.45 0.5 c∗ < ρ |c∗ > (e) 〈ρ|c∗〉 Figure 6.14: Conditional averages calculated in LES paradigm, ob- tained by using the modified PDF and regularization with laminar flame solution. —-: DNS; −−: CSE; · · · · ·: laminar flame solution. 162 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 c∗ < ω̇ i |c∗ > (a) 〈ω̇i|c∗〉 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 c∗ < ω̇ ii |c∗ > (b) 〈ω̇ii|c∗〉 Figure 6.15: Conditionally averaged reaction rates 〈ω̇|c∗〉, calculated in LES paradigm, conditional means obtained by using the modified PDF and regularization with laminar flame solution. —-: DNS, −−: CSE. 163 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ω̇DNS ω̇ i (a) 〈ω̇i〉 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ωDNS ω̇ ii (b) 〈ω̇ii〉 Figure 6.16: Reaction rates 〈ω̇〉 calculated in LES paradigm. —-: ω̇ = ω̇DNS , −−: ω̇i = 1.0454 ∗ ω̇DNS,i and ω̇ii = 1.0799 ∗ ω̇DNS,ii. 164 CSE closure method with the proposed PDF in premixed flames needs further tests in combustion simulations. 165 Chapter 7 Conclusions and Future Work 7.1 Conclusions In this thesis work, the Conditional Source-term Estimation method has been applied to study turbulent reacting flows. Firstly, CSE was combined with Laminar Flamelet Decomposition and then with Tra- jectory Generated Low Dimensional Manifold method to study non- premixed autoigniting jets of methane and methane mixed with ad- ditives. Both methods were implemented using OpenFOAM and ap- plied to study jet flameswith the same initial and boundary conditions, using the same grid and numerical scheme for comparison. The pre- dictions were compared to shock-tube experiments. The study shows that the CSE-TGLDMmethod is superior to the CSE-LFD method; this might be attributed to the wider range of solutions of detailed chemis- try included in the TGLDMmanifolds, which also allows CSE-TGLDM to model transition from ignition to combustion, as opposed to the CSE-LFD method which only works up to the moment of ignition. The effects of different injection parameters on the ignition delay time was discussed based on the predictions of the CSE-TGLDM me- thod. The fuel composition of the jet affects the ignition delay by chang- ing the chemical kinetics and fluid dynamics of the flame. Ethane ac- 166 celerates the autoignition of methane and nitrogen delays autoignition; however, the effects of modest quantities of these additives become negligible at high pre-combustion temperatures. Increasing injection pressure ratio reduces ignition delay and the effects become negligi- ble at high pressure ratios. The injection duration is irrelevant when ignition starts before the end of injection. The predictions of ignition delay time are generally consistent with the experimental data, which indicates the effectiveness of the CSE methods and the chemicalmechanism used in the jet simulations. How- ever, the scatter observed in the experimental data could not be cap- tured in the simulation. The scatter in ignition delaymight be attributed to the different realizations of chemical reaction paths, as well as the fluctuations in turbulent jets, which cannot be captured by continuum manifold methods and Reynolds Averaged Navier-Stokes simulation. To simulate these effects, a TGLDMmethod combining continuum and stochastic particle model generated trajectories might be used to simu- late the randomness of chemical reaction paths; Large Eddy Simulation is suggested to capture the fluctuations in turbulent flows. The recommended CSE-TGLDMmethod was applied to obtain clo- sure in LES of a piloted, non-premixedmethane jet flame, Sandia Flame D. In the thesis work, the adiabatic CSE-TGLDMmethodwas extended to model radiation in the flame. It is found that the CSE method can be used to model the turbulence/radiation interaction in a similar fash- ion to how it models the turbulence/chemistry interaction; the ma- jor sources of radiation in Sandia D flame are hot CO2 and H2O; the strongest radiation in the flame is around the stoichiometric flame height; 167 the radiation heat loss of the flame is insufficient to make a big differ- ence in the predictions of temperature and major species; however, the predictions of NO are improved by using the extended CSE-TGLDM model. The application of the CSE method in the study of premixed turbu- lent flames was explored. A conditioning variable based on the prod- uct mass fraction was proposed. A priori tests show that the presumed PDF of this progress variable may be well approximated with an ex- tended version of a laminar flame-based model. The precision of the PDF is sensitive to the chemical mechanism; it is suggested that the same mechanism be used for the solutions of the prototype laminar flame and the turbulent flame. Using the conditioning variable and its presumed PDF model, the CSE method was adapted to evaluate conditional averages and close chemical source terms in premixed flames. A priori tests were perform- ed in both the RANS and LES paradigms. The results of the tests are promising. However, the effectiveness and precision of chemical clo- sure with the CSE method in turbulent premixed flames need further tests in flames of different configurations and ultimately in turbulent premixed combustion simulations. 7.2 Future Work In the study of unsteady autoigniting jets, the CSE-LFD method ap- pears to be inferior to the CSE-TGLDM method. However, it is possi- ble to improve CSE-LFD predictions by enhancing the flamelet library to include solutions of extinction and re-ignition phenomena. More- 168 over, the a priori tests indicate that the LFD method might work well for premixed flames. The potential of the CSE-LFD method is worth exploring further. It is also found that the current RANS simulation is incapable of capturing the distribution of ignition delay caused by the fluctuations of chemical reactions and turbulence. LES of autoigniting jets with a TGLDM method combining the continuum approach and stochastic particle model to create trajectories is suggested to model these fluctu- ations. In the study of steady jet flame Sandia flame D, the CSE-TGLDM method has been extended to model radiation. The storage require- ments and retrieval of multiple manifolds turn out to be demanding. These requirements are still manageable in the current work; however, it will be necessary to use more manifolds to model strongly radiating flames, which poses a challenge for implementation. Another issue is the precision of the manifold in the low temperature range where the physical time scales are comparable to the time scales of the manifold. In future work, these issues should be addressed. In the study of premixed turbulent flames, the CSE method seems to be a promising tool to obtain chemical closure. In future study, other candidates of the presumed probability density function for the pro- posed conditioning variable might be studied; improvements to the extended laminar flame based PDF might be made by filtering. Ul- timately, the CSE with presumed PDF method should be applied to simulate turbulent premixed flames. 169 Bibliography [1] Information available at: http://www.centreforenergy.com/ET.asp, accessed 12 Dec. 2007. [2] Code of federal regulations (cfr) 40 parts 69, 80, 86: Control of air pollution from new motor vehicles: heavy duty engine and ve- hicle standards and highway diesel fuel sulfur control require- ments; final rule, US Federal Register 66(12) (2001) 5002–5050. [3] On road vehicle and engine emissions regulations, Canada Gazette Part II 137(1) (2003) 21–27. [4] G. P. McTaggart-Cowan, C. C. O. Reynolds, W. K. Bushe, Natu- ral gas fuelling for heavy-duty on-road use: Current trends and future direction, International Journal of Environmental Studies 63(4) (2006) 421–440. [5] Information available at: http://www.cngva.org/index.htm, accessed 12 Dec. 2007. [6] Information available at: http://www.ngvontario.com, accessed 12 Dec. 2007. [7] W. K. Bushe, H. Steiner, Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows, Physics of Fluids 11 (1999) 1896–1906. [8] D. Veynante, L. Vervisch, Turbulent combustion modeling, Progress in Energy and Combustion Science 28 (2002) 193–266. [9] P. Moin, K. Mahesh, Direct numerical simulation: a tool in tur- bulence research, Annual Review of Fluid Mechanics 30 (1998) 539–578. [10] V. Eswaran, S. B. Pope, An examination of forcing in direct numerical simulations of turbulence, Computational Fluids 16 (1988) 257–278. [11] J. Kim, P.Moin, R. Moser, Turbulence statistics in fully developed channel flow, Journal of Fluid Mechanics 177 (1987) 133–166. 170 [12] P. R. Spalart, Direct numerical study of leading-edge contamina- tion, in: Fluid Dynamics of Three-Dimensional Turbulent Shear Flows and Transition, 1989. [13] D. C. Wilcox, Turbulence modelling for CFD, DCW Industries, La Canada, 1998. [14] S. B. Pope, Turbulent flows, Cambridge University Press, 2001. [15] D. Veynante, T. Poinsot, Reynolds averaged and large eddy simulation modeling for turbulent combustion, in: O. Métais, J. Ferziger (Eds.), New tools in turbulence modelling, Berlin: Les Editions de Physique, Springer, 1997, pp. 105–140. [16] J. Smagorinsky, General circulation experiments with the primi- tive equations, Monthly Weather Review 91 (3) (1963) 99–164. [17] M. Germano, U. Piomelli, P. Moin, W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Physics of Fluids 3 (1991) 1760–1765. [18] P. Moin, K. Sqires, W. Cabot, S. Lee, A dynamic subgrid-scale model for compressible turbulence and scalar transport, Physics of Fluids 3 (1991) 2746–2757. [19] A. E. Tejada-Martı́z, K. E. Jansen, A dynamic smagorinskymodel with dynamic determination of the filter width ratio, Physics of Fluids 16 (2004) 2514–2528. [20] H. Pitsch, Large-eddy simulation of turbulent combustion, An- nual Review of Fluid Mechanics 38 (2006) 453–482. [21] N. Peters, Turbulent combustion, 1st Edition, Cambridge Uni- versity Press, 2000. [22] R. W. Bilger, Turbulent jet diffusion flames, Progress of Energy Combustion and Science (1976) 87–109. [23] P. A. Libby, F. A. Williams, Turbulent reacting flows, Academic press London, 1994. [24] R. W. Bilger, Turbulent flows with non-premixed reactants, Top- ics Applied Physics 44 (1980) 65–113. 171 [25] S. Burke, T. Schumann, Diffusion flames, Industrial and Engi- neering Chemistry 20 (1928) 998–1004. [26] A. W. Cook, J. J. Riley, A subgrid model for equilibrium chemis- try in turbulent flows, Physics of Fluids 6 (1994) 2868–2870. [27] J. Coupland, C. H. Priddin, Modelling of flow and combustion in a production gas turbine combustor, Turbulent Shear Flows 5 (1987) 310–323. [28] N. Peters, Laminar diffusion flamelet models in non-premixed turbulent combustion, Progress in Energy and Combustion Sci- ence 10 (3) (1984) 319–339. [29] H. Pitsch, N. Peters, A consistent flamelet formulation for non- premixed combustion considering differential diffusion effects, Combustion and Flame 114 (1998) 26–40. [30] N. Peters, Comment and reply on the ”Assessment of combus- tion and submodels for turbulent nonpremixed hydrocarbon flames”, Combustion and Flame 116 (1999) 675–676. [31] N. Swaminathan, R. Bilger, Assessment of combustion and sub- models for turbulent nonpremixed hydrocarbon flames, Com- bustion and Flame 116 (1999) 519–545. [32] N. Swaminathan, Flamelet regime in non-premixed combustion, Combustion and Flame 129 (2002) 217–219. [33] F. Mauss, D. Keller, N. Peters, A Lagrangian simulation of flamelet extinction and re-ignition in turbulent jet diffusion flames, Proceedings of the Combustion Institute 23 (1990) 693– 698. [34] H. Barths, N. Peters, N. Brehm, et al., Simulation of pollutant for- mation in a gas turbine combustor using undteady flamelets, in: Proceedings of the Combustion Institute, The Combustion Insti- tute, 1998, pp. 1841–1847. [35] H. Barths, C. Hasse, Simulation of combustion in direct injection diesel engines using an eulerian particle flamelet model, Pro- ceedings of the Combustion Institute 28 (2000) 1161–1168. 172 [36] H. Pitsch, Unsteady flamelet modeling of differential diffusion in turbulent jet diffusion flames, Combustion and Flame 123 (3) (2000) 358–374. [37] H. Pitsch, H. Steiner, Large-eddy simulation of a turbulent pi- loted methane/air diffusion flame (Sandia flame D), Physics of Fluids 12 (2000) 2541–2554. [38] S. P. Nandula, T. M. Brown, R. W. Pitz, Measurements of scalar dissipation in the reaction zones of turbulent nonpremixed h2-air flames., Combustion and Flame 99 (1994) 775–783. [39] R. W. Bilger, Conditional moment closure for turbulent reacting flows, Physics of Fluids 5 (1993) 436–444. [40] A. Y. Klimenko, Multicomponent diffusion of various admix- tures in turbulent flow, Fluid Dynamics 25 (3) (1990) 327–334. [41] Information available at: http://www.ca.sandia.gov/TNF/DataArch/FlameD.html, ac- cessed 12 Dec. 2007. [42] A. Klimenko, R. Bilger, Conditional moment closure for tur- bulent combustion, Progress in energy and combustion science 25 (6) (1999) 595–687. [43] S. Kim, K. Y. Huh, Use of the conditional moment closure model to predict no formation in a turbulent CH4/H2 flame over a bluff body, Combustion and Flame 130 (2002) 94–111. [44] S. H. Kim, K. Y. Huh, R. A. Fraser, Modeling autoignition of a turbulent methane jet by the conditional moment closure model, Proceedings of the Combustion Institute 28 (2000) 185–191. [45] C. Devaud, K. N. C. Bray, Assessment of the applicability of con- ditional moment closure model to a lifted turbulent flame: first order model, Combustion and Flame 132 (2003) 102–114. [46] S. H. Kim, K. Y. Huh, R. W. Bilger, Second-order conditional moment closure modeling of local extinction and reignition in turbulent non-premixed hydrocarbon flames, Proceedings of the Combustion Institute 29 (2002) 2131–2137. 173 [47] M. Roomina, R. Bilger, Conditional moment closure modeling of turbulent methanol jet flames, Combustion Theory and Mod- elling 3 (1999) 698–708. [48] M. Roomina, R. W. Bilger, Conditional moment closure (CMC) prediction of a turbulent methane-air flames, Combustion and Flame 125 (2001) 1176–1195. [49] M. Fairweather, R. M. Wooley, First-order conditional moment closure modeling of turbulent, non-premixed methane flames, Combustion and Flame 138 (2004) 3–19. [50] S. Kim, K. Y. Huh, L. Tao, Application of the elliptic con- ditional moment closure model to a two-dimensional non- premixed methanol bluff-body flame, Combustion and Flame 120 (75) (2000) 75–90. [51] N. Swaminathan, R. W. Bilger, Conditional variance equation and its analysis, Proceedings of the Combustion Institute (1998) 1191–1198. [52] E. Mastorakos, R. W. Bilger, Second-order conditional moment closure for the autoignition of turbulent flows, Physics of Fluids 10 (1998) 1246–1248. [53] C. M. Cha, G. Kosaly, H. Pitsch, Modeling extinction and reig- nition in turbulent nonpremixed combustion using a doubly- conditional moment closure approach, Physics of Fluids 13 (12) (2001) 3824–3834. [54] A. Kronenburg, Double conditioning of reactive scalar transport equations in turbulent nonpremixed flames, Physics of Fluids 16(7) (2004) 2640–2648. [55] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Nu- merical Recipes in Fortran, Cambridge University Press, Cam- bridge, UK, 1992. [56] 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. 174 [57] R. Grout, W. K. Bushe, C. Blair, Predicting the ignition delay of turbulent methane jets using conditional source-term estimation, To appear in Combustion Theory and Modelling. [58] R. Grout, Combustion modelling with conditional source-term estimation and laminar flamelet decomposition, Master’s thesis, University of British Columbia (2004). [59] M. Wang, J. Huang, W. K. Bushe, Simulation of a turbulent non- premixed flame using conditional source-term estimation with trajectory generated low-dimensional manifold, in: Thirty-First Symposium (International) on Combustion, The Combustion In- stitute, 2006, pp. 1701–1709. [60] J. Huang, W. K. Bushe, Simulation of an igniting methane jet us- ing conditional source-term estimation with a trajectory gener- ated low-dimensional manifold, To appear in Combustion The- ory and Modelling. [61] W. K. Bushe, H. Steiner, Laminar flamelet decomposition for con- ditional source-term estimation, Physics of Fluids 15 (2003) 1564– 1575. [62] M. Wang, Combustion modeling using conditional source-term estimation with flamelet decompositon and low dimensional manifolds, Ph.D. thesis, University of British Columbia (2006). [63] M.Wang,W. K. Bushe, Large eddy simulation of a turbulent non- premixed flame using conditional source-term estimation with laminar flamelet decomposition, submitted to Physics of Fluids. [64] R. Borghi, Turbulent combustion modelling, Progress in Energy and Combustion Science 14 (1988) 245–292. [65] D. B. Spalding, Mixing and chemical reaction in steady con- fined turbulent flames, Thirteenth Symposium (International) on Combusion, 1971, pp. 649–657. [66] B. F.Magnussen, B. H. Hjertager, Onmathematical models of tur- bulent combustion with special emphasis on soot formation and combustion, in: Sixteenth Symposium (International) on Com- bustion, The Combustion Institute, 1977, pp. 719–729. 175 [67] J. B. Moss, K. N. C. Bray, A unified statistical model of the pre- mixed turbulent flame, Acta Astronaut 4 (1977) 291–319. [68] K. N. Bray, P. A. Libby, Recent developments in the BML mo- del of premixed turbulent combustion, Academic Press, London, 1994. [69] T. Mantel, R. Borghi, A new model of premixed wrinkled flame propagation based on a scalar dissipation equation, Combustion and Flame 96(4) (1994) 443–457. [70] A. Trouvé, T. J. Poinsot, The evolution equation for the flame sur- face density in turbulent premixed combustion, Journal of Fluid Mechanics 278 (1994) 1–31. [71] F. A.Williams, Turbulent Combustion, SIAM, Philadelphia, 1985, the mathematics of combustion. [72] P. Nilsson, X. S. Bai, Effect of flame stretch and wrinkling on co formation in turbulent premixed combustion, 29th Symposium (International) on Combusion, 2002, pp. 1871–1879. [73] M. Z. Bodenstein, Eine theorie der photochemischen reaktions- geschwindigkeiten, Physics of Chemistry 85 (1913) 329–397. [74] N. Peters, R. Kee, The computation of stretched laminar methane-air diffusion flames using a reduced four-step mech- anism, Combustion and Flame 68 (1) (1987) 17–30. [75] J. Warnatz, Concentration, pressure, and temperature- dependence of the flame velocity in hydrogen-oxygen-nitrogen mixtures, Combustion Science and Technology 26 (5-6) (1981) 203–213. [76] M. D. Smooke, Reduced kinetic mechanisms and asymptotic ap- proximations for methane-air flames, in: M. D. Smooke (Ed.), Lecture Notes in Physics, Springer-Verlag, 1991, pp. 1–28. [77] N. Peters, B. Rogg (Eds.), Reduced kinetic mechanisms for appli- cations in combustion systems, Springer-Verlag, 1993. 176 [78] U. Maas, S. B. Pope, Simplifying chemical kinetics: Intrinsic low- dimensional manifolds in composition space, Combustion and Flame 88 (1992) 239–246. [79] Z. Ren, S. B. Pope, A. Vladimirsky, J. M. Guckenheimer, Applica- tion of the ICE-PICmethod for the dimension reduction of chem- ical kinetics, Journal of Chemical Physics 124 (2006) 114111. [80] J. A. van Oijen, L. P. H. de Goey, Modelling of premixed coun- terflow flames using the flamelet-generated manifold method, Combustion Theory and Modelling 6 (2002) 463–478. [81] J. A. van Oijen, L. P. H. de Goey, Modeling of premixed laminar flame using flame-generated low-dimensional manifolds, Seven- teenth International Colloquium on the Dynamics of Explosions and Reactive Systems (1999) 25–30. [82] U. Maas, S. B. Pope, Implementation of simplified chemical ki- netics based on intrinsic low-dimensional manifolds, Twenty- Fourth Symposium (International) on Combustion (1992) 103– 112. [83] J. Huang, Natural gas combustion under engine-relevant condi- tions, Ph.D. thesis, University of British Columbia (2006). [84] J. Nafe, U. Maas, A general algorithm for improving ildms, Com- bustion Theory and Modelling 6 (2002) 697–709. [85] O. Gicquel, N. Darabiha, D. Thevenin, Laminar premixed hydro- gen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion, in: Twenty-Eighth Sympo- sium (International) on Combustion, 2000, pp. 1901–1908. [86] H. Bongers, J. A. V. Oijen, L. P. H. de Goey, Intrinsic low- dimensional manifold method extended with diffusion, Pro- ceedings of the Combustion Institute 29 (2002) 1371–1378. [87] S. Pope, Computationally efficient implementation of com- bustion chemistry using in situ adaptive tabulation, Combus- tion Theory and Modelling 1 (1997) 41–63. 177 [88] S. B. Pope, U. Maas, Simplifying chemical kinetics: Trajec- tory generated low-dimensional manifolds, Mechanical and Aerospace Engineering Report: FDA 93-11, Cornell University. [89] S. Gordon, B. J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311. [90] 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. [91] S. B. Pope, Ten questions concerning the large-eddy simulation of turbulent flows, New Journal of Physics 6 (2004) 35–59. [92] R. Renka, Tripack: A constrained two-dimensional delaunay tri- angulation package, ACM Transactions on Mathematical soft- ware 22 (1996) 1–8. [93] M. Wang, A. Frisque, W. K. Bushe, Trajectory-generated low- dimensional manifolds generated using the stochastic particle method, to appear in Combustion Theory and Modelling. [94] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, Journal of Physics of Chemistry 81 (1977) 2340–2361. [95] D. T. Gillespie, A rigorous derivation of the chemical master equation, Physica A 188 (1992) 404–425. [96] A. Frisque, J. Schnakenberg, J. Huang, W. K. Bushe, Stochastic simulation of the variations in the autoignition delay time of premixed methane and air, Combustion Theory and Modelling 10 (2) (2006) 241–256. [97] Information available at: http://www.opencfd.co.uk/openfoam/, accessed 12 Dec. 2007. [98] G. D. Sullivan, J. Huang, T. X. Wang, W. K. Bushe, S. N. Rogak, Emissions variability in gaseous fuel direct injection compres- sion ignition combustion, Society of Automotive Engineers 114 (2005) 2005–01–0917. 178 [99] P. Ouellette, Direct injection of natural gas for diesel engine fuel- ing, Ph.D. thesis, University of British Columbia, Canada (1996). [100] F. P. Ricou, D. B.Spaiding, Measurement of entrainment by axial symmetric turbulent jet, Journal of FluidMechanics 11 (1961) 21– 32. [101] Information available at: http://www.fluent.com/, accessed 12 Dec. 2007. [102] P. Ouellette, P. G. Hill, Turbulent transient gas injections, Journal of Fluids Engineering 122 (2000) 743–753. [103] N. W. and M. H. Davy and W. K. Bushe, Ignition and emission characteristics of methane direct injection compression ignition combustion, 2006 CICS Spring Technical Meeting in Waterloo, Canada, 2006. [104] G. P. McTaggart-Cowan, N. Wu, S. N. Rogak, M. H. Davy, W. K. Bushe, Effects of fuel dilution on high-pressure non-premixed natural gas combustion, submitted to Combustion Science and Technology. [105] 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 144 (2006) 74–88. [106] K. Akselvoll, P. Moin, Large-eddy simulation of turbulent con- fined coannular jets, Journal of Fluid Mechanics 315 (1996) 387– 411. [107] P. M. Gresho, Incompressible fluid dynamics: some fundamental formulation issues, Annual Review of Fluid Mechanics 23 (1991) 413–453. [108] B. J. Boersma, G. Brethouwer, F. T. M. Nieuwstadt, A numerical investigation on the effect of the inflow conditions on the self- similar region of a round jet, Physics of Fluids 10(4) (1998) 899– 910. [109] Information available at: http://www-unix.mcs.anl.gov/mpi/, accessed 12 Dec. 2007. 179 [110] H. R. N. Jones, Radiation heat transfer, Oxford University Press, New York: Oxford University Press, 2000. [111] S. Chandrasekar, Radiative Transfer, Dover Publications, New York, 1960. [112] F. C. Lockwood, N. G. Shah, A new radiation solution method forincorporation in general combustion prediction procedures, Eighteenth Symposium (International) on Combusion, 1981, pp. 1405–1414. [113] P. J. Coelho, O. J. T. and D. Roekaerts, Spectral radiative ef- fects and turbulence-radiation-interaction in Sandia Flame D, in: TNF6 Proceedings, 2002. [114] P. J. Coelho, O. J. Teerling, D. Roekaerts, Spectral radiative effects and turbulence/radiation interaction in a non-luminous turbu- lent jet diffusion flame, Combustion and Flame 133 (April 2003) 75–91(17). [115] W. L. Grosshandler, Radcal: A narrow-band model for radiation calculations in a combustion environment, NIST Technical Note 1402. [116] Information available at: http://public.ca.sandia.gov/TNF/radiation.html, accessed 12 Dec. 2007. [117] B. Jin, M. Wang, W. K. Bushe, Modelling of mean spectral ra- diation intensities in sandia flame d, Proceedings of 2007 CICS Spring Technical Meeting, Banff, Canada. [118] M. Hossain, J. C. Jones, W. Malalasekera, Modelling of a bluff- body nonpremixed flame using a coupled radiation/flamelet combustion model, Flow, Turbulence and Combustion 67 (2001) 217–234. [119] Information available at: http://www.me.berkeley.edu/gri-mech/, accessed 12 Dec. 2007. 180 [120] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eit- eneer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, W. C. Gardiner, V. V. Lissianski, Z. Qin, Information available at: http://www.me.berkeley.edu/gri mech/, accessed 12 Dec. 2007. [121] M. E. Kounalakis, Y. R. Sivathanu, G. M. Faeth, Infrared radia- tion statistics of nonluminous turbulent diffusion flames, ASME Journal of Heat Transfer 113 (1992) 437–445. [122] Y. Zheng, R. S. Barlow, J. P. Gore, Measurements and calculations of spectral radiation intensities for turbulent non-premixed and partially premixed flames, ASME Journal of Heat Transfer 125 (2003) 678–686. [123] S. C. Kong, R. D. Reitz, Numerical study of premixed hcci en- gine combustion and its sensitivity to computational mesh and model uncertainties, Combustion Theory and Modelling (2003) 417–433. [124] R. W. Bilger, Turbulent flows with non-premixed reactants, Top- ics Applied Physics 44 (1980) 65–113. [125] P. A. Libby, F. A. Williams, A presumed pdf analysis of partially premixed turbulent combustion, Combustion Science and Tech- nology 161 (2001) 351–390. [126] F. A. Jaberi, R. S. Miller, P. Givi, Conditional statistics in turbulent scalar mixing and reaction, AlChE Journal 42 (1996) 1149–1152. [127] F. A. Jaberi, R. S. Miller, P. Givi, Non-gaussian scalar statistics in homogeneous turbulence, Journal of Fluid Mechanics 313 (1996) 241–282. [128] G. Ribert, M. Champion, P. Plion, Modelling turbulent reactive flows with variable equivalence ratio: application to the calcula- tion of a reactive shear layer, Combustion Science and Technol- ogy 176(5/6) (2004) 907–924. [129] P. Domingo, S. P. L. Vervisch, R. Hauguel, DNS of a premixed turbulent v flame and les of a ducted flame using a fsd-pdf sub- grid scale closure with fpi-tabulated chemistry, Combustion and Flame 143(4) (2005) 566–586. 181 [130] N. Swaminathan, R. W. Bilger, Analysis of conditional moment closure for turbulent premixed flames, Combustion Theory and Modelling 5 (2001) 241–260. [131] S. M. Martin, J. C. Kramlich, G. Kosaly, J. J. Riley, The premixed conditional moment closure method applied to idealized lean premixed gas turbine combustors, Journal of Engineering For Gas Turbines and Power-Transactions of the ASME 125 (2003) 895–900. [132] K. N. C. Bray, M. Champion, P. A. Libby, N. Swaminathan, Finite rate chemistry and presumed pdfmodels for premixed turbulent combustion, Combustion and Flame 146 (2006) 665–673. [133] R. Grout, An age extended progress variable for conditioning re- action rates, Submitted to Physics of Fluids. [134] K. W. Jenkins, R. S. Cant, Flame kernel interaction in turbulent environment, Third AFOSR Conference on DNS and LES, 1999, pp. 605–612. 182
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Conditional source-term estimation methods for turbulent...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Conditional source-term estimation methods for turbulent reacting flows Jin, Bei 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Conditional source-term estimation methods for turbulent reacting flows |
Creator |
Jin, Bei |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | Conditional Source-term Estimation (CSE) methods are used to obtain chemical closure in turbulent combustion simulation. A Laminar Flamelet Decomposition (LFD) and then a Trajectory Generated Low-Dimensional Manifold (TGLDM) method are combined with CSE in Reynolds-Averaged Navier Stokes (RANS) simulation of non-premixed autoigniting jets. Despite the scatter observed in the experimental data, the predictions of ignition delay from both methods agree reasonably well with the measurements. The discrepancy between predictions of these two methods can be attributed to different ways of generating libraries that contain information of detailed chemical mechanism. The CSE-TGLDM method is recommended for its seemingly better performance and its ability to transition from autoignition to combustion. The effects of fuel composition and injection parameters on ignition delay are studied using the CSE-TGLDM method. The CSE-TGLDM method is then applied in Large Eddy Simulation of a non-premixed, piloted jet flame, Sandia Flame D. The adiabatic CSE-TGLDM method is extended to include radiation by introducing a variable enthalpy defect to parameterize TGLDM manifolds. The results are compared to the adiabatic computation and the experimental data. The prediction of NO formation is improved, though the predictions of temperature and major products show no significant difference from the adiabatic computation due to the weak radiation of the flame. The scalar fields are then extracted and used to predict the mean spectral radiation intensities of the flame. Finally, the application of CSE in turbulent premixed combustion is explored. A product-based progress variable is chosen for conditioning. Presumed Probability Density Function (PDF) models for the progress variable are studied. A modified version of a laminar flame-based PDF model is proposed, which best captures the distribution of the conditional variable among all PDFs under study. A priori tests are performed with the CSE and presumed PDF models. Reaction rates of turbulent premixed flames are closed and compared to the DNS data. The results are promising, suggesting that chemical closure can be achieved in premixed combustion using the CSE method. |
Extent | 2569871 bytes |
Subject |
turbulent reacting flows |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2007-12-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066161 |
URI | http://hdl.handle.net/2429/232 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2008_spring_jin_bei.pdf [ 2.45MB ]
- Metadata
- JSON: 24-1.0066161.json
- JSON-LD: 24-1.0066161-ld.json
- RDF/XML (Pretty): 24-1.0066161-rdf.xml
- RDF/JSON: 24-1.0066161-rdf.json
- Turtle: 24-1.0066161-turtle.txt
- N-Triples: 24-1.0066161-rdf-ntriples.txt
- Original Record: 24-1.0066161-source.json
- Full Text
- 24-1.0066161-fulltext.txt
- Citation
- 24-1.0066161.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066161/manifest