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 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 flamebased 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. Ignition Delay of Methane Jets with Conditional Sourceterm 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 Turbulent Non-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 Density 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 Premixed Combustion with CSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.2 The CSE Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 iv 6.3 Results and Discussion : RANS tests . . . . . . . . . . . . . . . . . . . . . 6.3.1 Conditional Averages: Mechanism I . . . . . . . . . . . . . . . 6.3.2 Reaction Rate Predictions: Mechanism I . . . . . . . . . . . 6.3.3 Comparison with DNS: Mechanism II . . . . . . . . . . . . . 6.4 Results and Discussion : LES tests . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 Comparison with DNS: Mechanism I . . . . . . . . . . . . . . 6.4.2 Comparison with DNS: Mechanism II . . . . . . . . . . . . . 6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 142 146 147 151 151 161 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 2.2 2.3 2.4 2.5 2.6 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.9 4.10 Regimes in non-premixed turbulent combustion . . . . . . . . . . . . 19 Experiment measurement of scattered and conditional average temperature of Sandia D-flame at x/d = 30 . . . . . . . . . . . . . . . . . . 25 Borghi diagram for premixed combustion . . . . . . . . . . . . . . . . . . . 35 Regimes in premixed turbulent combustion . . . . . . . . . . . . . . . . . 37 Bray-Moss-Libby model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Comparing one dimensional manifold generated by ILDM and TGLDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Penetration length of a cold jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 TGLDM for a methane/ethane jet . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 LFD library of a methane/ethane jet . . . . . . . . . . . . . . . . . . . . . . . . 64 Maximum increase of T |ζ for methane/ethane jets . . . . . . . . 65 Comparison of CSE-TGLDM and CSE-LFD predictions . . . . . 67 Comparison of TGLDM and laminar flamelet library . . . . . . . . 68 Ignition Delay τ vs T0 of methane/nitrogen jets . . . . . . . . . . . . . 69 Ignition Delay τ vs T0 of methane and methane with additives jets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Conditional averaged temperature evaluated using TGLDMs of different fuel composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Comparison of Zst contours in different jets . . . . . . . . . . . . . . . . . 73 Temperature fields with Z contours at t = 0.9 ms . . . . . . . . . . . . 74 Temperature fields with Z contours at t = 1.1 ms . . . . . . . . . . . . 75 Ignition delay τ and pressure ratio . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Ignition delay τ and injection duration ti . . . . . . . . . . . . . . . . . . . . 78 Centerline profiles of temperature and species . . . . . . . . . . . . . . 92 Temperature radial profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Radial profiles of YCO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Radial profiles of YH2 O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Radial profiles of YNO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Conditional averaged CO2 at different locations . . . . . . . . . . . . . 97 Conditional averaged H2 O at different locations . . . . . . . . . . . . . 98 Conditional averaged NO at different locations . . . . . . . . . . . . . . 99 Spectral radiation intensities in Flame D . . . . . . . . . . . . . . . . . . . 104 Spectral radiation intensities in Flame D (cont.) . . . . . . . . . . . . 105 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 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 Probability Density Function in a turbulent premixed flame 111 Progress variable in an unstrained laminar flame . . . . . . . . . . . 114 Probability Density Function (PDF) and Cumulative Distribution Function (CDF) in a premixed turbulent flame . . . . . . . . . . . . . 120 PDF on planes in a premixed turbulent flame . . . . . . . . . . . . . . 121 CDF on planes in a premixed turbulent flame . . . . . . . . . . . . . . 122 PDF in different cubes in a premixed turbulent flame . . . . . . 124 CDF in different cubes in a premixed turbulent flame . . . . . . 125 PDF and CDF in the whole domain . . . . . . . . . . . . . . . . . . . . . . . . 127 Test II: PDF on different planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Test II: CDF on different planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Test II: PDF in different cubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 Test II: CDF in different cubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 PDF and CDF in the whole domain using the right and wrong chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 PDF on planes using the right and wrong chemistry . . . . . . . 134 CDF on planes using the right and wrong chemistry . . . . . . . 135 PDF in different cubes using the right and wrong chemistry 136 CDF in different cubes using the right and wrong chemistry 137 Comparison of laminar flames with different chemistry . . . . 138 Conditional averages calculated in RANS paradigm . . . . . . . 144 Conditional averages calculated using β-PDF . . . . . . . . . . . . . . 145 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 6.1 6.2 6.3 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 c deq D Da E F h I J Jh Jk Ka lF lref p P Pr Re sL Sc τc tF tt tk , tη uref vη v V T x Y Meaning Progress variable Equivalent diameter Diffusivity ¨ Damkohler number Activation Energy Body force Specific Enthalpy Identity Matrix Jacobian Matrix Enthalpy diffusivity Species molecular diffusivity Karlovitz Number Premixed flame length scale Non-premixed reference length scale Pressure Probability Density Function Prandtl number Reynolds number Laminar flame velocity Schmidt number Chemical time scale Premixed flame time scale Turbulent time scale Kolmogorov time scale Non-premixed reference velocity Kolmogorov velocity scale velocity vector Eigen(Schur) Basis Matrix Temperature Cartesian coordinates Mass Fraction x Units m m2 · s−1 m2 · s−1 cal·mol−1 J·kg−1 J/(m2 ·s) kg/(m2 ·s) m m bar m/s µs,ms or s µs,ms or s µs,ms or s µs,ms or s m/s m/s m/s K Symbol Z Y |c = c∗ Y |Z = ζ Z ′′2 BML CDF CMC CSE DNS ILDM LES LFD PDF QSSA RANS RTE SFL SPM TGLDM TRI UFL η Λ ρ τ ω˙ χ ζ Meaning Mixture Fraction Conditional Average of Y with Condition c = c∗ Conditional Average of Y with Condition Z = ζ Variance of Mixture Fraction Bray-Moss-Libby Model Cumulative Distribution Function Conditional Moment Closure Conditional Source-term Estimation Direct Numerical Simulation Intrinsically Low Dimensional Manifold Large Eddy Simulation Laminar Flamelet Decomposition Probability Density Function Quasi-steady State Approximation Reynolds Averaged Navier-Stokes Radiative Transfer Equations Steady Flamelet Library Stochastic Particle Model Trajectory-Generated Low Dimensional Manifold Turbulent-Radiation Interaction Unsteady Flamelet Library Kolmogorov length scale Diagonal Matrix of Eigenvalues Density Ignition Delay or Time Duration Chemical Source Term Scalar Dissipation Rate Sample space for mixture fraction xi Units m s−1 kg·m−3 µs,ms or s s−1 s−1 Acknowledgments The thesis work described here would not be possible without the support 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. Rogak, 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 precious 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, increasing concerns have been raised with respect to the environmental impacts of energy consumption. For centuries, the combustion of fossil fuels has been a major form of energy consumption that has been widely used in stationery power generation and transportation industries. Transportation in particular has been identified to be one of the largest sources of energy consumption, greenhouse gases and air pollutants 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 solution 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 compression ignition engines because of its relative abundance, lower emissions 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% reductions 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 handled 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 engines involve mixing low-pressure natural gas with air before the mixture 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 combustion of natural gas. This has motivated recent studies of the combustion of methane and methane mixed with additives. Among different research approaches such as experimental and theoretical studies, numerical simulation provides an efficient way of analyzing combustion phenomena of natural gas-like mixtures. Turbulent 2 combustion modelling 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 different mechanisms such as turbulence/chemistry and turbulence/radiation interactions. It also involves computational challenges such as the numerical solution of stiff ODE systems and the incorporation of complicated 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 combustion of natural gas-like mixtures. Although the cases under study all include methane as a major component in fuel, the modelling technique 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 introduced and their applications in non-premixed and premixed combustion 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) simulations of transient jet flames are performed to understand the autoignition of methane under engine-relevant conditions. Two different methods are combined with CSE to incorporate full chemical mechanisms 3 and their performances are compared head-to-head. The CSE with Trajectory Generated Low Dimensional Manifold (TGLDM) method is recommended 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 methane 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 include 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 postprocessing 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 premixed 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 extended laminar flame based PDF is proposed. The distribution prescribed by different PDFs are compared to Direct Numerical Simulation (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 models 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 turbulent premixed combustion. 5 Chapter 2 Literature Review 2.1 Introduction Turbulent reacting flows are widely observed in nature and in industrial devices. Numerical simulation of turbulent combustion is attractive 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 reacting 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 problem of turbulence modelling. • It is important to use detailed chemical kinetics to describe combustion 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 between different phases. For example, to study the liquid fuel injection and ignition in diesel engines, the breakdown and vaporization 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 turbulent 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 Computational Fluid Dynamics (CFD) methods, then introduces some combustion 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 basic balance equations such as the Navier-Stokes, continuity, species and 7 energy equations. Given the chemical mechanism and transport models, such equations may be solved directly without modelling turbulence or interactions between turbulence and other mechanisms. Alternatively, the equations may be averaged or filtered, resulting in unclosed 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): ∂ρ ∂ρuj + = 0, ∂t ∂xj (2.1) where ρ and uj are density and velocity respectively. • Momentum: ∂p ∂τij ∂ρui ∂ρuj ui + =− + + Fi , ∂t ∂xj ∂xi ∂xj (2.2) where τij is the viscous stress tensor and Fi is a body force. • Species transport equation: ∂ρYk ∂ρuj Yk ∂Jj k + =− + ω˙ k , ∂t ∂xj ∂xj (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: ∂p ∂ ∂ρht ∂ρuj ht + = + (Jj h + ui τij ) + uj Fj , ∂t ∂xj ∂t ∂xj (2.4) where ht is the total specific enthalpy defined by ht = h + ui ui /2, Jj h is the enthalpy diffusion, ui τij and ui Fj denote the effects of viscous and body forces respectively. For Newtonian fluids, τij is defined by: τij = µl ( ∂ui ∂uj ∂uk 2 + ) − µl δij ( ), ∂xj ∂xi 3 ∂xk (2.5) where µl is the molecular viscosity. Species molecular diffusivity Jj k can be approximated by Fick’s law: Jj k = − µl ∂Yk . Sck ∂xj Here Sck is the Schmidt number of species k defined by Sck = (2.6) µl ; ρDk it characterizes the simultaneous momentum and mass diffusion convection 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: N Jj h µl ∂h ∂Yk Pr , =− + − 1)hk ( P r ∂xj k=1 Sck ∂xj (2.7) where N is the number of species. The Prandtl number, P r, compares the diffusive transport of momentum and temperature; it is defined using the thermal diffusivity λ and constant pressure specific heat Cp as P r = µCp . λ The Lewis number of species k compares thermal and mass diffusivities and is defined by: Lek = λ Sck = . Pr ρCp Dk 9 (2.8) Eq. 2.7 can be simplified under the assumption of unity Lewis number; Eq. 2.3 and 2.4 are formally identical under the assumption of low Mach number. These assumptions may be applied to simplify turbulent 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 turbulent transport with viscous forces. It indicates the ratio of the large scales of turbulent motions to the Kolmogorov microscales. ¨ The turbulent Damkohler number is defined by: Da = τt , τc (2.10) where τt and τc are the turbulent and chemical time scale respectively. ¨ The Damkohler number measures how important the turbulence/chemistry 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 = where τk ≡ ν 1/2 ǫ τc , τk (2.11) 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 number 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 hundreds of species and reactions. Lastly, the boundary and initial conditions of a practical system are often difficult to define in DNS. Consequently, 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-averaged Navier-Stokes (RANS) approaches obtain the mean 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 for Q. 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 nonreacting turbulent flows. In turbulent combustion, because of the fluctuations in density caused 11 by heat release, Reynolds averaging generates unclosed terms including fluctuation correlations involving density. To avoid explicitly modelling such correlations, a Favre (mass weighted) averaging of a quantity Q is introduced: Q = Q + Q′′ , where Q = ρQ ρ and Q′′ = e ρ(Q−Q) ρ (2.13) = 0. The Favre averaged equations are: • Mass conservation (continuity): ∂ρ ∂ρuj + =0 ∂t ∂xj (2.14) • Momentum: ∂ρu′′i u′′j ∂ρui ∂ρuj ui ∂p ∂τij + =− − + + Fi ∂t ∂xj ∂xj ∂xi ∂xj (2.15) • Species transport equation: ∂ρu′′j Yk′′ ∂Jj k ∂ρYk ∂ρuj Yk + =− − + ω˙ k ∂t ∂xj ∂xj ∂xj (2.16) • Energy equation: ∂ρu′′j h′′t ∂p ∂ρht ∂ρuj ht ∂ + + =− + (Jj h + ui τij ) + uj Fj . (2.17) ∂t ∂xj ∂xj ∂t ∂xj Unclosed terms are evaluated by modelling or deriving balance equations and solving them. The Reynolds stresses u′′i u′′j are often modelled using turbulence models developed for non-reacting flows, such as the k − ǫ model [13], without explicitly including heat release effects. Species (u′′j Yk′′ ) and temperature (u′′j T ) turbulent fluxes are usually modelled under a gradient transport hypothesis: ρu′′j Yk′′ = − µt ∂ Yk , Sckt ∂xj 12 (2.18) 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 Jhj 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 associated with 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 expensive DNS and (perhaps overly) simple RANS methods. It explicitly computes the large scales that are larger than the mesh size, while modelling the small scales. Although its application in turbulent combustion [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 between heat release, hydrodynamic flow and acoustic waves. 13 In LES, a quantity Q is filtered in the spectral space where components 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 )dx1 dx2 dx3 = 1. −∞ −∞ (2.20) −∞ Typical examples include a cut-off filter in the spectral space, a box filter 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): ∂ρ ∂ρuj + =0 ∂t ∂xj (2.22) • Momentum: ∂p ∂ ∂τij ∂ρui ∂ρuj ui + =− [ρ (ui uj − ui uj )] − + + Fi (2.23) ∂t ∂xj ∂xj ∂xi ∂xj • Species transport equation: ∂ ∂ρYk ∂ρuj Yk ρ uj Yk − uj Yk + ω˙ k + =− ∂t ∂xj ∂xj 14 (2.24) • Energy equation: ∂ρht ∂ρuj ht ∂ ρ uj ht − uj ht + =− ∂t ∂xj ∂xj + ∂p ∂ Jj h + ui τij +uj Fj . + ∂t ∂xj (2.25) The unclosed terms include Reynolds stresses, the species and enthalpy fluxes, filtered laminar diffusion fluxes and filtered chemical reaction rates. Reynolds stresses are usually evaluated using a subgrid scale (SGS) turbulence model. Common SGS models include Smagorinsky models [16] and dynamic Smagorinsky models [17–19]. In a standard 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 turbulent 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 combustion models may be modified 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 in many industrial applications, where fuel and oxidizer are initially separated; both mixing and burning processes occur in the combustion chamber. Examples of non-premixed combustion include diesel engines and furnaces where fuel is supplied separately. In non-premixed combustion, the time for turbulent 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 the mixture 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 ∂ρui Z ∂ + = ∂t ∂xi ∂xi ρDZ ∂Z ∂xi . (2.27) The diffusion coefficient DZ can be approximated with the thermal diffusivity. The mixture fraction Z is attractive because there is no chemical 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 ∂ρuj Z ∂ + = ∂t ∂xj ∂xj ∂ρZ ′′2 ∂ρuj Z ′′2 ∂ + = ∂t ∂xj ∂xj ρD ∂Z − ρu′′i Z ′′ , ∂xj ρu′′i Z ′′ 2 − 2ρu′′i Z ′′ − ρχ. (2.28) (2.29) Here the flux u′′i Z ′′ might be modelled as: u′′i Z ′′ = −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 dissipation rate χ is used: χ = 2D| ▽ Z|2 (2.32) A diffusion time scale τχ may be defined using the inverse of the dissipation rate χ and this is often used to characterize extinction of nonpremixed flames. To study non-premixed turbulent combustion phenomena, different regime diagrams have been developed using dimensionless numbers. Some diagrams use Da and Re [23] to characterize different regimes. At large Da numbers the chemistry is sufficiently fast and the flame has a structure similar to the laminar flame; at small Da numbers where the chemical time scale is relatively large, extinction might occur. 17 Peters proposed a regime diagram to include the effects of turbulent mixing [21]. The diagram is defined using the ratio χq /χst and ′ Zst /(△Z)F . Here, χq is the extinction scalar dissipation rate, χst is the conditional Favre mean scalar dissipation rate at the stoichiometric ′ mixture fraction, Zst 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|st lD . (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) 1/4 where ε is a scaling factor. It is found [21] to be proportional to χst in 1/3 a methane-air diffusion flame with a four-step mechanism and χ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 according 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 fraction field is almost homogeneous. The contour of the mean stoichiometric mixture fraction in a jet diffusion flame is imposed on the plot, which shows the relevance between different zones and the local conditions in the jet. 2.3.1 Infinitely Fast Chemistry In non-premixed turbulent combustion, the time needed for convection and diffusion is usually much longer than that of chemical reactions. 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 presumed 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 = YFIF CM (Z); YO = YOIF CM (Z); T = T IF CM (Z). (2.37) Given a presumed PDF of mixture fraction, the averages of reactive scalars can be evaluated by integration: 1 YFIF CM = 0 YFIF CM (ζ)P (ζ; x, t)dζ (2.38) YOIF CM (ζ)P (ζ; x, t)dζ (2.39) T IF CM (ζ)P (ζ; x, t)dζ, (2.40) 1 YOIF CM = 0 1 T IF CM = 0 where ζ is the sample space of mixture fraction, P (ζ; x, t) is the presumed 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 because 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 maximum possible heat release. However, in turbulent flows where the local diffusion time scales are not so large, the fast chemistry assumption 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 equilibrium 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 solutions of such a prototype laminar flame can be obtained by solving the flamelet equations. The flamelet equations for non-premixed laminar flames are derived by applying the following transformation: ∂ ∂ ∂Z ∂ = + , ∂t ∂τ ∂t ∂Z ∂ ∂Z ∂ , = ∂x1 ∂x1 ∂Z ∂ ∂ ∂Z ∂ = + (k = 2, 3), ∂xk ∂Zk ∂xk ∂Z (2.41) (2.42) (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 being tangential to the stoichiometric isosurface, and τ = t. Assuming unity Lewis numbers for all species, the balance equations of the laminar flame are rewritten in the mixture fraction space as below: ρ ∂Yi ρχ ∂ 2 Yi = + ω˙ i − R(Yi ), ∂τ 2 ∂Z 2 ρχ ∂ 2 T ∂T = − ρ ∂τ 2 ∂Z 2 n i=1 1 ∂p hi ω˙i + { + qR } − R(T ), cp cp ∂t 21 (2.44) (2.45) where qR is the radiation heat loss. The flamelet equations are simplified 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 compared to the other terms. The simplified flamelet equations are: ρ ∂Yi ρχ ∂ 2 Yi = + ω˙ i , ∂τ 2 ∂Z 2 ρχ ∂ 2 T ∂T = ρ − ∂τ 2 ∂Z 2 n i=1 hi ω˙ i . cp (2.46) (2.47) The coordinate system of equations is reduced to (Z, τ ) after the transformation. 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: exp{−2[erfc−1 (2Z)]2 } χ (Z) = . exp{−2[erfc−1 (2Z0 )]2 } ∗ (2.48) Pitsch and Peters [29] obtained an alternate functional form for χ∗ using a symmetric one-dimensional mixing layer configuration: χ∗ (Z) = Z 2 lnZ . Z02 lnZ0 (2.49) Solutions of the flamelet equations are usually calculated in a preprocessing fashion and stored as flamelet libraries. Using such libraries, it is feasible to incorporate detailed chemical mechanisms in CFD computations. 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 turbulent flame can be obtained by using a SLFM library combined with 22 a joint PDF model of Z and χ: Yi = Z∗ χ∗ YiSLF M (Z ∗ , χ∗ )P (Z ∗ , χ∗ ; x, t)dχ∗ dZ ∗ , (2.50) where YiSLF M (Z ∗ , χ∗ ) is the SLFM solution and P (Z ∗ , χ∗ ; x, t) is the presumed joint PDF. SLFM models have been successfully used to study finite-rate chemistry problems. However, neglecting time-dependence on the dissipation rate makes it insufficient in the study of transient phenomena such as extinction and reignition. The effects of unsteadiness are included in Unsteady Laminar Flamelet Models (ULFM) [30–32]. Two different approaches, i.e. the Largrangian Flamelet Model (LFM) and the Eulerian Particle Flamelet Model (EPFM) have been used with ULFM. 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 method, flamelets were attached to these particles, which were used to trace different histories of scalar dissipation rates in the flow; a modeled 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 nonpremixed 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 nonlinearly 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-premixed combustion. In the absence of extinction and reignition, conditioning with Z 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 realizations 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 2500 Temperature(k) 2000 1500 1000 500 0 0 0.2 0.4 0.6 Mixture fraction 0.8 1 Figure 2.2: Experiment measurement of scattered and conditional average 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 ∂Q ∂Q ∂ζ ∂Y ′′ = + + , ∂t ∂t ∂ζ ∂t ∂t (2.53) ∂Q ∇ζ + ∇Y ′′ . ∂ζ (2.54) ∇Y = ∇Q + The molecular diffusion term is given by: ∂Q · (ρD∇Z) ∂ζ ∂Q ∂2Q + ρD(∇Z)2 2 + ρD∇Z · ∇ ∂ζ ∂ζ ∇ · (ρD∇Y ) = ∇ · (ρD∇Q) + (2.55) + ∇ · (ρD∇Y ′′ ). The CMC equation can be derived by substituting the above derivatives into the balance equations of Y (x, t); leading to: ρ|ζ ∂Q χ|ζ ∂ 2 Q ∂Q + ρ|ζ uj |ζ − ρ|ζ ∂t ∂xj 2 ∂ζ 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 ′′ + ρu · ∇Y′′ − ∇ · (Dρ∇Y′′ )|ζ . ∂t (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: ∇ · ( ρ|ζ u′′ Y′′ |ζ P(ζ)) ∇ · ( ρ|ζ (−Dt ∇Q)P (ζ)) eY ≈ ≈ . 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 χ|ζ ∂ 2 Q − = ω|ζ ˙ , ∂t 2 ∂ζ 2 (2.61) which resembles the flamelet equations. This shows the close relation between CMC and the flamelet model where the same one-dimensional mixing assumption is often adopted for species diffusion. The chemical source terms are often closed with the first order CMC hypothesis – that is, the conditional averages of the chemical reaction rate is approximated by evaluating the rate expressions with the conditional 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 fluctuations around the conditional averages are small. In CMC, conditional means of reactive scalars are obtained by solving the above CMC transport 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 hypothesis 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 local extinction and reignition. Although the predictions of the intermediates and fuel reaction rates are improved with the second-order closure method, the requirement of modelling additional unclosed terms make it less desirable. Another approach is to introduce a second conditioning variable. Cha et al. [53] suggested using the scalar dissipation rate as the second conditioning variable. The method predicts the extinction reasonably well, but predicts the onset of reignition too early. Kronenbug [54] proposed to use a normalized enthalpy parameter together with mixture fraction for double conditioning. The local extinction 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 conditioning 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 proposed by Bushe and Steiner [7]. It adopts the CMC first order closure hypothesis, but reduces the computational cost by eliminating the additional 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: 1 f (x, t) = 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 probability 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 presumed 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 assumed in many flows. For example, in an axisymmetric jet, variation in f |ζ on planes normal to the jet axis is usually very small. Generally, 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 |ζ where A j ∈ A, (2.65) denotes the average of ensemble A. Eq. 2.63 is then rewritten as: 1 f≈ 0 f |ζ A P (ζ; xj , t) dζ, 29 j ∈ A. (2.66) Such an equation may be written for each point in the ensemble. Using a numerical quadrature of N points in Z space, one can discretize the integral above and construct a system of linear equations. The conditional 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 solution; the inversion problem is reformulated as a minimization problem. 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 characteristics 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 matrix M ∗ 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 decreasing λ by a factor of 10 makes little difference in the solution of conditional 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 solutions. 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 appropriate 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 computational cost and precision of closure; however, to directly include detailed chemistry in CSE is still too expensive. The system of inversion 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 (TGLDM) method using a purpose-built high-order finite-volume solver. 31 Details of the TGLDM method will be discussed later in this chapter. The Laminar Flamelet Decomposition(LFD) method [61] is an alternative to using reduced chemical kinetics, which incorporates detailed chemical mechanisms in combustion simulation using flamelet libraries. In the decomposition algorithm, conditional means are approximated as a linear combination of basis functions: Nf f |ζ (t; A) ≈ ai Θi (ζ), (2.69) i=1 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: ρ ρ χ ∂ 2 ψi ∂ψi = + ω˙ i , ∂t Lei 2 ∂Z 2 (2.70) with Z being the mixture fraction and χ(Z, t) being the scalar dissipation 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 combination 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 dependence on mixture fraction χ(Z) can be defined using the solution of a one-dimensional mixing layer prescribed by: χ(Z) = e−2(erf c 32 −1 (2Z) 2 ). (2.71) The temporal function χ0 (t) can be defined arbitrarily as: χ0 (t) = c1 e−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 reactive flow under study; c3 is given different constant values to incorporate flamelets of high dissipation rates, which correspond to temporary fluctuations in scalar dissipation observed in turbulent combustion. Now, Eq. 2.66 can be rewritten as: 1 f (x, t) ≈ 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 solving 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 temperature and the mass fractions of CO and CH3 OH [57]. Other scalars such as the reaction rates can be evaluated with Eq. 2.73, using solutions of the coefficient vector and the basis functions. This makes CSE with laminar flamelet decomposition more similar to flamelet methods, 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 simulation 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 unsteady 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 include lean-burn gas turbines and homogeneous charge spark ignition engines. Different from non-premixed combustion, premixed combustion phenomena have a characteristic velocity scale- the laminar burning velocity- and a characteristic length scale- the flame thickness. These new parameters are used to define different regimes in premixed combustion. Regime diagrams have been proposed by numerous researchers [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 different non-dimensional parameters. For example, the Borghi diagram [64] uses a velocity ratio u′ /uref 34 = 1 To rn Pe Fl rtu am rb eF ed ro Fl nt am s, ele ts Well Stirred Reactor, Thickened Flame = 1 Da eT ′ R = 1 log( uuref ) 1 Ka Island Formation and Wrinkled Laminar Flame, Laminar Flame Fronts Flamelet Regime -1 1 ll log( lref ) Figure 2.3: Borghi diagram for premixed combustion and a length ratio ll /lref . In a classic Borghi diagram, u′ is the turbulent 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 turbulence 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 . sL 2 (2.75) ¨ The turbulent Damkohler and Karlovitz number are defined as: sL l , υ ′ lF (2.76) υη2 tF l2 = F2 = 2 , tη η sL (2.77) Da = Ka = 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δ = lδ2 = δ 2 Ka, η2 (2.78) where lδ is the thickness of the inner layer between the chemically inert 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 defined with the above characteristic scales as shown in Fig. 2.4. The ratios l/lF and υ ′ /sL parameterize the regime diagram, which are related as: υ′ = Re sL l lF −1 = Ka 2/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 3 10 K aδ = 1 broken reaction zones 2 10 υ′/s L η = lδ η = lF thin reaction zones 1 10 Ka = 1 Re = 1 corrugated flamelets 0 10 laminar flamelets wrinkled flamelets −1 10 −1 10 0 10 1 2 10 10 3 10 4 10 l/lF Figure 2.4: Regimes in premixed turbulent combustion. flames into wrinkled 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 consequently 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 result, 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 by Re > 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 inner layer. The corrugated flamelet and thin reaction zones regimes are separated by the line Ka = 1. Under such circumstances, the thickness 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 analysis 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 evaluated using the turbulent mixing time τt : ′′ 2 ω˙ F = −ρCEBU (YF )1/2 , τt (2.80) ′′ 2 where YF is the variance of the fuel mass fraction; CEBU is the EddyBreak-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 38 (2.81) 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 products in the local and fully burnt mixture. The reaction rate may be calculated 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 = ρ(c2 − c2 ) = ρ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 Dissipation Model (EDM). In EDM, the mean mass fraction of the deficient reactant is used and the minimum of the reaction 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 calculates 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 β(x, t) α(x, t) 0 1 Figure 2.5: Presumed probability density function in premixed turbulent 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 according 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, the model 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 mixtures at (x, t). The Dirac delta functions δ(c∗ ) and δ(1 − c∗ ) correspond to fresh mixture at c = 0 and complete combustion at c = 1 respectively. The functional form f (c∗ , x, t) for the probability of burning gas can be defined such that: 1 f (c∗ , x, t)dc∗ = 1 (2.87) 0 with f (0) = f (1) = 0. The coefficients are subject to the normalization of the PDF: 1 P (c∗ , x, t)dc∗ = 1, (2.88) α + β + γ = 1. (2.89) 0 from which one can derive: 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 transport. 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 normalization of PDF prescribed by Eq. 2.88 and averages of the progress variable c: 1 ρc∗ P (c∗ , x, t)dc∗ , ρc = 0 41 (2.91) from which one can derive: β= ρc , ρb α = 1 − β. (2.92) (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 ρu −1= − 1, Tu ρb ρu = (1 + τ )ρb = ρ(1 + τ c). (2.94) (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 mechanism and is therefore fixed. The PDF P (c) depends solely on c and τ. The governing equation of the progress variable is: ∂ρc + ▽ · (ρuc) = ▽ · (ρD ▽ c) + ω˙ ∂t (2.98) Under the assumptions of the BML model, 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)ω, ˙ 42 (2.100) where cm = R1 R01 0 cωf ˙ (c)dc ω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 ∂c . ∂xi ∂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 infinitely thin (i.e. c = 0 or c = 1), c2 = c is valid and c′′ 2 is given as: ρc′′ 2 = ρ(c − c)2 = ρ(c2 − c2 ) = ρc(1 − c). (2.104) Eq. 2.101 is rewritten as: ω˙ = 2 ρc(1 − c) , 2cm − 1 τ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: 1 ρQ = ρQ = ρQP (c)dc = αρu Qu + βρb Qb , (2.106) 0 where Qu and Qb 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 model [68]: ω c = ρu s0L I0 Σ, (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) , ˆy L (2.109) ˆ y is the crossing length scale; alternatively, it can be obtained where L by solving a transport equation. For example, an equation was proposed by Trouv´e and Poinsot [70]: ∂Σ ε Σ2 + ▽ · (υΣ) = ▽ · (Dt ▽ Σ) + C1 Σ − C2 sL . ∂t k 1−c 44 (2.110) 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 represent turbulent diffusion, production by flame stretch and flame surface 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 of G0 is fixed for a particular combustion event. A governing equation for G was proposed by Williams [71] as: ∂G + v · ▽G = sL | ▽ G| ∂t (2.111) This G-equation is applicable to thin flames that propagate at a welldefined burning velocity and is therefore suited for describing premixed turbulent combustion in the corrugated flamelets regime. Peters [21] formulated a G-equation in the thin reaction zones regime: ∂G ▽ · (ρD ▽ T ) + ωT | ▽ G|. + v · ▽G = ∂t ρ| ▽ T | (2.112) A common level set equation was formulated for both regimes: sL,s D ∂G + v ∗ · ▽∗ G = | ▽∗ G| − κ∗ | ▽∗ G|, ∗ ∂t νη ν (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η , κ∗ = ηκ, ▽∗ = η ▽ . 45 (2.114) 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 example [72], a flamelet library may be generated where the species, temperature 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 affordable computational cost, it is usually necessary to reduce the full mechanisms. The most commonly used methods include the Quasisteady 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 certain species; in partial equilibrium methods, they are derived by assuming small net reaction rates of fast reactions. These methods have been applied to study premixed and non-premixed flames with success [75–77]. Recently, new techniques have been developed and applied in turbulent 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-Dimensional Manifold (ILDM) method was first proposed 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 relevant to mixing in a turbulent non-premixed flame. The ILDM method 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 by means of a time scale analysis. The detailed chemical kinetic system is described with: ∂ψ = F(ψ) + Ξ(ψ, ∇ψ, ∆ψ), ∂t (2.115) where ψ = (h, p, ψ1 , ψ2 , ..., ψns )T is the scalar vector consisting of enthalpy, 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, diffusion, etc. For a simple homogeneous, adiabatic and isobaric system, Eq. 2.115 reduces to: ∂ψ = F(ψ). ∂t (2.116) The time scales of the system can be separated through a Schur decomposition analysis: J = VΛV, 47 (2.117) 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 eigenvalue representing the fastest process is at the bottom; V and V 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 manifold defined by: Vf F(ψ) = 0, (2.119) where Vf consists of the last n − nc rows taken from V, corresponding to the fast processes. Such a subspace can be found using Eq.2.119, together with additional equations representing system constraints and those of selected progress variables. The final form of the ILDM equation is: Vf F (ψ) = 0, G(ψ, τ ) = P (ψ, τ ) (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 chemical 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: ∂θ = S(θ) + Γ(ψ(θ), ∇ψ(θ), ∆ψ(θ)). ∂t (2.121) To apply ILDM in combustion simulations, the low-dimensional manifolds are usually calculated beforehand and stored as look-up tables. These tables can be retrieved during the combustion simulation and the themokinetic properties are evaluated using the progress variables 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, pressures and enthalpies are a subset of the conditions defined in ILDM computation. Another attractive feature of ILDM is that the mathematical 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 and Maas [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 temperatures, the diffusion and chemical time scales are of the same order. Therefore, the diffusion processes need to be included in the computation of manifolds. 49 Some researchers have combined ILDM with different ways of manifold generation to solve this problem. For example, a related approach called Flame Prolongation of ILDM (FPI) was studied by Gicquel et al. [85]. In this method, the low temperature regime in ILDM was constructed by using the solutions of laminar premixed flames. Oijen et al. [81] introduced the Flamelet-Generated Manifold (FGM) method 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-ILDM method, the ILDM strategy 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 tabulation of manifolds has been studied to reduce the requirements of storage. For instance, a method called in situ adaptive tabulation (ISAT) was combined with ILDM methods by Pope [87], which built the manifold 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) method, originally proposed by Pope and Maas [88]. The TGLDM method involves solving for physically realizable initial states of the system and integrating Eq. 2.115 to evolve the reaction trajectories. Each trajectory consists of solutions in the composition space starting from a given initial state of the system and evolving toward the chemical equilibrium 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 algorithm used in TGLDM guarantees a converged solution. The global optimization problem in ILDM is reduced to a local optimization problem of defining initial states for trajectories. As in ILDM, TGLDM manifolds 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 extremevalue-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 elements can be written for these species: ns Yc;j,i Yi = Ye;j j = 1, ..., ne , (2.122) i=1 where ns is the number of selected species, ne is the number of elements, 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 composition space. Huang and Bushe [60] used a slightly different approach. To construct an np dimensional manifold, constraints on progress variables are applied in addition to the element conservation equations: Pk,i Yi = 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 equations, 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 generate a one dimensional manifold which was compared to an ILDM and good agreement was observed. His analysis also showed that the TGLDM created with constrained equilibrium better approximated detailed chemistry than ILDM methods under moderate perturbation. The TGLDM generated was applied in three different methane/air reaction systems with success including an unstrained, premixed laminar flame, a non-premixed laminar flame and a perfectly stirred reactor. 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 YH2 O was constructed. The manifolds were stored using Delaunay triangulation [92] and retrieved by surface interpolation over the two-dimensional unstructured grid [92]. Five major species (O2 , CH4 , CO, CO2 and H2 O) were used to construct the equations for the initial states, including element conservations of C, H and O, and constrained equilibrium equations of the two progress variables. 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 becomes 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 system was reduced. The resulting trajectory chemistry 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 coupling of mechanisms such as fluid dynamics, chemical reactions, radiation 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 Estimation provides a closure method with a good balance between computational cost and precision. In previous research, CSE-TGLDM and CSELFD 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 compared head-to-head. Such a study would help to clarify the underlying 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 methane 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-Dimensional Manifold (TGLDM) method and then with a Laminar Flamelet Decomposition (LFD) method to reduce the complexity of detailed chemistry. The predictions of ignition delay time are compared with measurements of shock-tube experiments. 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 using 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 the momentum equations 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 significantly 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, Section 2.3.5 and 2.5.2. 3.3 Results and Discussion To validate and compare the CSE-TGLDM and CSE-LFD methods described 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 detected by identifying the first visible flame kernel that subsequently grows into a full jet flame. The experimental techniques used are described by Sullivan et al. [98]. 3.3.1 Simulation of a Non-reactive Jet A cold jet was simulated to validate solvers developed using OpenFOAM. The jet under discussion was studied experimentally and numerically 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 averaged over three independent tests. A scaling model was proposed by Ouellette [99] based on the vortexquasi-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 =Γ U0 t dn ρi ρ0 0.5 , (3.1) where Zt 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 penetration 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 under 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.494 MPa. Simulations of the non-reactive jet have been performed in an axisymmetric 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 computations 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 error because of the injection delay between issuing the command to inject 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 include 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 consistent with the measurements. The profiles of the OpenFOAM and Fluent 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 90 80 70 Zt /dn 60 50 40 30 20 10 0 0 0.5 1 1.5 t(s) 2 2.5 3 −3 x 10 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%C2 H6 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 Pi : injection pressure; ti : injection duration; T0 : pre-combustion a 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 different CSE approaches are implemented using the same code, numerical scheme, grid as well as the same initial and boundary conditions. A detailed chemical mechanism including 71 species and 379 reactions [105] is used in generating the TGLDM and LFD flamelet 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 YH2 O are selected as progress variables because of their relatively long formation time scales. As an example, Fig. 3.2 shows the TGLDM manifold and its Delaunay triangulation for a methane/ethane jet in YCO2 and YH2 O plane, generated at the stoichiometric mixture fraction, 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 solutions 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 pressure of 30 bar. The flamelets shown are all taken prior to ignition in the laminar flamelet time. It is observed that before ignition, the temperature change is not sufficient to differentiate between flamelets; while the combination of T , YCO and YCH3 OH might be used to detect different stages of the autoigniting process. Therefore, temperature and the mass fractions of CO and CH3 OH 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 delay time is largely insensitive to the definition used: the state changes so abruptly when ignition occurs that many different variables all change dramatically and almost simultaneously, such that it is possible to detect ignition by tracking any one of several different variables. For convenience, the ignition delay is defined based on a rise in the conditional mean temperature. Fig. 3.4 shows the maximum increase in the conditional means of temperature in the reaction field plotted against time 62 0.14 0.12 YH2 O 0.1 0.08 0.06 0.04 0.02 0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.12 0.14 0.16 YCO2 (a) TGLDM manifold 0.14 0.12 YH2 O 0.1 0.08 0.06 0.04 0.02 0 0 0.02 0.04 0.06 0.08 0.1 YCO2 (b) Delaunay triangulation Figure 3.2: TGLDM for methane/ethane jet at stoichiometric mixture fraction, T0 = 1200 K, P0 = 30 bar. 63 6 1200 5 1000 4 ωT (erg/g) T (K ) 11 1400 800 3 600 2 400 1 200 0 0.2 0.4 0.6 0.8 x 10 0 0 1 0.2 0.4 Z −5 4 4.5 1 0.6 0.8 1 x 10 4 3.5 3.5 3 3 YC H3 OH 2.5 YC O 0.8 −6 x 10 2 1.5 2.5 2 1.5 1 1 0.5 0 0 0.6 Z 0.5 0.2 0.4 0.6 0.8 0 0 1 0.2 0.4 Z Z Figure 3.3: LFD library of methane/ethane jet, T0 = 1200 K, P0 = 30 bar. —-: flamelet 5; · · · · ·: flamelet 10; − · − line: flamelet 20; −−: flamelet 25. 64 40 Max. Increase in T|ζ 35 30 25 20 15 10 5 0 0 0.2 0.4 0.6 0.8 t(s) 1 1.2 1.4 −3 x 10 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 dramatically. 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 simulated 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 fluctuations 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]. Although 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 discrepancy is due to the different techniques of incorporating detailed chemistry in TGLDM and LFD. For comparison, the TGLDM and laminar 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 (T˙ ) 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 consequence 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 interface 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 1.5 τ (ms) 1.0 0.5 0.7 0.75 0.8 0.85 1000/T (1/K ) (a) Jet A 1.5 τ (ms) 1.0 0.5 0.7 0.75 0.8 0.85 1000/T (1/K ) (b) Jet B Figure 3.5: Comparison of CSE-TGLDM and CSE-LFD predictions to experimental results. +: experimental data; : calculation with CSETGLDM; △ : calculation with CSE-LFD. 67 10 10 8 10 T˙ (K/s) 6 10 4 10 2 10 0 10 1000 1500 2000 2500 3000 T (k) 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 covers a much wider range of conditions and is able to transition from igniting to burning conditions seamlessly. This, combined with the (arguably) better predictions of ignition delay, leads to the conclusion that the CSE-TGLDM method 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, methane was diluted with nitrogen and autoignition was studied in shocktube experiments. The characteristics of the diluted jet are listed as Jet C in Table 3.1. The experimental results are discussed by McTaggart- 68 τ (ms) 1.5 1.0 0.5 0.7 0.75 0.8 0.85 1000/T (1/K ) 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 predictions 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 process, 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 obvious in the experimental data because of the scatter in the measurements. There is as yet insufficient data from the methane/ethane experiments 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 reduction in ignition delay time of the methane/ethane jet – in fact, the re- 69 τ (ms) 1.5 1.0 0.5 0.7 0.75 0.8 0.85 1000/T (1/K ) Figure 3.8: Ignition Delay τ vs T0 , predictions of CSE-TGLDM. ◦: CH4 ; : 90%CH4 + 10%C2 H6 ; △ : 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 reduced 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 nitrogen 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 temperatures at the stoichiometric mixture fraction is presented. Given the same conditional means of progress variables, temperatures interpolated from the TGLDM vary with the fuel composition. The temperatures for methane/ethane jets are the highest in most cases, followed by pure methane jets and then methane/nitrogen jets. The most obvious difference is observed at small values of YCO2 and YH2 O , corresponding 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 predictions 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 C2 H6 or N2 is added because of the larger molecular weight. The gaseous jet momentum is increased and mixing is enhanced. As a demonstration, the contours of stoichiometry Zst in the pure methane and methane with additives 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 1750 T (K) T (K) 1400 1350 1300 1650 1550 0 0.01 0 0.02 0.02 0.04 YCO2 YCO2 2100 T (K) T (K) 2400 2000 1900 2200 1800 0 0.02 0.04 0.06 0.08 0 YCO2 0.05 0.1 YCO2 Figure 3.9: Conditional averaged temperature evaluated using TGLDMs of different fuel composition, at stoichiometric mixture fraction and pre-combustion temperature 1200 K. −·−: 90%CH4 +10%C2 H6 , —-: CH4 , −−: 80%CH4 + 20%N2 . Top left: YH2 O = 0.02, top right: YH2 O = 0.045, bottom left: YH2 O = 0.07 , bottom right: YH2 O = 0.095. 72 r/deq 10 5 0 0 10 20 30 40 50 60 70 80 50 60 70 80 x/deq r/deq (a) t=0.9ms 10 5 0 0 10 20 30 40 x/deq (b) t=1.1ms Figure 3.10: Comparison of Zst contours in different jets. —-: pure methane; −−:80%CH4 + 20%N2 ; − · −: 90%CH4 + 10%C2 H6 . 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 methane/ethane and then methane jets. The temperature fields of the different jets are compared during injection and after injection stops. Contours of mixture fraction are imposed 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 additives. This comparison is of interest because of the possible application in compression ignition engines, where pilot fuel injection can be 73 1200 1000 r/deq 10 800 5 0 600 10 20 30 40 50 60 x/deq 400 (a) CH4 1200 1000 r/deq 10 800 5 0 600 10 20 30 40 50 60 x/deq 400 (b) 80%CH4 + 20%N2 1200 1000 r/deq 10 800 5 0 600 10 20 30 40 x/deq 50 60 400 (c) 90%CH4 + 10%C2 H6 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 r/deq 1200 1000 10 5 0 800 10 20 30 40 50 60 70 80 x/deq 600 (a) CH4 r/deq 1200 1000 10 5 0 800 10 20 30 40 50 60 70 80 x/deq 600 (b) 80%CH4 + 20%N2 r/deq 1200 1000 10 5 0 800 10 20 30 40 50 60 x/deq 70 80 600 (c) 90%CH4 + 10%C2 H6 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 temperatures, ignition delay seems to be mainly affected by changes in chemical 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 significant enough to make 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 delay 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 τ (ms) 1.0 0.5 1.5 2 2.5 3 3.5 4 4.5 5 5.5 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-combustion 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 experimental data. Despite the scatter in the measurements, the numerical predictions show that at the specified high pre-combustion temperature, ignition delay is insensitive to the injection duration. In fact, the profiles of maximum T |ζ for all injection durations coincide with one 77 1.2 τ (ms) 1.0 0.5 1 1.5 2.0 2.5 ti (ms) Figure 3.14: Ignition delay τ and injection duration ti . Cross: experimental 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 methane 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 generating 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 significantly smaller than the scatter in the experimental results. Diluting methane with nitrogen delays autoignition at low pre-combustion temperatures and the effects become less obvious at high temperatures. Increasing pressure ratio reduces ignition delay and the effects become negligible at high pressure ratios. The injection duration seems irrelevant when autoignition starts before the end of injection; its significance 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 turbulent flames. Conditional Source-term Estimation (CSE) has been applied 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 modeled. The models to be used should be a compromise among precision, complexity and computational cost. • The possibility of using CSE to describe Turbulence/Radiation Interaction (TRI) is to be tested. Similar to the Turbulence/Chemistry 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 analyzed. 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 premixture of C2 H2 , H2 , CO2 , N2 and air having nominally the same equilibrium 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]. The modelling work 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 subgridscale convective fluxes are modeled using the eddy viscosity and eddy 81 diffusivity models [56, 59]: 2 τ = −¯ ρ(uu − uu)) = ρ¯νt ((∇v) + (∇v)T − (∇ · v)I), 3 (4.1) −¯ ρ(uYj − uYj ) = ρ¯Dt,j ▽ Yj , (4.2) −¯ ρ(uh − uh) = kt ▽ T , (4.3) where νt is the eddy viscosity, Dt,j is the eddy diffusivity of species j, kt is the thermal eddy conductivity and I is an identity matrix. The turbulent transport coefficients νt , Dt,j and kt 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 volume method with second order accuracy in space is applied. A predictor-corrector-projection scheme is used to integrate in time, with a second-order AdamsBashford scheme in the predictor step and a semi-implicit Crank-Nicholson 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 kinetic mechanism. 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 the Message-Passing Interface (MPI) [109] for parallel computation on clusters. The chemical source-terms are closed with the extended CSE-TGLDM with radiation model. 82 4.3 CSE-TGLDM with 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 radiation 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 models, gas radiation property models and TRI modelling methods were discussed. In terms of flame radiation, it was suggested that the optically thin approximation (where flame absorption is neglected) might be adopted in the simulation of non-luminous flames; more accurate calculations might be done by solving radiative transfer equations with the discrete ordinates method (DOM), where the absorption of radiation was described as a function of local conditions. As for gas radiation properties, the spectral line-based weighted-sum-of-gray-gases (SLW) model was suggested as an alternative to the weighted-sumof-gray-gases (WSGG) model. Comparison was made among calculations 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 coefficients are adopted. Accurate predictions are obtained with the combination of complete RTE and SLW models. Overall, the temperature difference between the adiabatic calculation and those with various radiation 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-TGLDM model, several problems need to be solved. To begin with, 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 accuracy 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 TGLDM method including 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 RADCAL model by Grosshandler [115]. It is assumed that the flame is optically thin and each radiation point source has an unimpeded isotropic view of the cold surroundings. The radiation heat loss is calculated with [116]: nr Q = Q(T, Yi ) = 4σ 1 pi · ap,i 4 · T 4 − T∞ , (4.4) where Q[W/m3 ] is the radiation heat loss rated per unit volume, σ = 5.669×10−8 [W/m2 K 4 ] is the Steffan-Boltzmann constant, nr is the number of the species included in the calculation, pi is the partial pressure of species i in atmospheres (equal to mole fraction times local pressure), ap,i [m−1 atm−1 ] is the Plank mean 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 that CO2 and H2 O 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 1000 2 1000 3 1000 4 1000 5 )+c2 ∗( ) +c3 ∗( ) +c4 ∗( ) +c5 ∗( ), T T T T T (4.5) with the constant coefficients listed in Table 4.1 [116]. The CSE method is proposed to account for TRI where the nonlinearity 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 H2 O -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) 1 Q(x, t) ≈ Q|ζ P (ζ; x, t) dζ, (4.7) 0 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 energy 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 simplification might be acceptable since the overall fraction of radiative heat loss in Flame D is only around 5%. The formulation of such a method 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 simulation. 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: ns h= ns T h0i Yi hi = 1 cp,i dT + 1 , (4.8) T0 where ns is the number of species in the mixture and h0i 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 of multiple 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 generated by changing the initial temperature of fuel and oxidant. Multiple 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 the maximum radiation heat loss observed in the flame under study. These TGLDMs are then tabulated and retrieved later during combustion simulation. The modified CSE-TGLDM method can now evaluate reacting scalars using manifolds with the same enthalpy defects as those observed locally in the flame. The local enthalpy defect of the flame can be calculated directly using the instantaneous solutions of species and temperature, 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 specifically, 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 temperature Tj1 . Interpolation along the enthalpy defect dimension can be performed implicitly, assuming that reactive scalars can be retrieved from the TGLDM using the progress variables. At the same time, T |ζ can be evaluated using the CSE inversion method: 1 T ≈ 0 T |ζ A P˜Z, j (ζ) dζ, 88 j ∈ A. (4.13) Supposing the local enthalpy defect in the flame corresponds to a manifold 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 chemical source terms. In practice, TGLDMs are generated and stored for selected initial temperatures T0 , T1 , ..., Tj , which might not include the desired temperature 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) T |ζ (Tj1 ) = g( YCO2 |ζ , YH2 O |ζ ), (4.15) Here, is the conditional averaged temperature recovered by simple interpolation 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 TGLDM 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 ) 89 , (4.17) where ζst is the stoichiometric mixture fraction; or by solving the following problem: nz min 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 extended 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 library because of the relatively slow process of NO formation. The simulation 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 temperature and major products such as CO2 and H2 O shows minor differences between these two methods, which indicates that the approximation of adiabatic combustion is likely acceptable for the weakly radiating Flame D. The predictions of CO2 and H2 O are slightly improved in the far field, but the improvement is slight. On the other hand, the predictions of NO are quite different. As expected, 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 experimental data well. Slight over-prediction is observed at x/D = 15 and under-prediction at other stations. The over-prediction at x/D = 15 might be a result of accumulated error in the near field where the effects 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 accurate 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 H2 O as well as NO. As observed in the unconditional means, the predictions of YCO2 |ζ and YH2 O |ζ 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 obvious in the conditional means. Overall, the adiabatic computation and that including radiation both capture the trends observed in the experiment. Compared to the adiabatic CSE-TGLDM method, the requirement of storage for CSE-TGLDM with radiation increases dramatically because of the use of multiple TGLDMs. In terms of speed, the computation is slowed down because interpolation over each TGLDM is needed in order to evaluate the enthalpy defect level and to close the 91 0.2 6 0.15 YCO2 T/T0 8 4 2 0 0.1 0.05 0 20 40 0 60 0 20 x/D 40 60 40 60 x/D −4 0.2 1.5 x 10 0.15 YNO Y H2o 1 0.1 0.5 0.05 0 0 20 40 60 x/D 0 0 20 x/D Figure 4.1: Centerline profiles of Flame D. ∆: experimental data; — -: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSETGLDM with radiation. 92 x/D=7.5 x/D=15 4 4 T/T0 6 T/T0 6 2 0 2 0 2 4 0 6 0 8 8 6 6 T/T0 T/T0 r/D x/D=30 4 2 0 2 4 r/D x/D=45 6 8 4 2 0 2 4 r/D 6 8 0 0 5 r/D Figure 4.2: Temperature radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSETGLDM; − · −: predictions of CSE-TGLDM with radiation. 93 10 x/D=15 0.08 0.08 0.06 0.06 2 0.1 YCO YCO 2 x/D=7.5 0.1 0.04 0.02 0 0.04 0.02 0 2 4 0 6 0 r/D x/D=30 0.15 0.15 YCO YCO 4 r/D x/D=45 6 8 2 0.2 2 0.2 2 0.1 0.05 0 0.1 0.05 0 2 4 r/D 6 8 0 0 5 r/D Figure 4.3: Mass fraction of CO2 radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSETGLDM; − · −: predictions of CSE-TGLDM with radiation. 94 10 x/D=15 0.08 0.08 0.06 0.06 YH YH 2 O 0.1 2 O x/D=7.5 0.1 0.04 0.02 0 0.04 0.02 0 2 4 0 6 0 r/D x/D=30 0.15 0.15 2 0.1 YH 2 YH 4 r/D x/D=45 6 8 O 0.2 O 0.2 2 0.05 0 0.1 0.05 0 2 4 r/D 6 8 0 0 5 r/D Figure 4.4: Mass fraction of H2 O radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSETGLDM; − · −: predictions of CSE-TGLDM with radiation. 95 10 −5 6 −5 x/D=7.5 x 10 8 x/D=15 x 10 6 YNO Y NO 4 4 2 2 0 0 2 0 6 r/D x/D=30 −4 1 4 x 10 0 2 −4 1.5 x 10 4 r/D x/D=45 6 8 1 0.6 YNO YNO 0.8 0.4 0.5 0.2 0 0 2 4 r/D 6 8 0 0 5 r/D Figure 4.5: Mass fraction of NO radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSETGLDM; − · −: predictions of CSE-TGLDM with radiation. 96 10 <YCO2|ξ> 0.1 0.05 0 <YCO2|ξ> x/D=15 0 0.5 z x/D=30 0.1 0.05 0 0 0.5 z 1 0.1 0.05 0 1 <YCO2|ξ> <YCO2|ξ> x/D=7.5 0 0.5 z x/D=45 1 0 0.5 z 1 0.1 0.05 0 Figure 4.6: Conditional averages of CO2 radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 97 x/D=15 0.1 <YH2o|ξ> <YH2o|ξ> x/D=7.5 0.05 0 0 0.5 z x/D=30 0.1 0.05 0 1 0 0.5 z x/D=45 1 0 0.5 z 1 0.1 <YH2o|ξ> <YH2o|ξ> 0.15 0.05 0 0 0.5 z 1 0.1 0.05 0 Figure 4.7: Conditional averages of H2 O radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 98 −5 6 x 10 −5 x/D=7.5 8 x 10 x/D=15 0 0.5 z x/D=45 1 0.5 z 1 <YNO|ξ> <Y NO |ξ> 6 4 2 4 2 0 0 −5 8 x 10 0.5 z x/D=30 0 1 −4 x 10 1 <YNO|ξ> 4 <Y NO |ξ> 6 0.5 2 0 0 0.5 z 1 0 0 Figure 4.8: Conditional averages of NO radial profiles at different stations downstream. ∆: experimental data; —-: predictions of adiabatic CSE-TGLDM; − · −: predictions of CSE-TGLDM with radiation. 99 reaction rates accordingly. The additional dimension of enthalpy defect 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 different 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-TGLDM with radiation needs to be optimized. One suggestion is to combine the parallel computation technique with the extended CSE-TGLDM method, where a segment of the flame simulated on each individual processor has a relatively 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 assumptions 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 transport equations are solved) reasonably well. The radiation effects can be incorporated in CSE-TGLDM by introducing enthalpy defect as an additional 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-TGLDM with radiation should be further 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 transfer 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 RADCAL model, 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′λ,w e−κλ (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 the mean 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: 1 i′λ (l∗ , l) = 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 contributions from all elements. 4.4.2 Results and Discussion The radiation characteristics of the piloted methane-air diffusion flame, Sandia Flame D are studied. Large Eddy Simulation of the Flame D has been done with the CSE-TGLDM model 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 different 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% underprediction. 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 complicated 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 evaluated [117] with the CSE method using scalar fields extracted from the adiabatic calculation by Wang [59]. As shown in Fig. 4.10, the discrepancy is negligible, which might be explained by the small difference in temperature, CO2 and H2 O 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 intensities 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 7000 Iλmean,W/m2_sr_µm 6000 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 4 4.5 5 4 4.5 5 (a) x/D = 30 8000 Iλmean,W/m2_sr_µm 7000 6000 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 (b) x/D = 45 Figure 4.9: Spectral radiation intensities in Flame D. ◦: measured; —-: CSE method; −−: mean properties method. 104 6000 Iλmean,W/m2_sr_µm 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 4 4.5 5 (c) x/D = 60 Figure 4.9: Spectral radiation intensities in Flame D (cont.). ◦: measured; —-: 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 predictions of NO are improved using the extended CSE-TGLDM method. The mean spectral radiation intensities in the flame have been analyzed by solving radiative transfer equations and modelling TRI with CSE. The results agree well with experimental results. 105 7000 Iλmean,W/m2_sr_µm 6000 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 4 4.5 5 4 4.5 5 (a) x/D = 30 8000 7000 Iλmean,W/m2_sr_µm 6000 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 (b) x/D = 45 Figure 4.10: Spectral radiation intensities in Flame D, ◦: measured, —-: CSE-TGLDM with radiation, −−: adiabatic CSE-TGLDM. 106 6000 Iλmean,W/m2_sr_µm 5000 4000 3000 2000 1000 0 1 1.5 2 2.5 3 λ,µm 3.5 4 4.5 5 (c) x/D = 60 Figure 4.10: Spectral radiation intensities in Flame D (cont.). ◦: measured; —-: 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 combustion 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 flamelets regime and can be used to account for finite rate chemistry. Detailed 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 premixed 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 introduces an additional dimension while solving the transport equations for conditional means. In previous chapters, Conditional Source-term Estimation (CSE) has been shown to be a promising closure method for non-premixed flames [7, 57, 59, 60]. In CSE, conditional means are obtained by inverting integral 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 homogeneity, 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 variable, 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 using 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 function 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 laminar 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. Different presumed PDF models are studied and a modification is suggested to improve the precision of the laminar flame-based PDF model 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 definition, 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 reaction 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 illustrated in Fig. 5.1, showing the distribution of c at different locations in a turbulent premixed flame. Near the leading edge, where combustion has barely started, and the trailing edge, where the reactions ap- 110 35 x=10 x=30 x=60 x=90 x=120 30 25 p(c) 20 15 10 5 0 0 0.1 0.2 0.3 0.4 0.5 c 0.6 0.7 0.8 0.9 1 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 intermediate 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) −1 (c2 − c2 ) 111 (5.4) 1 0 where β(a, b) = b= a(1 − c) c (5.5) C= 1 , β(a, b) (5.6) (c∗ )a−1 (1 − c∗ )b−1 dc∗ 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 P (c∗ )dc∗ = A + B 1= 0 0 1 1 c∗ P (c∗ )dc∗ = B c= 0 0 1 c2 ∗2 1 ∗ ∗ c P (c )dc = B = 0 0 112 1 dc∗ + C ▽c∗ (5.9) c∗ dc∗ + C ▽c∗ (5.10) c∗ 2 ∗ dc + C. ▽c∗ (5.11) From these equations the coefficients can be easily obtained: c − c2 I1 − I2 (5.12) c2 I1 − cI2 I1 − I2 (5.13) B= C= A = 1 − B ∗ I0 − C, (5.14) where I0 , I1 and I2 are defined as: 1 I0 = 0 1 I1 = 0 1 I2 = 0 1 dc∗ ▽c∗ (5.15) c∗ dc∗ ▽c∗ (5.16) (c∗ )2 ∗ dc . ▽c∗ (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 laminar 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 laminar flame which is independent of the local conditions in the turbulent flame. Consequently, the PDF always rises dramatically towards infinity as c → 0 and c → 1, even if c is not near 0 or 1 and c2 is very small. 113 1 0.9 0.8 0.7 c 0.6 0.5 0.4 0.3 0.2 0.1 0 −Inf xmin x1 x2 xmax +Inf 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 modifications 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 ci +∆c/2 (P dc)i = P dc (5.18) ci −∆c/2 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 ∆c/2 (P dc)0 = P dc (5.19) P dc (5.20) 0 and 1 (P dc)n = 1−∆c/2 where (P dc)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 (P dc)0 and (P dc)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 located in the laminar flame such that the means and variance of c sampled 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 0 if c∗ < c1 or c∗ > c2 P (c∗ ; x, t) = B f (c∗ ) if c ≤ c∗ ≤ c 1 1 2 • Case 2: x1 < xmin < x2 < xmax A2 δ(c∗ ) + B2 f (c∗ ) if c∗ ≤ c2 ∗ P (c ; x, t) = 0 if c∗ > c2 • Case 3: xmin < x1 < xmax < x2 0 if c∗ < c1 ∗ P (c ; x, t) = B f (c∗ ) + C δ(1 − c∗ ) if c∗ ≥ c 3 3 1 • Case 4: x1 < xmin < xmax < x2 P (c∗ ; x, t) = A4 δ(c∗ ) + B4 f (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 − B2 I0 116 (5.24) c−1 I1 − I 0 B3 = C3 = 1 − B3 I0 (5.26) c − c2 I1 − I 2 (5.27) c2 I1 − cI2 I1 − I2 (5.28) B4 = C4 = (5.25) A4 = 1 − B4 ∗ I0 − C4 . (5.29) I0 , I1 and I2 are redefined as: xu 1 dc∗ ▽c∗ (5.30) c∗ dc∗ ▽c∗ (5.31) c∗ 2 ∗ dc , ▽c∗ (5.32) I0 = xl xu I1 = xl xu I2 = xl 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 using 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 flame was 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 and H corresponding to the full mechanism; P is a combination of H2 O and CO2 . The rates of the two reactions, in [mol/cc/s] are calculated as: ω˙ R1 = 7.8(1014 )e −7640 T [F ][R] (5.33) ω˙ R2 = 8.9(1016 )T −1.8 [O][R] (5.34) where [R] is the radical concentration given by: 4.5(105 )e [R] = [P ] 2300 T √ 15 |F | [O]0.5 [I]1.5 e−λ 4 |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 turbulence 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 bimodal 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 laminar PDF more closely matches the true PDF, although it over-predicts the magnitude of the local maximum at c∗ ≈ 0.97. The difference between PDFs from the original and the modified versions is more obvious in the Cumulative Distribution Function (CDF): the CDF calculated 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 locations. Overall, the β-PDF and the original laminar flame PDF predictions 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 20 18 16 14 p(c∗ ) 12 10 8 6 4 2 0 0 0.1 0.2 0.3 0.4 0.5 ∗ 0.6 0.7 0.8 0.9 1 c (a) Probability Density Function 1 0.9 0.8 p(c∗ )dc∗ 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 ∗ 0.6 0.7 0.8 0.9 1 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 15 10 8 10 p(c∗) p(c∗) 6 4 5 2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 0.1 0.2 (a) x/L = 0.234 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.8 0.9 1 (b) x/L = 0.391 25 40 20 30 p(c∗) p(c∗) 15 20 10 10 5 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0.99 (c) x/L = 0.625 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 (d) x/L = 1 Figure 5.4: Probability Density Function on planes in a premixed turbulent flame. —-: DNS; −−: modified PDF model; − · − line: original PDF by Bray et al.; · · · · ·: β-PDF. 121 1 0.8 0.8 0.6 0.6 p(c∗)dc∗ p(c∗)dc∗ 1 0.4 0.4 0.2 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 0.1 0.2 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.8 0.9 1 (b) x/L = 0.391 p(c∗)dc∗ p(c∗)dc∗ (a) x/L = 0.234 0.3 0.8 0.9 1 (c) x/L = 0.625 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 (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 laminar 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 interpolated 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 laminar 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 the modified PDF and the DNS seems less significant in the CDF calculation. In CSE, the presumed PDFs are used 123 30 20 25 15 p(c∗) p(c∗) 20 15 10 10 5 5 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 (a) c = 0.089, c2 = 0.048 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.9 1 (b) c = 0.406, c2 = 0.298 15 10 10 p(c∗) p(c∗) 15 5 0 0 5 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 (c) c = 0.565, c2 = 0.451 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 (d) c = 0.800, c2 = 0.674 Figure 5.6: Probability Density Function in different cubes in a premixed turbulent flame. —-: DNS; −−: modified PDF model; − · −: βPDF. 124 1 0.8 0.8 0.6 0.6 p(c∗)dc∗ p(c∗)dc∗ 1 0.4 0.4 0.2 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.4 c∗ 0.6 0.8 1 (b) c = 0.406, c2 = 0.298 p(c∗)dc∗ p(c∗)dc∗ (a) c = 0.089, c2 = 0.048 0.2 0.9 1 (c) c = 0.565, c2 = 0.451 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 (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: J f≈ 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 variable 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 precision 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 twostep simplified chemical mechanism, where the coefficients in the Arrhenius 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 20 p(c∗ ) 15 10 5 0 0 0.2 0.4 0.6 0.8 1 c∗ (a) Probability Density Function 1 p(c∗ )dc∗ 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 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 15 20 15 p(c∗) p(c∗) 10 10 5 5 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 0.2 (a) x/L = 0.234 0.4 c∗ 0.6 0.8 1 0.8 1 (b) x/L = 0.391 30 70 25 60 50 p(c∗) p(c∗) 20 15 10 30 20 5 0 0 40 10 0.2 0.4 c∗ 0.6 0.8 0 0 1 (c) x/L = 0.625 0.2 0.4 c∗ 0.6 (d) x/L = 1 Figure 5.9: Probability Density Function on different planes. —-: DNS; − · −: original laminar flame-based model; −−: modified PDF model; · · · · ·: β-PDF. 128 1 0.8 0.8 p(c∗)dc∗ p(c∗)dc∗ 1 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.2 0.4 c∗ 0.6 0.8 0 1 0 1 1 0.8 0.8 0.6 0.4 0.2 0.2 0 0.2 0.4 c∗ 0.6 c∗ 0.6 0.8 1 0.8 1 0.6 0.4 0 0.4 (b) x/L = 0.391 p(c∗)dc∗ p(c∗)dc∗ (a) x/L = 0.234 0.2 0.8 0 0 1 (c) x/L = 0.625 0.2 0.4 c∗ 0.6 (d) x/L = 1 Figure 5.10: Cumulative Distribution Function on different planes. —-: DNS; − · −: original laminar flame-based model; −−: modified PDF model; · · · · ·: β-PDF. 129 20 100 80 15 p(c∗) p(c∗) 60 10 40 5 20 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 18 9 16 8 14 7 12 6 10 5 8 3 4 2 2 1 0.2 0.4 c∗ 0.6 0.8 c∗ 0.6 0.8 1 4 6 0 0 0.4 (b) c¯ = 0.413,c¯2 = 0.282 p(c∗) p(c∗) (a) c¯ = 0.089,c¯2 = 0.046 0.2 0 0 1 (c) c¯ = 0.461,c¯2 = 0.311 0.2 0.4 c∗ 0.6 0.8 1 (d) c¯ = 0.689,c¯2 = 0.505 Figure 5.11: Probability Density Function in different cubes. —-: DNS; −−: modified PDF model; − · −: β-PDF. 130 1 0.8 0.8 p(c∗)dc∗ p(c∗)dc∗ 1 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 1 1 0.8 0.8 0.6 0.4 0.2 0.2 0.2 0.4 c∗ 0.6 0.8 c∗ 0.6 0.8 1 0.6 0.4 0 0 0.4 (b) c¯ = 0.413,c¯2 = 0.282 p(c∗)dc∗ p(c∗)dc∗ (a) c¯ = 0.089,c¯2 = 0.046 0.2 0 0 1 (c) c¯ = 0.461,c¯2 = 0.311 0.2 0.4 c∗ 0.6 0.8 1 (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 evaluated with the modified PDF model; the predictions using the ‘right’ chemistry (based on laminar flame solution II) and the ‘wrong’ chemistry (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 precision of the modified laminar flamelet PDF is strongly affected by the chemical kinetics. When the laminar flame solution is based on a different 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 solutions 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 chemistry, 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 15 p(c∗ ) 10 5 0 0 0.2 0.4 0.6 0.8 1 c∗ (a) Probability Density Function 1 p(c∗ )dc∗ 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 c∗ (b) Cumulative Distribution Function Figure 5.13: PDF and CDF in the whole domain. —-: DNS; − · −: modified PDF model with laminar flame solution I; −−: modified PDF model. 133 15 10 10 p(c∗) p(c∗) 15 5 5 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 0.2 (a) x/L = 0.234 0.4 c∗ 0.6 0.8 1 0.8 1 (b) x/L = 0.391 50 70 60 40 50 p(c∗) p(c∗) 30 20 40 30 20 10 10 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 (c) x/L = 0.625 0.2 0.4 c∗ 0.6 (d) x/L = 1 Figure 5.14: Probability Density Function on different planes. —-: DNS; −·−: modified PDF model with laminar flame solution I; −−: modified PDF model. 134 1 0.8 0.8 p(c∗)dc∗ p(c∗)dc∗ 1 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.2 0.4 c∗ 0.6 0.8 0 1 0 1 1 0.8 0.8 0.6 0.4 0.2 0.2 0 0.2 0.4 c∗ 0.6 c∗ 0.6 0.8 1 0.8 1 0.6 0.4 0 0.4 (b) x/L = 0.391 p(c∗)dc∗ p(c∗)dc∗ (a) x/L = 0.234 0.2 0.8 0 1 (c) x/L = 0.625 0 0.2 0.4 c∗ 0.6 (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 70 30 60 25 50 p(c∗) p(c∗) 20 40 30 15 10 20 5 10 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 18 9 16 8 14 7 12 6 10 5 8 3 4 2 2 1 0.2 0.4 c∗ 0.6 0.8 c∗ 0.6 0.8 1 4 6 0 0 0.4 (b) c¯ = 0.413,c¯2 = 0.282 p(c∗) p(c∗) (a) c¯ = 0.089,c¯2 = 0.046 0.2 0 0 1 (c) c¯ = 0.461,c¯2 = 0.311 0.2 0.4 c∗ 0.6 0.8 1 (d) c¯ = 0.689,c¯2 = 0.505 Figure 5.16: Probability Density Function in different cubes. —-: DNS; −·−: modified PDF model with laminar flame solution I; −−: modified PDF model. 136 1 0.8 0.8 p(c∗)dc∗ p(c∗)dc∗ 1 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0.2 0.4 c∗ 0.6 0.8 0 0 1 1 1 0.8 0.8 0.6 0.4 0.2 0.2 0.2 0.4 c∗ 0.6 0.8 c∗ 0.6 0.8 1 0.6 0.4 0 0 0.4 (b) c¯ = 0.413,c¯2 = 0.282 p(c∗)dc∗ p(c∗)dc∗ (a) c¯ = 0.089,c¯2 = 0.046 0.2 0 0 1 (c) c¯ = 0.461,c¯2 = 0.311 0.2 0.4 c∗ 0.6 0.8 1 (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 1 chem−1 chem−2 max(f )/f 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 c∗ Figure 5.18: Comparison of laminar flames with different chemistry ences in f lead to the error in the PDF model when the ‘wrong’ chemistry 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 laminar 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 progress 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 kinetics. 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 subsequently 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 quantities conditional on the progress variable c: ω|c ˙ ∗ ≈ ω(T ˙ |c∗ , Yi |c∗ , ρ|c∗ ), 140 (6.1) 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 conditional averages of different quantities by inverting integral equations, taking advantage of the statistical homogeneity of conditional averages on pre-determined surfaces in the flow field. Similar to non-premixed 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: 1 f (x, t) ≈ 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 wellapproximated by the PDF models described previously, Eq. 6.3 becomes a relatively simple linear system that can be inverted so long as some kind of regularization is provided. Thus, the conditional averages of different scalars can be approximated. Chemical closure can then be achieved by invoking the CMC hypothesis, Eq. 6.1. The chemical 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 averages can be assumed and conditional means obtained at the previous time step used for regularization [57]. In an a priori test, such information is unavailable. Alternately, two methods 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 similarity is expected between the mean flow and the laminar flame. • Solve the inversion equation iteratively from some arbitrary initial condition and use the solutions of a previous iteration for regularization. The computation stops when the solutions become converged or when some maximum number of iterations are reached. 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 orthogonal to the direction of the main flow are extracted from the DNS data base. The reaction rates are also evaluated with scalars conditional 142 on c, using the presumed PDFs and the CSE formulation described previously. An assumption of statistical homogeneity in the conditional averages over the flow field is adopted to construct the linear equation system for inversion. As discussed previously, a suitable regularization method is important for evaluating conditional means in CSE. The effectiveness of the two proposed methods 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 solutions 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 method are extremely poor, and are therefore not presented. Regularization with the laminar flame solution yields an acceptable approximation 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.25 0.04 0.2 0.03 0.15 < YF |c∗ > < YO |c∗ > 0.05 0.02 0.01 0 0 0.1 0.05 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 0.1 0.2 (a) YF |c∗ 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.7 0.8 0.9 1 (b) YO |c∗ 3000 0.03 0.025 2500 0.02 < T |c∗ > < YI |c∗ > 2000 0.015 1500 0.01 1000 0.005 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 500 0 1 0.1 0.2 (c) YI |c∗ 0.3 0.4 0.5 c∗ 0.6 (d) T |c∗ 0.45 0.4 0.35 < ρ|c∗ > 0.3 0.25 0.2 0.15 0.1 0.05 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 (e) ρ|c∗ Figure 6.1: Conditional averages calculated in RANS paradigm, obtained by using modified laminar flame PDF and different regularization methods. —-: DNS; −−: regularization with laminar flame solution (CSE-1); − · −: regularization with iterations up to 2000 times (CSE-2); · · · · ·: laminar flame solution. 144 0.06 0.25 0.05 0.2 0.04 < YO |c∗ > < YF |c∗ > 0.15 0.03 0.1 0.02 0.05 0.01 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 0 0 1 0.1 0.2 (a) YF |c∗ 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.7 0.8 0.9 1 (b) YO |c∗ 0.03 3000 0.025 2500 0.02 < T |c∗ > < YI |c∗ > 2000 0.015 1500 0.01 1000 0.005 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 500 0 1 0.1 (c) YI |c∗ 0.2 0.3 0.4 0.5 c∗ 0.6 (d) T |c∗ 0.45 0.4 0.35 < ρ|c∗ > 0.3 0.25 0.2 0.15 0.1 0.05 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 01 (e) ρ|c∗ Figure 6.2: Conditional averages calculated in RANS paradigm, obtained 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 simulations. 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 evaluated using the modified laminar flame PDF and the β-PDF. The seemingly 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 unacceptable. The chemical reactions rates on different planes are closed by integration. 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 second 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 agree with the DNS data, indicating the effectiveness of chemical 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 regularization (method CSE-1); the conditional means obtained using the CSE inversion method are closer to the DNS than the a priori information 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 introduced for regularization contains error. Altogether, the predictions 147 1.4 1.2 < ω˙ i |c∗ > 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 (a) ω˙ i |c∗ 1.2 1 < ω˙ ii |c∗ > 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 ∗ 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 3 2.5 < ω˙ i |c∗ > 2 1.5 1 0.5 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 (a) ω˙ i |c∗ 3 2.5 < ω˙ ii |c∗ > 2 1.5 1 0.5 0 0 0.1 0.2 0.3 0.4 0.5 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 1 ω˙ i 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 x/L (a) ω˙ i 1 ω˙ ii 0.8 0.6 0.4 0.2 0 0 0.2 0.4 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 regularization with the laminar flame solution. —-: DNS; −−: CSE. 150 of the CSE model 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 reactive scalars. The filtered chemical reactions rates in different cubes are evaluated by integration and plotted against the filtered DNS reaction rates. A constrained least squares analysis is performed to evaluate the precision 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.05 0.22 0.18 0.035 0.16 0.03 0.14 < YO |c∗ > 0.2 0.04 < YF |c∗ > 0.045 0.025 0.12 0.02 0.1 0.015 0.08 0.01 0.06 0.005 0 0 0.04 0.2 0.4 0.6 0.8 0.02 0 1 0.2 0.4 c∗ 0.6 0.8 1 0.8 1 c∗ (a) YF |c∗ (b) YO |c∗ 0.04 2000 0.035 1800 0.03 1600 < T |c∗ > < YI |c∗ > 0.025 1400 0.02 1200 0.015 1000 0.01 800 0.005 0 0 0.2 0.4 0.6 0.8 600 0 1 0.2 0.4 c∗ 0.6 c∗ (c) YI |c∗ (d) T |c∗ 0.5 0.45 < ρ|c∗ > 0.4 0.35 0.3 0.25 0.2 0 0.2 0.4 0.6 0.8 1 c∗ (e) ρ|c∗ Figure 6.6: Conditional averages calculated in RANS paradigm, obtained by using modified laminar flame PDF and different regularization methods. —-: DNS; −−: regularization with laminar flame solution (CSE-1); − · −: regularization with iterations up to 2000 times (CSE-2); · · · · ·: laminar flame solution. 152 1.4 1.2 < ω˙ i |c∗ > 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 c∗ (a) ω˙ i |c∗ 1.4 1.2 < ω˙ ii |c∗ > 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 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 1 0.9 0.8 0.7 ω˙ i 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 x/L (a) ω˙ i 1 0.9 0.8 0.7 ω˙ ii 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 x/L (b) ω˙ ii Figure 6.8: Reaction rates ω˙ on planes in a premixed turbulent flame. Conditional means obtained by using the modified PDF and regularization with the laminar flame solution. —-: DNS; −−: CSE. 154 0.045 0.25 0.04 0.2 0.035 0.15 < YO |c∗ > < YF |c∗ > 0.03 0.025 0.02 0.1 0.015 0.01 0.05 0.005 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 c∗ 0 0 1 0.1 0.2 (a) YF |c∗ 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.7 0.8 0.9 1 (b) YO |c∗ 0.03 3000 0.025 2500 0.02 < T |c∗ > < YI |c∗ > 2000 0.015 1500 0.01 1000 0.005 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 500 0 1 c∗ 0.1 0.2 (c) YI |c∗ 0.3 0.4 0.5 c∗ 0.6 (d) T |c∗ 0.45 0.4 0.35 < ρ|c∗ > 0.3 0.25 0.2 0.15 0.1 0.05 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 (e) ρ|c∗ Figure 6.9: Conditional averages calculated in LES paradigm, obtained by using the modified PDF and regularization with laminar flame solution. —-: DNS; −−: CSE; · · · · ·: laminar flame solution. 155 1.2 1 < ω˙ i |c∗ > 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 c∗ 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 (a) ω˙ i |c∗ 1.2 1 < ω˙ ii |c∗ > 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 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 promising, although it is somewhat disappointing that there might be a positive 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 variable 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 repeated. 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 1.2 1 ω˙ i 0.8 0.6 0.4 0.2 0 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 ω˙ DNS (a) ω˙ i 1.6 1.4 1.2 ω˙ ii 1 0.8 0.6 0.4 0.2 0 −0.5 0 0.5 1 1.5 ω˙ DNS (b) ω˙ ii Figure 6.11: Reaction rates ω˙ calculated in LES paradigm. —-: ω˙ = ω˙ DN S ; −−: ω˙ i = 1.222 ∗ ω˙ DN S,i and ω˙ ii = 1.649 ∗ ω˙ DN S,ii . 158 1 ω˙ i 0.8 0.6 0.4 0.2 0 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 ω˙ DNS (a) ω˙ i 1.4 1.2 ω˙ ii 1 0.8 0.6 0.4 0.2 0 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ω˙ DNS (b) ω˙ ii Figure 6.12: Reaction rates ω˙ calculated in LES paradigm, with denser grid in c-space. —-: ω˙ = ω˙ DN S ; −−: ω˙ i = 1.1345 ∗ ω˙ DN S,i and ω˙ ii = 1.4912 ∗ ω˙ DN S,ii . 159 1 0.9 0.8 0.7 ω˙ i 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 ω˙ DNS (a) ω˙ i 1 0.9 0.8 0.7 ω˙ ii 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 ω˙ DNS (b) ω˙ ii Figure 6.13: Reaction rates ω˙ calculated in LES paradigm, evaluated by using the modified PDF and conditional mean reaction rates extracted directly from the DNS. —-: ω˙ = ω˙ DN S ; −−: ω˙ i = 0.9710 ∗ ω˙ DN S,i and ω˙ ii = 0.9684 ∗ ω˙ DN S,ii . 160 6.4.2 Comparison with DNS: Mechanism II The LES tests described in the previous section is repeated for the turbulent 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 inversion 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 respectively. 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 A product-based progress variable is introduced as a conditioning variable for the CSE model. The mean reaction rates in premixed turbulent flames are modelled using the CSE method and the proposed PDF model. The closure method is tested in both RANS and LES paradigms. In the a priori tests, the conditional means of reactive scalars are obtained 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 the modified 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.05 0.25 0.045 0.04 0.2 0.03 0.15 < YO |c∗ > < YF |c∗ > 0.035 0.025 0.02 0.1 0.015 0.01 0.05 0.005 0 0 0.2 0.4 0.6 0.8 0 0 1 0.2 0.4 c∗ 0.6 0.8 1 0.8 1 c∗ (a) YF |c∗ (b) YO |c∗ 2000 0.04 0.035 1800 0.03 1600 < T |c∗ > < YI |c∗ > 0.025 1400 0.02 0.015 1200 0.01 1000 0.005 0 0 0.1 0.2 0.3 0.4 0.5 ∗ 0.6 0.7 0.8 0.9 800 0 1 0.2 0.4 c 0.6 c∗ (c) YI |c∗ (d) T |c∗ 0.5 0.45 < ρ|c∗ > 0.4 0.35 0.3 0.25 0.2 0 0.2 0.4 0.6 0.8 1 c∗ (e) ρ|c∗ Figure 6.14: Conditional averages calculated in LES paradigm, obtained by using the modified PDF and regularization with laminar flame solution. —-: DNS; −−: CSE; · · · · ·: laminar flame solution. 162 1.4 1.2 < ω˙ i |c∗ > 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 c∗ (a) ω˙ i |c∗ 1.4 1.2 < ω˙ ii |c∗ > 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 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 1 0.9 0.8 0.7 ω˙ i 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 ω˙ DNS (a) ω˙ i 1 0.9 0.8 0.7 ω˙ ii 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 ω DNS (b) ω˙ ii Figure 6.16: Reaction rates ω˙ calculated in LES paradigm. —-: ω˙ = ω˙ DN S , −−: ω˙ i = 1.0454 ∗ ω˙ DN S,i and ω˙ ii = 1.0799 ∗ ω˙ DN S,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 Trajectory Generated Low Dimensional Manifold method to study nonpremixed autoigniting jets of methane and methane mixed with additives. Both methods were implemented using OpenFOAM and applied to study jet flames with the same initial and boundary conditions, using the same grid and numerical scheme for comparison. The predictions were compared to shock-tube experiments. The study shows that the CSE-TGLDM method is superior to the CSE-LFD method; this might be attributed to the wider range of solutions of detailed chemistry included in the TGLDM manifolds, 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 method. The fuel composition of the jet affects the ignition delay by changing 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 negligible 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 chemical mechanism used in the jet simulations. However, the scatter observed in the experimental data could not be captured in the simulation. The scatter in ignition delay might 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 TGLDM method combining continuum and stochastic particle model generated trajectories might be used to simulate the randomness of chemical reaction paths; Large Eddy Simulation is suggested to capture the fluctuations in turbulent flows. The recommended CSE-TGLDM method was applied to obtain closure in LES of a piloted, non-premixed methane jet flame, Sandia Flame D. In the thesis work, the adiabatic CSE-TGLDM method was 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 fashion to how it models the turbulence/chemistry interaction; the major sources of radiation in Sandia D flame are hot CO2 and H2 O; 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 difference 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 turbulent flames was explored. A conditioning variable based on the product mass fraction was proposed. A priori tests show that the presumed PDF of this progress variable may be well approximated with an extended 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 performed in both the RANS and LES paradigms. The results of the tests are promising. However, the effectiveness and precision of chemical closure 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 appears to be inferior to the CSE-TGLDM method. However, it is possible 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 fluctuations. In the study of steady jet flame Sandia flame D, the CSE-TGLDM method has been extended to model radiation. The storage requirements 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 proposed conditioning variable might be studied; improvements to the extended laminar flame based PDF might be made by filtering. Ultimately, 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 vehicle standards and highway diesel fuel sulfur control requirements; 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, Natural 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 turbulence 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 contamination, 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´etais, 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 primitive 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 smagorinsky model with dynamic determination of the filter width ratio, Physics of Fluids 16 (2004) 2514–2528. [20] H. Pitsch, Large-eddy simulation of turbulent combustion, Annual Review of Fluid Mechanics 38 (2006) 453–482. [21] N. Peters, Turbulent combustion, 1st Edition, Cambridge University 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, Topics Applied Physics 44 (1980) 65–113. 171 [25] S. Burke, T. Schumann, Diffusion flames, Industrial and Engineering Chemistry 20 (1928) 998–1004. [26] A. W. Cook, J. J. Riley, A subgrid model for equilibrium chemistry 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 Science 10 (3) (1984) 319–339. [29] H. Pitsch, N. Peters, A consistent flamelet formulation for nonpremixed combustion considering differential diffusion effects, Combustion and Flame 114 (1998) 26–40. [30] N. Peters, Comment and reply on the ”Assessment of combustion and submodels for turbulent nonpremixed hydrocarbon flames”, Combustion and Flame 116 (1999) 675–676. [31] N. Swaminathan, R. Bilger, Assessment of combustion and submodels for turbulent nonpremixed hydrocarbon flames, Combustion 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 formation in a gas turbine combustor using undteady flamelets, in: Proceedings of the Combustion Institute, The Combustion Institute, 1998, pp. 1841–1847. [35] H. Barths, C. Hasse, Simulation of combustion in direct injection diesel engines using an eulerian particle flamelet model, Proceedings 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 piloted 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 admixtures in turbulent flow, Fluid Dynamics 25 (3) (1990) 327–334. [41] Information available at: http://www.ca.sandia.gov/TNF/DataArch/FlameD.html, accessed 12 Dec. 2007. [42] A. Klimenko, R. Bilger, Conditional moment closure for turbulent 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 conditional 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 Modelling 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 conditional moment closure model to a two-dimensional nonpremixed 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 reignition in turbulent nonpremixed combustion using a doublyconditional 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, Numerical Recipes in Fortran, Cambridge University Press, Cambridge, 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 nonpremixed flame using conditional source-term estimation with trajectory generated low-dimensional manifold, in: Thirty-First Symposium (International) on Combustion, The Combustion Institute, 2006, pp. 1701–1709. [60] J. Huang, W. K. Bushe, Simulation of an igniting methane jet using conditional source-term estimation with a trajectory generated low-dimensional manifold, To appear in Combustion Theory and Modelling. [61] W. K. Bushe, H. Steiner, Laminar flamelet decomposition for conditional 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 nonpremixed 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 confined turbulent flames, Thirteenth Symposium (International) on Combusion, 1971, pp. 649–657. [66] B. F. Magnussen, B. H. Hjertager, On mathematical models of turbulent combustion with special emphasis on soot formation and combustion, in: Sixteenth Symposium (International) on Combustion, The Combustion Institute, 1977, pp. 719–729. 175 [67] J. B. Moss, K. N. C. Bray, A unified statistical model of the premixed turbulent flame, Acta Astronaut 4 (1977) 291–319. [68] K. N. Bray, P. A. Libby, Recent developments in the BML model 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´e, T. J. Poinsot, The evolution equation for the flame surface 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 reaktionsgeschwindigkeiten, 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 mechanism, Combustion and Flame 68 (1) (1987) 17–30. [75] J. Warnatz, Concentration, pressure, and temperaturedependence 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 approximations 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 applications in combustion systems, Springer-Verlag, 1993. 176 [78] U. Maas, S. B. Pope, Simplifying chemical kinetics: Intrinsic lowdimensional manifolds in composition space, Combustion and Flame 88 (1992) 239–246. [79] Z. Ren, S. B. Pope, A. Vladimirsky, J. M. Guckenheimer, Application of the ICE-PIC method for the dimension reduction of chemical kinetics, Journal of Chemical Physics 124 (2006) 114111. [80] J. A. van Oijen, L. P. H. de Goey, Modelling of premixed counterflow 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, Seventeenth International Colloquium on the Dynamics of Explosions and Reactive Systems (1999) 25–30. [82] U. Maas, S. B. Pope, Implementation of simplified chemical kinetics based on intrinsic low-dimensional manifolds, TwentyFourth Symposium (International) on Combustion (1992) 103– 112. [83] J. Huang, Natural gas combustion under engine-relevant conditions, Ph.D. thesis, University of British Columbia (2006). [84] J. Nafe, U. Maas, A general algorithm for improving ildms, Combustion Theory and Modelling 6 (2002) 697–709. [85] O. Gicquel, N. Darabiha, D. Thevenin, Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion, in: Twenty-Eighth Symposium (International) on Combustion, 2000, pp. 1901–1908. [86] H. Bongers, J. A. V. Oijen, L. P. H. de Goey, Intrinsic lowdimensional manifold method extended with diffusion, Proceedings of the Combustion Institute 29 (2002) 1371–1378. [87] S. Pope, Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation, Combustion Theory and Modelling 1 (1997) 41–63. 177 [88] S. B. Pope, U. Maas, Simplifying chemical kinetics: Trajectory 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 triangulation package, ACM Transactions on Mathematical software 22 (1996) 1–8. [93] M. Wang, A. Frisque, W. K. Bushe, Trajectory-generated lowdimensional 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 compression ignition combustion, Society of Automotive Engineers 114 (2005) 2005–01–0917. 178 [99] P. Ouellette, Direct injection of natural gas for diesel engine fueling, 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 Fluid Mechanics 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 autoignition 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 confined 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 selfsimilar 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 effects 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 turbulent 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 radiation 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 bluffbody 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. Eiteneer, 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 radiation 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 engine 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, Topics 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 Technology 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 calculation of a reactive shear layer, Combustion Science and Technology 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 subgrid 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 pdf models for premixed turbulent combustion, Combustion and Flame 146 (2006) 665–673. [133] R. Grout, An age extended progress variable for conditioning reaction 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-12-31
pdf
Page Metadata
Item Metadata
Title | Conditional source-term estimation methods for turbulent reacting flows |
Creator |
Jin, Bei |
Publisher | University of British Columbia |
Date | 2007 |
Date Issued | 2007-12-20T17:45:48Z |
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 |
File Format | application/pdf |
Language | eng |
Collection |
Electronic Theses and Dissertations (ETDs) 2008+ |
Date Available | 2007-12-20 |
Provider | Vancouver : University of British Columbia Library |
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 |
Graduation Date | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_2008_spring_jin_bei.pdf [ 2.45MB ]
- Metadata
- JSON: 1.0066161.json
- JSON-LD: 1.0066161+ld.json
- RDF/XML (Pretty): 1.0066161.xml
- RDF/JSON: 1.0066161+rdf.json
- Turtle: 1.0066161+rdf-turtle.txt
- N-Triples: 1.0066161+rdf-ntriples.txt
- Original Record: 1.0066161 +original-record.json
- Full Text
- 1.0066161.txt
- Citation
- 1.0066161.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 19 | 8 |
United States | 18 | 1 |
Canada | 12 | 1 |
Saudi Arabia | 10 | 0 |
Belgium | 4 | 0 |
Japan | 4 | 0 |
Russia | 3 | 0 |
United Kingdom | 3 | 1 |
France | 3 | 0 |
Germany | 3 | 4 |
Spain | 2 | 0 |
Vietnam | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 18 | 14 |
Shenzhen | 12 | 7 |
Ashburn | 9 | 0 |
Brussels | 4 | 0 |
University Park | 4 | 0 |
Beijing | 4 | 1 |
Toronto | 3 | 0 |
Vancouver | 3 | 0 |
Saint Petersburg | 3 | 0 |
Burnaby | 3 | 0 |
Abilene | 3 | 0 |
Oxford | 2 | 0 |
Tokyo | 2 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066161/manifest