Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Conditional source-term estimation methods for turbulent reacting flows Jin, Bei 2007

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


ubc_2008_spring_jin_bei.pdf [ 2.45MB ]
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

Full Text

CONDITIONAL SOURCE-TERM ESTIMATIONMETHODS FOR TURBULENT REACTING FLOWSbyBEI JINBE Tsinghua University, 2000ME Tsinghua University, 2002A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIES(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIADecember 2007? Bei Jin, 2007AbstractConditional Source-term Estimation (CSE) methods are used to obtainchemical closure in turbulent combustion simulation.A Laminar Flamelet Decomposition (LFD) and then a TrajectoryGenerated Low-Dimensional Manifold (TGLDM) method are combin-ed with CSE in Reynolds-Averaged Navier Stokes (RANS) simulationof non-premixed autoigniting jets. Despite the scatter observed in theexperimental data, the predictions of ignition delay from both meth-ods agree reasonably well with the measurements. The discrepancybetween predictions of these two methods can be attributed to differ-ent ways of generating libraries that contain information of detailedchemical mechanism. The CSE-TGLDM method is recommended forits seemingly better performance and its ability to transition from au-toignition to combustion. The effects of fuel composition and injectionparameters on ignition delay are studied using the CSE-TGLDM met-hod.The CSE-TGLDM method is then applied in Large Eddy Simulationof a non-premixed, piloted jet flame, Sandia Flame D. The adiabaticCSE-TGLDM method is extended to include radiation by introducinga variable enthalpy defect to parameterize TGLDM manifolds. The re-sults are compared to the adiabatic computation and the experimentaldata. The prediction of NO formation is improved, though the predic-tions of temperature and major products show no significant differencefrom the adiabatic computation due to the weak radiation of the flame.The scalar fields are then extracted and used to predict the mean spec-tral radiation intensities of the flame.Finally, the application of CSE in turbulent premixed combustionis explored. A product-based progress variable is chosen for condi-tioning. Presumed Probability Density Function (PDF) models for theprogress variable are studied. A modified version of a laminar flame-based PDF model is proposed, which best captures the distribution ofthe conditional variable among all PDFs under study. A priori tests areperformed with the CSE and presumed PDF models. Reaction ratesof turbulent premixed flames are closed and compared to the DNSdata. The results are promising, suggesting that chemical closure canbe achieved in premixed combustion using the CSE method.iiTable of ContentsAbstract ............................................................ iiTable of Contents .................................................. iiiList of Tables ...................................................... viList of Figures ..................................................... viiNomenclature ..................................................... xAcknowledgments ................................................ xiiChapter 1. Introduction and Thesis Outline .................... 11.1 Introduction ............................................... 11.2 Outline .................................................... 3Chapter 2. Literature Review .................................... 62.1 Introduction ............................................... 62.2 CFD Tools for Turbulent Combustion ...................... 72.2.1 Direct Numerical Simulation ........................ 82.2.2 Reynolds Averaged Navier Stokes .................. 112.2.3 Large Eddy Simulation ............................. 132.3 Non-premixed Combustion ............................... 152.3.1 Infinitely Fast Chemistry ........................... 192.3.2 Laminar Flamelet Model ........................... 212.3.3 Conditional Moment Closure ...................... 242.3.4 Conditional Source-term Estimation ................ 282.3.5 Including Chemistry in CSE ........................ 312.4 Premixed Combustion .................................... 342.4.1 Eddy-Break-Up .................................... 382.4.2 The Bray-Moss-Libby Model ....................... 402.4.3 Models Based on Surface Area Estimation .......... 442.4.4 Models Based on the G-equation ................... 452.5 Reduction of Detailed Reaction Mechanism .............. 462.5.1 Intrinsic Low-Dimensional Manifold ............... 472.5.2 Trajectory Generated Low-Dimensional Manifold .. 502.6 Summary ................................................. 54iiiChapter 3. Ignition Delay of Methane Jets with Conditional Source-term Estimation ...................................... 563.1 Introduction .............................................. 563.2 Model Formulation ....................................... 563.2.1 Numerical Scheme ................................. 563.3 Results and Discussion ................................... 573.3.1 Simulation of a Non-reactive Jet .................... 573.3.2 Simulations with CSE-TGLDM and CSE-LFD ...... 613.3.3 Effects of Different Fuel Composition .............. 683.3.4 Effects of Different Injection Parameters ............ 763.4 Conclusions .............................................. 78Chapter 4. Modelling Radiative Effects of a Turbulent Non-premixedFlame using CSE-TGLDM ........................... 804.1 Introduction .............................................. 804.2 LES of Sandia Flame D using CSE-TGLDM ............... 814.3 CSE-TGLDM with Radiation ............................. 834.3.1 Radiation Models .................................. 834.3.2 Formulation ....................................... 844.3.3 Results and Discussion ............................. 904.4 Radiation Intensities Analysis using CSE ................ 1014.4.1 Formulation ...................................... 1014.4.2 Results and Discussion ............................ 1024.5 Conclusion .............................................. 105Chapter 5. A Progress Variable and Presumed Probability Den-sity Function for Premixed Combustion ............ 1085.1 Introduction ............................................. 1085.2 The PDF Models ........................................ 1105.2.1 beta-PDF ............................................ 1115.2.2 A Laminar Flame-based PDF ...................... 1125.2.3 Modified Laminar Flame-based PDF .............. 1145.3 Comparison with DNS .................................. 1175.3.1 DNS Configuration ............................... 1175.3.2 Comparisons of PDF Models ...................... 1185.3.3 Effects of Chemical Kinetics ....................... 1325.4 Conclusions ............................................. 138Chapter 6. Closure of Reaction Rates in Turbulent Premixed Com-bustion with CSE .................................... 1406.1 Introduction ............................................. 1406.2 The CSE Formulation .................................... 140iv6.3 Results and Discussion : RANS tests ..................... 1426.3.1 Conditional Averages: Mechanism I ............... 1426.3.2 Reaction Rate Predictions: Mechanism I ........... 1466.3.3 Comparison with DNS: Mechanism II ............. 1476.4 Results and Discussion : LES tests ....................... 1516.4.1 Comparison with DNS: Mechanism I .............. 1516.4.2 Comparison with DNS: Mechanism II ............. 1616.5 Conclusions ............................................. 161Chapter 7. Conclusions and Future Work ....................... 1667.1 Conclusions ............................................. 1667.2 Future Work ............................................. 168Bibliography ...................................................... 170vList of Tables3.1 Autoigniting jets with various fuel composition .............. 614.1 Coefficients for Planck mean absorption coefficients. ......... 86viList of Figures2.1 Regimes in non-premixed turbulent combustion ............ 192.2 Experiment measurement of scattered and conditional averagetemperature of Sandia D-flame at x/d = 30 .................. 252.3 Borghi diagram for premixed combustion ................... 352.4 Regimes in premixed turbulent combustion ................. 372.5 Bray-Moss-Libby model .................................... 402.6 Comparing one dimensional manifold generated by ILDM andTGLDM .................................................... 523.1 Penetration length of a cold jet .............................. 603.2 TGLDM for a methane/ethane jet ........................... 633.3 LFD library of a methane/ethane jet ........................ 643.4 Maximum increase of angbracketleftT|zetaangbracketright for methane/ethane jets ........ 653.5 Comparison of CSE-TGLDM and CSE-LFD predictions ..... 673.6 Comparison of TGLDM and laminar flamelet library ........ 683.7 Ignition Delay tau vs T0 of methane/nitrogen jets ............. 693.8 Ignition Delay tau vs T0 of methane and methane with additivesjets ......................................................... 703.9 Conditional averaged temperature evaluated using TGLDMs ofdifferent fuel composition ................................... 723.10 Comparison of Zst contours in different jets ................. 733.11 Temperature fields with Z contours at t = 0.9 ms ............ 743.12 Temperature fields with Z contours at t = 1.1 ms ............ 753.13 Ignition delay tau and pressure ratio .......................... 773.14 Ignition delay tau and injection duration ti .................... 784.1 Centerline profiles of temperature and species .............. 924.2 Temperature radial profiles ................................. 934.3 Radial profiles of YCO2 ...................................... 944.4 Radial profiles of YH2O ...................................... 954.5 Radial profiles of YNO ....................................... 964.6 Conditional averaged CO2 at different locations ............. 974.7 Conditional averaged H2O at different locations ............. 984.8 Conditional averaged NO at different locations .............. 994.9 Spectral radiation intensities in Flame D ................... 1044.9 Spectral radiation intensities in Flame D (cont.) ............ 1054.10 Spectral radiation intensities in Flame D based on adiabatic andnon-adiabatic computations ............................... 106vii4.10 Spectral radiation intensities in Flame D based on adiabatic andnon-adiabatic computations (cont.) ........................ 1075.1 Probability Density Function in a turbulent premixed flame 1115.2 Progress variable in an unstrained laminar flame ........... 1145.3 Probability Density Function (PDF) and Cumulative DistributionFunction (CDF) in a premixed turbulent flame ............. 1205.4 PDF on planes in a premixed turbulent flame .............. 1215.5 CDF on planes in a premixed turbulent flame .............. 1225.6 PDF in different cubes in a premixed turbulent flame ...... 1245.7 CDF in different cubes in a premixed turbulent flame ...... 1255.8 PDF and CDF in the whole domain ........................ 1275.9 Test II: PDF on different planes ............................ 1285.10 Test II: CDF on different planes ............................ 1295.11 Test II: PDF in different cubes .............................. 1305.12 Test II: CDF in different cubes .............................. 1315.13 PDF and CDF in the whole domain using the right and wrongchemistry .................................................. 1335.14 PDF on planes using the right and wrong chemistry ....... 1345.15 CDF on planes using the right and wrong chemistry ....... 1355.16 PDF in different cubes using the right and wrong chemistry 1365.17 CDF in different cubes using the right and wrong chemistry 1375.18 Comparison of laminar flames with different chemistry .... 1386.1 Conditional averages calculated in RANS paradigm ....... 1446.2 Conditional averages calculated using beta-PDF .............. 1456.3 Conditionally averaged chemical reaction rates calculated in RANSparadigm .................................................. 1486.4 Conditionally averaged chemical reaction rates calculated in RANSparadigm using beta-PDF .................................... 1496.5 Reaction rates on planes in a premixed turbulent flame .... 1506.6 Test II: Conditional averages calculated in RANS paradigm 1526.7 Test II: Conditionally averaged chemical reaction rates calculatedin RANS paradigm ........................................ 1536.8 Test II: Chemical reaction rates on planes in a premixed turbulentflame ...................................................... 1546.9 Conditional averages calculated in LES paradigm .......... 1556.10 Conditionally averaged chemical reaction rates calculated in LESparadigm .................................................. 1566.11 Chemical reaction rates calculated in LES paradigm ........ 1586.12 Chemical reaction rates calculated in LES paradigm, with densergrid for conditioning ...................................... 159viii6.13 Chemical reaction rates calculated in LES paradigm, conditionalmeans extracted from DNS ................................ 1606.14 Test II: Conditional averages calculated in LES paradigm .. 1626.15 Test II: Conditionally averaged chemical reaction rates calculatedin LES paradigm ........................................... 1636.16 Test II: Chemical reaction rates calculated in LES paradigm 164ixNomenclatureSymbol Meaning Unitsc Progress variabledeq Equivalent diameter mD Diffusivity m2? s-1Da Damk?hler number m2? s-1E Activation Energy cal?mol-1F Body forceh Specific Enthalpy J?kg-1I Identity MatrixJ Jacobian MatrixJh Enthalpy diffusivity J/(m2?s)Jk Species molecular diffusivity kg/(m2?s)Ka Karlovitz NumberlF Premixed flame length scale mlref Non-premixed reference length scale mp Pressure barP Probability Density FunctionPr Prandtl numberRe Reynolds numbersL Laminar flame velocity m/sSc Schmidt numbertauc Chemical time scale ?s,ms or stF Premixed flame time scale ?s,ms or stt Turbulent time scale ?s,ms or stk, teta Kolmogorov time scale ?s,ms or suref Non-premixed reference velocity m/sveta Kolmogorov velocity scale m/sv velocity vector m/sV Eigen(Schur) Basis MatrixT Temperature Kx Cartesian coordinatesY Mass FractionxSymbol Meaning UnitsZ Mixture FractionangbracketleftY |c = casteriskmathangbracketright Conditional Average of Y with Condi-tion c = casteriskmathangbracketleftY |Z = zetaangbracketright Conditional Average of Y with Condi-tion Z = zetaZprimeprime2 Variance of Mixture FractionBML Bray-Moss-Libby ModelCDF Cumulative Distribution FunctionCMC Conditional Moment ClosureCSE Conditional Source-term EstimationDNS Direct Numerical SimulationILDM Intrinsically Low Dimensional Manif-oldLES Large Eddy SimulationLFD Laminar Flamelet DecompositionPDF Probability Density FunctionQSSA Quasi-steady State ApproximationRANS Reynolds Averaged Navier-StokesRTE Radiative Transfer EquationsSFL Steady Flamelet LibrarySPM Stochastic Particle ModelTGLDM Trajectory-Generated Low Dimensio-nal ManifoldTRI Turbulent-Radiation InteractionUFL Unsteady Flamelet Libraryeta Kolmogorov length scale m? Diagonal Matrix of Eigenvalues s-1rho Density kg?m-3tau Ignition Delay or Time Duration ?s,ms or s? Chemical Source Term s-1chi Scalar Dissipation Rate s-1zeta Sample space for mixture fractionxiAcknowledgmentsThe thesis work described here would not be possible without the sup-port and help of many people. I would like to thank all my teachers,old and new friends in China and Canada who have seen me throughthis long journey.I would like to thank Drs. M. H. Davy, J. J. Feng and S. N. Ro-gak, who kindly agree to be my research committee despite their busyschedules.I would like to thank my colleagues and friends in my researchgroup. Thanks to Drs. Jian Huang and Mei Wang for their help on myproject. Thanks to Drs. Gord McTaggart-Cowan and Andrea Frisqueand Heather Jones for being supportive and encouraging all along.I would like to thank my supervisor, Dr. W. Kendal Bushe, for hisadvice, patience and continuous support in every way throughout myproject. My PhD study has been a great experience thanks to the pre-cious opportunities and guidance he gave me.Special thanks to my family - my parents, sister, brother-in-law andnephew, who always believe in me and give me their unconditionallove.Lastly, I would like to dedicate this work to my father, who inspiredand motivated me every step of the way. I wish my efforts would havemade him proud and I will always cherish the memory of him.xiiChapter 1Introduction and Thesis Outline1.1 IntroductionThe last century has witnessed soaring gas prices, deteriorating airquality and alarming global climate changes [1]. In recent years, in-creasing concerns have been raised with respect to the environmentalimpacts of energy consumption. For centuries, the combustion of fos-sil fuels has been a major form of energy consumption that has beenwidely used in stationery power generation and transportation indus-tries. Transportation in particular has been identified to be one of thelargest sources of energy consumption, greenhouse gases and air pol-lutants emissions. As a result, governments now set more and morestringent standards [2,3] for vehicles powered by internal combustionengines. Great efforts have been made to improve fuel efficiency andreduce harmful emissions of traditional IC engines. One promising so-lution is to substitute conventional fuels such as gasoline and dieselwith clean-burning fuels.Natural gas is a viable alternative fuel to replace diesel in compres-sion ignition engines because of its relative abundance, lower emis-sions and lower cost [4]. In terms of availability, Canada is one of theworld?s largest producers of natural gas [4]. As for emission, natural1gas as a transportation fuel has the potential of significantly reducingharmful emissions. For example, it has been demonstrated [5, 6] thatmedium- and heavy-duty natural gas engines have over 90% reduc-tions of carbon monoxide (CO) and particulate matter and more than50% reduction in nitrogen oxides (NOx) relative to commercial dieselengines. Natural gas has high antiknock performance, eliminating theneed to add antiknock ingredients to the fuel. It is safe to use if han-dled properly. Conventional diesel engines can be converted to run onnatural gas with relatively minor changes while maintaining similarmaximum horsepower.Despite these benefits, using natural gas does not automaticallyguarantee lower emissions. The performance of the engine dependson the type of engine technology deployed and the implementation.Currently, the most common use of natural gas as a transportationfuel is in the form of compressed gas (CNG). Spark-ignited CNG en-gines involve mixing low-pressure natural gas with air before the mix-ture enters the combustion chamber. Alternatively, pilot diesel may beinjected for ignition. Recently, High Pressure Direct Injection (HPDI)technology enables direct injection of CNG into the cylinder, allowingdiesel engines to operate on clean-burning CNG. To implement thesetechnologies properly, it is crucial to understand the ignition and com-bustion of natural gas. This has motivated recent studies of the com-bustion of methane and methane mixed with additives.Among different research approaches such as experimental and the-oretical studies, numerical simulation provides an efficient way of ana-lyzing combustion phenomena of natural gas-like mixtures. Turbulent2combustion modelling is known to be challenging. It involves a varietyof closure problems including the traditional problems of turbulencemodelling and unclosed terms caused by the coupling between differ-ent mechanisms such as turbulence/chemistry and turbulence/radia-tion interactions. It also involves computational challenges such as thenumerical solution of stiff ODE systems and the incorporation of com-plicated chemical kinetics in CFD.This thesis focuses on the numerical study of turbulent combustion.The objective is to apply the Conditional Source-term Estimation (CSE)method [7] to study properties of non-premixed and premixed com-bustion of natural gas-like mixtures. Although the cases under studyall include methane as a major component in fuel, the modelling tech-nique can be applied to other gaseous fuels with proper adaption.1.2 OutlineThe thesis starts with a comprehensive review in Chapter 2 to discusscombustion simulation subjects relevant to the thesis work. Basic toolsof CFD are described first. Different combustion models are then in-troduced and their applications in non-premixed and premixed com-bustion are discussed. Some commonly used techniques for includingdetailed chemistry in CFD are introduced and the chapter ends with alist of questions to be studied in the thesis work.In Chapter 3, Reynolds Averaged Navier-Stokes (RANS) simula-tions of transient jet flames are performed to understand the autoigni-tion of methane under engine-relevant conditions. Two different meth-ods are combined with CSE to incorporate full chemical mechanisms3and their performances are compared head-to-head. The CSE with Tra-jectory Generated Low Dimensional Manifold (TGLDM) method is rec-ommended and used thereafter to study non-premixed jet flames. Thepredictions of ignition delay are compared to experimental data; theeffects of additives and injection parameters on the autoignition of me-thane are analyzed.In Chapter 4, Large Eddy Simulation (LES) of a non-premixed steadyjet flame is performed with the recommended CSE-TGLDM method.An extended version of the CSE-TGLDM method is proposed to in-clude a radiation model in the simulation. The effects of radiation heatloss on the temperature field and species formation are analyzed. Themean spectral radiation characteristics are then analyzed in a post-processing fashion where the effectiveness of using the CSE methodto model turbulence/radiation interaction is further proved.Chapters 5 and 6 discuss the application of the CSE method in pre-mixed flames. In Chapter 5, a conditioning variable is introduced andthe presumed probability density function (PDF) of this variable isstudied. PDF models from previous studies are reviewed and an ex-tended laminar flame based PDF is proposed. The distribution pre-scribed by different PDFs are compared to Direct Numerical Simula-tion (DNS) results. The precision of PDF models is analyzed and theimportance of chemical kinetics in the PDF computation is discussed.In Chapter 6, a priori tests of CSE with various presumed PDF mod-els are performed for turbulent premixed flames. Reactive scalars suchas conditional averages of species are evaluated and chemical reactionrates are closed using the CSE method; comparison to the DNS data4shows the potential of the CSE method as a tool to get closure in tur-bulent premixed combustion.5Chapter 2Literature Review2.1 IntroductionTurbulent reacting flows are widely observed in nature and in indus-trial devices. Numerical simulation of turbulent combustion is attrac-tive because it is affordable and informative. Great progress has beenmade thanks to the development of modern computer technology thatenables the simulation of complex turbulent combustion phenomena.Despite this progress, it is still challenging to simulate turbulent re-acting flows. A wide variety of problems are coupled and need to besolved:? It is computationally expensive to directly solve for all the scalesin a turbulent flow, which introduces the traditional closure prob-lem of turbulence modelling.? It is important to use detailed chemical kinetics to describe com-bustion processes, especially in the study of ignition, extinctionand pollutant formation. Such mechanisms usually include largenumber 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 the6system, which requires modelling of complicated interactions be-tween different phases. For example, to study the liquid fuel in-jection and ignition in diesel engines, the breakdown and vapor-ization of the liquid droplets, turbulent mixing and combustionneed to be modelled.? Radiation as an inherent heat transfer mechanism influences theflame structure, soot and NOX formation. Radiation modelling isknown to be complicated; yet accurate calculation of radiationeffects is necessary especially in luminous flames.? The problems listed above are coupled in turbulent combustionsimulation due to the interactions between turbulence and othermechanisms. For example, turbulence/chemistry and turbulence/radiation interactions lead to typical closure problems in turbu-lent combustion modelling.Clearly, turbulent combustion modeling is a complicated and broadsubject. This chapter will address a few aspects of this subject thatare relevant to the thesis work. The review starts with basic Compu-tational Fluid Dynamics (CFD) methods, then introduces some com-bustion models for non-premixed and premixed turbulent combustionand ends with a discussion of different methods used to incorporatedetailed chemical kinetics in CFD.2.2 CFD Tools for Turbulent CombustionNumerical study of turbulent combustion involves the solution of ba-sic balance equations such as the Navier-Stokes, continuity, species and7energy equations. Given the chemical mechanism and transport mod-els, such equations may be solved directly without modelling turbu-lence or interactions between turbulence and other mechanisms. Al-ternatively, the equations may be averaged or filtered, resulting in un-closed terms that need to be modelled. Major CFD tools for turbulentcombustion simulation include Direct Numerical Simulation (DNS),Reynolds-Averaged Navier-Stokes (RANS) methods and Large EddySimulation (LES).2.2.1 Direct Numerical SimulationAs indicated by the name, DNS method directly solves the governingequations of turbulent reactive flows [8]:? Mass conservation (continuity):partialdiffrhopartialdifft +partialdiffrhoujpartialdiffxj = 0, (2.1)where rho and uj are density and velocity respectively.? Momentum:partialdiffrhouipartialdifft +partialdiffrhoujuipartialdiffxj = -partialdiffppartialdiffxi +partialdifftauijpartialdiffxj +Fi, (2.2)where tauij is the viscous stress tensor and Fi is a body force.? Species transport equation:partialdiffrhoYkpartialdifft +partialdiffrhoujYkpartialdiffxj = -partialdiffJjkpartialdiffxj + ?omegak, (2.3)where Yk is the mass fraction of the species k, Jjk is the speciesdiffusive flux and ? k is the species chemical reaction rate.8? Energy equation:partialdiffrhohtpartialdifft +partialdiffrhoujhtpartialdiffxj =partialdiffppartialdifft +partialdiffpartialdiffxj (Jjh +uitauij) +ujFj, (2.4)where ht is the total specific enthalpy defined by ht = h +uiui/2,Jjh is the enthalpy diffusion, uitauij and uiFj denote the effects ofviscous and body forces respectively.For Newtonian fluids, tauij is defined by:tauij = ?l( partialdiffuipartialdiffxj+ partialdiffujpartialdiffxi) - 23?ldeltaij(partialdiffukpartialdiffxk), (2.5)where ?l is the molecular viscosity. Species molecular diffusivity Jjkcan be approximated by Fick?s law:Jjk = - ?lSckpartialdiffYkpartialdiffxj . (2.6)Here Sck is the Schmidt number of species k defined by Sck = ?lrhoDk;it characterizes the simultaneous momentum and mass diffusion con-vection processes in turbulent flows. Dk is the molecular diffusivity ofspecies k relative to the major species. The enthalpy diffusion Jjh canbe approximated by Fourier?s law:Jjh = - ?lPrbracketleftBiggpartialdiffhpartialdiffxj +Nsummationdisplayk=1( PrSck-1)hk partialdiffYkpartialdiffxjbracketrightBigg, (2.7)where N is the number of species. The Prandtl number, Pr, comparesthe diffusive transport of momentum and temperature; it is definedusing the thermal diffusivity lambda and constant pressure specific heat Cpas Pr = ?Cplambda .The Lewis number of species k compares thermal and mass diffu-sivities and is defined by:Lek = SckPr = lambdarhoCpDk. (2.8)9Eq. 2.7 can be simplified under the assumption of unity Lewis num-ber; Eq. 2.3 and 2.4 are formally identical under the assumption of lowMach number. These assumptions may be applied to simplify turbu-lent combustion modelling under proper circumstances.Besides those numbers already discussed, listed below are someother dimensionless numbers that are commonly used to characterizeturbulent reactive flows.The turbulent Reynolds number is defined as:Re = uprimeltnu , (2.9)where uprime is the velocity root mean square related to turbulent kineticenergy k by uprime = radicalbig2/3k, lt is the turbulence integral length scale andnu is the kinematic viscosity. The Reynolds number, Re, compares tur-bulent transport with viscous forces. It indicates the ratio of the largescales of turbulent motions to the Kolmogorov microscales.The turbulent Damk?hler number is defined by:Da = tauttauc, (2.10)where taut and tauc are the turbulent and chemical time scale respectively.The Damk?hler number measures how important the turbulence/ch-emistry interaction is in comparison to other processes.The Karlovitz number is used to characterize interactions betweenturbulent motions and flame structure; it is defined as:Ka = tauctauk, (2.11)where tauk equivalence parenleftbignuepsilon1 parenrightbig1/2 is the Kolmogorov time scale, epsilon1 is the average rate ofenergy dissipation per unit mass.10In DNS, the balance equations Eq. 2.1- 2.7 are solved directly afterdiscretization. Meshes used in DNS need to be fine enough to resolvethe smallest eddies in the turbulent flow. DNS of turbulent combustionis often prohibitively expensive. Firstly and mostly, the computationalcost is high to resolve all time and length scales which might differ byorders of magnitude. To simulate a real flow with large Re, the num-ber of cells needed can be overwhelming [9]. Moreover, the number ofscalars that need to be solved for in turbulent combustion is often verylarge. A detailed chemical kinetic system might include tens or hun-dreds of species and reactions. Lastly, the boundary and initial condi-tions of a practical system are often difficult to define in DNS. Con-sequently, DNS is mostly used to analyze turbulent flames in simpleconfigurations [10?12] such as homogeneous isotropic turbulent flowsand mixing layers.2.2.2 Reynolds Averaged Navier StokesReynolds-averaged Navier-Stokes (RANS) approaches obtain the meanflow field by solving the averaged governing equations. Each quantityQ is decomposed into a mean value Q and a deviation Qprime away fromthe mean:Q = Q +Qprime, (2.12)where Qprime = 0. The instantaneous balance equations may be ensembleor time averaged to derive the transport equations for Q. This Reynoldsaveraging procedure leads to unclosed terms such as uprimeQprime that need tobe modelled. This method has been widely used in simulations of non-reacting turbulent flows.In turbulent combustion, because of the fluctuations in density caused11by heat release, Reynolds averaging generates unclosed terms includ-ing fluctuation correlations involving density. To avoid explicitly mod-elling such correlations, a Favre (mass weighted) averaging of a quan-tity Q is introduced:Q = tildewide +Qprimeprime, (2.13)where tildewide = rhoQrho and tildewiderprimeprime = rho(Q-eQ)rho = 0.The Favre averaged equations are:? Mass conservation (continuity):partialdiffrhopartialdifft +partialdiffrho tildewidejpartialdiffxj = 0 (2.14)? Momentum:partialdiffrhotildewideipartialdifft +partialdiffrho tildewidej tildewideipartialdiffxj = -partialdiffrho tildewidestprimeprimei uprimeprimejpartialdiffxj -partialdiffppartialdiffxi +partialdifftildewiderijpartialdiffxj +tildewideFi (2.15)? Species transport equation:partialdiffrho tildewidekpartialdifft +partialdiffrho tildewidej tildewidekpartialdiffxj = -partialdiffrhotildewideprimeprimej Y primeprimekpartialdiffxj -partialdiffJjkpartialdiffxj + ?omegak (2.16)? Energy equation:partialdiffrhotildewidetpartialdifft +partialdiffrho tildewidej tildewidetpartialdiffxj = -partialdiffrho tildewidestprimeprimej hprimeprimetpartialdiffxj +partialdiffppartialdifft +partialdiffpartialdiffxj (Jjh +uitauij) +ujFj. (2.17)Unclosed terms are evaluated by modelling or deriving balanceequations and solving them. The Reynolds stresses tildewidestprimeprimei uprimeprimej are often mod-elled using turbulence models developed for non-reacting flows, suchas the k - epsilon1 model [13], without explicitly including heat release ef-fects. Species (tildewideprimeprimej Y primeprimek ) and temperature (tildewidestprimeprimej T) turbulent fluxes are usu-ally modelled under a gradient transport hypothesis:rhotildewideprimeprimej Y primeprimek = - ?tScktpartialdifftildewidekpartialdiffxj , (2.18)12where ?t is the turbulent viscosity estimated with turbulence modelsand Sckt is a turbulent Schmidt number of species k. At high Reynoldsnumbers, laminar diffusive fluxes such as Jkj and Jhj are usually smallcompared to turbulent transport. Mean species chemical reaction rates? k need to be closed with turbulent combustion models, which will bediscussed in greater detail later in this chapter.RANS simulation solves for the averaged quantities, providing thestatistical means of scalars associated with turbulent flames. Because ofthe 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 SimulationLarge Eddy Simulation [14] is a compromise between extremely ex-pensive DNS and (perhaps overly) simple RANS methods. It explic-itly computes the large scales that are larger than the mesh size, whilemodelling the small scales. Although its application in turbulent com-bustion [15] is still at an early stage, it is a promising method becauseof the following features:? Compared to DNS, LES only solves for large structures, which iscomputationally more affordable.? The unsteady large scale mixing is simulated, instead of beingaveraged as in RANS.? LES might capture the instabilities caused by the interactions be-tween heat release, hydrodynamic flow and acoustic waves.13In LES, a quantity Q is filtered in the spectral space where compo-nents larger than a given cut-off frequency are suppressed, or in thephysical space where weighted averaging of the quantity is performedin a given volume. The filtering operation is defined by:Q(x) =integraldisplayQ(xasteriskmath)F(x -xasteriskmath)dxasteriskmath (2.19)where F is the LES filter, normalized such that:integraldisplay +infinity-infinityintegraldisplay +infinity-infinityintegraldisplay +infinity-infinityF(x1,x2,x3)dx1dx2dx3 = 1. (2.20)Typical examples include a cut-off filter in the spectral space, a box fil-ter in the physical space and a Gaussian filter in the physical space [14].In turbulent combustion, a mass weighted, Favre filtering is usuallyused:rho tildewide(x) =integraldisplayrhoQ(xasteriskmath)F(x -xasteriskmath)dxasteriskmath (2.21)Any quantity Q may be decomposed into a filtered component Q anda fluctuating component Qprime such that Q = Q +Qprime. The filtered balanceequations are:? Mass conservation (continuity):partialdiffrhopartialdifft +partialdiffrhotildewidejpartialdiffxj = 0 (2.22)? Momentum:partialdiffrhotildewideipartialdifft +partialdiffrho tildewidej tildewideipartialdiffxj = -partialdiffpartialdiffxj [rho( tildewidestuiuj - tildewideui tildewideuj)] -partialdiffppartialdiffxi +partialdifftauijpartialdiffxj +Fi (2.23)? Species transport equation:partialdiffrho tildewidekpartialdifft +partialdiffrho tildewidej tildewidekpartialdiffxj = -partialdiffpartialdiffxjbracketleftBigrhoparenleftBig tildewidestujYk - tildewidejparenrightBig tildewideYkbracketrightBig+ ? k (2.24)14? Energy equation:partialdiffrhotildewidetpartialdifft +partialdiffrho tildewidej tildewidetpartialdiffxj = -partialdiffpartialdiffxjbracketleftBigrhoparenleftBigtildewidestujht - tildewidej tildewidetparenrightBigbracketrightBig+partialdiffppartialdifft+ partialdiffpartialdiffxjparenleftBigJjh +uitauijparenrightBig+ujFj.(2.25)The unclosed terms include Reynolds stresses, the species and en-thalpy fluxes, filtered laminar diffusion fluxes and filtered chemical re-action rates. Reynolds stresses are usually evaluated using a subgridscale (SGS) turbulence model. Common SGS models include Smagorin-sky models [16] and dynamic Smagorinsky models [17?19]. In a stan-dard dynamic procedure, the coefficients of the SGS model are solvedby assuming that they are the same at the grid and test-filter level. Thedynamic subgrid models usually have better performance.LES has been successfully used to study a wide variety of turbu-lent reacting flows. Since neither RANS nor LES resolves the smallestscales, at which the rate-controlling molecular mixing processes andchemical reactions occur, the formulations of these two methods bothrequire modelling of these scales. Therefore, most of the RANS com-bustion models may be modified and applied in LES. A comprehensivereview of LES for turbulent combustion is presented by Pitsch [20].2.3 Non-premixed CombustionNon-premixed combustion is observed in many industrial applications,where fuel and oxidizer are initially separated; both mixing and burn-ing processes occur in the combustion chamber. Examples of non-pre-mixed combustion include diesel engines and furnaces where fuel issupplied separately. In non-premixed combustion, the time for turbu-lent mixing is usually much longer than the chemical time scale [21].15As a result, the mixing process is often the rate-controlling process.It is observed that the reaction rates of non-premixed combustion areclosely related to the state of mixing [8].A conserved scalar, the mixture fraction Z, is often used to describethe state of mixing in non-premixed flames. The mixture fraction canbe defined in various ways. Different definitions of Z [21] essentiallyindicate 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 inthe fuel stream. One definition is given by Bilger [22]:Z = Yi -Yi,OYi,F -Yi,O, (2.26)where Yi is the mass fraction of element i in the local mixture, andYi,F and Yi,O are the the mass fractions of element i in the fuel andoxidizer 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 bychemical reactions. The balance equation of Z is:partialdiffrhoZpartialdifft +partialdiffrhouiZpartialdiffxi =partialdiffpartialdiffxiparenleftbiggrhoDZ partialdiffZpartialdiffxiparenrightbigg. (2.27)The diffusion coefficient DZ can be approximated with the thermal dif-fusivity. The mixture fraction Z is attractive because there is no chemi-cal source term in the balance equation. It is found that in non-premixedflames, 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 innon-premixed reacting flows. These two variables tildewide and tildewidestprimeprime2 might bemodelled or solved for using their transport equations. For example,16in RANS simulation, the following transport equations can be derived:partialdiffrhotildewidepartialdifft +partialdiffrho tildewidej tildewidepartialdiffxj =partialdiffpartialdiffxjparenleftbiggrhoD partialdiffZpartialdiffxj-rhotildewideprimeprimei Zprimeprimeparenrightbigg, (2.28)partialdiffrhotildewidestprimeprime2partialdifft +partialdiffrho tildewidejtildewidestprimeprime2partialdiffxj =partialdiffpartialdiffxjparenleftbiggrhotildewiderprimeprimei Zprimeprime2parenrightbigg-2rhotildewideprimeprimei Zprimeprime -rhotildewide. (2.29)Here the flux tildewideprimeprimei Zprimeprime might be modelled as:tildewideuprimeprimei Zprimeprime = -Dt partialdiff tildewideZpartialdiffxi ; (2.30)The mean scalar dissipation rate defined by:tildewide = tildewiderZ|2 (2.31)needs to be modelled.Scalar dissipation rate is important in the study of non-premixedcombustion. It also appears in the laminar flamelet equations that willbe discussed in detail later. In laminar flows, the instantaneous dissi-pation rate chi is used:chi = 2D| triangleinvZ|2 (2.32)A diffusion time scale tauchi may be defined using the inverse of the dis-sipation rate chi and this is often used to characterize extinction of non-premixed flames.To study non-premixed turbulent combustion phenomena, differ-ent regime diagrams have been developed using dimensionless num-bers. Some diagrams use Da and Re [23] to characterize different regimes.At large Da numbers the chemistry is sufficiently fast and the flame hasa structure similar to the laminar flame; at small Da numbers where thechemical time scale is relatively large, extinction might occur.17Peters proposed a regime diagram to include the effects of turbu-lent mixing [21]. The diagram is defined using the ratio chiq/tildewidest andZprimest/(triangleZ)F . Here, chiq is the extinction scalar dissipation rate, tildewidest is theconditional Favre mean scalar dissipation rate at the stoichiometricmixture fraction, Zprimest is the fluctuation around the mean value at thestoichiometric mixture fraction, (triangleZ)F is the diffusion thickness in theZ space defined as:(triangleZ)F = | triangleinvZ|stlD. (2.33)The diffusion length lD in physical space is defined by:lD =parenleftbiggDstalphaparenrightbigg1/2, (2.34)where Dst is the diffusion coefficient and alpha is the strain rate. (triangleZ)R isa reaction zone thickness in mixture fraction space defined as:(triangleZ)R = epsilon(triangleZ)F , (2.35)where epsilon is a scaling factor. It is found [21] to be proportional to chi1/4st ina methane-air diffusion flame with a four-step mechanism and chi1/3st ina diffusion flame with a one-step mechanism.As shown in Fig.2.1, extinction occurs at high scalar dissipationrates where chiq/tildewidest < 1. Beyond this, three regimes are defined accord-ing to the fluctuations in mixture fraction. In the separated flameletsregime, the fluctuations extend to sufficiently lean and rich mixtureand the diffusion layers around the reaction zone are separated. Inthe connected flame zones regime, the fluctuation is relatively smalldue to intense mixing or partial premixing of the fuel stream. In theconnected reaction zones regime, the mixture fraction fluctuations are18Figure 2.1: Regimes in non-premixed turbulent combustion.smaller than the thickness of the reaction zones and the mixture frac-tion field is almost homogeneous. The contour of the mean stoichio-metric mixture fraction in a jet diffusion flame is imposed on the plot,which shows the relevance between different zones and the local con-ditions in the jet.2.3.1 Infinitely Fast ChemistryIn non-premixed turbulent combustion, the time needed for convec-tion and diffusion is usually much longer than that of chemical reac-tions. Therefore, it is appropriate to assume infinitely fast chemistry inthe study of global properties [24]. With this assumption, the reactivescalars in the flame can be modeled as a function of Z.For example, non-premixed flames can be modelled using a pre-sumed Probability Density Function (PDF) and the infinitely fast che-19mistry model (IFCM). Assuming an infinitely fast single step chemicalreaction [25]:F +O arrowrightP, (2.36)the reactive scalars are related to the mixture fraction by:YF = Y IFCMF (Z); YO = Y IFCMO (Z); T = TIFCM(Z). (2.37)Given a presumed PDF of mixture fraction, the averages of reactivescalars can be evaluated by integration:tildewideY IFCMF =integraldisplay 10Y IFCMF (zeta) tildewide(zeta; x,t)dzeta (2.38)tildewideY IFCMO =integraldisplay 10Y IFCMO (zeta) tildewide(zeta; x,t)dzeta (2.39)tildewideTIFCM =integraldisplay 10TIFCM(zeta) tildewide(zeta; x,t)dzeta, (2.40)where zeta is the sample space of mixture fraction, tildewide(zeta; x,t) is the pre-sumed PDF at (x,t). A beta-PDF prescribed by tildewide and tildewidestprimeprime2 is usually usedto model the air/fuel mixing in non-premixed combustion [26]. Themeans and variances of mixture fraction can be obtained by solvingtheir transport equations.The infinitely fast chemistry with a presumed PDF is attractive be-cause it calculates the averages of reactive scalars without solving theirtransport equations. A chemical equilibrium condition can be adoptedin case of multiple-step chemistry [27]. The IFCM method providesvaluable information about the global flame structure and the maxi-mum possible heat release. However, in turbulent flows where the lo-cal diffusion time scales are not so large, the fast chemistry assump-tion is no longer appropriate and non-equilibrium effects need to beincluded. Moreover, relatively slow chemistry is poorly predicted in20the IFCM method; therefore, it should not be used to study formationof particulate matter, CO or NOx.2.3.2 Laminar Flamelet ModelFlamelet models [28] introduce the scalar dissipation rate chi as a secondparameter besides Z to account for the departure from IFCM equilib-rium states. It is assumed that the local balance between diffusion andreaction in a turbulent reactive flow is similar to that of a laminar flamewith the state of mixing characterized by the same Z and chi. The solu-tions of such a prototype laminar flame can be obtained by solving theflamelet equations.The flamelet equations for non-premixed laminar flames are de-rived by applying the following transformation:partialdiffpartialdifft =partialdiffpartialdifftau +partialdiffZpartialdifftpartialdiffpartialdiffZ, (2.41)partialdiffpartialdiffx1 =partialdiffZpartialdiffx1partialdiffpartialdiffZ, (2.42)partialdiffpartialdiffxk =partialdiffpartialdiffZk +partialdiffZpartialdiffxkpartialdiffpartialdiffZ (k = 2,3), (2.43)where (x1,x2,x3,t) is the original coordinate system in the physicalspace; (Z1,Z2,Z3,tau) is the new coordinate system, defined as Z1 = triangleinvZ|triangleinvZ|being locally normal to the stoichiometric flame surface, Z2 and Z3 be-ing tangential to the stoichiometric isosurface, and tau = t. Assumingunity Lewis numbers for all species, the balance equations of the lami-nar flame are rewritten in the mixture fraction space as below:rhopartialdiffYipartialdifftau = rhochi2 partialdiff2YipartialdiffZ2 + ?omegai -R(Yi), (2.44)rhopartialdiffTpartialdifftau = rhochi2 partialdiff2TpartialdiffZ2 -nsummationdisplayi=1hicp ?omegai +1cp {partialdiffppartialdifft +qR} -R(T), (2.45)21where qR is the radiation heat loss. The flamelet equations are sim-plified by neglecting the radiation and pressure fluctuation terms inEq. 2.45, as well as R(Yi) and R(T), which are of a lower order com-pared to the other terms. The simplified flamelet equations are:rhopartialdiffYipartialdifftau = rhochi2 partialdiff2YipartialdiffZ2 + ?omegai, (2.46)rhopartialdiffTpartialdifftau = rhochi2 partialdiff2TpartialdiffZ2 -nsummationdisplayi=1hicp ?omegai. (2.47)The coordinate system of equations is reduced to (Z,tau) after the trans-formation. The scalar dissipation rate chi is used to relate the laminarflame solutions to the turbulent flow. A typical model for chi takes theform chi = chi0chiasteriskmath(Z), where chiasteriskmath(Z) is some presumed function describingthe relationship between Z and chi. Peters [28] developed an analyticalsolution of chiasteriskmath(Z) using the configuration of an infinite one dimensionalmixing layer:chiasteriskmath(Z) = exp{-2[erfc-1(2Z)]2}exp{-2[erfc-1(2Z0)]2}. (2.48)Pitsch and Peters [29] obtained an alternate functional form for chiasteriskmath us-ing a symmetric one-dimensional mixing layer configuration:chiasteriskmath(Z) = Z2lnZZ20 lnZ0 . (2.49)Solutions of the flamelet equations are usually calculated in a pre-processing fashion and stored as flamelet libraries. Using such libraries,it is feasible to incorporate detailed chemical mechanisms in CFD com-putations. In flamelet libraries, reactive scalars are often tabulated forgiven values of the two control parameters as Yi(Z,chi). The solutions offlamelet equations without the time dependent terms are used in theSteady Laminar Flamelet Model (SLFM). The mean properties of a tur-bulent flame can be obtained by using a SLFM library combined with22a joint PDF model of Z and chi:tildewideYi =integraldisplayZasteriskmathintegraldisplaychiasteriskmathY SLFMi (Zasteriskmath,chiasteriskmath) tildewide(Zasteriskmath,chiasteriskmath; x,t)dchiasteriskmathdZasteriskmath, (2.50)where Y SLFMi (Zasteriskmath,chiasteriskmath) is the SLFM solution and tildewide(Zasteriskmath,chiasteriskmath; x,t) is the pre-sumed joint PDF.SLFM models have been successfully used to study finite-rate che-mistry problems. However, neglecting time-dependence on the dissi-pation rate makes it insufficient in the study of transient phenomenasuch as extinction and reignition. The effects of unsteadiness are in-cluded in Unsteady Laminar Flamelet Models (ULFM) [30?32]. Twodifferent approaches, i.e. the Largrangian Flamelet Model (LFM) andthe Eulerian Particle Flamelet Model (EPFM) have been used with ULF-M. For example, Mauss et al. [33] studied extinction and re-ignition in aturbulent jet flame using unsteady flamelets. In his work, a Lagrangiantime measured along the stoichiometric interface was used to accountfor the effects of time-evolving chi. Barths et al. [34, 35] used unsteadyflamelet solutions by tracing imaginary marker particles. In his me-thod, flamelets were attached to these particles, which were used totrace different histories of scalar dissipation rates in the flow; a mod-eled Eulerian convective-diffusive equation was solved to evaluate theprobability of finding these particles at specific locations.The flamelet models considerably relax the hypothesis required inthe infinite fast chemistry model and have been widely used in non-premixed combustion modelling. Efforts have been made to improveflamelet models, such as including differential diffusion effects [36,37]and including diffusion in the direction tangential to the isosurface ofZst [38], etc.232.3.3 Conditional Moment ClosureConditional Moment Closure (CMC) methods were first proposed byBilger [39] and Klimenko [40] independently. CMC is based on theobservation that in non-premixed systems, the reaction rates are non-linearly and crucially dependent on the stoichiometry; the fluctuationsaround mean properties are significantly reduced after conditioningwith a stoichiometry variable. Therefore, accurate closure might beachieved in non-premixed combustion by using conditional moments.Mixture fraction Z is used as the conditioning variable in non-pre-mixed combustion. In the absence of extinction and reignition, condi-tioning with Z is usually sufficient to reduce the impact of nonlinearity.For illustration, temperature at the station x/d = 30 in Sandia flameD [41] is plotted against mixture fraction Z in Fig. 2.2. The scattereddata points all lie within close vicinity of the conditional average oftemperature.Bilger [39,40,42] developed CMC equations with a decompositionapproach. In his work, the conditional average of a scalar Y (x,t) isdefined as:Q(zeta,x,t) equivalence angbracketleftY (x,t)|Z(x,t) = zetaangbracketright equivalence angbracketleftY |zetaangbracketright, (2.51)where the angular brackets is an ensemble average over different real-izations of the flow and the vertical bar means ?given the condition?.The basic decomposition of Y (x,t) can be expressed as:Y (x,t) = Q(Z(x,t),x,t) +Y primeprime(x,t), (2.52)where Q is the conditional average of Y given Z(x,t) = zeta and Y primeprime(x,t)is the fluctuation around the conditional means. The derivatives of Y240 0.2 0.4 0.6 0.8 105001000150020002500Mixture fractionTemperature(k)Figure 2.2: Experiment measurement of scattered and conditional av-erage temperature of Sandia D-flame at x/d = 30 (grey dots: scattereddata; solid black line: conditional average temperature).25can be rewritten as:partialdiffYpartialdifft =partialdiffQpartialdifft +partialdiffQpartialdiffzetapartialdiffzetapartialdifft +partialdiffY primeprimepartialdifft , (2.53)nablaY = nablaQ + partialdiffQpartialdiffzeta nablazeta +nablaY primeprime. (2.54)The molecular diffusion term is given by:nabla? (rhoDnablaY ) = nabla? (rhoDnablaQ) + partialdiffQpartialdiffzeta ? (rhoDnablaZ) (2.55)+rhoD(nablaZ)2 partialdiff2Qpartialdiffzeta2 +rhoDnablaZ ? nablapartialdiffQpartialdiffzeta+nabla? (rhoDnablaY primeprime).The CMC equation can be derived by substituting the above deriva-tives into the balance equations of Y (x,t); leading to:angbracketleftrho|zetaangbracketrightpartialdiffQpartialdifft +angbracketleftrho|zetaangbracketrightangbracketleftuj|zetaangbracketright partialdiffQpartialdiffxj-angbracketleftrho|zetaangbracketrightangbracketleftchi|zetaangbracketright2 partialdiff2Qpartialdiffzeta2 (2.56)= angbracketleftrho|zetaangbracketrightangbracketleft ? |zetaangbracketright +eQ +eY ,whereangbracketleftchi|zetaangbracketright = 2DangbracketleftnablaZ ? nablaZ|zetaangbracketright (2.57)is the conditional scalar dissipation rate that can be modeled [43?46];other unclosed terms are:eQ equivalence angbracketleftnabla(rhoDnablaQ) +rhoDnablaZ ? nablapartialdiffQpartialdiffzeta |zetaangbracketright, (2.58)andeY equivalence -angbracketleftrhopartialdiffYprimeprimepartialdifft +rhou ? nablaYprimeprime -nabla? (DrhonablaYprimeprime)|zetaangbracketright. (2.59)The term eQ represents contributions from molecular diffusion downthe mean gradients of Q and Z; it is negligible at large Re numbers.26The other term eY needs to be modelled. Klimenko [42] proposed toclose eY with:eY approxequal nabla? (angbracketleftrho|zetaangbracketrightangbracketleftuprimeprimeYprimeprime|zetaangbracketrightP(zeta))P(zeta)angbracketleftrho|zetaangbracketright approxequalnabla? (angbracketleftrho|zetaangbracketright(-DtnablaQ)P(zeta))P(zeta)angbracketleftrho|zetaangbracketright . (2.60)Similar to Eq. 2.56, a CMC equation can be derived for enthalpy. Ina homogeneous field, the CMC equations reduce to:partialdiffQpartialdifft -angbracketleftchi|zetaangbracketright2partialdiff2Qpartialdiffzeta2 = angbracketleft ?omega|zetaangbracketright, (2.61)which resembles the flamelet equations. This shows the close relationbetween CMC and the flamelet model where the same one-dimensionalmixing assumption is often adopted for species diffusion.The chemical source terms are often closed with the first order CMChypothesis ? that is, the conditional averages of the chemical reactionrate is approximated by evaluating the rate expressions with the con-ditional means of the species mass fraction vector, temperature anddensity:? |zeta approxequal ?parenleftBigT|zeta,Yi|zeta,rho|zetaparenrightBig. (2.62)The first order CMC hypothesis is justifiable as long as the fluctua-tions around the conditional averages are small. In CMC, conditionalmeans of reactive scalars are obtained by solving the above CMC trans-port equations. It has been successfully applied to study a wide varietyof problems [43,45,47?50].In cases with extinction and reignition, the first order CMC hypoth-esis might be imprecise and refinements to the closure hypothesis needto be applied. One solution is to use high order closure [51,52]. Kim etal. [46] used the second order conditional moment closure of CMC to27study turbulent non-premixed hydrocarbon flames with significant lo-cal extinction and reignition. Although the predictions of the interme-diates and fuel reaction rates are improved with the second-order clo-sure method, the requirement of modelling additional unclosed termsmake it less desirable. Another approach is to introduce a second con-ditioning variable. Cha et al. [53] suggested using the scalar dissipationrate as the second conditioning variable. The method predicts the ex-tinction reasonably well, but predicts the onset of reignition too early.Kronenbug [54] proposed to use a normalized enthalpy parameter to-gether with mixture fraction for double conditioning. The local extinc-tion and reignition can both be captured successfully.The traditional CMC methods are computationally expensive notonly because of the unclosed terms in CMC transport equations, butalso because of the need to directly solve CMC equations ? each con-ditioning variable introduces a new dimension to the system. Someresearchers have experimented with using coarse grids in the spatialdimensions to fix this problem.2.3.4 Conditional Source-term EstimationThe Conditional Source-term Estimation (CSE) method was first pro-posed by Bushe and Steiner [7]. It adopts the CMC first order closurehypothesis, but reduces the computational cost by eliminating the ad-ditional dimension introduced in CMC. Instead of directly solving theconditional averaged transport equations, the conditional means areobtained by inversion of an integral equation.28The integral equations for inversion are built based on the identity:tildewidef(x,t) =integraldisplay 10f|zeta tildewide(zeta; x,t) dzeta, (2.63)where f represents a reactive scalar such as density, temperature andspecies concentration, tildewide(zeta; x,t) is the density-weighted or Favre prob-ability density function (FPDF) defined (for LES implementation) by:tildewideP(zeta; x,t) = 1rho(x,t)integraldisplayDrho(xprime,t)delta[zeta -Z(xprime,t)]G(x -xprime)dxprime, (2.64)where zeta is the sample space of mixture fraction, delta is the delta functionand Z(xprime,t) is the instantaneous mixture fraction in the flow field. Thescalar tildewide(x,t) is obtained by solving a transport equation; the FPDF canbe obtained by solving a FPDF-transport equation, or by using a pre-sumed function form such as the beta-PDF.It is observed in experiments and DNS that the spatial gradientof f|zeta is small [7]. Therefore, spatial homogeneity in f|zeta might be as-sumed in many flows. For example, in an axisymmetric jet, variationin f|zeta on planes normal to the jet axis is usually very small. Gener-ally, it is assumed that some ensemble A may be selected in the flowfield with statistical homogeneity in the conditional averages, i.e., theconditional means at an individual point j can be approximated by anensemble conditional average:(f|zeta)j approxequalangbracketleftBigf|zetaangbracketrightBigAj element A, (2.65)where angbracketleftangbracketright denotes the average of ensemble A. Eq. 2.63 is then rewrittenas:tildewidef approxequalintegraldisplay 10angbracketleftBigf|zetaangbracketrightBigAtildewideP(zeta; xj,t) dzeta, j element A. (2.66)29Such an equation may be written for each point in the ensemble. Us-ing a numerical quadrature of N points in Z space, one can discretizethe integral above and construct a system of linear equations. The con-ditional means can then be obtained by inverting Eq. 2.66, which is aFredholm equation of the first kind [55].Inverting for the conditional means greatly reduces computationalcost; however, Eq. 2.66 is usually ill-posed. To address this issue, Steinerand Bushe [56] proposed to simultaneously minimize the residual ofEq. 2.66 and the derivative of the conditional means with respect tomixture fraction. Later Grout et al. [57] proposed to use the Tikhonovregularization method that introduces a priori knowledge of the solu-tion; the inversion problem is reformulated as a minimization prob-lem. In this approach, the target is to minimize the residual of Eq. 2.66as well as the deviation from the a priori information:min{||MangbracketleftY |zetaangbracketright -Y || +lambda||angbracketleftY |zetaangbracketright -angbracketleftY |zetaangbracketright0||}, (2.67)where M is the coefficient matrix of the discretized integral equationsrepresented by Eq. 2.66, lambda is a coefficient evaluated based on the char-acteristics of the inversion problem and angbracketleftY |zetaangbracketright0 is the a priori informa-tion introduced. In Grout?s formulation [58], lambda is evaluated with:lambda = tr(MasteriskmathM)/n, (2.68)where tr(MasteriskmathM) is the trace of the n by n matrix MasteriskmathM. Using Eq. 2.68for lambda leads to a well-behaved solution angbracketleftY |zetaangbracketright. Grout?s study also showsthat the CSE inversion method is not sensitive to lambda; increasing or de-creasing lambda by a factor of 10 makes little difference in the solution of con-ditional means. angbracketleftY |zetaangbracketright0 is chosen to be the solution from a previous time30step or iteration, which is proved to stabilize the inversion process.Truncation and re-scaling is then applied to eliminate non-physical so-lutions. The Tikhonov regularization method has been applied in CSEwith success [57?59]. It has also been found that optimizing the meshesin conditioning variable sample space and normalizing the inversionequations can improve the solutions.Compared to CMC, CSE is more computationally efficient. It doesnot require the constraining assumptions of the fast chemistry andlaminar flamelet models. Therefore, it can be used to study variousnon-premixed flames and possibly premixed flames as long as an ap-propriate conditioning variable and an effective inversion algorithmcan be selected.2.3.5 Including Chemistry in CSEThe closure method of CSE is desirable considering the reduced com-putational cost and precision of closure; however, to directly includedetailed chemistry in CSE is still too expensive. The system of inver-sion equations will be too large and the species transport equationsmight be too stiff to be handled by available ODE solvers. To addressthis issue, different methods can be combined with CSE model.The CSE method was first implemented in LES to study a pilotedmethane/air diffusion flame with a simplified two-step chemistry [7].The predictions of scalars such as temperature, mixture fraction andmajor species were consistent with the experimental data.Huang and Bushe [60] combined CSE with a chemistry reductiontechnique, the Trajectory Generated Low-Dimensional Manifold (TG-LDM) method using a purpose-built high-order finite-volume solver.31Details of the TGLDM method will be discussed later in this chapter.The Laminar Flamelet Decomposition(LFD) method [61] is an al-ternative to using reduced chemical kinetics, which incorporates de-tailed chemical mechanisms in combustion simulation using flameletlibraries. In the decomposition algorithm, conditional means are ap-proximated as a linear combination of basis functions:angbracketleftf|zetaangbracketright(t; A) approxequalNfsummationdisplayi=1aiThetai(zeta), (2.69)where Thetai(zeta) is the basis function.The basis functions should capture the relationship between theconditional means and the conditioning variable. For non-premixedflames, one option is to use the laminar flame solutions, where a libraryof basis functions (a flamelet library) can be generated by solving theunsteady laminar flamelet equation:rhopartialdiffpsiipartialdifft = rhoLeichi2partialdiff2psiipartialdiffZ2 + ?omegai, (2.70)with Z being the mixture fraction and chi(Z,t) being the scalar dis-sipation rate. Solutions at different flamelet times can be stored andused as the basis functions for decomposition. Alternatively, one canuse the solutions of steady laminar flamelet equations or a combina-tion of both [62, 63]. The choices of scalar dissipation rate chi includeusing constant values or using a presumed functional form such aschi(Z,t) = chi0(t)chi(Z). According to Peters [28], the functional depen-dence on mixture fraction chi(Z) can be defined using the solution of aone-dimensional mixing layer prescribed by:chi(Z) = e-2(erfc-1(2Z))2. (2.71)32The temporal function chi0(t) can be defined arbitrarily as:chi0(t) = c1e-t/c2 +c3, (2.72)where coefficients c1 and c2 are evaluated to approximate the behaviorof mean scalar dissipation rate on the stoichiometric interface in the re-active flow under study; c3 is given different constant values to incor-porate flamelets of high dissipation rates, which correspond to tempo-rary fluctuations in scalar dissipation observed in turbulent combus-tion.Now, Eq. 2.66 can be rewritten as:tildewidef(x,t) approxequalintegraldisplay 10parenlefttpparenleftbtNfsummationdisplayi=1aiThetai(zeta)parenrighttpparenrightbt tildewideP(x,t; zeta)dzeta. (2.73)With Eq. 2.73, one solves for the coefficient vector ai rather than solv-ing for the conditional averages directly. In Grout?s implementation,the equations of inversion were constructed using selected scalars thatwere found to best differentiate between flamelet solutions at variousstages of the ignition process in flamelet time, which included tem-perature and the mass fractions of CO and CH3OH [57]. Other scalarssuch as the reaction rates can be evaluated with Eq. 2.73, using solu-tions of the coefficient vector and the basis functions. This makes CSEwith laminar flamelet decomposition more similar to flamelet meth-ods, rather than CMC methods.The CSE-LFD method was first studied in an a priori test by Busheand Steiner [61]. Grout [58] implemented CSE-LFD in Fluent to studydifferent flames. He used a steady flamelet library in a RANS simu-lation of a piloted steady jet flame, Sandia flame D and an unsteadyflamelet library for autoigniting jets under engine-relevant conditions.33It was argued that the unsteady flamelet library better captured thetransient characteristics of autoignition processes. Wang and Bushe [62,63] implemented CSE-LFD in LES to study Sandia flame D. It wasfound that a mixed flamelet library consisting of both steady and un-steady flamelet solutions led to better predictions of the flame. Bothstudies showed that by optimizing the library of basis functions andthe combination of scalars for inversion, the precision of closure withCSE-LFD could be improved. However, Wang and Bushe [63] foundthat the LES results were sensitive to the composition of the flameletlibrary, which suggests that one must be very careful in constructingthe flamelet library if one is to use the CSE-LFD method.2.4 Premixed CombustionFor premixed combustion, the fuel and oxidizer are completely mixedbefore the combustion starts. Examples of premixed combustion in-clude lean-burn gas turbines and homogeneous charge spark ignitionengines. Different from non-premixed combustion, premixed combus-tion phenomena have a characteristic velocity scale- the laminar burn-ing velocity- and a characteristic length scale- the flame thickness. Thesenew parameters are used to define different regimes in premixed com-bustion. Regime diagrams have been proposed by numerous research-ers [21]. The diagrams are usually based on the comparison betweenthe turbulent and chemical length scales and time scales- as was thecase in non-premixed combustion. Regimes are defined based on dif-ferent non-dimensional parameters.For example, the Borghi diagram [64] uses a velocity ratio uprime/uref34Da = 11-1log( lllref)Well Stirred Reactor,log(uprimeuref)Island Formation andWrinkled Laminar Flame,Flamelet RegimeTorn Flame Fronts,Perturbed FlameletsThickened FlameReT = 1Laminar Flame FrontsKa = 11Figure 2.3: Borghi diagram for premixed combustionand a length ratio ll/lref. In a classic Borghi diagram, uprime is the turbu-lent intensity, uref = sL is the laminar flame velocity, ll is the integralscale and lref = deltaL is the laminar flame thickness. As shown in Fig. 2.3,at Ka < 1, where island formation and wrinkled laminar flames areobserved, the flamelet regime is defined. The laminar flamelet model(to be discussed later) works well for flames within this regime. AtDa < 1, the thickened flame regime is defined. Here small scale eddiespenetrate the reaction zones and the flamelet concept is inapplicable.In between is the torn flame front regime. Within this regime, the tur-bulence interacts with the flame without completely overwhelming it.Peters [30] proposed a regime diagram for premixed flames where,for scaling purposes, it is assumed that the diffusivities of all scalarsare equal and the Schmidt number is unity. The length and time scales35of the flame are defined by:lF = DsL, (2.74)tF = DsL2. (2.75)The turbulent Damk?hler and Karlovitz number are defined as:Da = sLlupsilonprimelF, (2.76)Ka = tFteta= l2Feta2 =upsilon2etas2L , (2.77)where teta, eta and veta are the Kolmogorov time, length and velocity scalesrespectively; upsilonprime is the turbulent intensity. A second Karlovitz numberis defined by:Kadelta = l2deltaeta2 = delta2Ka, (2.78)where ldelta is the thickness of the inner layer between the chemically in-ert preheat zone and the oxidation layer. The inner layer is where thereaction process of a premixed flame is kept alive. The thickness ldelta isevaluated using ldelta = deltalF , where the value of delta is typically around onetenth.Different regimes can be defined with the above characteristic scalesas shown in Fig. 2.4. The ratios l/lF and upsilonprime/sL parameterize the regimediagram, which are related as:upsilonprimesL = Reparenleftbigg llFparenrightbigg-1= Ka2/3parenleftbigg llFparenrightbigg1/3. (2.79)The line Re = 1 separates turbulent flames from laminar flames.The lines upsilonprime = sL, Ka = 1 and Kadelta = 1 divide premixed turbulent3610-1 100 101 102 103 10410-1100101102103l/lFupsilonprime/sLRe = 1eta = ldeltaKadelta = 1eta = lFKa = 1brokenreactionzonesthinreactionzonescorrugated flameletswrinkled flameletslaminarflameletsFigure 2.4: Regimes in premixed turbulent combustion.flames into wrinkled and corrugated flamelet regimes, the thin reactionzones and the broken reaction zones regimes.In the wrinkled flamelets regime, even large eddies are not largeenough to compete with the laminar flame propagation, which con-sequently dominates the combustion process. In the broken reactionzones, the inner layer is disturbed by the smaller eddies. Enhancedheat loss may cause temperature decrease and loss of radicals; as a re-sult, extinction might occur. For practical reasons, these two regimesare of less interest in industrial applications.The corrugated flamelets regime features Ka < 1 and upsilonprime greaterequal sL greaterequal upsiloneta.In this regime, small eddies interact with the advancing laminar flame.The thin reaction zones regime is characterized by Re > 1, Kaeta < 1 andKa > 1. The upper limit of this regime is defined by Kadelta = 1 and in37most cases Ka = 100 can be used as an approximation. In this regime,small eddies may increase scalar mixing, but cannot penetrate the in-ner layer. The corrugated flamelet and thin reaction zones regimes areseparated by the line Ka = 1. Under such circumstances, the thick-ness of the flame is equal to the Kolmogorov length scale; tF = teta andupsiloneta = sL can be derived from Eq. 2.77. Many premixed combustionmodels are developed based on the characteristics of different regimes.2.4.1 Eddy-Break-UpThe Eddy-Break-Up (EBU) model is based on a phenomenological anal-ysis and is primarily used in premixed combustion with high Re andDa numbers. It is argued [65] that in cases where mixing determinesthe rate of chemical reactions, the cascade process from the integral tomolecular scales of turbulent mixing also controls the reaction process.In case of excessive oxidizer, the mean reaction rate of fuel can be eval-uated using the turbulent mixing time taut:? F = -rhoCEBU (tildewidestY primeprimeF2)1/2taut , (2.80)where tildewidestprimeprimeF 2 is the variance of the fuel mass fraction; CEBU is the Eddy-Break-Up constant (on the order of unity); taut is the turbulent mixingtime estimated from the turbulent kinetic energy k and its dissipationrate epsilon by using taut = k/epsilon.A progress variable c for the chemical reaction is defined such thatc = 0 in unburnt gas and c = 1 in fully burnt mixture. Examples of cinclude a normalized temperature or a normalized mass fraction:c = T -TuTb -Tu, (2.81)38c = YPYP,b, (2.82)where T, Tu and Tb are the temperatures of the local mixture, unburntgas and fully burnt mixture, YP and YP,b are the mass fractions of prod-ucts in the local and fully burnt mixture. The reaction rate may be cal-culated using a progress variable c as:? = -rhoCEBUradicalBigtildewidercprimeprime2taut . (2.83)Eq. 2.83 leads to inconsistencies at tildewide = 0 and tildewide = 1 when the flamefront is infinitely thin (i.e. c = 0 or c = 1). Under such conditions, c2 = cis valid and tildewiderprimeprime2 is given by:rhotildewiderprimeprime2 = rho(c -tildewide)2 = rho( tildewide2 -tildewide2) = rhotildewide(1 -tildewide). (2.84)Substituting Eq. 2.84 in the reaction rate equation, one can find d ? /dtildewideis infinite at tildewide = 0 and tildewide = 1. A corrected version was proposed byBorghi [8]:? = CEBUrhotildewidec(1 -tildewidec)taut. (2.85)The EBU model was modified by Magnussen and Hjertager [66],leading to the Eddy Dissipation Model (EDM). In EDM, the mean massfraction of the deficient reactant is used and the minimum of the reac-tion rates for fuel, oxidizer and products is used to evaluate the meanchemical source terms.The EBU model replaces the chemical time scale of an assumedone-step chemical reaction with the turbulent time scale taut and calcu-lates the reaction rates simply as a function of some known quantities.Because of its simplicity, it has been widely used in commercial CFDpackages. However, because the chemical kinetics are eliminated from390 1alpha(x, t)beta(x, t)Figure 2.5: Presumed probability density function in premixed turbu-lent combustion, prescribed by the Bray-Moss-Libby model.the modelled reaction rates, the model represents the fast chemistrylimit only. The model constant CEBU sometimes needs to be tuned ac-cording to specific problems. It is also observed that the model tendsto over-estimate the reaction rates.2.4.2 The Bray-Moss-Libby ModelThe Bray-Moss-Libby (BML) [67] model implements the classical flameletconcept 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 flameletconcept. The presumed probability density function of the progressvariable c at a given time and location (x,t) is modeled with:P(casteriskmath,x,t) = alpha(x,t)delta(casteriskmath) +beta(x,t)delta(1 -casteriskmath) +gamma(x,t)f(casteriskmath,x,t). (2.86)The presumed PDF prescribed by Eq. 2.86 is presented in Fig. 2.5. Itdescribes the mixture at (x,t) as a sum of unburnt, fully burnt andburning gases, whose contributions are represented by the three terms40on the r.h.s. of Eq. 2.86 respectively. The coefficients alpha, beta and gamma denotethe probabilities of having fresh, completely burnt and reacting mix-tures at (x,t). The Dirac delta functions delta(casteriskmath) and delta(1 -casteriskmath) correspondto fresh mixture at c = 0 and complete combustion at c = 1 respec-tively. The functional form f(casteriskmath,x,t) for the probability of burning gascan be defined such that:integraldisplay 10f(casteriskmath,x,t)dcasteriskmath = 1 (2.87)with f(0) = f(1) = 0.The coefficients are subject to the normalization of the PDF:integraldisplay 10P(casteriskmath,x,t)dcasteriskmath = 1, (2.88)from which one can derive:alpha +beta +gamma = 1. (2.89)The model is further simplified in cases where Re >> Da >> 1.Under such circumstances, the premixed turbulent flame is extremelythin and the combustion is mainly controlled by the turbulent trans-port. Consequently, gamma << 1, gamma << alpha and gamma << beta are valid. The PDF isreduced to a quasi-bimodal form:P(casteriskmath,x,t) = alpha(x,t)delta(casteriskmath) +beta(x,t)delta(1 -casteriskmath). (2.90)The coefficients alpha and beta can be easily evaluated using the normaliza-tion of PDF prescribed by Eq. 2.88 and averages of the progress vari-able tildewide:rhotildewide=integraldisplay 10rhocasteriskmathP(casteriskmath,x,t)dcasteriskmath, (2.91)41from which one can derive:beta = rhotildewidecrhob, (2.92)alpha = 1 -beta. (2.93)Here rhob is the density of fully burnt mixture.The reaction heat release factor tau = Tb/Tu -1 is often introduced inthe BML analysis. At constant pressure P, the following equations canbe derived by assuming perfect gas behavior:tau = TbTu-1 = rhourhob-1, (2.94)rhou = (1 +tau)rhob = rho(1 +tautildewide). (2.95)Substituting these in Eq. 2.93 and Eq. 2.92 and one obtains:alpha = 1 -tildewidec1 +tautildewidec, (2.96)beta = (1 +tau)tildewidec1 +tautildewidec . (2.97)The heat release factor tau is determined by the specific chemical mec-hanism and is therefore fixed. The PDF P(c) depends solely on tildewide andtau.The governing equation of the progress variable is:partialdiffrhocpartialdifft +triangleinv? (rhouc) = triangleinv? (rhoD triangleinvc) + ?omega (2.98)Under the assumptions of the BML model, where the progress variablec is equal to zero or unity, an equation can be derived for reaction rates:2rhoD triangleinvc ? triangleinvc = 2c ? - ? (2.99)After averaging, Eq. 2.99 can be rewritten as:2rhoD triangleinvc ? triangleinvc = (2cm -1) ? , (2.100)42where cm =R 10 c ?omegaf(c)dcR 10 ?omegaf(c)dcis a progress variable characterizing the chemicalkinetics. The mean reaction rate is:? = 2 rhochi2cm -1, (2.101)where rhochi is the scalar dissipation rate of the progress variable c:rhochi = rhotildewide = rhoD triangleinvc ? triangleinvc = rhoD partialdiffcpartialdiffxipartialdiffcpartialdiffxi . (2.102)The scalar dissipation rate tildewide can be evaluated by solving a transportequation [69] or by modelling:rhochi approxequal rhocprimeprime2taut , (2.103)where taut is the turbulent time scale. Assuming the flame front is in-finitely thin (i.e. c = 0 or c = 1), c2 = c is valid and tildewiderprimeprime2 is given as:rhotildewiderprimeprime2 = rho(c -tildewide)2 = rho( tildewide2 -tildewide2) = rhotildewide(1 -tildewide). (2.104)Eq. 2.101 is rewritten as:? = 22cm -1rhotildewide(1 -tildewide)taut , (2.105)which is consistent with the EBU model represented by Eq. 2.85. Thelink between the BML and EBU models can be also established usingflamelet analysis [8].Using the progress variable, the average of any quantity may beevaluated by conditional averaging. For example, in the bimodal formof BML, the Favre average of a scalar Q can be calculated as:rho tildewide = rhoQ =integraldisplay 10rhoQP(c)dc = alpharhouQu +betarhobQb, (2.106)where Qu and Qb are the conditional means in the fresh gases and fullyburnt mixture respectively.432.4.3 Models Based on Surface Area EstimationAnother main category of turbulent premixed combustion models isbased on the concept of flame surface density ?. By definition, ? is theflame surface area per unit volume. The mean reaction rate of species iis modeled with:? i = ? i?, (2.107)where ? i is the local reaction rate per unit of flame area integrated inthe direction normal to the flame surface. It is related to the propertiesof the local flame structure and can be estimated using a prototypelaminar flame in a pre-processing fashion. This method decouples thechemical kinetics described by ? i and the turbulence/flame interactiondescribed by ?.Bray et al. suggested that the mean chemical source term must beproportional to the flamelet crossing frequency and proposed a mo-del [68]:omegac = rhous0LI0?, (2.108)where s0L is the laminar burning velocity of the unstretched flame andI0 is the stretch factor.The flame surface density ? can be calculated using an algebraicexpression such as [68]:? approxequal c(1 -c)?Ly, (2.109)where ?y is the crossing length scale; alternatively, it can be obtainedby solving a transport equation. For example, an equation was pro-posed by Trouv? and Poinsot [70]:partialdiff?partialdifft +triangleinv? (tildewideupsilon?) = triangleinv? (Dt triangleinv?) +C1epsilonk? -C2sL?21 -c. (2.110)44The terms on the left hand side of the equation represent the local rateof change and convection, and the terms on the right hand side rep-resent turbulent diffusion, production by flame stretch and flame sur-face annihilation respectively. Other forms of the balance equation for? have been developed in various premixed combustion models [8]where the source and sink terms are modelled differently.2.4.4 Models Based on the G-equationThe premixed combustion models discussed previously involve theuse of a progress variable. An alternative approach is the G-equationmethod, where a non-reacting scalar G is introduced. An iso-surfaceat a critical value G0 divides the flow field into regions of burnt gasesand unburnt gases. The value of G0 is fixed for a particular combustionevent. A governing equation for G was proposed by Williams [71] as:partialdiffGpartialdifft +v ? triangleinvG = sL| triangleinvG| (2.111)This G-equation is applicable to thin flames that propagate at a well-defined burning velocity and is therefore suited for describing pre-mixed turbulent combustion in the corrugated flamelets regime. Pe-ters [21] formulated a G-equation in the thin reaction zones regime:partialdiffGpartialdifft +v ? triangleinvG =bracketleftbiggtriangleinv? (rhoD triangleinvT) +omegaTrho| triangleinvT|bracketrightbigg| triangleinvG|. (2.112)A common level set equation was formulated for both regimes:partialdiffGpartialdifftasteriskmath +vasteriskmath ? triangleinvasteriskmathG = sL,snueta | triangleinvasteriskmath G| - Dnu kappaasteriskmath| triangleinvasteriskmath G|, (2.113)where the independent variables and curvature are normalized by theKolmogorov length, time and velocity scales:tasteriskmath = t/teta, xasteriskmath = s/eta, vasteriskmath = v/veta, kappaasteriskmath = etakappa, triangleinvasteriskmath = eta triangleinv. (2.114)45Analogous to the role of the conserved scalar Z in non-premixedflames, G and/or its variance are often used in premixed combustionmodels to evaluate various properties of premixed flames. For exam-ple [72], a flamelet library may be generated where the species, tem-perature and density of the prototype flames are stored as functions ofdistance from the G = 0 surface and strain rate; it is then incorporatedin the combustion simulation using a presumed probability densityfunction.2.5 Reduction of Detailed Reaction MechanismA detailed chemical reaction system generally includes many speciesand elementary reactions. To simulate turbulent reactive flows at af-fordable computational cost, it is usually necessary to reduce the fullmechanisms. The most commonly used methods include the Quasi-steady State Approximation [73] (QSSA) and the partial equilibriummethods [74]. In QSSA, analytical approximations of reaction rates arederived by assuming equal consumption and formation rates for cer-tain species; in partial equilibrium methods, they are derived by as-suming small net reaction rates of fast reactions. These methods havebeen applied to study premixed and non-premixed flames with suc-cess [75?77].Recently, new techniques have been developed and applied in tur-bulent combustion simulation, which typically involve computationwith full chemistry in a pre-processing step and tabulation. Examplesinclude the Intrinsic Low-Dimensional Manifold (ILDM) [78] method,Invariant Constrained-equilibrium Edge Pre-Image Curve (ICE-PIC)46method [79] and methods based on flamelet concepts [80,81].2.5.1 Intrinsic Low-Dimensional ManifoldThe Intrinsic Low-Dimensional Manifold (ILDM) method was first pro-posed by Maas and Pope [78]. The method is based on the fact that awide 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 ontomanifolds of fewer dimensions on longer time-scales such as those rel-evant to mixing in a turbulent non-premixed flame. The ILDM methodis used to find such a manifold, which is parameterized by progressvariables and stored as a look-up table in combustion simulation.The first step is to define the low-dimensional manifold by means ofa time scale analysis. The detailed chemical kinetic system is describedwith:partialdiffpsipartialdifft = F(psi) +Xi(psi,nablapsi,deltapsi), (2.115)where psi = (h,p,psi1,psi2,...,psins)T is the scalar vector consisting of en-thalpy, pressure and all ns species in the detailed chemistry, psik = YkWkisthe specific mole number of species k ([mol/kg]), F(Psi) is the chemicalsource-term and Xi denotes physical processes such as convection, dif-fusion, etc. For a simple homogeneous, adiabatic and isobaric system,Eq. 2.115 reduces to:partialdiffpsipartialdifft = F(psi). (2.116)The time scales of the system can be separated through a Schur decom-position analysis:J = V?tildewide, (2.117)47where J is the Jacobian matrix of the chemical source term F(psi) withJij = partialdiffFipartialdiffpsij , 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. Theeigenvalues lambdak in ?are sorted in descending order such that the eigen-value representing the fastest process is at the bottom; V and tildewide aresorted accordingly. The characteristic chemical time scales of differentprocesses are evaluated with:tauk = lambda-1k . (2.118)Different time scales are grouped into nf fast time scales and nc slowtime scales, where the fast processes quickly attenuate to quasi-steadystates and the reaction system converges onto an nc-dimensional man-ifold defined by:tildewiderVfF(psi) = 0, (2.119)where tildewiderf consists of the last n -nc rows taken from tildewide, correspondingto the fast processes. Such a subspace can be found using Eq.2.119, to-gether with additional equations representing system constraints andthose of selected progress variables. The final form of the ILDM equa-tion is:G(psi,tau) =parenlefttpparenleftbttildewiderVfF(psi)P(psi,tau)parenrighttpparenrightbt = 0, (2.120)where P contains the parametric equations such as the adiabaticity,constant pressure and element conservation.Once the low dimensional manifold is defined, the detailed chemi-cal kinetic system can be parametrized by a small number of variablestheta = (theta1,theta2,...,thetaN) (N < n). The n-dimensional system can be projected48onto the new manifold with the governing equations:partialdiffthetapartialdifft = S(theta) +G(psi(theta),nablapsi(theta),deltapsi(theta)). (2.121)To apply ILDM in combustion simulations, the low-dimensionalmanifolds are usually calculated beforehand and stored as look-up ta-bles. These tables can be retrieved during the combustion simulationand the themokinetic properties are evaluated using the progress vari-ables solved in the flow computation. Since the calculation of ILDM isonly based on the analysis of the reaction system, they can be used invarious flame configurations as long as the element composition, pres-sures and enthalpies are a subset of the conditions defined in ILDMcomputation. Another attractive feature of ILDM is that the mathe-matical formulation of ILDM provides an automatic error control ofthe simplification scheme.The ILDM method has been used to study laminar and turbulentflames of various configurations [82,83]. However, the ILDM methodis based on the assumption that the fast chemical time scales can bewell separated from the slow ones, which is not necessarily true for anychemical system. To solve this problem, Nafe and Maas [84] introducedthe Slow Invariant Manifolds (SIM) methods, where the original ILDMcomputation was used as an initial guess and a post-processing stepwas applied such that only movements perpendicular to the manifoldwere allowed. The resultant manifold was found to be more accuratethan the original ILDM. Another issue is that at relatively low temper-atures, the diffusion and chemical time scales are of the same order.Therefore, the diffusion processes need to be included in the computa-tion of manifolds.49Some researchers have combined ILDM with different ways of ma-nifold generation to solve this problem. For example, a related ap-proach called Flame Prolongation of ILDM (FPI) was studied by Gic-quel et al. [85]. In this method, the low temperature regime in ILDMwas constructed by using the solutions of laminar premixed flames. Oi-jen et al. [81] introduced the Flamelet-Generated Manifold (FGM) me-thod based on the similarity of flame structure between a propagatingthree-dimensional flame and a one-dimensional flame. A combinationof ILDM and FGM methods called Phase-space ILDM (PS-ILDM) wasstudied by Bongers et al. [86]. In the PS-ILDM method, the ILDM strat-egy was used to decouple the fast and slow processes in the solutionof flamelet equations.Besides the efforts to improve the quality of manifolds, the tabula-tion of manifolds has been studied to reduce the requirements of stor-age. For instance, a method called in situ adaptive tabulation (ISAT)was combined with ILDM methods by Pope [87], which built the ma-nifold as needed during simulation.2.5.2 Trajectory Generated Low-Dimensional ManifoldAn alternative technique of generating low dimensional manifolds isthe Trajectory Generated Low-Dimensional Manifold (TGLDM) met-hod, originally proposed by Pope and Maas [88]. The TGLDM methodinvolves solving for physically realizable initial states of the systemand integrating Eq. 2.115 to evolve the reaction trajectories. Each tra-jectory consists of solutions in the composition space starting from agiven initial state of the system and evolving toward the chemical equi-librium state. Theoretically, if a viable low dimension manifold does50exist for a chemical system, the trajectories with initial states on themanifold will stay on the manifold. In this sense, such trajectories areequivalent to a SIM. An important advantage over ILDM is that the al-gorithm used in TGLDM guarantees a converged solution. The globaloptimization problem in ILDM is reduced to a local optimization prob-lem of defining initial states for trajectories. As in ILDM, TGLDM man-ifolds can be calculated and tabulated in pre-processing and retrievedlater during a combustion simulation.To select the initial states, Pope and Maas [88] adopted the extreme-value-of-major-species method, in which linear equations were solvedbased on the conservation of elements. For example, major species inthe reaction systems may be selected and the mass conservation of el-ements can be written for these species:nssummationdisplayi=1Yc;j,iYi = Ye;j j = 1,...,ne, (2.122)where ns is the number of selected species, ne is the number of ele-ments, Yc;j,i is the mass fraction of element j in species i, Yi is the massfraction of species i and Ye;j is the mass fraction of element j in thefresh mixture. By solving such equations, initial states can be foundwith extreme values of mass fractions at the boundary of the composi-tion space.Huang and Bushe [60] used a slightly different approach. To con-struct an np dimensional manifold, constraints on progress variablesare applied in addition to the element conservation equations:Pk,iYi = Yp;k; i = 1...ns,k = (2.123)where np is the number of progress variables that parameterize the low51Figure 2.6: Comparing one dimensional manifold generated by ILDMand TGLDM.dimensional manifold, Pk,i is the coefficient for the parametric equa-tions, Yp;k is the assigned value of progress variable k. A natural choiceof the constraints is to use the equilibrium state in the subspace, whichcan be calculated with the method of Lagrange multipliers [89,90] orGibbs function continuation [91]. Huang [83] used this method to gen-erate a one dimensional manifold which was compared to an ILDMand good agreement was observed. His analysis also showed that theTGLDM created with constrained equilibrium better approximated de-tailed chemistry than ILDM methods under moderate perturbation.The TGLDM generated was applied in three different methane/air re-action systems with success including an unstrained, premixed lami-nar flame, a non-premixed laminar flame and a perfectly stirred reac-tor.52Later Huang and Bushe [60] combined TGLDM with CSE to studya methane jet flame where fuel was injected into a shock-tube underengine-relevant conditions. In this study, a two dimensional manifoldparameterized by YCO2 and YH2O was constructed. The manifolds werestored using Delaunay triangulation [92] and retrieved by surface in-terpolation over the two-dimensional unstructured grid [92]. Five ma-jor species (O2, CH4, CO, CO2 and H2O) were used to construct theequations for the initial states, including element conservations of C, Hand O, and constrained equilibrium equations of the two progress vari-ables. The predictions of autoignition characteristics of the jets werefound to be promising.In another work on autoigniting jets, Wang et al. [62,93] discussedthat at low pre-combustion temperatures, the ignition delay time be-comes exceedingly long such that it is extremely difficult for the ODEsolver to obtain a solution for the trajectory. In previous studies, thisproblem was circumvented by mixing the unburned mixture with asmall amount of mixture at equilibrium. The initial state point wastherefore shifted slightly into a state of a higher equivalent temperatureand the stiffness of the system was reduced. The resulting trajectory ch-emistry does not describe autoignition correctly. To solve this problem,the Stochastic Particle Model (SPM) [94?96] was used to generate theautoignition trajectory at low temperatures. In Wang?s study, the SPMwas used to calculate several realizations of trajectories, which werefiltered, averaged, and combined with TGLDM trajectories generatedusing a traditional continuum approach. The SPM-extended TGLDMwas then used (in the same high-order code as Huang and Bushe [60])53to study a methane jet and the fluctuations in ignition delays due torandomness 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 approachas discussed in Huang and Bushe?s work [60] was adopted to closethe chemical source terms in LES. The predictions of major species,temperature field as well as the conditional averages agreed well withthe experiment.2.6 SummaryTurbulent combustion simulation is challenging because of the cou-pling of mechanisms such as fluid dynamics, chemical reactions, radi-ation and phase changes, etc. Depending on the characteristics of thecombustion phenomena, different methods can be used to approachthis problem. Among these methods, Conditional Source-term Estima-tion provides a closure method with a good balance between computa-tional cost and precision. In previous research, CSE-TGLDM and CSE-LFD have been successfully implemented in different studies. Morework needs to be done to explore the applications of CSE:? Although CSE-TGLDM and CSE-LFD have been implementedin individual studies, these two methods have never been com-pared head-to-head. Such a study would help to clarify the un-derlying differences between these two methods and comparethe performance of them on a fair basis.? So far, the CSE applications have never accounted for the effectsof radiation heat loss, which is important especially in the pre-54diction of NO formation and might not be neglected in luminousflames.? In previous work, the application of CSE is limited to modellingnon-premixed flames. The possibility of using CSE in premixedflames is worth exploring.These research problems are studied in this thesis work and will bediscussed in the following chapters.55Chapter 3Ignition Delay of Methane Jets withConditional Source-term Estimation3.1 IntroductionConditional source-term estimation (CSE) has been used to study theautoignition process in non-premixed jet flames of methane and me-thane mixed with additives, where conditional means of temperatureand species are used to close chemical source-terms. CSE is combinedwith a Trajectory Generated Low-Dimensional Manifold (TGLDM) me-thod and then with a Laminar Flamelet Decomposition (LFD) methodto reduce the complexity of detailed chemistry. The predictions of ig-nition delay time are compared with measurements of shock-tube ex-periments.3.2 Model Formulation3.2.1 Numerical SchemeReynolds-Averaged Navier-Stokes (RANS) simulation of autoignitingjets has been performed in an axisymmetric computational domain us-ing the standard k - epsilon1 [14] turbulence model. Besides continuity andmomentum equations, Reynolds averaged transport equations [8] of56energy, selected species, the mean and variance of mixture fraction aresolved in a cylindrical coordinate.In the current study, body force is neglected in the momentum equa-tions and radiation heat loss is neglected in the energy equation. It isassumed that the Schmidt number Sc is 0.7; the Lewis number Le isunity and the diffusivities of individual gaseous species are not signif-icantly different. The chemical source terms in the Reynolds averagedenergy and species equations are closed with CSE. For details of theCSE-LFD and CSE-TGLDM methods, please refer to Chapter 2, Sec-tion 2.3.5 and Results and DiscussionTo validate and compare the CSE-TGLDM and CSE-LFD methods de-scribed above, these methods have been implemented in OpenFOAM[97] to simulate a variety of shock-tube experiments where the shockis used to compress air and fuel is injected along the centerline of theshock-tube with an injector mounted on the end-plate. The resultingjet ignites and the ignition delay time is measured using a high-speedcamera, imaging the entire jet in the visual spectrum. Ignition is de-tected by identifying the first visible flame kernel that subsequentlygrows into a full jet flame. The experimental techniques used are de-scribed by Sullivan et al. [98].3.3.1 Simulation of a Non-reactive JetA cold jet was simulated to validate solvers developed using Open-FOAM. The jet under discussion was studied experimentally and nu-merically by Ouellette [99] . In his experiments, cold methane jets were57injected at high Reynolds number through a nozzle and mixed withcold air in the shock-tube. Because of the high pressure ratio upstreamof the chamber and the high exit pressure at the nozzle, underexpandedjets were observed. The penetration length was measured and aver-aged over three independent tests.A scaling model was proposed by Ouellette [99] based on the vortex-quasi-steady-state assumption and the mass/momentum entrainmentmeasurement from Ricou and Spalding [100]. According to the scalingmodel, the penetration of a self-similar transient jet can be modeled as:parenleftBiggZtdnradicalBigrhoirho0parenrightBigg= Gradicalbiggpi4parenleftBiggU0tdnradicalBigrhoirho0parenrightBigg0.5, (3.1)where Zt is the penetration length of the jet, dn is the nozzle diameter, rhoiand rho0 are the densities of the gas at the nozzle exit and in the ambientair respectively, U0 is the jet exit velocity, and G is an empirical penetra-tion number with a value around 3. An equivalent diameter definedby deq = dnradicalBigrhoirho0 is introduced to non-dimensionalize the penetrationlength Zt in post-processing.To test the CFD solvers developed with OpenFOAM, the predictedpenetration length of the jet is compared to experimental data as wellas predictions of Fluent [101] and the scaling models. The cold jet un-der study is injected from a 0.5 mm nozzle at the injection pressure2.285 MPa and temperature 300 K into the shock-tube filled with coldair at the pressure of 1.494 MPa.Simulations of the non-reactive jet have been performed in an axi-symmetric computational domain over a structured mesh. The meshis 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 in58the solution. The inflow boundary conditions are defined by assuminga polytropic expansion of the fuel jet at the injector tip. The penetrationlength is evaluated using the method proposed by Ouellette [102]. Thecold jet is also simulated using the commercial software Fluent. Forcomparison, the configurations of the OpenFOAM and Fluent compu-tations are set as close to each other as possible.Fig. 3.1 shows the non-dimensionalized penetration length of themethane jet as a function of time. The predictions of the OpenFOAMand Fluent simulations are compared with the experimental data andthe prediction of the scaling model. The experimental data contain er-ror because of the injection delay between issuing the command to in-ject and the actual onset of injection. Therefore, it is more important tocompare the simulation results to the scaling model prediction, whichis calculated in two steps. Firstly, the scaling model is rewritten to in-clude an injection delay time and G is estimated by curve fitting thenew function form to the experimental data; secondly, Eq. 3.1 is usedto calculate the penetration length as a function of time (without theinjection delay).The figure shows that the simulation results are generally consis-tent with the measurements. The profiles of the OpenFOAM and Flu-ent predictions are close to each other, while the OpenFOAM and thescaling model calculation are barely distinguishable. The comparisonshows that for the two dimensional RANS simulation of this cold jet,the performance of the solver developed with OpenFOAM is at leastas good as that of Fluent. On the other hand, the OpenFOAM solvercompletes the computation of an injection duration of around 3 ms in590 0.5 1 1.5 2 2.5 3x 10-30102030405060708090t(s)Zt/dnFigure 3.1: Penetration length of a cold jet. square: experimental data; ?-: simulation with OpenFOAM, --: simulation with Fluent; ? ? ? ? ?:prediction of the scaling model.60Table 3.1: Autoigniting jets with various fuel compositionParametera Pi (bar) ti (ms) Fuel T0 (K)Jet A 120 1.0 100%CH4 1200-1400Jet B 120 1.0 90%CH4 +10%C2H6 1200-1400Jet C 120 1.0 80%CH4 +20%N2 1200-1400Jet D 60-120 1.0 100%CH4 1300Jet E 120 1.0-2.5 100%CH4 1300a Pi: injection pressure; ti: injection duration; T0: pre-combustiontemperature; fuel composition is given by volume.about half the time that Fluent takes. The simulations of autoignitingjets hereafter are performed using OpenFOAM.3.3.2 Simulations with CSE-TGLDM and CSE-LFDThe autoignition of jets of pure methane and of methane mixed withadditives has been simulated with OpenFOAM. For all the cases understudy, the nozzle diameter is 0.28 mm and cold fuel jets at around 300K are injected into the shock-tube with a pre-combustion pressure of 30bar. Other important characteristics of the investigated jets are listed inTable 3.1 The experimental results are provided by Wu et al. [103,104].To compare the CSE-TGLDM and CSE-LFD methods, the two dif-ferent CSE approaches are implemented using the same code, numer-ical scheme, grid as well as the same initial and boundary conditions.A detailed chemical mechanism including 71 species and 379 reac-tions [105] is used in generating the TGLDM and LFD flamelet libraries.This mechanism is tuned to premixed shock-tube experiments.In the simulation with the CSE-TGLDM model, the mass fractions61of YCO2 and YH2O are selected as progress variables because of their rel-atively long formation time scales. As an example, Fig. 3.2 shows theTGLDM manifold and its Delaunay triangulation for a methane/ethanejet in YCO2 and YH2O plane, generated at the stoichiometric mixture frac-tion, with a pre-combustion temperature of 1200 K and at a constantpressure of 30 bar.In the simulation with the CSE-LFD model, unsteady flamelet so-lutions are used as basis functions for decomposition. For example,Fig. 3.3 shows flamelet solutions of a methane/ethane jet extracted atvarious flamelet times (corresponding to different flamelet numbers),with a pre-combustion temperature of 1200 K and at a constant pres-sure of 30 bar. The flamelets shown are all taken prior to ignition in thelaminar flamelet time. It is observed that before ignition, the temper-ature change is not sufficient to differentiate between flamelets; whilethe combination of T, YCO and YCH3OH might be used to detect differ-ent stages of the autoigniting process. Therefore, temperature and themass fractions of CO and CH3OH are selected as scalars to constructthe inversion equations of CSE in this study.It has been found in earlier work [57] that the simulated ignition de-lay time is largely insensitive to the definition used: the state changesso abruptly when ignition occurs that many different variables all changedramatically and almost simultaneously, such that it is possible to de-tect ignition by tracking any one of several different variables. For con-venience, the ignition delay is defined based on a rise in the conditionalmean temperature. Fig. 3.4 shows the maximum increase in the condi-tional means of temperature in the reaction field plotted against time620 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.1600. TGLDM manifold0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.1600. Delaunay triangulationFigure 3.2: TGLDM for methane/ethane jet at stoichiometric mixturefraction, T0 = 1200 K, P0 = 30 bar.630 0.2 0.4 0.6 0.8 1200400600800100012001400ZT(K)0 0.2 0.4 0.6 0.8 10123456 x 1011ZomegaT(erg/g)0 0.2 0.4 0.6 0.8 100.511.522.533.54 x 10-5ZYCO0 0.2 0.4 0.6 0.8 100.511.522.533.544.5 x 10-6ZYCH3OHFigure 3.3: LFD library of methane/ethane jet, T0 = 1200 K, P0 = 30 bar.?-: flamelet 5; ? ? ? ? ?: flamelet 10; - ? - line: flamelet 20; --: flamelet25.640 0.2 0.4 0.6 0.8 1 1.2 1.4x 10-30510152025303540t(s)Max.IncreaseinangbracketleftT|zetaangbracketrightFigure 3.4: Maximum increase of angbracketleftT|zetaangbracketright predicted with CSE-TGLDM forseveral different initial temperatures with Jet B. ?-: T0 = 1400 K; --:T0 = 1350 K; ? ? ? ? ?: T0 = 1300 K; -? -: T0 = 1250 K; openbullet: T0 = different pre-combustion temperatures for the methane/ethane jetcalculated by CSE-TGLDM model. It can be seen that this maximumvalue is almost zero early in the simulation and later increases dramat-ically. Similar observations are found with CSE-LFD calculation and inthe simulations of pure methane jets. In this study, the ignition delaytime is calculated as the time when the dramatic change in conditionalmean temperature begins.The autoignition of pure methane and methane/ethane jets is sim-ulated using both CSE-TGLDM and CSE-LFD methods. Predictions ofignition delay are compared with experiments, as shown in Fig. 3.5.Predictions with both methods appear to agree with the experimentalresults despite (if not because of) the scatter in the experimental data.65The scatter observed in the measurements can be attributed to fluctua-tions in turbulent mixing as well as to different realizations of chemicalreaction paths - even in homogeneous, quiescent ignition, it has beenfound that the ignition delay time is a random variable [93, 96]. Al-though the present models capture the trend of ignition delay, furtherwork is needed to account for these fluctuations.The ignition delay predicted with CSE-TGLDM is somewhat betterthan that predicted with with CSE-LFD ? the predictions of CSE-LFDbeing universally shorter than those with CSE-TGLDM. This discrep-ancy is due to the different techniques of incorporating detailed ch-emistry in TGLDM and LFD. For comparison, the TGLDM and lam-inar flamelet library generated at the stoichiometric mixture fractionfor the 1200 K methane/ethane jet are plotted in Fig. 3.6. For a giventemperature, the rate of change of temperature ( ? ) in the unsteadyflamelet library is generally higher than that of the TGLDM, especiallyin the temperature range of 1200?1500 K within which autoignitionoccurs. The TGLDM also includes more data points with lower changerates that are missing from the flamelet library. This difference betweenthe predictions obtained with the two libraries appears to be a conse-quence of the modelling technique. For the unsteady flamelet library,the scalar dissipation rate is modelled according to the approximatebehavior of the mean scalar dissipation along the stoichiometric inter-face in the jets. That the modelled scalar dissipation rate cannot exactlyreproduce all possible realizations of the scalar dissipation rate in thereal flame could be a source of the discrepancy.Another issue is that the solutions in the flamelet library are ex-660.7 0.75 0.8 0.85      0.5            1.0            1.51000/T (1/K)tau(ms)(a) Jet A0.7  0.75 0.8 0.85      0.5            1.0            1.51000/T (1/K)tau(ms)(b) Jet BFigure 3.5: Comparison of CSE-TGLDM and CSE-LFD predictions toexperimental results. +: experimental data; square: calculation with CSE-TGLDM; triangle: calculation with CSE-LFD.671000 1500 2000 2500 30001001021041061081010T(k)?T(K/s)Figure 3.6: Comparison of TGLDM and laminar flamelet library. +:laminar flamelet library; ?: TGLDM library.clusively igniting flamelets: the flamelet method is not able to accountfor the transition from ignition to burning. The TGLDM library cov-ers a much wider range of conditions and is able to transition fromigniting to burning conditions seamlessly. This, combined with the (ar-guably) better predictions of ignition delay, leads to the conclusion thatthe CSE-TGLDM method is the better of the two approaches available.3.3.3 Effects of Different Fuel CompositionThe CSE-TGLDM method is adopted to simulate jets of different fuelcomposition. To study the effects of dilution on the ignition delay, met-hane was diluted with nitrogen and autoignition was studied in shock-tube experiments. The characteristics of the diluted jet are listed as JetC in Table 3.1. The experimental results are discussed by McTaggart-680.7 0.75 0.8 0.85   0.5            1.0            1.5         1000/T (1/K)tau(ms)Figure 3.7: Ignition Delay tau vs T0. +: experimental data; square: predictionsof CSE-TGLDM.Cowan et al. [104]. This series of experiments is simulated and the pre-dictions are compared to the experiments in Fig. 3.7. The predictionsof ignition delay time capture the trend observed in the experimentaldata.To analyze the effects of different additives on the autoignition pro-cess, the measured and predicted ignition delay of three series of jets(from jet A to jet C) are compared. The effects of additives are not ob-vious in the experimental data because of the scatter in the measure-ments. There is as yet insufficient data from the methane/ethane ex-periments to be able to draw statistically significant conclusions aboutthe effects of ethane addition on the ignition these jets. However, asis shown in Fig. 3.8, the CSE-TGLDM model predicts a modest reduc-tion in ignition delay time of the methane/ethane jet ? in fact, the re-690.7 0.75 0.8 0.850.5            1.0            1.51000/T (1/K)tau(ms)Figure 3.8: Ignition Delay tau vs T0, predictions of CSE-TGLDM. openbullet: CH4;square: 90%CH4 +10%C2H6; triangle: 80%CH4 +20%N2.duction is considerably smaller than the scatter in either the methaneor methane/ethane data. At high pre-combustion temperatures, eventhis modest reduction becomes negligible. Further work is necessary topredict pollutant formation so as to study the trade-off between the re-duced ignition delay and increased emissions that result from addingethane to the fuel.Similarly, it is hard to draw conclusions on the effect of dilutionfrom the experiment because of the scatter in the measurements andthe amount of available experimental data points. The comparison ofignition delay predictions in Fig. 3.8 shows that the addition of nitro-gen delays methane autoignition. The effect is most significant at thelow pre-combustion temperature of 1200 K and becomes negligible athigher temperatures. Further work is necessary to predict the change70in NOx emission caused by the dilution.Mixing methane with additives can change the chemical kineticsof the flame. Not only does it change the stoichiometry of the flame,but it also affects the initial combustion rate by changing the energydensity of the fuel. These effects are reflected in the TGLDM libraries.For example, in Fig. 3.9, a comparison of conditional mean temper-atures at the stoichiometric mixture fraction is presented. Given thesame conditional means of progress variables, temperatures interpo-lated from the TGLDM vary with the fuel composition. The tempera-tures for methane/ethane jets are the highest in most cases, followedby pure methane jets and then methane/nitrogen jets. The most ob-vious difference is observed at small values of YCO2 and YH2O, corre-sponding to conditions at the start of ignition. For large values of theprogress variables, the difference is small; angbracketleftT|zetaangbracketright for pure methane andmethane/nitrogen jets are almost the same, which corresponds to theconditions approaching equilibrium. As a result, it is expected that pre-dictions of ignition delay change with fuel composition. In this study,the CSE-TGLDM model predicts ignition delay being the shortest formethane/ethane jets and the longest for jets diluted with nitrogen.Additives also affect the mixing of fuel and oxidizer. Given thesame injection parameters, more fuel is injected when C2H6 or N2 isadded because of the larger molecular weight. The gaseous jet momen-tum is increased and mixing is enhanced. As a demonstration, the con-tours of stoichiometry Zst in the pure methane and methane with addi-tives jets are studied at injection pressure 120 bar and pre-combustiontemperature 1200 K. Contours of Zst shown in Fig. 3.10 are recorded710 0.01 0.02130013501400YCO2T(K)0 0.02 0.04155016501750YCO2T(K)0 0.02 0.04 0.06 0.081800190020002100YCO2T(K)0 0.05 0.122002400YCO2T(K)Figure 3.9: Conditional averaged temperature evaluated usingTGLDMs of different fuel composition, at stoichiometric mixture frac-tion and pre-combustion temperature 1200 K. -?-: 90%CH4+10%C2H6,?-: CH4, --: 80%CH4 + 20%N2. Top left: YH2O = 0.02, top right:YH2O = 0.045, bottom left: YH2O = 0.07 , bottom right: YH2O = 0.095.720 10 20 30 40 50 60 70 800510x/deqr/deq(a) t=0.9ms0 10 20 30 40 50 60 70 800510x/deqr/deq(b) t=1.1msFigure 3.10: Comparison of Zst contours in different jets. ?-: pure me-thane; --:80%CH4 +20%N2; -? -: 90%CH4 t = 0.9ms (during injection), t = 1.1ms (after injection stops). Thefigures show that the methane/additive jets penetrate further than thepure methane jet. The Zst of methane/nitrogen jets is larger than thoseof methane and methane/ethane jets. When the contours at the samemixture fraction are compared for the three jets, the methane/nitrogenjet always penetrates furthest and expands widest, followed by meth-ane/ethane and then methane jets.The temperature fields of the different jets are compared during in-jection and after injection stops. Contours of mixture fraction are im-posed on temperature fields to compare the mixing of fuel and air. Asimilar observation as discussed previously is found in Fig. 3.11- 3.12,where enhanced mixing is exhibited in jets of methane mixed with ad-ditives. This comparison is of interest because of the possible applica-tion in compression ignition engines, where pilot fuel injection can be7310 20 30 40 50 600510x/deqr/deq40060080010001200(a) CH410 20 30 40 50 600510x/deqr/deq40060080010001200(b) 80%CH4 + 20%N210 20 30 40 50 600510x/deqr/deq40060080010001200(c) 90%CH4 + 10%C2H6Figure 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.7410 20 30 40 50 60 70 800510x/deqr/deq60080010001200(a) CH410 20 30 40 50 60 70 800510x/deqr/deq60080010001200(b) 80%CH4 + 20%N210 20 30 40 50 60 70 800510x/deqr/deq60080010001200(c) 90%CH4 + 10%C2H6Figure 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.75used to ignite the gas near the tip of the main jet and the shape of thejet thus affects the engine performance.Altogether, the effects of additives are a result of the competitionbetween chemical kinetics and mixing. The changes in ignition delayare related to the percentage and properties of the additives, as wellas the pre-combustion conditions. At low pre-combustion tempera-tures, ignition delay seems to be mainly affected by changes in chem-ical kinetics, where ethane addition accelerates ignition and nitrogenretards ignition. At high temperatures, the formation of the ignitionkernel seems to be mainly limited by the mixing of fuel and air, wherethe chemical kinetics of methane and methane with additives are fastenough and the changes in mixing caused by additives are not signifi-cant enough to make any noticeable difference because of the relativelyshort time for mixing prior to ignition.3.3.4 Effects of Different Injection ParametersThe effects of different injection parameters are also studied with theCSE-TGLDM method. Simulations are performed for pure methanejets at various pressure ratios. In this series of simulations, injectionpressure Pi is varied and predictions of ignition delay are compared tothe measurements of Wu [103]. The characteristics of the jets are listedin Table 3.1 as Jet D.The predicted and measured ignition delays are plotted against thepressure ratio in Figure 3.13. The results show that the autoignition de-lay is reduced with increased pressure ratio. For this series of jets, thepre-combustion temperature of 1300 K is relatively high and turbulentmixing seems to be the rate-limiting process of autoignition. Increas-761.5 2 2.5 3 3.5 4 4.5 5 5.5   0.5            1.0   tau(ms)Figure 3.13: Ignition delay tau and pressure ratio. +: experimental data;square: predictions of injection pressure helps to enhance mixing and therefore reducesignition delay time. This decrease becomes negligible at relatively highpressure ratios, when the effect of enhanced mixing is no longer thedominant factor and the autoignition process seems to be limited bychemical kinetics.The effect of varying injection duration is studied at the pre-com-bustion temperature of 1300 K using the CSE-TGLDM method. Thecharacteristics of the jet are listed in Table 3.1 as Jet E.The predicted ignition delays are compared with the experimen-tal data. Despite the scatter in the measurements, the numerical pre-dictions show that at the specified high pre-combustion temperature,ignition delay is insensitive to the injection duration. In fact, the pro-files of maximum angbracketleftT|zetaangbracketright for all injection durations coincide with one771   1.5 2.0 2.50.5            1.0   1.2ti (ms)tau(ms)Figure 3.14: Ignition delay tau and injection duration ti. Cross: experi-mental data; square: predictions of CSE-TGLDM.other during the autoignition process and only slight differences areobserved later in the burning stage.The irrelevance of injection duration to ignition delay, however, isinconclusive. For the specified temperature 1300 K, the ignition delayis relatively short and autoignition starts before the end of injection inall cases under study. At lower temperatures, ignition might start afterinjection stops and the ignition delay might be affected by injectionduration because of the changes in mixing behavior.3.4 ConclusionsThis chapter has presented simulations of autoigniting jets of meth-ane and methane mixed with additives. The model predictions of the78ignition delay are consistent with experimental results. However, thescatter observed in the experiments is not accounted for in the presentmodels. The discrepancy between the predictions of the two methodsappears to be due to the different modeling techniques used for gen-erating the libraries of full chemistry. Ultimately, it is found that theCSE-TGLDM approach is superior to the CSE-LFD approach. Addingethane to the methane jet results in a modest reduction in the ignitiondelay time at lower initial temperatures. However, this reduction is sig-nificantly smaller than the scatter in the experimental results. Dilutingmethane with nitrogen delays autoignition at low pre-combustion tem-peratures and the effects become less obvious at high temperatures. In-creasing pressure ratio reduces ignition delay and the effects becomenegligible at high pressure ratios. The injection duration seems irrel-evant when autoignition starts before the end of injection; its signifi-cance under other circumstances needs further study.79Chapter 4Modelling Radiative Effects of aTurbulent Non-premixed Flame usingCSE-TGLDM4.1 IntroductionUnlike RANS, Large Eddy Simulation of turbulent combustion solveslarge scale motions and might capture the transient behavior of tur-bulent flames. Conditional Source-term Estimation (CSE) has been ap-plied to study non-premixed turbulent flames in LES. As discussed inChapter 2, CSE-LFD and later CSE-TGLDM [62] has been used in LESto model the Sandia Flame D ? a piloted, non-premixed turbulent jetflame that has been widely studied numerically and experimentally. Inspite of the general good agreement between the predictions of CSEand the experimental data, there are some unanswered questions:? The radiation properties of Sandia Flame D remain to be mod-eled. The models to be used should be a compromise among pre-cision, complexity and computational cost.? The possibility of using CSE to describe Turbulence/RadiationInteraction (TRI) is to be tested. Similar to the Turbulence/Che-mistry Interaction, TRI is complicated by the nonlinearity of the80radiation process. Theoretically, TRI might be modeled by CSE ina similar way to that with which it closes chemical source terms.? The effects of neglecting radiation heat loss on the predictionsof temperature and species formation in Flame D need to be ana-lyzed. Previous studies with CSE are all based on the assumptionof adiabatic combustion. The CSE methods are to be extended tomodel radiation.In this chapter, the CSE-TGLDM method is extended and applied inLES of Sandia Flame D to answer these questions.4.2 LES of Sandia Flame D using CSE-TGLDMThe Sandia D-flame is a piloted, non-premixed CH4/Air flame. Themain jet is composed of 25% CH4 and 75% air. The nozzle diameter is7.2 mm and the Reynolds number is 22400. The pilot flame burns a pre-mixture of C2H2, H2, CO2, N2 and air having nominally the same equi-librium composition and enthalpy as the main jet. The jet is enclosedby an air co-flow with the velocity of 0.9 m/s. A detailed descriptionfor the configuration of this flame and experiment documentation areprovided 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 filteredbalance equations of continuity, momentum, energy and species in LEShave been discussed in Chapter 2. In this work, the unresolved subgrid-scale convective fluxes are modeled using the eddy viscosity and eddy81diffusivity models [56,59]:tildewide = -?( tildewider - tildewidetildewide)) = ?tildewidet((nablatildewide) +(nablatildewide)T - 23(nabla? tildewide)I), (4.1)-?(tildewidestj - tildewide tildewidej) = ?tildewidet,j triangleinv tildewidej, (4.2)-?(tildewider - tildewidetildewide) = tildewidet triangleinv tildewide, (4.3)where tildewidet is the eddy viscosity, tildewidet,j is the eddy diffusivity of speciesj, tildewidet is the thermal eddy conductivity and I is an identity matrix. Theturbulent transport coefficients tildewidet, tildewidet,j and tildewidet are computed using thedynamic model for compressible flows by Moin et al. [17,18].LES is performed in spherical coordinates on a computational gridof 192 ? 84 ? 48 control volumes in the streamwise, cross-stream andazimuthal directions respectively [56,59]. A finite volume method withsecond order accuracy in space is applied. A predictor-corrector-projec-tion scheme is used to integrate in time, with a second-order Adams-Bashford scheme in the predictor step and a semi-implicit Crank-Nichol-son scheme in the corrector step. A timestep splitting technique is usedin the computation of conditional averaged source terms to overcomethe stiffness of the chemical kinetic mechanism. The inflow boundary isdefined according to the experimental data, with velocity fluctuationssuperimposed on the measured mean velocity profile. A convectiveboundary condition [106] is used for the outflow boundary to stabilizethe simulation. A traction-free condition [107,108] is used at the lateralboundary to allow for a flux of ambient fluid across the boundary ofthe spreading jet. The code is developed using the Message-Passing In-terface (MPI) [109] for parallel computation on clusters. The chemicalsource-terms are closed with the extended CSE-TGLDM with radiationmodel.824.3 CSE-TGLDM with Radiation4.3.1 Radiation ModelsDifferent 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 radiationproperties might be evaluated using the filtered values of species andtemperature. However, similar to Turbulence/Chemistry interactions,the fluctuations around the means of relevant scalars affect the radia-tion quantities. This effect, or Turbulence/Radiation interaction (TRI),leads to closure problems similar to chemical closure.Coelho [113,114] studied Flame D using various models to evaluateradiation characteristics of gaseous species, combined with differentmethods to account for TRI. The pros and cons of a few radiation mod-els, gas radiation property models and TRI modelling methods werediscussed. In terms of flame radiation, it was suggested that the opti-cally thin approximation (where flame absorption is neglected) mightbe adopted in the simulation of non-luminous flames; more accuratecalculations might be done by solving radiative transfer equations withthe discrete ordinates method (DOM), where the absorption of radia-tion was described as a function of local conditions. As for gas radi-ation properties, the spectral line-based weighted-sum-of-gray-gases(SLW) model was suggested as an alternative to the weighted-sum-of-gray-gases (WSGG) model. Comparison was made among calcula-tions based on a variety of combinations: assuming an adiabatic flameor optically thin flame, or solving radiative transfer equations (RTE)combined with a SLW radiation property model or a model where ab-83sorption coefficients are approximated by curve fitting; combined withmethods that fully or partially account for TRI effects. It was foundthat the optically thin approach tends to over-predict the radiation heatloss; a complete RTE method only slightly reduces the overpredictionwhen the medium is treated as gray and Planck mean absorption coef-ficients are adopted. Accurate predictions are obtained with the com-bination of complete RTE and SLW models. Overall, the temperaturedifference between the adiabatic calculation and those with various ra-diation models was found to be less than 150 K in all calculations.4.3.2 FormulationTo include radiation effects in the combustion simulation using theCSE-TGLDM model, several problems need to be solved. To begin with,the radiative properties of gaseous species need to be modelled withprecision at an affordable cost. Secondly, a radiation model must bechosen. To use complicated models such as those solving the RadiativeTransfer Equation (RTE) may significantly increase the computationalcost of LES. A relatively simple radiation model with acceptable accu-racy is needed to evaluate the radiation heat loss term in the energyequation. Thirdly, Turbulence/Radiation Interaction (TRI) needs to beaccounted for in the calculation. Similar to the Turbulence/ChemistryInteraction, TRI is caused by turbulent fluctuations. Lastly, the TGLDMmethod for reducing detailed chemical kinetics has so far been basedon the assumption of adiabatic combustion. A TGLDM method includ-ing the effects of radiation needs to be developed.Considering the performance of different models and additionalcomputational cost, an optically thin flame model is used in this work,84combined with a gas radiation property model derived from the RAD-CAL model by Grosshandler [115]. It is assumed that the flame is opti-cally thin and each radiation point source has an unimpeded isotropicview of the cold surroundings. The radiation heat loss is calculatedwith [116]:Q = Q(T,Yi) = 4sigmaparenleftBigg nrsummationdisplay1pi ? ap,iparenrightBigg? parenleftbigT4 -T4infinityparenrightbig , (4.4)where Q[W/m3] is the radiation heat loss rated per unit volume, sigma =5.669?10-8[W/m2K4] is the Steffan-Boltzmann constant, nr is the num-ber of the species included in the calculation, pi is the partial pressureof species i in atmospheres (equal to mole fraction times local pres-sure), ap,i[m-1atm-1] is the Plank mean absorption coefficient of speciesi and Tinfinity is the background temperature.In the context of the optically thin model, the inclusion of speciessuch as CH4 and CO is not essential for methane flames and may beneglected. Moreover, previous study [117] shows that CO2 and H2Oarethe major radiating species in Sandia Flame D. Only these two speciesare considered in the computation of radiation heat loss. The Planckmean absorption coefficients of the gaseous species are evaluated as afunction of temperature by curve fitting within the temperature rangeof 300 K to 2500 K based on the RADCAL model:ap = c0+c1asteriskmath(1000T )+c2asteriskmath(1000T )2+c3asteriskmath(1000T )3+c4asteriskmath(1000T )4+c5asteriskmath(1000T )5,(4.5)with the constant coefficients listed in Table 4.1 [116].The CSE method is proposed to account for TRI where the nonlin-earity is partially accounted for with conditional averaging. The radia-85Table 4.1: Coefficients for Planck mean absorption coefficients.species c0 c1 c2 c3 c4 c5H2O -0.2309 -1.1239 9.4153 -2.9988 0.51382 -1.8684 ?10-5CO2 18.741 -121.31 273.50 -194.05 56.310 -5.8169tion heat loss in the turbulent flame is calculated with:angbracketleftQ|zetaangbracketright = Q(angbracketleftT|zetaangbracketright,angbracketleftYi|zetaangbracketright), (4.6)tildewideQ(x,t) approxequalintegraldisplay 10angbracketleftQ|zetaangbracketright tildewide(zeta; x,t) dzeta, (4.7)where the radiation properties of gaseous species are evaluated usingconditional means. The radiative heat loss is presented in the energyequation as a sink term.The effects of radiation heat loss also need to be considered in theTGLDM calculation. If the radiation term is directly included in the en-ergy equation in the TGLDM computation, the trajectories will deviatefrom the desired equilibrium states and tail off to extinction. Moreover,radiation quantities at a particular location in a flame is related to localconditions as well as global conditions at distant locations. Therefore,it is physically unreasonable to develop trajectories with an arbitraryheat loss term or one evaluated using only the instantaneous solutionsof species and temperature.A simple solution is to assume that radiation of the flame is notstrong enough to significantly affect the trajectories. This simplifica-tion might be acceptable since the overall fraction of radiative heatloss in Flame D is only around 5%. The formulation of such a me-thod would be almost the same as the adiabatic CSE-TGLDM, where86an adiabatic TGLDM library will be used and the only difference isthe sink term added to the energy equation in the combustion simu-lation. Alternatively, an additional parameter, the enthalpy defect canbe introduced in the TGLDM calculation to describe the relationshipbetween the manifolds and the local radiation conditions in the flame.The concept of enthalpy defect has been used in laminar flameletmodels to account for non-adiabatic effects [118]. It is the differencebetween the actual and the adiabatic enthalpy. The specific enthalpy ofa gaseous mixture can be calculated with:h =nssummationdisplay1Yihi =nssummationdisplay1parenleftbiggh0i +integraldisplay TT0cp,idTparenrightbigg, (4.8)where ns is the number of species in the mixture and h0i is the enthalpyof formation of species i at the standard reference temperature T0. Theadiabatic enthalpy is:had = ho +Z(hf -ho), (4.9)where ho and hf are the specific enthalpies of the oxidant and fuel inthe fresh mixture respectively. The enthalpy defect triangleh is defined as:triangleh = h -had. (4.10)The calculation of TGLDM with enthalpy defect as an additionalparameter involves the computation of multiple adiabatic TGLDMs formixtures at various pre-combustion temperatures. Firstly, an adiabaticTGLDM can be solved at the given pre-combustion temperature T0 offuel and oxidant for a sampled mixture fraction, where h0 = had(T0)and triangleh0 = 0. Another adiabatic TGLDM can be generated for the samemixture fraction at a lower initial temperature T1, where h1 = had(T1).87For this TGLDM, the enthalpy defect is:triangleh1 = h1 -had(T0) = had(T1) -had(T0). (4.11)Several TGLDMs with different enthalpy defect values can be gener-ated by changing the initial temperature of fuel and oxidant. Multi-ple TGLDMs can be generated at temperatures T0 > T1 > T2 > ... >Tj, where j is the number of TGLDM libraries and Tj should be lowenough to account for the maximum radiation heat loss observed in theflame under study. These TGLDMs are then tabulated and retrievedlater during combustion simulation.The modified CSE-TGLDM method can now evaluate reacting scal-ars using manifolds with the same enthalpy defects as those observedlocally in the flame. The local enthalpy defect of the flame can be calcu-lated directly using the instantaneous solutions of species and temper-ature, as shown in Eq. 4.8 and Eq. 4.10. Once the enthalpy defect trianglehis calculated, interpolation needs to be performed based on progressvariable (as in the adiabatic CSE-TGLDM) as well as triangleh. More specifi-cally, TGLDMs of two initial temperatures can be found such that:trianglehj1 <=< triangleh <=< trianglehj1+1, (4.12)where trianglehj1 is the enthalpy defect for the manifold of an initial temper-ature Tj1.Interpolation along the enthalpy defect dimension can be performedimplicitly, assuming that reactive scalars can be retrieved from the TG-LDM using the progress variables. At the same time, angbracketleftT|zetaangbracketright can be eval-uated using the CSE inversion method:tildewideT approxequalintegraldisplay 10angbracketleftBigT|zetaangbracketrightBigA?PZ,j (zeta) dzeta, j element A. (4.13)88Supposing the local enthalpy defect in the flame corresponds to a man-ifold of an initial temperature Tr, ideally, the temperature interpolatedusing this manifold angbracketleftT|zetaangbracketright(Tr) should be equal to the inverted one atconvergence: angbracketleftT|zetaangbracketright approxequal angbracketleftT|zetaangbracketright(Tr). By comparing the interpolated andinverted angbracketleftT|zetaangbracketright, one can decide which manifold has the same enthalpydefect as that of the local mixture and should be used to evaluate chem-ical source terms.In practice, TGLDMs are generated and stored for selected initialtemperatures T0,T1,...,Tj, which might not include the desired tem-perature Tr. In this case, interpolation in the dimension of enthalpydefect is needed. Specifically, two TGLDM libraries are selected suchthat:angbracketleftT|zetaangbracketright(Tj1) <=< angbracketleftT|zetaangbracketright <=< angbracketleftT|zetaangbracketright(Tj1+1), (4.14)Here,angbracketleftT|zetaangbracketright(Tj1) = g(angbracketleftYCO2|zetaangbracketright,angbracketleftYH2O|zetaangbracketright), (4.15)is the conditional averaged temperature recovered by simple interpo-lation over the manifold generated at the initial temperature of Tj1.Eq. 4.14 indicates Tj1 <=< Tr <=< Tj1+1. The conditional mean of a reactivescalar angbracketleftf|zetaangbracketright can then be calculated with:angbracketleftf|zetaangbracketright = angbracketleftf|zetaangbracketright(Tj1) asteriskmath (1 -psi) +angbracketleftf|zetaangbracketright(Tj1+1) asteriskmath psi, (4.16)where angbracketleftf|zetaangbracketright(Tj1) is the conditional mean interpolated using the TG-LDM manifold generated at the initial temperature Tj1. Interpolationin the direction of enthalpy defect is performed between manifolds ofTj1 and Tj1+1, where the interpolation factor is evaluated with:psi = angbracketleftT|zetastangbracketright -angbracketleftT|zetastangbracketright(Tj1)angbracketleftT|zetastangbracketright(Tj1+1) -angbracketleftT|zetastangbracketright(Tj1), (4.17)89where zetast is the stoichiometric mixture fraction; or by solving the fol-lowing problem:minbraceleftBigg nzsummationdisplay1parenleftBigangbracketleftT|zetaangbracketright(Tj1) asteriskmath (1 -psi) +angbracketleftT|zetaangbracketright(Tj1+1) asteriskmath psi -angbracketleftT|zetaangbracketrightparenrightBig2bracerightBigg, (4.18)where nz is the number of sample mixture fractions.4.3.3 Results and DiscussionLarge Eddy Simulation of Sandia Flame D is performed using the ex-tended CSE-TGLDM method. GRI-Mech 2.11 [119] is used instead ofGRI-Mech 3.0 [120] because it gave more accurate predictions of NO inprevious study. A transport equation is solved for NO, with chemicalsource terms retrieved from the TGLDM library. This leads to a moreaccurate prediction compared to extracting NO from the TGLDM li-brary because of the relatively slow process of NO formation. The sim-ulation results are compared to experimental data as well as Wang etal.?s simulation [59] where adiabatic combustion was assumed.Figs. 4.1 to 4.5 show the centerline profiles and radial profiles oftemperature and species. Generally speaking, the predictions of bothadiabatic CSE-TGLDM [59] and CSE-TGLDM with radiation are ingood agreement with the experimental data. The comparison of tem-perature and major products such as CO2 and H2Oshows minor differ-ences between these two methods, which indicates that the approxima-tion of adiabatic combustion is likely acceptable for the weakly radiat-ing Flame D. The predictions of CO2 and H2O are slightly improved inthe far field, but the improvement is slight.On the other hand, the predictions of NO are quite different. As ex-pected, the adiabatic CSE-TGLDM over-predicts NO formation, espe-90cially in the far field where radiation is relatively strong. CSE-TGLDMwith radiation gives better predictions of NO, which agree with the ex-perimental data well. Slight over-prediction is observed at x/D = 15and 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 ef-fects of the pilot flame are not precisely described. The under-predictionof NO at other stations is expected because of the optically thin modelused for radiation, which tends to over-predict heat loss in radiativeCH4 flames due to the strong absorption by CO2 around the 4.3 ?mwave length. This under-prediction might be improved if a more accu-rate radiation model including gaseous species absorption is used toevaluate the energy sink term in combustion simulation.Figs. 4.6 to 4.8 show the conditional means of progress variablesCO2 and H2O as well as NO. As observed in the unconditional means,the predictions of angbracketleftYCO2|zetaangbracketright and angbracketleftYH2O|zetaangbracketright are slightly improved in themid- and far-field using CSE-TGLDM with radiation; however, thedifference is small. Fig. 4.8 shows the conditional means of NO. Theunder-prediction of NO caused by the optically thin model is more ob-vious in the conditional means. Overall, the adiabatic computation andthat including radiation both capture the trends observed in the exper-iment.Compared to the adiabatic CSE-TGLDM method, the requirementof storage for CSE-TGLDM with radiation increases dramatically be-cause of the use of multiple TGLDMs. In terms of speed, the com-putation is slowed down because interpolation over each TGLDM isneeded in order to evaluate the enthalpy defect level and to close the910 20 40 6002468x/DT/T00 20 40 6000. 20 40 6000. 20 40 6000.511.5 x 10-4x/DYNOFigure 4.1: Centerline profiles of Flame D. delta: experimental data; ?-: predictions of adiabatic CSE-TGLDM; - ? -: predictions of CSE-TGLDM with radiation.920 2 4 60246r/DT/T0x/D=7.50 2 4 6 80246r/DT/T0x/D=150 2 4 6 802468x/D=30r/DT/T00 5 1002468x/D=45r/DT/T0Figure 4.2: Temperature radial profiles at different stations down-stream. delta: experimental data; ?-: predictions of adiabatic CSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.930 2 4 600. 2 4 6 800. 2 4 6 800. 5 1000. 4.3: Mass fraction of CO2 radial profiles at different stationsdownstream. delta: experimental data; ?-: predictions of adiabatic CSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.940 2 4 600. 2 4 6 800. 2 4 6 800. 5 1000. 4.4: Mass fraction of H2O radial profiles at different stationsdownstream. delta: experimental data; ?-: predictions of adiabatic CSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.950 2 4 60246 x 10-5r/DYNOx/D=7.50 2 4 6 802468 x 10-5r/DYNOx/D=150 2 4 6 800. x 10-4 x/D=30r/DYNO0 5 1000.511.5 x 10-4 x/D=45r/DYNOFigure 4.5: Mass fraction of NO radial profiles at different stationsdownstream. delta: experimental data; ?-: predictions of adiabatic CSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.960 0.5 100.050.1z<YCO2|xi>x/D=7.50 0.5 100.050.1z<YCO2|xi>x/D=150 0.5 100.050.1z<YCO2|xi>x/D=300 0.5 100.050.1z<YCO2|xi>x/D=45Figure 4.6: Conditional averages of CO2 radial profiles at different sta-tions downstream. delta: experimental data; ?-: predictions of adiabaticCSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.970 0.5 100.050.1z<YH2o|xi>x/D=7.50 0.5 100.050.1z<YH2o|xi>x/D=150 0.5 100.050.1z<YH2o|xi>x/D=300 0.5<YH2o|xi>x/D=45Figure 4.7: Conditional averages of H2O radial profiles at different sta-tions downstream. delta: experimental data; ?-: predictions of adiabaticCSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.980 0.5 10246 x 10-5z<YNO|xi>x/D=7.50 0.5 102468 x 10-5z<YNO|xi>x/D=150 0.5 102468 x 10-5z<YNO|xi>x/D=300 0.5 100.51x 10-4z<YNO|xi>x/D=45Figure 4.8: Conditional averages of NO radial profiles at different sta-tions downstream. delta: experimental data; ?-: predictions of adiabaticCSE-TGLDM; -? -: predictions of CSE-TGLDM with radiation.99reaction rates accordingly. The additional dimension of enthalpy de-fect also leads to more unsmoothness in the interpolated reaction rates,which might cause the ODE solver to choose smaller sub-timesteps toconverge solutions in the time-splitting. In this study of Flame D, theradiation is not strong and only three TGLDMs generated at differ-ent enthalpy defect levels are used; the increase of computational costis still manageable. However, a more dense grid in enthalpy defect isnecessary for highly radiating flames and CSE-TGLDM with radiationneeds to be optimized. One suggestion is to combine the parallel com-putation technique with the extended CSE-TGLDM method, where asegment of the flame simulated on each individual processor has a rel-atively small range of enthalpy defect and accordingly only a smallnumber of TGLDMs (that might differ from processor to processor)need to be loaded and used on each processor.Altogether, the adiabatic combustion and optically thin flame as-sumptions are both acceptable in the study of Sandia Flame D. Theadiabatic CSE-TGLDM method can be used to study non-premixedturbulent flames without strong radiation, such as Sandia Flame D;the method predicts temperature and major species (for which trans-port equations are solved) reasonably well. The radiation effects can beincorporated in CSE-TGLDM by introducing enthalpy defect as an ad-ditional variable to parameterize the manifold and interpolating basedon the progress variables as well as the enthalpy defect. The extendedCSE-TGLDM with radiation method improves the prediction of NOformation. Optimization of CSE-TGLDM with radiation should be fur-ther studied to deal with the increase of storage requirement and com-100putational cost.4.4 Radiation Intensities Analysis using CSE4.4.1 FormulationThe scalar fields predicted with CSE-TGLDM are used to evaluate themean spectral radiation intensities in the Flame D. The radiation trans-fer equations are solved for various radiation paths with RADCAL, aprogram based on a narrow band radiation model [115]. To account forTRI in the analysis of radiation intensities, a stochastic time and spacetime series simulation (TASS) method might be adopted [121]. TASSanalyzes statistics of radiation properties using many realizations ofthe flame under study. In this work, the CSE model is combined withthe RADCAL model to simulate TRI.In the RADCAL model, the radiative transfer equation (RTE) for anabsorbing and emitting medium is solved by breaking the line-of-sightinto uniform elements. The solution of the spectral intensity iprimelambda is [115]:iprimelambda(l) = iprimelambda,we-kappalambda(l) +integraldisplay kappalambda(l)0ib,lambda(lasteriskmath)exp[-(kappalambda(l) -kappalambda(lasteriskmath))]dkappalambda(lasteriskmath) (4.19)where ib,lambda is the Planck blackbody function, kappalambda is the optical thicknessdefined as kappalambda = integraltext l0 alambda(lasteriskmath)dlasteriskmath, lambda is the wavelength and the subscript wrefers to the boundary condition of a wall.A rough estimate of the mean spectral intensities can be obtained byusing the means of relevant scalars in the RADCAL model. In such acalculation, due to the nonlinearity of Eq. 4.19, the effects of fluctuatingscalar fields on the radiation intensities are neglected. To correct forthis, the CSE model can be combined with the RADCAL model in away similar to the closure for the reaction rates. The contributions of101uniform elements to the overall radiation intensity can be evaluatedusing conditional means:iprimelambda(lasteriskmath,l) =integraldisplay 10angbracketleftiprimelambda(lasteriskmath,l)|zetaangbracketrightP(zeta)dzeta, (4.20)where angbracketleftiprimelambda(lasteriskmath,l)|zetaangbracketright is evaluated with conditional averages of reactivescalars. The mean radiation intensity iprimelambda(l) is the summation of the con-tributions from all elements.4.4.2 Results and DiscussionThe radiation characteristics of the piloted methane-air diffusion flame,Sandia Flame D are studied. Large Eddy Simulation of the Flame Dhas been done with the CSE-TGLDM model discussed previously. Thecode has been revised and re-run to extract the reactive scalars neededin the RADCAL model. An example of the conditional means at dif-ferent stations downstream of the nozzle has been shown in Figs. 4.6and 4.7. The mean spectral radiation intensities along diametric pathsat x/D = 30, x/D = 45 and x/D = 60 are calculated using the meanproperties method and the CSE method as described in the previoussection. The results are compared with the experimental data of Zhenget al. [122] in Fig. 4.9.Predictions of both methods appear to capture the trends in themeasurements. The spectra show that the main sources of radiation inFlame D are water vapor and CO2. The strongest radiation intensity isobserved at x/D = 45, which is close to the stoichiometric flame heightof x/D = 47. The predictions of the CSE method are generally higherthan the mean value calculation as a result of accounting for TRI.The radiation quantities predicted with the CSE method seem to be102closer to the experimental data only at x/D = 60. However, this doesnot necessarily mean that using the mean values is the better choice. Inthe work of Zheng et al. [122], the experimental data was used in themean value method, and this lead to as much as a 50% underpredic-tion. The seemingly good performance of the mean properties methodhere might be as much attributable to a slight discrepancy between thepredicted mean scalar fields and those of the real flame.Zheng et al. [122] also provided predictions of TRI using the com-plicated TASS method with the measured scalars as input data. TheCSE predictions here are 15% to 22% higher than those of the meanproperties method around peak values; this magnitude of discrepancyis similar to that observed with the TASS method. Ultimately, the CSEmethod seems promising for modelling TRI.For comparison, the spectral radiation intensities are also evalu-ated [117] with the CSE method using scalar fields extracted from theadiabatic calculation by Wang [59]. As shown in Fig. 4.10, the discrep-ancy is negligible, which might be explained by the small difference intemperature, CO2 and H2O predictions of the adiabatic CSE-TGLDMand CSE-TGLDM with radiation methods. It again proves that FlameD is not strongly radiating and an adiabatic assumption is acceptablefor such a flame.This section has presented the analysis of spectral radiation inten-sities in the turbulent, piloted, Sandia Flame D. The CSE method hasbeen used to model TRI, predicting higher radiation intensities thanthe mean value method. The model predictions are consistent with theexperimental data.1031 1.5 2 2.5 3 3.5 4 4.5 501000200030004000500060007000lambda,?mIlambdamean,W/m2_sr_?m(a) x/D = 301 1.5 2 2.5 3 3.5 4 4.5 5010002000300040005000600070008000lambda,?mIlambdamean,W/m2_sr_?m(b) x/D = 45Figure 4.9: Spectral radiation intensities in Flame D. openbullet: measured; ?-:CSE method; --: mean properties method.1041 1.5 2 2.5 3 3.5 4 4.5 50100020003000400050006000lambda,?mIlambdamean,W/m2_sr_?m(c) x/D = 60Figure 4.9: Spectral radiation intensities in Flame D (cont.). openbullet: mea-sured; ?-: CSE method; --: mean properties method.4.5 ConclusionThe CSE method has been used to analyze the radiation properties ofthe Sandia Flame D. The study shows that CSE can be used to modelthe Turbulence/Radiation interaction (TRI) effects. The CSE-TGLDMmethod has been extended with the enthalpy defect as an additionaldimension in the manifolds. It has been combined with an opticallythin radiation model and applied in the LES of Flame D. The predic-tions of NO are improved using the extended CSE-TGLDM method.The mean spectral radiation intensities in the flame have been ana-lyzed by solving radiative transfer equations and modelling TRI withCSE. The results agree well with experimental results.1051 1.5 2 2.5 3 3.5 4 4.5 501000200030004000500060007000lambda,?mIlambdamean,W/m2_sr_?m(a) x/D = 301 1.5 2 2.5 3 3.5 4 4.5 5010002000300040005000600070008000lambda,?mIlambdamean,W/m2_sr_?m(b) x/D = 45Figure 4.10: Spectral radiation intensities in Flame D, openbullet: measured, ?-:CSE-TGLDM with radiation, --: adiabatic CSE-TGLDM.1061 1.5 2 2.5 3 3.5 4 4.5 50100020003000400050006000lambda,?mIlambdamean,W/m2_sr_?m(c) x/D = 60Figure 4.10: Spectral radiation intensities in Flame D (cont.). openbullet: mea-sured; ?-: CSE-TGLDM with radiation; --: adiabatic CSE-TGLDM.107Chapter 5A Progress Variable and PresumedProbability Density Function forPremixed Combustion5.1 IntroductionAs discussed in Chapter 2, partially premixed and premixed combus-tion also have wide applications in industry. Numerical simulationpromises to assist in the design of premixed combustion devices [34,123]. One major category in the models for turbulent premixed flamesis statistical models that use a presumed probability density function(PDF) [124?128]. The PDF approaches are not limited to the thin flame-lets regime and can be used to account for finite rate chemistry. De-tailed chemistry can be included by using tabulated data. Because ofthese advantages, PDF models have been applied in RANS and LES ofpremixed flames [129].Recently, efforts have been made to introduce the well-establishedConditional Moment Closure (CMC) model [39, 40], which has beenwidely used to model non-premixed flames [43, 45, 46, 48, 50], to pre-mixed flame simulations. Klimenko and Bilger developed the theoryof using CMC for premixed flames [42]. The method was applied inthe study of premixed hydrogen-air flames by Swaminathan and Bil-108ger [130] and in modelling lean premixed gas turbines by Martin [131].Although effective, CMC is computationally expensive because it in-troduces an additional dimension while solving the transport equa-tions for conditional means.In previous chapters, Conditional Source-term Estimation (CSE) hasbeen shown to be a promising closure method for non-premixed flames[7,57,59,60]. In CSE, conditional means are obtained by inverting in-tegral equations and chemical source terms are closed by invoking thefirst conditional moment hypothesis from CMC. It only requires theCMC closure hypothesis and the assumption of statistical homogene-ity, which do not exclude its application in premixed flames.To apply CSE to premixed flames, a suitable conditioning variableis needed. Since the fluctuations in temperature and mass fractionsin premixed flames are mostly related to a combustion progress vari-able, this seems to be a logical choice for this purpose. Bray et al. [132]studied three different PDFs of a species-based progress variable c us-ing DNS with single step chemistry. The definition of c was such thatc = 0 in reactants and c = 1 in products. A beta-PDF, a twin delta func-tion PDF and a PDF based on unstrained laminar flame propertieswere adopted to predict the mean reaction rates in turbulent premixedflames with finite rate chemistry. Although the beta-PDF and the lam-inar flamelet based PDF agreed with the DNS when the normalizedvariance g = cprime2/[c(1 - c)] was around 1, the precision of modelledPDFs under other conditions was not tested. It was also suggestedthat the laminar flamelet based PDF be used only in the thin flameletregime, where the Kolmogorov scales of the turbulence are larger than109the characteristic scales of the laminar flame.In this study, a product-based progress variable is introduced. Dif-ferent presumed PDF models are studied and a modification is sug-gested to improve the precision of the laminar flame-based PDF modelof Bray et al. [132]. The DNS data from Grout?s work [133] is used toevaluate the precision of different PDF models in the RANS and LESparadigms.5.2 The PDF ModelsA combustion progress variable can be defined with the mass fractionof products:c(x,t) = YP(x,t)/YP,eq, (5.1)where YP,eq is the mass fraction of products at equilibrium. By defi-nition, c = 0 in pure reactants and c = 1 in fully burnt mixture atequilibrium. The presumed PDF takes the form:P(casteriskmath; x,t) approxequal p(casteriskmath; c,c2), (5.2)where, presumably, a transport equation would be solved for c andpossibly also for c2 in a real implementation in a CFD code. The re-action rates are zero at c = 0 and c = 1, which is similar to what isseen in non-premixed flames: in non-premixed flames, reaction ratesare typically also zero at a mixture fraction Z = 0 and Z = 1.An example of P(casteriskmath) taken from the DNS data of Grout [133] is il-lustrated in Fig. 5.1, showing the distribution of c at different locationsin a turbulent premixed flame. Near the leading edge, where combus-tion has barely started, and the trailing edge, where the reactions ap-1100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 105101520253035cp(c)x=10x=30x=60x=90x=120Figure 5.1: Probability Density Function in a turbulent premixed flame.proach equilibrium, a bimodal distribution is observed and might becaptured by a beta-PDF. In the interior of the flame, because of the inter-mediate reactions, the formation of products might be incomplete evenif the fuel is mostly consumed. A flame-based PDF might capture thedistribution under such conditions if the premixed turbulent flame isextremely thin.5.2.1 beta-PDFAlthough the beta-PDF approximates the real distribution only in bimodaldistributions, it is tested here because of its wide use in both premixedand non-premixed flames. A beta-PDF of c is defined as:P(casteriskmath; x,t) = C (casteriskmath)a-1 (1 -casteriskmath)b-1 , (5.3)where a, b and C are functions of local statistical moments of c:a = cparenleftbigg c(1 -c)(c2 -c2) -1parenrightbigg(5.4)111b = a(1 -c)c (5.5)C = 1beta(a,b), (5.6)where beta(a,b) = integraltext 10 (casteriskmath)a-1 (1 -casteriskmath)b-1dcasteriskmath is the usual beta function.5.2.2 A Laminar Flame-based PDFA laminar flame-based PDF model has been proved to work well inthe thin flamelet regime [132]. The model is based on the observationthat in the limit of a thin flame, the iso-surfaces of c are all parallel. Theinterior of the PDF is assumed to be inversely proportional to | triangleinv c|.The PDF is defined as:P(casteriskmath; x,t) = Adelta(casteriskmath) +Bf(casteriskmath) +Cdelta(1 -casteriskmath), (5.7)where f(casteriskmath) is calculated using the solution of an unstrained laminarflame:f(casteriskmath) = 1triangleinvcasteriskmath . (5.8)The coefficients A, B and C are obtained based on the first threemoment equations of c:1 =integraldisplay 10P(casteriskmath)dcasteriskmath = A +Bintegraldisplay 101triangleinvcasteriskmath dcasteriskmath +C (5.9)c =integraldisplay 10casteriskmathP(casteriskmath)dcasteriskmath = Bintegraldisplay 10casteriskmathtriangleinvcasteriskmath dcasteriskmath +C (5.10)c2 =integraldisplay 10casteriskmath2P(casteriskmath)dcasteriskmath = Bintegraldisplay 10casteriskmath2triangleinvcasteriskmath dcasteriskmath +C. (5.11)112From these equations the coefficients can be easily obtained:B = c -c2I1 -I2 (5.12)C = c2I1 -cI2I1 -I2 (5.13)A = 1 -B asteriskmath I0 -C, (5.14)where I0, I1 and I2 are defined as:I0 =integraldisplay 101triangleinvcasteriskmath dcasteriskmath (5.15)I1 =integraldisplay 10casteriskmathtriangleinvcasteriskmath dcasteriskmath (5.16)I2 =integraldisplay 10(casteriskmath)2triangleinvcasteriskmath dcasteriskmath. (5.17)This PDF model has an unsolved problem. The PDF is constructedusing the solution of an unstrained laminar flame. Ideally, the domainfor the planar laminar flame would extend from negative infinity topositive infinity. If one were to have an ?ideal? simulation of the lam-inar flame with a domain extending from negative infinity to positiveinfinity, the PDF one would obtain from this would just be two deltasat 0 and 1. In reality, the computational domain of the laminar flameis finite and this affects the precision of the PDF as c arrowright 0 and c arrowright 1.Moreover, by definition, the profile of this PDF is determined by a lam-inar flame which is independent of the local conditions in the turbulentflame. Consequently, the PDF always rises dramatically towards infin-ity as c arrowright0 and c arrowright1, even if c is not near 0 or 1 and c2 is very small.113-Inf xmin x1 x2 xmax +Inf00. 5.2: Progress variable c in an unstrained laminar flame.5.2.3 Modified Laminar Flame-based PDFTo address the issues in the laminar flame-based PDF, some modifi-cations to the original model are proposed. The idea is to find a partof the laminar flame that approximates the conditions in the turbulentflame characterized by c and c2. Fig. 5.2 shows the profile of c in anunstrained laminar flame. Recognizing that the PDF will be stored atdiscrete points in c-space, a domain of interest for the laminar flamecan be determined by the discretization scheme.A common approach in implementing presumed PDFs is to storethe PDF in a look-up table; it is usually easier to store the integral(Pdc)i =integraldisplay ci+deltac/2ci-deltac/2Pdc (5.18)at discrete points i in c-space. It is possible to refine the spacing locallyin c-space to better resolve the reaction zone. However, for this simpledescription (and in the later tests), it is assumed that there will be a114constant spacing deltac in c-space. The first and last intervals in c-spacemust be treated differently, since the space is bounded by 0 and 1, suchthat(Pdc)0 =integraldisplay deltac/20Pdc (5.19)and(Pdc)n =integraldisplay 11-deltac/2Pdc (5.20)where (Pdc)n is the last interval in c-space; this way Eq. 5.18 is usedover the interval 1 < i < n -1.This discretization can be used to provide bounds for using thelaminar flame solution in assembling a PDF for a turbulent flame. Alower bound can be defined at cmin = deltac/2 and an upper bound atcmax = 1 - deltac/2. The positions xmin and xmax correspond to locationsin the laminar flame where c = cmin and c = cmax. The modelled PDFwould then comprise delta functions at 0 and 1 (at (Pdc)0 and (Pdc)n inthe table) and the function f(casteriskmath) over all of the other intervals, as wasdone in the original formulation of Bray et al. [132]. Further, however,it is possible to eliminate one or both of the delta functions from the PDFand truncate the function f(casteriskmath) as required in the event that the local cand c2 in the turbulent flame indicate that this is needed.Given the local c and c2 in a turbulent flame, a truncated area [x1,x2](with corresponding values of the progress variable [c1,c2]) can be lo-cated in the laminar flame such that the means and variance of c sam-pled within [x1,x2] match those in the turbulent flame. The functionalform of the modified PDF changes depending on the locations of x1and x2 relative to xmin and xmax. Four cases are possible:115? Case 1: xmin < x1 < x2 < xmaxP(casteriskmath; x,t) =bracelefttpbraceleftmidbraceleftbt0 if casteriskmath < c1 or casteriskmath > c2B1f(casteriskmath) if c1 <=< casteriskmath <=< c2? Case 2: x1 < xmin < x2 < xmaxP(casteriskmath; x,t) =bracelefttpbraceleftmidbraceleftbtA2delta(casteriskmath) +B2f(casteriskmath) if casteriskmath <=< c20 if casteriskmath > c2? Case 3: xmin < x1 < xmax < x2P(casteriskmath; x,t) =bracelefttpbraceleftmidbraceleftbt0 if casteriskmath < c1B3f(casteriskmath) +C3delta(1 -casteriskmath) if casteriskmath greaterequal c1? Case 4: x1 < xmin < xmax < x2P(casteriskmath; x,t) = A4delta(casteriskmath) +B4f(casteriskmath) +C4delta(1 -casteriskmath). (5.21)Under all four circumstances, the first three moment equations of ccan be used to locate the optimal area [x1,x2] within the laminar flameand evaluate the unknown coefficients. The unknowns are solved suchthat 1 = integraltext 10 P(casteriskmath)dcasteriskmath is satisfied and the deviation from the targeted cand c2 is minimized. The coefficients in the equations can be evaluatedas:B1 = 1I0(5.22)B2 = cI1(5.23)A2 = 1 -B2I0 (5.24)116B3 = c -1I1 -I0(5.25)C3 = 1 -B3I0 (5.26)B4 = c -c2I1 -I2 (5.27)C4 = c2I1 -cI2I1 -I2 (5.28)A4 = 1 -B4 asteriskmath I0 -C4. (5.29)I0, I1 and I2 are redefined as:I0 =integraldisplay xuxl1triangleinvcasteriskmath dcasteriskmath (5.30)I1 =integraldisplay xuxlcasteriskmathtriangleinvcasteriskmath dcasteriskmath (5.31)I2 =integraldisplay xuxlcasteriskmath2triangleinvcasteriskmath dcasteriskmath, (5.32)where xl = max(x1,xmin) and xu = min(x2,xmax).5.3 Comparison with DNS5.3.1 DNS ConfigurationThe DNS [133] was performed in an inflow-outflow configuration us-ing the code SENGA [134], where the fully compressible Navier-Stokesequations were solved explicitly with the 3rd order Runge-Kutta timeadvances. The unburnt mixture consisted of 5% methane, 20% oxygenand 75%nitrogen at 900K. The flame was solved in a 1283 domain, with117a physical edge domain of 6 mm, and the inflow velocity was forcedto provide a decaying, isotropic, homogeneous turbulent field throughwhich the flame propagates.The chemical mechanism used in the DNS of laminar and turbulentflames is a two-step mechanism:F +O arrowrightI +P (R1)I +O arrowright2P (R2)where F and O are model fuel and oxidizer which are analogous toCH4 and O2; I is a model intermediate which is a combination of COand H corresponding to the full mechanism; P is a combination of H2Oand CO2. The rates of the two reactions, in [mol/cc/s] are calculated as:? R1 = 7.8(1014)e-7640T [F][R] (5.33)? R2 = 8.9(1016)T-1.8[O][R] (5.34)where [R] is the radical concentration given by:[R] = 4.5(105)e2300T[P] [O]0.5[I]1.5e-lambdaradical154|F||O| (5.35)The parameter lambda is related to the ratio of the rate of the elementaryfuel consumption reaction to the chain branching reaction. Reducing lambdaresults in a disproportionate increase in the flame speed relative to itseffect on the local gradient and giving a more easily resolvable flamefront.5.3.2 Comparisons of PDF ModelsThe proposed PDFs are tested against the DNS of a premixed turbulentflame in the both the RANS and LES paradigms. For the RANS test,118the distributions of c are studied over the whole domain and on planesorthogonal to the direction of the mean flow ? for this flow, the tur-bulence on these planes should be statistically homogeneous. The PDFprofiles from the DNS are compared with the three presumed PDFs.Fig. 5.3 shows the distribution of c over the whole domain at 40evenly spaced discrete points in c-space. The beta-PDF captures the bi-modal shape of distribution, but the magnitude and location of thepeaks are nearly symmetrical, which fails to capture the shape of thetrue PDF. The laminar flame-based PDF does a better job than the beta-PDF, especially for large values of c. The original version of the laminarflame-based PDF requires that there be a delta at casteriskmath = 1, which does notmatch the behavior of the true PDF; the modified version of the lami-nar PDF more closely matches the true PDF, although it over-predictsthe magnitude of the local maximum at casteriskmath approxequal 0.97. The difference be-tween PDFs from the original and the modified versions is more ob-vious in the Cumulative Distribution Function (CDF): the CDF calcu-lated with the original version departs from the DNS as casteriskmath increases,while the CDF calculated with the modified version almost perfectlyreproduces the CDF from the DNS.Figs. 5.4 and 5.5 show the distributions of c on planes at different lo-cations. Overall, the beta-PDF and the original laminar flame PDF predic-tions contain significant errors. On the other hand, the modified PDFagrees well with the DNS, although it generally fails to capture theexact location and magnitude of the peak on the products side of thedistribution.For the LES tests of the PDFs, the distributions of c in cubes consist-1190 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 102468101214161820casteriskmathp(casteriskmath)(a) Probability Density Function0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Cumulative Distribution FunctionFigure 5.3: Presumed Probability Density Function and CumulativeDistribution Function in a premixed turbulent flame. ?-: DNS; --:modified PDF model; -? -: original PDF by Bray et al.; ? ? ? ? ?: beta-PDF.1200 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10246810casteriskmathp(casteriskmath)(a) x/L = 0.2340 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1051015casteriskmathp(casteriskmath)(b) x/L = 0.3910 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.990510152025casteriskmathp(casteriskmath)(c) x/L = 0.6250 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1010203040casteriskmathp(casteriskmath)(d) x/L = 1Figure 5.4: Probability Density Function on planes in a premixed tur-bulent flame. ?-: DNS; --: modified PDF model; - ? - line: originalPDF by Bray et al.; ? ? ? ? ?: beta-PDF.1210 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x/L = 0.2340 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x/L = 0.3910 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x/L = 0.6250 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x/L = 1Figure 5.5: Cumulative Distribution Function on planes in a premixedturbulent flame. ?-: DNS; --: modified PDF model; - ? -: originalPDF by Bray et al.; ? ? ? ? ?: beta-PDF.122ing of 4 ? 32 ? 32 grid cells are analyzed. The results of the originallaminar flame PDF are not presented since the improvements in themodified version are similar to those observed in the RANS tests. Forthis test, PDF and CDF tables were generated for 40 grid points in thespace of the conditioning variable c. The tables are calculated using themodified laminar flame model for given values of c and c2.As a demonstration, the distributions in four cubes are shown inFigs. 5.6 and 5.7. The distribution calculated with the modified lami-nar flame PDF is obtained by look-up and linear interpolation in thepre-calculated tables. Observations from different cubes show that thebeta-PDF model predictions usually contain major errors. The interpo-lated values from tables of the modified PDF model agree well withthe DNS, with the biggest discrepancy again being around the peak inthe distribution towards the product side.Among the three presumed PDFs, the modified version of the lam-inar flame PDF model best captures the c distribution in the premixedturbulent flame. In the PDF plots, the maximum of the PDF from theDNS on the products side of the distribution tend to appear roundedoff, likely due to turbulent mixing. This is typically the only featurethat the modified PDF fails to capture. One way to improve this mightbe to add a filtering operator to the PDF:?f(casteriskmath) =integraldisplay +infinity-infinityf(c+)g(c+,casteriskmath)dc+. (5.36)Further work would be necessary to define an appropriate filter kernelg.The discrepancy between the modified PDF and the DNS seems lesssignificant in the CDF calculation. In CSE, the presumed PDFs are used1230 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1051015202530casteriskmathp(casteriskmath)(a) c = 0.089, c2 = 0.0480 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 105101520casteriskmathp(casteriskmath)(b) c = 0.406, c2 = 0.2980 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1051015casteriskmathp(casteriskmath)(c) c = 0.565, c2 = 0.4510 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1051015casteriskmathp(casteriskmath)(d) c = 0.800, c2 = 0.674Figure 5.6: Probability Density Function in different cubes in a pre-mixed turbulent flame. ?-: DNS; --: modified PDF model; - ? -: beta-PDF.1240 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 c = 0.089, c2 = 0.0480 0.2 0.4 0.6 0.8 c = 0.406, c2 = 0.2980 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 c = 0.565, c2 = 0.4510 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 c = 0.800, c2 = 0.674Figure 5.7: Cumulative Distribution Function in different cubes in apremixed turbulent flame. ?-: DNS; --: modified PDF model; - ? -:beta-PDF.125to discretize integral equations:f approxequalJsummationdisplayi=1angbracketleftf|casteriskmathjangbracketrightP(casteriskmathj; c,c2)dcasteriskmathj, (5.37)where J denotes the number of grid points in the conditioning vari-able space. In these equations, the precision of P(casteriskmathj; c,c2)dcasteriskmathj is moreof concern. It is more related to the CDF. Considering the good CDFpredictions of the modified PDF, the discrepancy might be acceptablein the context of the CSE model. This argument will be explored in theCSE tests in Chapter 6.Another series of tests was performed to further analyze the preci-sion of the PDF models. The same RANS and LES tests are repeatedfor a second set of DNS data. This set of DNS data is of a turbulentpremixed flame with lower turbulence intensity with a different two-step simplified chemical mechanism, where the coefficients in the Ar-rhenius expressions are tuned to accelerate the DNS computation. Forclarity, the DNS data obtained using the original chemical kinetics willbe referred to as solution ?I?, and data obtained using the modifiedchemical kinetics as solution ?II? hereafter.The results of the a priori tests with data set II are shown in Fig. 5.8to Fig. 5.12. As observed in the first series of tests, the modified laminarflame-based PDF model generally captures the distribution of c in theturbulent flame; overall, its predictions are better than those of the beta-PDF and the original flamelet PDF model. The discrepancies betweenthe modified PDF and the DNS from data set II are in similar locationsto those found using set I. This suggests that this modified PDF maywork well regardless of the chemical kinetics, although considerablymore work would be required to state this conclusively.1260 0.2 0.4 0.6 0.8 105101520casteriskmathp(casteriskmath)(a) Probability Density Function0 0.2 0.4 0.6 0.8 Cumulative Distribution FunctionFigure 5.8: PDF and CDF in the whole domain. ?-: DNS; -? -: originallaminar flame-based model; --: modified PDF model; ? ? ? ? ?: beta-PDF.1270 0.2 0.4 0.6 0.8 1051015casteriskmathp(casteriskmath)(a) x/L = 0.2340 0.2 0.4 0.6 0.8 105101520casteriskmathp(casteriskmath)(b) x/L = 0.3910 0.2 0.4 0.6 0.8 1051015202530casteriskmathp(casteriskmath)(c) x/L = 0.6250 0.2 0.4 0.6 0.8 1010203040506070casteriskmathp(casteriskmath)(d) x/L = 1Figure 5.9: Probability Density Function on different planes. ?-: DNS;- ? -: original laminar flame-based model; --: modified PDF model;? ? ? ? ?: beta-PDF.1280 0.2 0.4 0.6 0.8 x/L = 0.2340 0.2 0.4 0.6 0.8 x/L = 0.3910 0.2 0.4 0.6 0.8 x/L = 0.6250 0.2 0.4 0.6 0.8 x/L = 1Figure 5.10: Cumulative Distribution Function on different planes. ?-:DNS; - ? -: original laminar flame-based model; --: modified PDFmodel; ? ? ? ? ?: beta-PDF.1290 0.2 0.4 0.6 0.8 1020406080100casteriskmathp(casteriskmath)(a) ? = 0.089, ?2 = 0.0460 0.2 0.4 0.6 0.8 105101520casteriskmathp(casteriskmath)(b) ? = 0.413, ?2 = 0.2820 0.2 0.4 0.6 0.8 1024681012141618casteriskmathp(casteriskmath)(c) ? = 0.461, ?2 = 0.3110 0.2 0.4 0.6 0.8 10123456789casteriskmathp(casteriskmath)(d) ? = 0.689, ?2 = 0.505Figure 5.11: Probability Density Function in different cubes. ?-: DNS;--: modified PDF model; -? -: beta-PDF.1300 0.2 0.4 0.6 0.8 ? = 0.089, ?2 = 0.0460 0.2 0.4 0.6 0.8 ? = 0.413, ?2 = 0.2820 0.2 0.4 0.6 0.8 ? = 0.461, ?2 = 0.3110 0.2 0.4 0.6 0.8 ? = 0.689, ?2 = 0.505Figure 5.12: Cumulative Distribution Function in different cubes. ?-:DNS; --: modified PDF model; -? -: beta-PDF.1315.3.3 Effects of Chemical KineticsTo study the sensitivity of the proposed PDF model to the chemicalkinetics, the distribution of c in the second turbulent flame is evalu-ated with the modified PDF model; the predictions using the ?right?chemistry (based on laminar flame solution II) and the ?wrong? chemi-stry (based on laminar flame solution I) are compared to the DNS. Theresults are shown in Figs. 5.13 to 5.17.The results of both the RANS and LES tests show that the preci-sion of the modified laminar flamelet PDF is strongly affected by thechemical kinetics. When the laminar flame solution is based on a dif-ferent chemical mechanism, the performance of the PDF model seemsto be far less satisfying, especially in the calculation of the CDF. Thepredictions of the beta-PDF are not shown for better visibility; however,the modified PDF predictions are still closer to the DNS than the beta-PDF even when the ?wrong? chemistry (i.e., laminar flame solution I) isused.The discrepancy between the PDF model predictions and the DNSis increased when different chemical kinetics are used for the solu-tions of the one dimensional laminar flame and the premixed turbulentflame. According to the modified PDF, the distribution of c between[0,1] is proportional to f, the inverse of product gradient as defined inEq. 5.8. In Fig. 5.18, comparison is made between normalized productgradient of the laminar flame solution I and II corresponding to twodifferent chemical mechanisms. It is shown that with different chemi-stry, the profiles of normalized f differ significantly, with a shift in thelocation of peak and changes in the slopes of the curve. Such differ-1320 0.2 0.4 0.6 0.8 1051015casteriskmathp(casteriskmath)(a) Probability Density Function0 0.2 0.4 0.6 0.8 Cumulative Distribution FunctionFigure 5.13: PDF and CDF in the whole domain. ?-: DNS; -? -: mod-ified PDF model with laminar flame solution I; --: modified PDF mo-del.1330 0.2 0.4 0.6 0.8 1051015casteriskmathp(casteriskmath)(a) x/L = 0.2340 0.2 0.4 0.6 0.8 1051015casteriskmathp(casteriskmath)(b) x/L = 0.3910 0.2 0.4 0.6 0.8 101020304050casteriskmathp(casteriskmath)(c) x/L = 0.6250 0.2 0.4 0.6 0.8 1010203040506070casteriskmathp(casteriskmath)(d) x/L = 1Figure 5.14: Probability Density Function on different planes. ?-: DNS;-?-: modified PDF model with laminar flame solution I; --: modifiedPDF model.1340 0.2 0.4 0.6 0.8 x/L = 0.2340 0.2 0.4 0.6 0.8 x/L = 0.3910 0.2 0.4 0.6 0.8 x/L = 0.6250 0.2 0.4 0.6 0.8 x/L = 1Figure 5.15: Cumulative Distribution Function on different planes. ?-:DNS; - ? -: modified PDF model with laminar flame solution I; --:modified PDF model.1350 0.2 0.4 0.6 0.8 1010203040506070casteriskmathp(casteriskmath)(a) ? = 0.089, ?2 = 0.0460 0.2 0.4 0.6 0.8 1051015202530casteriskmathp(casteriskmath)(b) ? = 0.413, ?2 = 0.2820 0.2 0.4 0.6 0.8 1024681012141618casteriskmathp(casteriskmath)(c) ? = 0.461, ?2 = 0.3110 0.2 0.4 0.6 0.8 10123456789casteriskmathp(casteriskmath)(d) ? = 0.689, ?2 = 0.505Figure 5.16: Probability Density Function in different cubes. ?-: DNS;-?-: modified PDF model with laminar flame solution I; --: modifiedPDF model.1360 0.2 0.4 0.6 0.8 ? = 0.089, ?2 = 0.0460 0.2 0.4 0.6 0.8 ? = 0.413, ?2 = 0.2820 0.2 0.4 0.6 0.8 ? = 0.461, ?2 = 0.3110 0.2 0.4 0.6 0.8 ? = 0.689, ?2 = 0.505Figure 5.17: Cumulative Distribution Function in different cubes. ?-:DNS; - ? -: modified PDF model with laminar flame solution I; --:modified PDF model.1370 0.2 0.4 0.6 0.8 5.18: Comparison of laminar flames with different chemistryences in f lead to the error in the PDF model when the ?wrong? chem-istry is used. It is also observed that the error introduced by differentkinetics is less significant at relatively small values of the variance cprimeprime2.Under such conditions, the modified PDF might predict a non-zerodistribution over a small range of casteriskmath around the mean value ? to matchcprimeprime2 in the turbulent flame. Only a relatively small segment of the lam-inar flame solutions (corresponding to casteriskmath within this range) needs tobe used. The error introduced by the laminar solution using a differentchemical mechanism therefore might be less significant.5.4 ConclusionsDifferent forms for the presumed PDF of a product-based reaction prog-ress variable have been tested against DNS of turbulent premixed flames.On the whole, the beta-PDF rarely captures the distribution of the vari-138able. A modification is suggested to extend the laminar flame basedPDF proposed by Bray et al.. [132]. Comparison shows improvementsin PDF and CDF predictions with the modified laminar flame PDF. Theprecision of the modified PDF model is affected by the chemical kinet-ics. The same mechanism used in the turbulent combustion simulationshould be adopted to solve the prototype laminar flame that would beused to calculate the presumed PDF if this is to best match the subse-quently modelled turbulent flame.139Chapter 6Closure of Reaction Rates inTurbulent Premixed Combustion withCSE6.1 IntroductionIn previous studies, closure with CSE was obtained for non-premixedflames using the mixture fraction Z for conditioning. For premixedflames, chemical closure might be obtained by using the product basedprogress variable c for conditioning. In this chapter, the mean reactionrates in turbulent premixed flames are evaluated using the CSE modeland the presumed PDFs proposed in Chapter 5. The DNS data fromGrout?s work [133] is used to evaluate the precision of closure withCSE.6.2 The CSE FormulationThe formulation of CSE with c as the conditional variable is similarto the one with Z. The first conditional moment closure hypothesis ofCMC is adopted, which evaluates chemical source terms using quanti-ties conditional on the progress variable c:? |casteriskmath approxequal ? (T|casteriskmath,Yi|casteriskmath,rho|casteriskmath), (6.1)140where casteriskmath is the sample space for the progress variable.As discussed in previous chapters, instead of solving for the CMCtransport equations for the conditional averages, CSE obtains condi-tional averages of different quantities by inverting integral equations,taking advantage of the statistical homogeneity of conditional aver-ages on pre-determined surfaces in the flow field. Similar to non-pre-mixed combustion simulation, ensembles of points with homogeneitymay be selected in premixed flames, where the conditional means atdifferent points are equal to those of the ensemble:f|casteriskmath(x,t) approxequal angbracketleftf|casteriskmathangbracketright(t; A). (6.2)Here, A is the ensemble of homogeneity and x is a point within thatensemble.The CSE inversion equations can be constructed in a way similar tothe non-premixed cases if the probability density function of the newconditioning variable c is known. For a selected ensemble of statisticalhomogeneity, such integrals may be written for different points:f(x,t) approxequalintegraldisplay 10angbracketleftf|casteriskmathangbracketright(t; A)P(x,t; casteriskmath)dcasteriskmath, (6.3)where P(x,t; casteriskmath) is the local probability density function at point x.Assuming that the probability density function of c can be well-approximated by the PDF models described previously, Eq. 6.3 be-comes a relatively simple linear system that can be inverted so longas some kind of regularization is provided. Thus, the conditional av-erages of different scalars can be approximated. Chemical closure canthen be achieved by invoking the CMC hypothesis, Eq. 6.1. The chem-ical source term in the transport equations of the unconditional meansis closed by integrating over the sample space of casteriskmath, using Eq. 6.3.141In the inversion process, a linear regularization method can be usedto overcome ill-conditioning of the problem. In the simulation of realflows, it has been found that temporal continuity in the conditional av-erages can be assumed and conditional means obtained at the previoustime step used for regularization [57]. In an a priori test, such informa-tion is unavailable. Alternately, two methods are suggested for the CSEinversion computation:? Use the laminar flame solutions for regularization. The argumentis that the premixed turbulent flame is extremely thin and simi-larity is expected between the mean flow and the laminar flame.? Solve the inversion equation iteratively from some arbitrary ini-tial condition and use the solutions of a previous iteration for reg-ularization. The computation stops when the solutions becomeconverged or when some maximum number of iterations are reach-ed. The idea is to mimic the time-advancement in a CSE-basedsimulation with this iteration process.6.3 Results and Discussion : RANS tests6.3.1 Conditional Averages: Mechanism IThe mean reaction rates are closed using the presumed PDF and theCSE method. Tests are performed using the turbulent flame solutionsI, with the two-step mechanism defined previously by Eqs. 5.33 to 5.35.In the RANS tests, the mean reaction rates on parallel planes or-thogonal to the direction of the main flow are extracted from the DNSdata base. The reaction rates are also evaluated with scalars conditional142on c, using the presumed PDFs and the CSE formulation described pre-viously.An assumption of statistical homogeneity in the conditional aver-ages over the flow field is adopted to construct the linear equation sys-tem for inversion. As discussed previously, a suitable regularizationmethod is important for evaluating conditional means in CSE. The ef-fectiveness of the two proposed methods are compared, using the samepresumed PDFs in the CSE inversion computation.The results presented in Fig. 6.1 are obtained using the modifiedlaminar flame PDF. The method of regularizing with the laminar flamesolution yields good estimates of the conditional means for all reactivescalars. The iteration method which started from an initial conditionof zero for each scalar either fails to obtain converged solutions within2000 iterations, or yields solutions less desirable than those of the firstmethod. The most significant errors are observed in the iteration solu-tions of angbracketleftYI|casteriskmathangbracketright. Since the second reaction in the two-step mechanism ishighly sensitive to the change in intermediates, inaccuracy here causeseven bigger errors in the reaction rate calculation.The inversion methods are also tested with the beta-PDF to explorethe sensitivity of CSE method to the presumed shape of the PDF. Theresults are presented in Fig. 6.2. The solutions with the iteration me-thod are extremely poor, and are therefore not presented. Regulariza-tion with the laminar flame solution yields an acceptable approxima-tion for the conditional means, but the deviation from the DNS is largerthan that observed in the solutions with the modified PDF.The laminar flame solutions used for regularization are also plotted1430 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<YF|casteriskmath>(a) angbracketleftYF |casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<YO|casteriskmath>(b) angbracketleftYO|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.0050.010.0150.020.0250.03casteriskmath<YI|casteriskmath>(c) angbracketleftYI|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 150010001500200025003000casteriskmath<T|casteriskmath>(d) angbracketleftT|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<rho|casteriskmath>(e) angbracketleftrho|casteriskmathangbracketrightFigure 6.1: Conditional averages calculated in RANS paradigm, ob-tained by using modified laminar flame PDF and different regular-ization methods. ?-: DNS; --: regularization with laminar flame so-lution (CSE-1); - ? -: regularization with iterations up to 2000 times(CSE-2); ? ? ? ? ?: laminar flame solution.1440 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<YF|casteriskmath>(a) angbracketleftYF |casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<YO|casteriskmath>(b) angbracketleftYO|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.0050.010.0150.020.0250.03casteriskmath<YI|casteriskmath>(c) angbracketleftYI|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 150010001500200025003000casteriskmath<T|casteriskmath>(d) angbracketleftT|casteriskmathangbracketright0   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<rho|casteriskmath>(e) angbracketleftrho|casteriskmathangbracketrightFigure 6.2: Conditional averages calculated in RANS paradigm, ob-tained by using beta-PDF and regularization with laminar flame solution.?-: DNS; --: CSE.145in Fig. 6.1, which happen to coincides with the DNS conditional meansin this test. This seems to cast questions on the effectiveness of the CSEinversion method, which will be clarified in the following discussionand later in another series of a priori tests.The problems with the iteration methods might be attributed to thelack of a stabilization mechanism that would be present in CSE simula-tions. In previous RANS and LES simulations using CSE, the changesin conditional means affect the reactive scalars through the chemicalsource terms, and changes in reactive scalars are directly reflected inthe conditional means. The interactions between unconditional andconditional solutions appear to stabilize the CSE inversion process.Such a mechanism does not exist in the a priori tests, which may beleading to the deviation and fluctuations observed in the results fromthe iteration method.Regularization with laminar flame solutions works relatively wellwith both the beta-PDF and the modified PDF. This observation indicatesthat a method similar to LFD [61] might work well in premixed flamesimulations, in which a linear combination of laminar flame solutionsat various strain rates might provide a better approximation to theflame.6.3.2 Reaction Rate Predictions: Mechanism IThe chemical reaction rates are modelled and compared with the DNS.Figs. 6.3 and 6.4 present the conditionally averaged reaction rates eval-uated using the modified laminar flame PDF and the beta-PDF. The seem-ingly small errors in the solutions of conditionally averaged reactivescalars are amplified in the calculation of reaction rates. The CSE mo-146del with the modified PDF over-predicts the reaction rate angbracketleft ? |casteriskmathangbracketright by upto 20%, while the model with beta-PDF predictions are obviously unac-ceptable.The chemical reactions rates on different planes are closed by in-tegration. Only predictions with the modified PDF are presented. Asshown in Fig. 6.5, overall, the model predictions agree reasonably wellwith the DNS when the modified PDF is used. The peak reaction ratesthat occur in the middle of the flow field are captured well for bothreactions.6.3.3 Comparison with DNS: Mechanism IIThe a priori test described in the previous section is repeated for a sec-ond set of data, turbulent flame solution II. The conditional means ofreactive scalars and mean reaction rates are evaluated using the CSEand presumed PDF method; the results of RANS tests are shown inFigs. 6.6 to 6.8.Similar observations are found in this series of tests. The resultsgenerally agree with the DNS data, indicating the effectiveness of chem-ical closure with the CSE and the proposed PDF model in premixedflames. It is noteworthy that in this test, the conditional means fromthe DNS are different from the laminar flame solutions used for reg-ularization (method CSE-1); the conditional means obtained using theCSE inversion method are closer to the DNS than the a priori informa-tion from the laminar solutions. This clarifies the question raised inthe first series of tests: good estimation of conditional means can beobtained using the CSE method, even when the a priori information in-troduced for regularization contains error. Altogether, the predictions1470 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<?omegai|casteriskmath>(a) angbracketleft ?i|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<?omegaii|casteriskmath>casteriskmath(b) angbracketleft ?ii|casteriskmathangbracketrightFigure 6.3: Conditionally averaged reaction rates angbracketleft ? |casteriskmathangbracketright calculated inRANS paradigm. Conditional means obtained by using the modifiedPDF and regularization with the laminar flame solution. ?-: DNS; --:CSE.1480 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.511.522.53casteriskmath<?omegai|casteriskmath>(a) angbracketleft ?i|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.511.522.53casteriskmath<?omegaii|casteriskmath>(b) angbracketleft ?ii|casteriskmathangbracketrightFigure 6.4: Conditionally averaged reaction rates angbracketleft ? |casteriskmathangbracketright calculated inRANS paradigm. Conditional means obtained by using the beta-PDF andregularization with the laminar flame solution. ?-: DNS; --: CSE.1490 0.2 0.4 0.6 0.8 angbracketleft ?iangbracketright0 0.2 0.4 0.6 0.8 angbracketleft ?iiangbracketrightFigure 6.5: Reaction rates angbracketleft ? angbracketright on planes in a premixed turbulent flame.Conditional means obtained by using the modified PDF and regular-ization with the laminar flame solution. ?-: DNS; --: CSE.150of the CSE model capture the trend observed in the DNS data, whereasthe reaction rates tend to be over-predicted, as in the first test.6.4 Results and Discussion : LES tests6.4.1 Comparison with DNS: Mechanism IA priori tests are performed in the LES paradigm using the turbulentflame solution I.In the LES tests, the mean reaction rates in 4 ? 32 ? 32 cubes areestimated with CSE and the modified PDF. Again, the predictions arecompared to the data extracted from the DNS. Again, the assumptionof statistical homogeneity in the conditional averages throughout theflow field is adopted and the linear regularization with an unstrainedlaminar flame solution is used to solve for the conditional means. Thesolutions are presented in Fig. 6.9. As was seen in the RANS test, thepredicted conditional means are close to the DNS; the laminar flamesolution is barely distinguishable from the DNS.The chemical reaction rates are then modelled and compared withthe DNS results. In Fig. 6.10, the conditionally averaged reaction ratesare compared to the DNS. As in the RANS test, the reaction rates areover-predicted because of the errors in the conditional means of reac-tive scalars.The filtered chemical reactions rates in different cubes are evalu-ated by integration and plotted against the filtered DNS reaction rates.A constrained least squares analysis is performed to evaluate the pre-cision of the closure. As shown in Fig. 6.11, the modelled reaction ratesare close to the DNS results. The slopes of the fitted lines for the two re-1510 0.2 0.4 0.6 0.8 100.0050.010.0150.020.0250.030.0350.040.0450.05casteriskmath<YF|casteriskmath>(a) angbracketleftYF |casteriskmathangbracketright0 0.2 0.4 0.6 0.8<YO|casteriskmath>(b) angbracketleftYO|casteriskmathangbracketright0 0.2 0.4 0.6 0.8 100.0050.010.0150.020.0250.030.0350.04casteriskmath<YI|casteriskmath>(c) angbracketleftYI|casteriskmathangbracketright0 0.2 0.4 0.6 0.8 1600800100012001400160018002000casteriskmath<T|casteriskmath>(d) angbracketleftT|casteriskmathangbracketright0 0.2 0.4 0.6 0.8<rho|casteriskmath>(e) angbracketleftrho|casteriskmathangbracketrightFigure 6.6: Conditional averages calculated in RANS paradigm, ob-tained by using modified laminar flame PDF and different regular-ization methods. ?-: DNS; --: regularization with laminar flame so-lution (CSE-1); - ? -: regularization with iterations up to 2000 times(CSE-2); ? ? ? ? ?: laminar flame solution.1520 0.2 0.4 0.6 0.8<?omegai|casteriskmath>(a) angbracketleft ?i|casteriskmathangbracketright0 0.2 0.4 0.6 0.8<?omegaii|casteriskmath>(b) angbracketleft ?ii|casteriskmathangbracketrightFigure 6.7: Conditionally averaged reaction rates angbracketleft ? |casteriskmathangbracketright calculated inRANS paradigm. Conditional means obtained by using the modifiedPDF and regularization with the laminar flame solution. ?-: DNS; --:CSE.1530 0.2 0.4 0.6 0.8 angbracketleft ?iangbracketright0 0.2 0.4 0.6 0.8 angbracketleft ?iiangbracketrightFigure 6.8: Reaction rates angbracketleft ? angbracketright on planes in a premixed turbulent flame.Conditional means obtained by using the modified PDF and regular-ization with the laminar flame solution. ?-: DNS; --: CSE.1540 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.0050.010.0150.020.0250.030.0350.040.045casteriskmath<YF|casteriskmath>(a) angbracketleftYF |casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<YO|casteriskmath>(b) angbracketleftYO|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.0050.010.0150.020.0250.03casteriskmath<YI|casteriskmath>(c) angbracketleftYI|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 150010001500200025003000casteriskmath<T|casteriskmath>(d) angbracketleftT|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<rho|casteriskmath>(e) angbracketleftrho|casteriskmathangbracketrightFigure 6.9: Conditional averages calculated in LES paradigm, obtainedby using the modified PDF and regularization with laminar flame so-lution. ?-: DNS; --: CSE; ? ? ? ? ?: laminar flame solution.1550 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<?omegai|casteriskmath>(a) angbracketleft ?i|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9<?omegaii|casteriskmath>(b) angbracketleft ?ii|casteriskmathangbracketrightFigure 6.10: Conditionally averaged reaction rates angbracketleft ? |casteriskmathangbracketright, calculated inLES paradigm, conditional means obtained by using the modified PDFand regularization with laminar flame solution. ?-: DNS; --: CSE.156actions are 1.222 and 1.649; the Pearson product moment coefficients ofcorrelation are 0.9508 and 0.6924 respectively. The results seem promis-ing, although it is somewhat disappointing that there might be a posi-tive bias in the CSE prediction of reaction rates.Optimization of the PDF table might further improve the precisionof closure. For illustration, a slightly denser grid in conditional vari-able space with 46 points was used. The same procedure of generatinga PDF table, inverting for conditional means, invoking first conditionalmoment closure and integrating to obtain filtered reaction rates is re-peated. As shown in Fig. 6.12, the slopes of the fitted lines are reducedto 1.135 and 1.491, and the coefficients of correlation are increased to0.9614 and 0.7167, indicating enhanced precision of closure.For further analysis, the mean reaction rates are calculated usingthe conditional means extracted from the DNS and the modified PDF.The results are plotted against the DNS reaction rates in Fig. 6.13. Theslopes 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. Thisanalysis shows the error observed in Fig. 6.11 is mainly attributable toerrors in the conditional means. As discussed previously, in real CSEmodeling, the interaction between conditional means and the reactivescalars is likely to yield better solutions in the CSE inversion, whichmight improve the precision of closure as well. Therefore, the resultsshown in Fig. 6.13 give more confidence in the closure method withCSE and the modified PDF.157-0.2 0 0.2 0.4 0.6 0.8 1 DNS?omegai(a) angbracketleft ?iangbracketright-0.5 0 0.5 1 1.500. DNS?omegaii(b) angbracketleft ?iiangbracketrightFigure 6.11: Reaction rates angbracketleft ? angbracketright calculated in LES paradigm. ?-: ? =? DNS; --: ? i = 1.222 asteriskmath ? DNS,i and ? ii = 1.649 asteriskmath ? DNS,ii.158-0.2 0 0.2 0.4 0.6 0.8 1 DNS?omegai(a) angbracketleft ?iangbracketright-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.400. DNS?omegaii(b) angbracketleft ?iiangbracketrightFigure 6.12: Reaction rates angbracketleft ? angbracketright calculated in LES paradigm, with densergrid in c-space. ?-: ? = ? DNS; --: ? i = 1.1345 asteriskmath ? DNS,i and ? ii =1.4912 asteriskmath ? DNS,ii.1590 0.2 0.4 0.6 0.8 DNS?omegai(a) angbracketleft ?iangbracketright0 0.2 0.4 0.6 0.8 DNS?omegaii(b) angbracketleft ?iiangbracketrightFigure 6.13: Reaction rates angbracketleft ? angbracketright calculated in LES paradigm, evaluatedby using the modified PDF and conditional mean reaction rates ex-tracted directly from the DNS. ?-: ? = ? DNS; --: ? i = 0.9710 asteriskmath ? DNS,iand ? ii = 0.9684 asteriskmath ? DNS,ii.1606.4.2 Comparison with DNS: Mechanism IIThe LES tests described in the previous section is repeated for the tur-bulent flame solution II. The conditional means of reactive scalars andunconditional means of reaction rates are compared to the DNS dataas shown in Figs. 6.14 to 6.16.As observed in the RANS test II, the conditional means of the DNSare different from the laminar flame solutions; despite this imprecise apriori information, the conditional means obtained with the CSE inver-sion method are quite close to the DNS results. As observed in the LEStest I, the predicted reaction rates are consistent with the DNS data.The slopes of the fitted lines in Fig. 6.16 are 1.0454 and 1.0799 respec-tively. The tendency of overprediction is still exhibited in this series oftest; however, it is less significant than that of LES test I.6.5 ConclusionsA product-based progress variable is introduced as a conditioning vari-able for the CSE model. The mean reaction rates in premixed turbulentflames are modelled using the CSE method and the proposed PDF mo-del. The closure method is tested in both RANS and LES paradigms.In the a priori tests, the conditional means of reactive scalars are ob-tained using the CSE inversion method. Good estimates of conditionalmeans are obtained when the inversion equations are regularized witha laminar flame solution. The a priori tests further validate the modifiedversion of the laminar flame based PDF model and the closure methodof CSE. The predictions of the mean reaction rates agree reasonablywell with the DNS in both RANS and LES tests. The application of the1610 0.2 0.4 0.6 0.8 100.0050.010.0150.020.0250.030.0350.040.0450.05casteriskmath<YF|casteriskmath>(a) angbracketleftYF |casteriskmathangbracketright0 0.2 0.4 0.6 0.8<YO|casteriskmath>(b) angbracketleftYO|casteriskmathangbracketright0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.0050.010.0150.020.0250.030.0350.04casteriskmath<YI|casteriskmath>(c) angbracketleftYI|casteriskmathangbracketright0 0.2 0.4 0.6 0.8 1800100012001400160018002000casteriskmath<T|casteriskmath>(d) angbracketleftT|casteriskmathangbracketright0 0.2 0.4 0.6 0.8<rho|casteriskmath>(e) angbracketleftrho|casteriskmathangbracketrightFigure 6.14: Conditional averages calculated in LES paradigm, ob-tained by using the modified PDF and regularization with laminarflame solution. ?-: DNS; --: CSE; ? ? ? ? ?: laminar flame solution.1620 0.2 0.4 0.6 0.8<?omegai|casteriskmath>(a) angbracketleft ?i|casteriskmathangbracketright0 0.2 0.4 0.6 0.8<?omegaii|casteriskmath>(b) angbracketleft ?ii|casteriskmathangbracketrightFigure 6.15: Conditionally averaged reaction rates angbracketleft ? |casteriskmathangbracketright, calculated inLES paradigm, conditional means obtained by using the modified PDFand regularization with laminar flame solution. ?-: DNS, --: CSE.1630 0.2 0.4 0.6 0.8 angbracketleft ?iangbracketright0 0.2 0.4 0.6 0.8 angbracketleft ?iiangbracketrightFigure 6.16: Reaction rates angbracketleft ? angbracketright calculated in LES paradigm. ?-: ? =? DNS, --: ? i = 1.0454 asteriskmath ? DNS,i and ? ii = 1.0799 asteriskmath ? DNS,ii.164CSE closure method with the proposed PDF in premixed flames needsfurther tests in combustion simulations.165Chapter 7Conclusions and Future Work7.1 ConclusionsIn this thesis work, the Conditional Source-term Estimation methodhas been applied to study turbulent reacting flows. Firstly, CSE wascombined with Laminar Flamelet Decomposition and then with Tra-jectory Generated Low Dimensional Manifold method to study non-premixed autoigniting jets of methane and methane mixed with ad-ditives. Both methods were implemented using OpenFOAM and ap-plied to study jet flames with the same initial and boundary conditions,using the same grid and numerical scheme for comparison. The pre-dictions were compared to shock-tube experiments. The study showsthat the CSE-TGLDM method is superior to the CSE-LFD method; thismight be attributed to the wider range of solutions of detailed chemis-try included in the TGLDM manifolds, which also allows CSE-TGLDMto model transition from ignition to combustion, as opposed to theCSE-LFD method which only works up to the moment of ignition.The effects of different injection parameters on the ignition delaytime was discussed based on the predictions of the CSE-TGLDM me-thod. The fuel composition of the jet affects the ignition delay by chang-ing the chemical kinetics and fluid dynamics of the flame. Ethane ac-166celerates the autoignition of methane and nitrogen delays autoignition;however, the effects of modest quantities of these additives becomenegligible at high pre-combustion temperatures. Increasing injectionpressure ratio reduces ignition delay and the effects become negligi-ble at high pressure ratios. The injection duration is irrelevant whenignition starts before the end of injection.The predictions of ignition delay time are generally consistent withthe experimental data, which indicates the effectiveness of the CSEmethods and the chemical mechanism used in the jet simulations. How-ever, the scatter observed in the experimental data could not be cap-tured in the simulation. The scatter in ignition delay might be attributedto the different realizations of chemical reaction paths, as well as thefluctuations in turbulent jets, which cannot be captured by continuummanifold methods and Reynolds Averaged Navier-Stokes simulation.To simulate these effects, a TGLDM method combining continuum andstochastic particle model generated trajectories might be used to simu-late the randomness of chemical reaction paths; Large Eddy Simulationis suggested to capture the fluctuations in turbulent flows.The recommended CSE-TGLDM method was applied to obtain clo-sure in LES of a piloted, non-premixed methane jet flame, Sandia FlameD. In the thesis work, the adiabatic CSE-TGLDM method was extendedto model radiation in the flame. It is found that the CSE method can beused to model the turbulence/radiation interaction in a similar fash-ion to how it models the turbulence/chemistry interaction; the ma-jor sources of radiation in Sandia D flame are hot CO2 and H2O; thestrongest radiation in the flame is around the stoichiometric flame height;167the radiation heat loss of the flame is insufficient to make a big differ-ence in the predictions of temperature and major species; however, thepredictions of NO are improved by using the extended CSE-TGLDMmodel.The application of the CSE method in the study of premixed turbu-lent flames was explored. A conditioning variable based on the prod-uct mass fraction was proposed. A priori tests show that the presumedPDF of this progress variable may be well approximated with an ex-tended version of a laminar flame-based model. The precision of thePDF is sensitive to the chemical mechanism; it is suggested that thesame mechanism be used for the solutions of the prototype laminarflame and the turbulent flame.Using the conditioning variable and its presumed PDF model, theCSE method was adapted to evaluate conditional averages and closechemical source terms in premixed flames. A priori tests were perform-ed in both the RANS and LES paradigms. The results of the tests arepromising. However, the effectiveness and precision of chemical clo-sure with the CSE method in turbulent premixed flames need furthertests in flames of different configurations and ultimately in turbulentpremixed combustion simulations.7.2 Future WorkIn the study of unsteady autoigniting jets, the CSE-LFD method ap-pears to be inferior to the CSE-TGLDM method. However, it is possi-ble to improve CSE-LFD predictions by enhancing the flamelet libraryto include solutions of extinction and re-ignition phenomena. More-168over, the a priori tests indicate that the LFD method might work wellfor premixed flames. The potential of the CSE-LFD method is worthexploring further.It is also found that the current RANS simulation is incapable ofcapturing the distribution of ignition delay caused by the fluctuationsof chemical reactions and turbulence. LES of autoigniting jets with aTGLDM method combining the continuum approach and stochasticparticle model to create trajectories is suggested to model these fluctu-ations.In the study of steady jet flame Sandia flame D, the CSE-TGLDMmethod has been extended to model radiation. The storage require-ments and retrieval of multiple manifolds turn out to be demanding.These requirements are still manageable in the current work; however,it will be necessary to use more manifolds to model strongly radiatingflames, which poses a challenge for implementation. Another issue isthe precision of the manifold in the low temperature range where thephysical 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 seemsto be a promising tool to obtain chemical closure. In future study, othercandidates of the presumed probability density function for the pro-posed conditioning variable might be studied; improvements to theextended laminar flame based PDF might be made by filtering. Ul-timately, the CSE with presumed PDF method should be applied tosimulate turbulent premixed flames.169Bibliography[1] Information available at:, accessed 12 Dec.2007.[2] Code of federal regulations (cfr) 40 parts 69, 80, 86: Control of airpollution from new motor vehicles: heavy duty engine and ve-hicle standards and highway diesel fuel sulfur control require-ments; final rule, US Federal Register 66(12) (2001) 5002?5050.[3] On road vehicle and engine emissions regulations, CanadaGazette Part II 137(1) (2003) 21?27.[4] G. P. McTaggart-Cowan, C. C. O. Reynolds, W. K. Bushe, Natu-ral gas fuelling for heavy-duty on-road use: Current trends andfuture direction, International Journal of Environmental Studies63(4) (2006) 421?440.[5] Information available at:, accessed 12 Dec. 2007.[6] Information available at:, accessed 12 Dec. 2007.[7] W. K. Bushe, H. Steiner, Conditional moment closure for largeeddy simulation of nonpremixed turbulent reacting flows,Physics of Fluids 11 (1999) 1896?1906.[8] D. Veynante, L. Vervisch, Turbulent combustion modeling,Progress in Energy and Combustion Science 28 (2002) 193?266.[9] P. Moin, K. Mahesh, Direct numerical simulation: a tool in tur-bulence research, Annual Review of Fluid Mechanics 30 (1998)539?578.[10] V. Eswaran, S. B. Pope, An examination of forcing in directnumerical simulations of turbulence, Computational Fluids 16(1988) 257?278.[11] J. Kim, P. Moin, R. Moser, Turbulence statistics in fully developedchannel flow, Journal of Fluid Mechanics 177 (1987) 133?166.170[12] P. R. Spalart, Direct numerical study of leading-edge contamina-tion, in: Fluid Dynamics of Three-Dimensional Turbulent ShearFlows 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 eddysimulation modeling for turbulent combustion, in: O. M?tais,J. Ferziger (Eds.), New tools in turbulence modelling, Berlin: LesEditions de Physique, Springer, 1997, pp. 105?140.[16] J. Smagorinsky, General circulation experiments with the primi-tive equations, Monthly Weather Review 91 (3) (1963) 99?164.[17] M. Germano, U. Piomelli, P. Moin, W. H. Cabot, A dynamicsubgrid-scale eddy viscosity model, Physics of Fluids 3 (1991)1760?1765.[18] P. Moin, K. Sqires, W. Cabot, S. Lee, A dynamic subgrid-scalemodel for compressible turbulence and scalar transport, Physicsof Fluids 3 (1991) 2746?2757.[19] A. E. Tejada-Mart?z, K. E. Jansen, A dynamic smagorinsky modelwith dynamic determination of the filter width ratio, Physics ofFluids 16 (2004) 2514?2528.[20] H. Pitsch, Large-eddy simulation of turbulent combustion, An-nual Review of Fluid Mechanics 38 (2006) 453?482.[21] N. Peters, Turbulent combustion, 1st Edition, Cambridge Uni-versity Press, 2000.[22] R. W. Bilger, Turbulent jet diffusion flames, Progress of EnergyCombustion and Science (1976) 87?109.[23] P. A. Libby, F. A. Williams, Turbulent reacting flows, Academicpress London, 1994.[24] R. W. Bilger, Turbulent flows with non-premixed reactants, Top-ics Applied Physics 44 (1980) 65?113.171[25] S. Burke, T. Schumann, Diffusion flames, Industrial and Engi-neering Chemistry 20 (1928) 998?1004.[26] A. W. Cook, J. J. Riley, A subgrid model for equilibrium chemis-try in turbulent flows, Physics of Fluids 6 (1994) 2868?2870.[27] J. Coupland, C. H. Priddin, Modelling of flow and combustionin a production gas turbine combustor, Turbulent Shear Flows 5(1987) 310?323.[28] N. Peters, Laminar diffusion flamelet models in non-premixedturbulent combustion, Progress in Energy and Combustion Sci-ence 10 (3) (1984) 319?339.[29] H. Pitsch, N. Peters, A consistent flamelet formulation for non-premixed combustion considering differential diffusion effects,Combustion and Flame 114 (1998) 26?40.[30] N. Peters, Comment and reply on the ?Assessment of combus-tion and submodels for turbulent nonpremixed hydrocarbonflames?, Combustion and Flame 116 (1999) 675?676.[31] N. Swaminathan, R. Bilger, Assessment of combustion and sub-models for turbulent nonpremixed hydrocarbon flames, Com-bustion and Flame 116 (1999) 519?545.[32] N. Swaminathan, Flamelet regime in non-premixed combustion,Combustion and Flame 129 (2002) 217?219.[33] F. Mauss, D. Keller, N. Peters, A Lagrangian simulation offlamelet extinction and re-ignition in turbulent jet diffusionflames, Proceedings of the Combustion Institute 23 (1990) 693?698.[34] H. Barths, N. Peters, N. Brehm, et al., Simulation of pollutant for-mation in a gas turbine combustor using undteady flamelets, in:Proceedings of the Combustion Institute, The Combustion Insti-tute, 1998, pp. 1841?1847.[35] H. Barths, C. Hasse, Simulation of combustion in direct injectiondiesel engines using an eulerian particle flamelet model, Pro-ceedings of the Combustion Institute 28 (2000) 1161?1168.172[36] H. Pitsch, Unsteady flamelet modeling of differential diffusionin turbulent jet diffusion flames, Combustion and Flame 123 (3)(2000) 358?374.[37] H. Pitsch, H. Steiner, Large-eddy simulation of a turbulent pi-loted methane/air diffusion flame (Sandia flame D), Physics ofFluids 12 (2000) 2541?2554.[38] S. P. Nandula, T. M. Brown, R. W. Pitz, Measurements of scalardissipation in the reaction zones of turbulent nonpremixed h2-airflames., Combustion and Flame 99 (1994) 775?783.[39] R. W. Bilger, Conditional moment closure for turbulent reactingflows, Physics of Fluids 5 (1993) 436?444.[40] A. Y. Klimenko, Multicomponent diffusion of various admix-tures in turbulent flow, Fluid Dynamics 25 (3) (1990) 327?334.[41] Information available at:, ac-cessed 12 Dec. 2007.[42] A. Klimenko, R. Bilger, Conditional moment closure for tur-bulent combustion, Progress in energy and combustion science25 (6) (1999) 595?687.[43] S. Kim, K. Y. Huh, Use of the conditional moment closure modelto predict no formation in a turbulent CH4/H2 flame over a bluffbody, Combustion and Flame 130 (2002) 94?111.[44] S. H. Kim, K. Y. Huh, R. A. Fraser, Modeling autoignition of aturbulent methane jet by the conditional moment closure model,Proceedings of the Combustion Institute 28 (2000) 185?191.[45] C. Devaud, K. N. C. Bray, Assessment of the applicability of con-ditional moment closure model to a lifted turbulent flame: firstorder model, Combustion and Flame 132 (2003) 102?114.[46] S. H. Kim, K. Y. Huh, R. W. Bilger, Second-order conditionalmoment closure modeling of local extinction and reignition inturbulent non-premixed hydrocarbon flames, Proceedings of theCombustion Institute 29 (2002) 2131?2137.173[47] M. Roomina, R. Bilger, Conditional moment closure modelingof turbulent methanol jet flames, Combustion Theory and Mod-elling 3 (1999) 698?708.[48] M. Roomina, R. W. Bilger, Conditional moment closure (CMC)prediction of a turbulent methane-air flames, Combustion andFlame 125 (2001) 1176?1195.[49] M. Fairweather, R. M. Wooley, First-order conditional momentclosure modeling of turbulent, non-premixed methane flames,Combustion and Flame 138 (2004) 3?19.[50] S. Kim, K. Y. Huh, L. Tao, Application of the elliptic con-ditional moment closure model to a two-dimensional non-premixed methanol bluff-body flame, Combustion and Flame120 (75) (2000) 75?90.[51] N. Swaminathan, R. W. Bilger, Conditional variance equationand its analysis, Proceedings of the Combustion Institute (1998)1191?1198.[52] E. Mastorakos, R. W. Bilger, Second-order conditional momentclosure for the autoignition of turbulent flows, Physics of Fluids10 (1998) 1246?1248.[53] C. M. Cha, G. Kosaly, H. Pitsch, Modeling extinction and reig-nition in turbulent nonpremixed combustion using a doubly-conditional moment closure approach, Physics of Fluids 13 (12)(2001) 3824?3834.[54] A. Kronenburg, Double conditioning of reactive scalar transportequations in turbulent nonpremixed flames, Physics of Fluids16(7) (2004) 2640?2648.[55] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Nu-merical Recipes in Fortran, Cambridge University Press, Cam-bridge, UK, 1992.[56] H. Steiner, W. K. Bushe, Large eddy simulation of a turbulentreacting jet with conditional source-term estimation, Physics ofFluids 13 (2001) 754?769.174[57] R. Grout, W. K. Bushe, C. Blair, Predicting the ignition delay ofturbulent methane jets using conditional source-term estimation,To appear in Combustion Theory and Modelling.[58] R. Grout, Combustion modelling with conditional source-termestimation and laminar flamelet decomposition, Master?s thesis,University of British Columbia (2004).[59] M. Wang, J. Huang, W. K. Bushe, Simulation of a turbulent non-premixed flame using conditional source-term estimation withtrajectory generated low-dimensional manifold, in: Thirty-FirstSymposium (International) on Combustion, The Combustion In-stitute, 2006, pp. 1701?1709.[60] J. Huang, W. K. Bushe, Simulation of an igniting methane jet us-ing conditional source-term estimation with a trajectory gener-ated low-dimensional manifold, To appear in Combustion The-ory and Modelling.[61] W. K. Bushe, H. Steiner, Laminar flamelet decomposition for con-ditional source-term estimation, Physics of Fluids 15 (2003) 1564?1575.[62] M. Wang, Combustion modeling using conditional source-termestimation with flamelet decompositon and low dimensionalmanifolds, Ph.D. thesis, University of British Columbia (2006).[63] M. Wang, W. K. Bushe, Large eddy simulation of a turbulent non-premixed flame using conditional source-term estimation withlaminar flamelet decomposition, submitted to Physics of Fluids.[64] R. Borghi, Turbulent combustion modelling, Progress in Energyand Combustion Science 14 (1988) 245?292.[65] D. B. Spalding, Mixing and chemical reaction in steady con-fined turbulent flames, Thirteenth Symposium (International) onCombusion, 1971, pp. 649?657.[66] B. F. Magnussen, B. H. Hjertager, On mathematical models of tur-bulent combustion with special emphasis on soot formation andcombustion, in: Sixteenth Symposium (International) on Com-bustion, The Combustion Institute, 1977, pp. 719?729.175[67] J. B. Moss, K. N. C. Bray, A unified statistical model of the pre-mixed turbulent flame, Acta Astronaut 4 (1977) 291?319.[68] K. N. Bray, P. A. Libby, Recent developments in the BML mo-del of premixed turbulent combustion, Academic Press, London,1994.[69] T. Mantel, R. Borghi, A new model of premixed wrinkled flamepropagation based on a scalar dissipation equation, Combustionand Flame 96(4) (1994) 443?457.[70] A. Trouv?, T. J. Poinsot, The evolution equation for the flame sur-face density in turbulent premixed combustion, Journal of FluidMechanics 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 coformation in turbulent premixed combustion, 29th Symposium(International) on Combusion, 2002, pp. 1871?1879.[73] M. Z. Bodenstein, Eine theorie der photochemischen reaktions-geschwindigkeiten, Physics of Chemistry 85 (1913) 329?397.[74] N. Peters, R. Kee, The computation of stretched laminarmethane-air diffusion flames using a reduced four-step mech-anism, Combustion and Flame 68 (1) (1987) 17?30.[75] J. Warnatz, Concentration, pressure, and temperature-dependence of the flame velocity in hydrogen-oxygen-nitrogenmixtures, Combustion Science and Technology 26 (5-6) (1981)203?213.[76] M. D. Smooke, Reduced kinetic mechanisms and asymptotic ap-proximations for methane-air flames, in: M. D. Smooke (Ed.),Lecture Notes in Physics, Springer-Verlag, 1991, pp. 1?28.[77] N. Peters, B. Rogg (Eds.), Reduced kinetic mechanisms for appli-cations in combustion systems, Springer-Verlag, 1993.176[78] U. Maas, S. B. Pope, Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space, Combustion andFlame 88 (1992) 239?246.[79] Z. Ren, S. B. Pope, A. Vladimirsky, J. M. Guckenheimer, Applica-tion of the ICE-PIC method for the dimension reduction of chem-ical kinetics, Journal of Chemical Physics 124 (2006) 114111.[80] J. A. van Oijen, L. P. H. de Goey, Modelling of premixed coun-terflow flames using the flamelet-generated manifold method,Combustion Theory and Modelling 6 (2002) 463?478.[81] J. A. van Oijen, L. P. H. de Goey, Modeling of premixed laminarflame using flame-generated low-dimensional manifolds, Seven-teenth International Colloquium on the Dynamics of Explosionsand Reactive Systems (1999) 25?30.[82] U. Maas, S. B. Pope, Implementation of simplified chemical ki-netics based on intrinsic low-dimensional manifolds, Twenty-Fourth Symposium (International) on Combustion (1992) 103?112.[83] J. Huang, Natural gas combustion under engine-relevant condi-tions, Ph.D. thesis, University of British Columbia (2006).[84] J. Nafe, U. Maas, A general algorithm for improving ildms, Com-bustion Theory and Modelling 6 (2002) 697?709.[85] O. Gicquel, N. Darabiha, D. Thevenin, Laminar premixed hydro-gen/air counterflow flame simulations using flame prolongationof ILDM with differential diffusion, in: Twenty-Eighth Sympo-sium (International) on Combustion, 2000, pp. 1901?1908.[86] H. Bongers, J. A. V. Oijen, L. P. H. de Goey, Intrinsic low-dimensional manifold method extended with diffusion, Pro-ceedings of the Combustion Institute 29 (2002) 1371?1378.[87] S. Pope, Computationally efficient implementation of com-bustion chemistry using in situ adaptive tabulation, Combus-tion Theory and Modelling 1 (1997) 41?63.177[88] S. B. Pope, U. Maas, Simplifying chemical kinetics: Trajec-tory generated low-dimensional manifolds, Mechanical andAerospace Engineering Report: FDA 93-11, Cornell University.[89] S. Gordon, B. J. McBride, Computer program for calculation ofcomplex chemical equilibrium compositions and applications,NASA Reference Publication 1311.[90] W. C. Reynolds, The element potential method for chemicalequilibrium analysis: implementation in the interactive programstanjan, Stanford University Report ME 270 HO No.7.[91] S. B. Pope, Ten questions concerning the large-eddy simulationof turbulent flows, New Journal of Physics 6 (2004) 35?59.[92] R. Renka, Tripack: A constrained two-dimensional delaunay tri-angulation package, ACM Transactions on Mathematical soft-ware 22 (1996) 1?8.[93] M. Wang, A. Frisque, W. K. Bushe, Trajectory-generated low-dimensional manifolds generated using the stochastic particlemethod, to appear in Combustion Theory and Modelling.[94] D. T. Gillespie, Exact stochastic simulation of coupled chemicalreactions, Journal of Physics of Chemistry 81 (1977) 2340?2361.[95] D. T. Gillespie, A rigorous derivation of the chemical masterequation, Physica A 188 (1992) 404?425.[96] A. Frisque, J. Schnakenberg, J. Huang, W. K. Bushe, Stochasticsimulation of the variations in the autoignition delay time ofpremixed methane and air, Combustion Theory and Modelling10 (2) (2006) 241?256.[97] Information available at:, accessed 12 Dec. 2007.[98] G. D. Sullivan, J. Huang, T. X. Wang, W. K. Bushe, S. N. Rogak,Emissions variability in gaseous fuel direct injection compres-sion ignition combustion, Society of Automotive Engineers 114(2005) 2005?01?0917.178[99] P. Ouellette, Direct injection of natural gas for diesel engine fuel-ing, Ph.D. thesis, University of British Columbia, Canada (1996).[100] F. P. Ricou, D. B.Spaiding, Measurement of entrainment by axialsymmetric turbulent jet, Journal of Fluid Mechanics 11 (1961) 21?32.[101] Information available at:, accessed 12 Dec. 2007.[102] P. Ouellette, P. G. Hill, Turbulent transient gas injections, Journalof Fluids Engineering 122 (2000) 743?753.[103] N. W. and M. H. Davy and W. K. Bushe, Ignition and emissioncharacteristics of methane direct injection compression ignitioncombustion, 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-premixednatural gas combustion, submitted to Combustion Science andTechnology.[105] J. Huang, W. K. Bushe, Experimental and kinetic study of au-toignition in methane/ethane/air and methane/propane/airmixtures under engine-relevant conditions, Combustion andFlame 144 (2006) 74?88.[106] K. Akselvoll, P. Moin, Large-eddy simulation of turbulent con-fined coannular jets, Journal of Fluid Mechanics 315 (1996) 387?411.[107] P. M. Gresho, Incompressible fluid dynamics: some fundamentalformulation issues, Annual Review of Fluid Mechanics 23 (1991)413?453.[108] B. J. Boersma, G. Brethouwer, F. T. M. Nieuwstadt, A numericalinvestigation on the effect of the inflow conditions on the self-similar region of a round jet, Physics of Fluids 10(4) (1998) 899?910.[109] Information available at:, 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, NewYork, 1960.[112] F. C. Lockwood, N. G. Shah, A new radiation solution methodforincorporation in general combustion prediction procedures,Eighteenth Symposium (International) on Combusion, 1981, pp.1405?1414.[113] P. J. Coelho, O. J. T. and D. Roekaerts, Spectral radiative ef-fects and turbulence-radiation-interaction in Sandia Flame D, in:TNF6 Proceedings, 2002.[114] P. J. Coelho, O. J. Teerling, D. Roekaerts, Spectral radiative effectsand turbulence/radiation interaction in a non-luminous turbu-lent jet diffusion flame, Combustion and Flame 133 (April 2003)75?91(17).[115] W. L. Grosshandler, Radcal: A narrow-band model for radiationcalculations in a combustion environment, NIST Technical Note1402.[116] Information available at:, accessed 12Dec. 2007.[117] B. Jin, M. Wang, W. K. Bushe, Modelling of mean spectral ra-diation intensities in sandia flame d, Proceedings of 2007 CICSSpring Technical Meeting, Banff, Canada.[118] M. Hossain, J. C. Jones, W. Malalasekera, Modelling of a bluff-body nonpremixed flame using a coupled radiation/flameletcombustion model, Flow, Turbulence and Combustion 67 (2001)217?234.[119] Information available at:, accessed 12 Dec.2007.180[120] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eit-eneer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song,W. C. Gardiner, V. V. Lissianski, Z. Qin, Information availableat: mech/, accessed 12 Dec.2007.[121] M. E. Kounalakis, Y. R. Sivathanu, G. M. Faeth, Infrared radia-tion statistics of nonluminous turbulent diffusion flames, ASMEJournal of Heat Transfer 113 (1992) 437?445.[122] Y. Zheng, R. S. Barlow, J. P. Gore, Measurements and calculationsof spectral radiation intensities for turbulent non-premixed andpartially premixed flames, ASME Journal of Heat Transfer 125(2003) 678?686.[123] S. C. Kong, R. D. Reitz, Numerical study of premixed hcci en-gine combustion and its sensitivity to computational mesh andmodel uncertainties, Combustion Theory and Modelling (2003)417?433.[124] R. W. Bilger, Turbulent flows with non-premixed reactants, Top-ics Applied Physics 44 (1980) 65?113.[125] P. A. Libby, F. A. Williams, A presumed pdf analysis of partiallypremixed turbulent combustion, Combustion Science and Tech-nology 161 (2001) 351?390.[126] F. A. Jaberi, R. S. Miller, P. Givi, Conditional statistics in turbulentscalar mixing and reaction, AlChE Journal 42 (1996) 1149?1152.[127] F. A. Jaberi, R. S. Miller, P. Givi, Non-gaussian scalar statistics inhomogeneous turbulence, Journal of Fluid Mechanics 313 (1996)241?282.[128] G. Ribert, M. Champion, P. Plion, Modelling turbulent reactiveflows with variable equivalence ratio: application to the calcula-tion of a reactive shear layer, Combustion Science and Technol-ogy 176(5/6) (2004) 907?924.[129] P. Domingo, S. P. L. Vervisch, R. Hauguel, DNS of a premixedturbulent v flame and les of a ducted flame using a fsd-pdf sub-grid scale closure with fpi-tabulated chemistry, Combustion andFlame 143(4) (2005) 566?586.181[130] N. Swaminathan, R. W. Bilger, Analysis of conditional momentclosure for turbulent premixed flames, Combustion Theory andModelling 5 (2001) 241?260.[131] S. M. Martin, J. C. Kramlich, G. Kosaly, J. J. Riley, The premixedconditional moment closure method applied to idealized leanpremixed gas turbine combustors, Journal of Engineering ForGas Turbines and Power-Transactions of the ASME 125 (2003)895?900.[132] K. N. C. Bray, M. Champion, P. A. Libby, N. Swaminathan, Finiterate chemistry and presumed pdf models for premixed turbulentcombustion, Combustion and Flame 146 (2006) 665?673.[133] R. Grout, An age extended progress variable for conditioning re-action rates, Submitted to Physics of Fluids.[134] K. W. Jenkins, R. S. Cant, Flame kernel interaction in turbulentenvironment, Third AFOSR Conference on DNS and LES, 1999,pp. 605?612.182


Citation Scheme:


Usage Statistics

Country Views Downloads
United States 7 1
China 4 1
United Kingdom 3 1
Germany 2 2
Canada 2 0
France 2 0
Japan 2 0
Vietnam 1 0
City Views Downloads
Beijing 4 1
Unknown 4 2
Abilene 3 0
Ashburn 2 0
Québec 2 0
Tokyo 2 0
Oxford 2 0
Hanoi 1 0
Sunnyvale 1 0
Washington 1 0
Stuttgart 1 1

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items

Admin Tools

To re-ingest this item use button below, on average re-ingesting will take 5 minutes per item.


To clear this item from the cache, please use the button below;

Clear Item cache