Numerical Simulation of Turbulent Premixed Flames with Conditional Source-Term Estimation by Mohammad Mahdi Salehi B.Sc. Sharif University of Technology, 2004 M.Sc. Sharif University of Technology, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2012 c© Mohammad Mahdi Salehi 2012 Abstract Conditional Source-term Estimation (CSE) is a closure model for turbulence- chemistry interactions. This model is based on the conditional moment closure hypothesis for the chemical reaction source terms. The conditional scalar field is estimated by solving an integral equation using inverse meth- ods. CSE was originally developed for - and has been used extensively in - non-premixed combustion. This work is the first application of this com- bustion model to predictive simulations of turbulent premixed flames. The underlying inverse problem is diagnosed with rigorous mathematical tools. CSE is coupled with a Trajectory Generated Low-Dimensional Manifold (TGLDM) model for chemistry. The CSE-TGLDM combustion model is used with both Reynolds-Averaged Navier-Stokes (RANS) and Large-Eddy Simulation (LES) turbulence models to simulate two different turbulent pre- mixed flames. Also in this work, the Presumed Conditional Moment (PCM) turbulent combustion model is employed. This is a simple flamelet model which is used with the Flame Prolongation of ILDM (FPI) chemistry re- duction technique. The PCM-FPI approach requires a presumption for the shape of the probability density function of reaction progress variable. Two shapes have been examined: the widely used β-function and the Modified Laminar Flamelet PDF (MLF-PDF). This model is used in both RANS and large-eddy simulation of a turbulent premixed Bunsen burner. Radial dis- tributions of the calculated temperature field, axial velocity and chemical species mass fraction have been compared with experimental data. This comparison shows that using the MLF-PDF leads to predictions that are similar, and often superior to those obtained using the β-PDF. Given that the new PDF is based on the actual chemistry – as opposed to the ad hoc nature of the β-PDF – these results suggest that it is a better choice for the statistical description of the reaction progress variable. ii Preface The contents of chapter 3 are published in: • M. M. Salehi and W. K. Bushe. Presumed PDF modeling for RANS simulation of turbulent premixed flames. Combustion Theory and Modelling, 14:381-403, 2010. In this publication, I was responsible for conducting all parts of the research and preparing the manuscript. Dr. Bushe supervised the research and pro- vided feedback and reviewed the manuscript. Part of the contents of chapter 4 are published in: • M. M. Salehi, W. K. Bushe and K. J. Daun. Application of Conditional Source-term Estimation model for Turbulence-Chemistry Interactions in a Premixed Flame. Combustion Theory and Modelling, 16:301-320, 2012. In this publication, I was responsible for conducting all parts of the research and preparing the manuscript. Dr. Bushe supervised the research and pro- vided feedback and reviewed the manuscript. Dr. Daun also provides some valuable feedback and reviewed the manuscript. In chapter 5 and 6 the CFFC computer program was used. This program was developed at the University of Toronto under the supervision of Dr. Groth. I modified this program in collaboration with Ms. Shahbazian for the specific case of this research. Ms. Shahbazian provided help with us- ing this code and also prepared the chemistry library. I was responsible for changing the simulation case to the Bunsen burner studied in this work and changing the boundary conditions, especially the turbulent inlet condition. Also, I modified the code to be able to deliver the final statistical analysis results. Ms. Shahbazian also prepared the TGLDM chemistry library and added the extra transport equation for CSE to the existing flamelet solver. I was responsible for implementation of the CSE routines and the presumed PDF model to the code and running the two different cases and comparing iii Preface the results. Part of the contents of chapter 5 are accepted for publication: • M. M. Salehi, W. K. Bushe, N. Shahbazian and C. P. T. Groth. Modi- fied Laminar Flamelet Presumed Probability Density Function for LES of Premixed Turbulent Combustion. Proceedings of the Combustion Institute, in press, 2012. In this publication, I was responsible for preparing the manuscript. Dr. Bushe supervised the research and provided feedback and reviewed the manuscript. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Symbols and Abbreviations . . . . . . . . . . . . . . . . . . . . . xvi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.1 Reynolds-Averaged Navier-Stokes . . . . . . . . . . . 6 2.1.2 Large-Eddy Simulation . . . . . . . . . . . . . . . . . 7 2.2 Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Reduced Chemistry . . . . . . . . . . . . . . . . . . . 11 2.2.2 Laminar Flames with Detailed Chemistry . . . . . . . 13 2.3 Turbulence-Chemistry Interactions . . . . . . . . . . . . . . . 17 2.3.1 Regimes in Premixed Turbulent Combustion . . . . . 17 2.3.2 Turbulent Combustion Modelling in RANS . . . . . . 22 2.3.3 Turbulent Combustion Modelling in LES . . . . . . . 26 3 PCM-FPI Model in RANS Simulation of Turbulent Pre- mixed Flames . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 v Table of Contents 3.1.1 Presumed PDF . . . . . . . . . . . . . . . . . . . . . 28 3.1.2 Transport Equations for c̃ and c̃′′2 . . . . . . . . . . . 32 3.2 Summary of the Governing Equations . . . . . . . . . . . . . 38 3.3 Validation Case . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.4 Numerical Simulation Code . . . . . . . . . . . . . . . . . . . 40 3.5 Computational Domain and Boundary Conditions . . . . . . 41 3.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.6.1 Effect of Different Scalar Dissipation Rate Models . . 43 3.6.2 Turbulence Model . . . . . . . . . . . . . . . . . . . . 44 3.6.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 48 3.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4 CSE in RANS Simulation of Turbulent Premixed Flames 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2.1 Conditional Source-term Estimation . . . . . . . . . . 61 4.2.2 Ill-posed Inverse Problem . . . . . . . . . . . . . . . . 63 4.2.3 Parameter-Choice Methods . . . . . . . . . . . . . . . 66 4.2.4 Nonlinearity . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.5 Chemistry Reduction . . . . . . . . . . . . . . . . . . 67 4.2.6 Governing Equations . . . . . . . . . . . . . . . . . . 69 4.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.3.1 Validation Case . . . . . . . . . . . . . . . . . . . . . 72 4.3.2 Simulation Code . . . . . . . . . . . . . . . . . . . . . 73 4.3.3 Grid Resolution Study . . . . . . . . . . . . . . . . . 74 4.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . 74 4.4.1 Regularization and Inversion Results . . . . . . . . . 75 4.4.2 Simulation vs. Experiment . . . . . . . . . . . . . . . 79 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5 PCM-FPI in LES Simulation of Turbulent Premixed Flames 91 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.2.1 Filtered Laminar Flame Speed . . . . . . . . . . . . . 93 5.2.2 LES SGS Closures . . . . . . . . . . . . . . . . . . . . 96 5.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.3.1 Turbulent Inlet Boundary Condition . . . . . . . . . . 98 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . 100 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 vi Table of Contents 6 CSE in LES Simulation of Turbulent Premixed Flames . . 116 6.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.1.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.3.1 Test Case –I . . . . . . . . . . . . . . . . . . . . . . . 119 6.3.2 Test Case – II . . . . . . . . . . . . . . . . . . . . . . 120 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . 128 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 7.1.1 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . 128 7.1.2 Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . 128 7.1.3 Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . 129 7.1.4 Chapter 6 . . . . . . . . . . . . . . . . . . . . . . . . 129 7.1.5 Summary of Accomplishments . . . . . . . . . . . . . 129 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Appendix A: Bayesian Inverse Problem . . . . . . . . . . . . . 148 Appendix B: Second-order CMC . . . . . . . . . . . . . . . . . . 151 vii List of Tables 3.1 Model coefficients . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Characteristics of the turbulent Bunsen burner of Chen et al.[28] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 viii List of Figures 2.1 A sample of range of time scales in a turbulent reacting flow. 11 2.2 Schematic representation of low-dimensional manifold tech- nique for chemistry reduction. . . . . . . . . . . . . . . . . . . 12 2.3 Schematic of two different types of reacting flows: non-premixed and premixed flames. . . . . . . . . . . . . . . . . . . . . . . . 14 2.4 Chemical structure of one-dimensional unstrained laminar pre- mixed flame for a lean methane-air mixture. . . . . . . . . . . 15 2.5 Schematic representation of laminar flame thickness definition. 16 2.6 Schematic of a stagnation point laminar premixed flame. . . . 17 2.7 Interaction between turbulence and premixed flame. . . . . . 18 2.8 Different regimes in turbulent premixed combustion [10]. . . . 19 2.9 Modified turbulent premixed combustion regime diagram [116]. 21 3.1 Realizable values of mean and variance of reaction progress variable. The shaded area shows the applicable range of the PDF of equation 3.2 for a stoichiometric methane-air flame. Symbols denote the experimental data for Flame F3 of the Bunsen burner of Chen et al. [28] at axial positions x/D = 8.5 (+), x/D = 6.5 (♦), x/D = 4.5 () and x/D = 2.5 (◦). . 30 3.2 Four possible cases of the modified laminar flamelet PDF . . 31 3.3 Range of the applicability of each of the four possible cases of the modified laminar flamelet PDF for a lean mixture of methane and air at equivalence ratio of 0.6. . . . . . . . . . . 33 3.4 Range of the applicability of each of the four possible cases of the modified laminar flamelet PDF for a stoichiometric methane-air mixture. . . . . . . . . . . . . . . . . . . . . . . . 34 3.5 Gradient vs. counter-gradient diffusion in turbulent premixed flames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.6 Schematic of the computational domain and boundary con- ditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.7 Computational grid in OpenFOAM . . . . . . . . . . . . . . . 42 ix List of Figures 3.8 Effect of changing the time scale ratio on axial distribution of reaction progress variable for flame F3. . . . . . . . . . . . 44 3.9 Radial profile of several quantities calculated by using the algebraic model of equation 3.17 for scalar dissipation rate. The results are for flame F3 at x/D = 8.5. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.10 Radial profile of several quantities calculated by using the algebraic model of equation 3.15 for scalar dissipation rate. The results are for flame F3 at x/D = 8.5. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.11 Radial profile of mean axial velocity and mean turbulent kinetic energy for cold and reacting flow of Flame F3 at x/D = 8.5. Symbols denote experimental data [28], “-·-” is Standard k−ε, “—” is Standard k−ε with Pope’s correction and “- -” is Standard k − ε with Cε1 = 1.6 . . . . . . . . . . . 47 3.12 The effect of using the model of Zhang and Rutland [170] in the Standard k − ε turbulence model with Cε1 = 1.6 in reacting flow simulation of flame F3 at three axial locations. Symbols denote experimental data [28]. . . . . . . . . . . . . 48 3.13 Temperature contours in flame F3. . . . . . . . . . . . . . . . 49 3.14 Production rate of reaction progress variable in flame F3. . . 50 3.15 Contours of variance of reaction progress variable in flame F3. 51 3.16 Radial profile of mean axial velocity. Symbols denote exper- imental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 52 3.17 Radial profile of mean turbulent kinetic energy. Symbols de- note experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.18 Radial profile of reduced temperature. Tb is the adiabatic flame temperature, Tb = 2248. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . 54 x List of Figures 3.19 Radial profile of methane mass fraction. Symbols denote ex- perimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.20 Radial profile of O2 mass fraction. Symbols denote experi- mental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 56 3.21 Radial profile of CO2 mass fraction. Symbols denote experi- mental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 57 3.22 Radial profile of water mass fraction. Symbols denote exper- imental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 58 3.23 Radial profile of CO mass fraction. Symbols denote experi- mental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 59 4.1 Filter functions for Tikhonov and TSVD regularization methods 66 4.2 Trajectory Generated Low-Dimensional Manifold. Different trajectories are generated using the solution of an unstrained laminar premixed flame. . . . . . . . . . . . . . . . . . . . . . 69 4.3 H2 mass fraction for different trajectories. . . . . . . . . . . . 70 4.4 CH4 reaction rates for different trajectories. . . . . . . . . . . 71 4.5 Interaction between different physical sub-models in a pre- mixed CSE-TGLDM reacting flow simulation code. . . . . . . 73 4.6 Grid convergence test . . . . . . . . . . . . . . . . . . . . . . 75 4.7 Singular Values of the kernel of the inverse problem for the initial and converged field. . . . . . . . . . . . . . . . . . . . . 76 4.8 L-curves for regularization parameters. The solid line repre- sents the Tikhonov regularization using the laminar flamelet solution for ~α0 and the symbols are for TSVD. . . . . . . . . 77 4.9 L-curves for regularization parameters. The solid line repre- sents the Tikhonov regularization without any guess for the solution, i.e. ~α0 = 0 and the symbols are for TSVD. . . . . . 78 4.10 Conditional average of water mass fraction. Symbols denote the unstrained laminar flamelet solution, “—” is the result of using Tikhonov regularization with the flamelet solution as an initial guess, “-·-” is obtained from using the Tikhonov regularization without any initial guess and “- -” depicts the result of TSVD approach. . . . . . . . . . . . . . . . . . . . . 79 xi List of Figures 4.11 Axial velocity in a strained laminar premixed flame. . . . . . 80 4.12 Effect of increasing aerodynamic strain on Conditional av- erage of H2O mass fraction in a laminar premixed flame of stoichiometric methane-air mixture in RtP configuration. . . 81 4.13 Radial profiles of mean axial velocity at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. U0 is the mean axial velocity at the inlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.14 Radial profiles of mean turbulent kinetic energy at four differ- ent axial locations for turbulent Bunsen flame F3 [28]. Sym- bols denote experimental data and solid lines are the simu- lation results using CSE combustion model. k0 is the mean turbulent kinetic energy at the inlet. . . . . . . . . . . . . . . 83 4.15 Radial profiles of mean temperature at four different axial locations for turbulent Bunsen flame F3 [28], Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. Tb is the temperature of the hot pilot products at the inlet. . . . . . . . . . . . . . . . . . 84 4.16 Radial profiles of mean CH4 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . 85 4.17 Radial profiles of mean O2 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . . . . . 86 4.18 Radial profiles of mean CO2 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . 87 4.19 Radial profiles of mean H2O mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . 88 4.20 Radial profiles of mean OH mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . 89 xii List of Figures 4.21 Radial profiles of mean CO mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. . . . . . . . . . . . . . . 90 5.1 Four possible shapes for the Modified Laminar Flamelet PDF for lean methane-air premixed flame. Sc is the segregation factor Sc = c̃v/(c̃(1− c̃)) . . . . . . . . . . . . . . . . . . . . 94 5.2 Schematic of laminar vs. filtered flame speed. . . . . . . . . . 95 5.3 Filtered laminar flame speed for different filter sizes using β- PDF and MLF-PDF. `F is the laminar flame thickness. . . . 96 5.4 Velocity and vorticity contours in flame F3. . . . . . . . . . . 99 5.5 The temporal variation of the L2-norm of the residual of the velocity field in LES of a periodic pipe. . . . . . . . . . . . . . 100 5.6 Contours of mean velocity and sub-grid scale turbulent kinetic energy in the LES simulation of a periodic pipe . . . . . . . . 101 5.7 Crosswise velocity vectors in LES of a periodic pipe. . . . . . 102 5.8 Mean velocity profiles of LES simulation of a periodic pipe. . 103 5.9 Instantaneous contour of density in kg/m3 for flame F3. . . . 105 5.10 Instantaneous velocity and vorticity contours in flame F3. . . 106 5.11 Instantaneous contours of two radical species in flame F3. . . 107 5.12 Radial profiles of mean axial velocity at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. U0 is the mean axial velocity at the inlet. . . . . . . . . . . . 108 5.13 Radial profiles of mean temperature at four different axial lo- cations for turbulent bunsen flame F3 [28]. Tb is the adiabatic flame temperature, Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . 109 5.14 Radial profiles of mean CH4 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β- PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.15 Radial profiles of mean O2 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 111 xiii List of Figures 5.16 Radial profiles of mean CO2 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β- PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.17 Radial profiles of mean H2O mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β- PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.18 Radial profiles of mean OH mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β- PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.19 Radial profiles of mean CO mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β- PDF and “—” is the result of using the modified laminar flamelet PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.1 Location of the Chen burner on LES turbulent premixed regime diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.2 Singular values of the design matrix A for Chen burner at several different axial locations. . . . . . . . . . . . . . . . . . 120 6.3 The conditional averges in Chen burner at different axial lo- cations. Symbolds denote the unstrained laminar flamelet results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.4 Location of the Gulder burner on the premixed LES regime diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.5 Singular values of the design matrix A for Gulder burner at several different axial locations. . . . . . . . . . . . . . . . . . 123 6.6 Filtered control variable c2 conditioned on c1 for the Gulder burner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.7 c2 vs. c1 in the counter-flow premixed flame configuration. . . 124 6.8 Instantaneous velocity and vorticity contours in Gulder burner.125 6.9 Instantaneous temperature and methane mass fraction con- tours in Gulder burner. . . . . . . . . . . . . . . . . . . . . . 126 xiv List of Figures 6.10 Radial profiles of the normalized temperature c̃T = (T̃ − Tu)/(Tb − Tu) at two axial locations. Solid line is the aver- aged LES simulation results and symbols denote experimental results [169]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 xv Symbols and Abbreviations Alphanumeric Symbols Symbol Meaning Unit A Design matrix A Area m2 aT Tangential strain rate 1/s ~b Unconditional average vector in CSE ensemble B Bray number c Progress variable cp Heat capacity J/(kg.K) c̃v Sub-grid scale variance of c D Central jet diameter m D Molecular diffusion constant m2/s Da Damköhler number G LES filter function, G-field in the G-equation h Total enthalpy J/kg J Molecular diffusion flux kg/(m2s) k Turbulent kinetic energy m2/s2 Ka Karlovitz number La Markstein length m ` Large eddy length scale m `F Laminar flame thickness m `δ Flame inner layer thickness m ~n Flame front normal vector p Pressure Pa P Probability density function Pr Prandtl number Re Turbulent Reynolds number Sc Schmidt number Sd Flame displacement speed m/s Sij Rate of strain tensor 1/s xvi Symbols and Abbreviations Symbol Meaning Unit SL Laminar flame speed m/s t Time s tT Large-scale turbulent time scale s tF Flame time scale s T Temperature K ui Velocity component in i-direction m/s Yk Mass fraction of the k-th species xi Spatial coordinate in i-th direction m ~x Spatial coordinate, ~x = (x1, x2, x3) m Greek Symbols Symbol Meaning Unit ~α The Conditional average in discrete form ~α0 a priori knowledge of the Conditional average in discrete form ε Rate of dissipation of turbulent kinetic energy m2/s3 c Scalar dissipation rate 1/s ρ Density kg/m3 τij Viscous stress tensor kg/(m.s2) τb Heat release index τc Unresolved turbulent scalar flux kg/(m2.s) µ Coefficient of laminar viscosity kg/(m.s) µt Coefficient of turbulent viscosity kg/(m.s) ν Coefficient of laminar kinematic viscosity m2/s νt Coefficient of turbulent kinematic viscosity m2/s η Kolmogorov length scale m ∆ LES filter scale m σij Sub-filter scale turbulent stress tensor kg/(m.s2) σi i-th Singular value Σ Flame surface density 1/s δij Kronecker delta function δ Dirac delta function ω̇ Chemical production rate kg/(m3.s) λ Thermal conductivity J/(m.s.k) xvii Symbols and Abbreviations Operators Symbol Meaning c Instantaneous value of c in a turbulent field c∗ Statistical random variable associated with c c Time-averaged (mean) value of c c′ Fluctuating component c′ = c− c 〈c〉 Spatially-filtered value of c c̃ Favre-averaged or Favre-filtered value of c c′′ Favre fluctuating component c′′ = c− c̃ Abbreviations Abbreviation Meaning BML Bray Moss Libby CSE Conditional Source-term Estimation CFD Computational Fluid Dynamics CMC Conditional Moment Closure DNS Direct Numerical Simulation FDF Filtered Density Function FGM Flame-Generated Manifold FPI Flame Prolongation of ILDM FSD Flame Surface Density ILDM Intrinsic Low-Dimensional Manifold LES Large-Eddy Simulation MLF Modified Laminar Flamelet MMC Multiple Mapping Conditioning PCM Presumed Conditional Moment PDF Probability Density Function RANS Reynolds-Averaged Navier-Stokes SGS Sub-Grid Scale SVD Singular Value Decomposition TGLDM Trajectory Generated Low-Dimensional Manifold TSVD Truncated Singular Value Decomposition xviii Acknowledgements The past five years have been one of the most memorable parts of my life. During these years, I have received help and support from many people to whom I would like to express my sincere appreciation. First of all, I would like to thank my supervisor Prof. Kendal Bushe who patiently introduced me the challenging field of turbulent combustion. I would like to acknowledge his continuous support, advice and patience. Without his help, this work would have not been accomplished. I would like to thank Dr. Carl Ollivier-Gooch, Dr. Martin Davy and Dr. James Feng, who were in my research committee and spent time reading my proposal and the first draft of this thesis. I would also like to acknowledge the financial support of NSERC, MPRIME and Rolls-Royce Canada. Computations were performed on different super- computers of SciNet and Westgrid consortia of Compute Canada. I would like to thank all my colleagues and friends at Clean Energy Research Center (CERC) for their valuable helps especially Amirabbas Ali- abadi, Dr. Chan, Hong Tsui, Eric Kastanis, Ehsan Faghani, Girish Nivarti, Dr. Saunders, Dr. Shield, Jean Logan and Arminta Chicka. Thanks to Gra- ham Hendra for proofreading part of the thesis. I would also like to thank the members of the CFD and Propulsion group of University of Toronto Institute for Aerospace Studies (UTIAS) especially Prof. Groth, Ms. Shah- bazian, Dr. Charest, Dr. Jha, Dr. Hernandez-Perez and Dr. Ivan. I would also like to thank my parents and all my friends at the University of British Columbia for their support and the great times that I had with them, in particular Sina Chiniforoosh, Mohammadsadegh Mashayekhi and Mohammad Mokmeli. Finally, I would like to thank my wonderful wife Saeedeh for her patience, understanding and support during these years. xix Chapter 1 Introduction Combustion of fossil fuels has been the primary source of energy for heating, cooking, transportation and industrial applications for many years. Coal, oil and natural gas are the main fossil fuel resources that are limited in supply. Also, uncontrolled and incomplete combustion of such fuels results in a high level of pollutants such as oxides of nitrogen (NOx), Carbon monoxide, soot and unburnt hydrocarbons. Moreover, combustion of fossil fuels produce carbon dioxide, which is a greenhouse gas that contributes to global warm- ing. The world has agreed to reduce the emissions of greenhouse gases in order to control global warming and limit climate change problems [153]. Although fossil fuel resources are limited, the explored reserves are not going to be depleted for at least a century at the current consuming rate [166]. Also, biofuel and biomass are renewable energy resources that can be used in combustion devices to replace fossil fuels. Furthermore, new technolo- gies are emerging for capture and storage of the carbon dioxide emissions from biofuel and fossil fuel combustion [1]. As a result, the environmental side effects of combustion are reduced and this source of energy becomes more sustainable. On the other hand, other sources of energy might not be effectively sustainable [8]. Pollutants especially NOx, CO and particulate matter (PM) are ma- jor side effects of combustion that can cause serious health problems [104]. Hence, governments set emission standards to limit the amount of pollu- tants produced by transportation vehicles and industry [132]. Combustion pollutions can be reduced by new designs and technologies. Premixed and partially-premixed burning is a promising way to accomplish this goal [164]. In premixed combustion, the combustion temperature can be controlled and reduced especially in lean burning. As a result, the NOx emissions reduce substantially. Also, CO and particulate matter are the result of incomplete combustion that can be reduced in efficient lean premixed combustors. These devices still involve major design challenges. For example, in lean- burn industrial gas turbines - used for stationary power - a common problem in designing the combustor is that there are thermo-acoustic instabilities that must be avoided to prevent catastrophic failures. 1 1.1. Objectives Numerical simulation has been proven to be a helpful tool in improv- ing designs of combustion devices. Development of computer models for turbulent combustion simulation is very challenging due to the inherent multi-physics and multi-scale nature of the problem. Both combustion and turbulence are multi-scale phenomena and it is not possible to resolve all the scales in practical applications due to limited computational resources. Hence, numerical simulation of turbulent reacting flows requires reduced mathematical modelling and the effect of unresolved scales needs to be mod- elled properly to obtain accurate - and hence useful - results. 1.1 Objectives The purpose of this work was to develop computer codes for simulation of turbulent premixed flames based on advanced modelling approaches for this multi-scale and muli-physics phenomena. These new modelling techniques can be summarized as follows: • Application of the modified laminar flamelet probability density func- tion (MLF-PDF) as a new presumd PDF model for the statistics of the reaction progress variable in turbulent premixed flames. This pre- sumed PDF model is a key sub-model in flamelet combustion models and the conditional source-term estimation (CSE) model. The results of using this model in both RANS and LES context are presented. • Application of the conditional source-term estimation (CSE) turbulent combustion model to turbulent premixed combustion. It is shown that this model is stable and converges to meaningful results using RANS and LES turbulence models. 1.2 Outline An introduction about turbulence, chemistry and the interaction between these two multi-scale phenomena is presented in chapter 2. This introduc- tory chapter mainly focuses on the physical aspects of these phenomena and different modelling approaches in the literature. Chapter 3 is about numeri- cal simulation of three premixed flames in a Bunsen burner. In this chapter, a simple RANS approach is used for turbulence modelling and a flamelet model is used for the turbulent combustion model. This flamelet model is based on a presumed PDF where two different PDFs are used and the re- sults are compared. In chapter 4, the conditional source-term estimation 2 1.2. Outline (CSE) turbulent combustion model is used with a RANS turbulence model. The theory section in this chapter deals with the mathematical challenges of using CSE in turbulent premixed combustion. Chapter 5 is similar to chap- ter 3 with an LES turbulence model instead of RANS. The challenges in us- ing LES such as turbulent inlet boundary conditions and high-performance computing are explained in this chapter. Chapter 6 is about application of the CSE turbulent combustion model in large-eddy simulation of turbulent premixed flames. Finally, chapter 7 provides conclusions and discusses op- portunities for future work that have remained out of the research presented in this thesis. 3 Chapter 2 Background Simulation of turbulent reacting flows requires understanding and modelling turbulence, chemistry and the interaction between these two phenomena at different scales. In this chapter, the physics of these phenomena are briefly explained and different modelling strategies for turbulence, chemistry and their correlation are explained. The emphasis is on the interaction between turbulence and chemistry, for which a more comprehensive literature review is presented. 2.1 Turbulence Turbulence is a highly non-linear and multi-scale phenomena. The complex physical aspects of turbulent flows are not fully understood. As a result, modelling turbulent flows is still challenging. In this section, first some im- portant physical characteristics of turbulent flows are explained and then the main modelling approaches are briefly explained. The governing equa- tions of a fluid flow consists of conservation of mass and momentum. The conservation of mass principle reads: ∂ρ ∂t + ∂ ∂xj (ρuj) = 0 (2.1) where ρ is density and uj for j = 1, 2, 3 are the three components of velocity. The conservation of momentum is a form of representing Newton’s second law: ∂ ∂t (ρui) + ∂ ∂xj (ρuiuj) = − ∂p ∂xi + ∂τij ∂xj (2.2) here the effect of body force is ignored. In this equation, p is pressure and τij is the stress tensor. For a Newtonian fluid, assuming Stokes’s hypothesis, this tensor is [48]: τij = µ ( ∂ui ∂xj + ∂uj ∂xi − 2 3 ∂ul ∂xl δij ) (2.3) 4 2.1. Turbulence where µ is the coefficient of viscosity. The last term in the above equa- tion vanishes for incompressible flows. Equation 2.2 is called Navier-Stokes equation. This equation is highly non-linear and the solution can be time de- pendent and chaotic which is called turbulence. Turbulence is the last great unresolved problem of classical physics as the great physicist Richard Feyn- man mentioned [21]. A turbulent flow is a multi-scale phenomena which can be considered as a flow of different eddies with different sizes [130]. These sizes range from large eddies, which are of the same order of the flow length scales, to the smallest scale, which is called the Kolmogorov length scale. The energy cascade is mainly from the large eddies to the smaller eddies. Basically, the large eddies are not stable and break down to smaller eddies through the vortex stretching phenomenon. This phenomenon continues un- til the eddy size is small enough that the turbulent kinetic energy of the eddy is dissipated to heat by molecular viscosity. According to Kolmogorov’s first similarity hypothesis for high Reynolds number turbulent flows, the small- est eddies have a universal statistical characteristics which depends only on viscosity and the rate of dissipation of turbulent kinetic energy. Hence, the Kolmogorov length and time scale can be defined as: η ≡ ( ν3 ε )1/4 (2.4a) tη ≡ (ν ε )1/2 (2.4b) where ν is the kinematic viscosity and ε is the rate of dissipation of turbulent kinetic energy. For the turbulent scales that are sufficiently bigger than the Kolmogorov scales but still smaller than the large scales, the turbulence statistics are only dependent on ε. This is based on the Kolmogorov’s second similarity hypothesis [130]. As a result, the turbulent velocity U and time scale τ of these scales can be related: U2/τ = ε (2.5) In Direct Numerical Simulation (DNS) of turbulent flows all the turbulent scales are resolved. The ratio of the large length scales to small scales can be approximated [130]: `/η ∼ Re3/4 (2.6) where ` is a representative length scale for large eddies which can be the integral length scale. Re is the turbulent Reynolds number. Based on this equation, direct numerical simulation of actual turbulent flows which have 5 2.1. Turbulence high Reynold number requires computational resources that exceeds the capacity of available supercomputers. While DNS is not a practical tool, it is a very useful research approach for studying the physics of a turbulent flow and development of reduced models for practical applications [106]. Apart from DNS, there are two main approaches to simulate turbulent flows: Reynolds-Averaged Navier-Stokes (RANS) and Large-Eddy Simula- tion (LES). In RANS, the averaged flow field is obtained and in LES the spatially filtered turbulent flow is captured. 2.1.1 Reynolds-Averaged Navier-Stokes In RANS, the governing equations of the averaged flow field are solved. These governing equations are obtained via Reynolds decomposition: ui = ui + u′i (2.7) here, the instantaneous velocity field u is decomposed into a mean velocity field ui and a fluctuating part u′i. Using this decomposition, one can obtain the mean field governing equations from the Navier-Stokes and continuity equations. It is easy to show that the non-linear terms result in correlation terms that are unknown. This is called the turbulent closure problem. For an incompressible flow, the averaged Navier-Stokes equation reads: ρ ( ∂ui ∂t + uj ∂ui ∂xj ) = − ∂p ∂xi + ∂ ∂xj ( µ ∂ui ∂xj − ρu′iu′j ) (2.8) where ρu′iu ′ j is the unclosed term and is called Reynolds stress tensor. µ ∂ui ∂xj is the viscous stress which is much smaller than the Reynolds stress and is usually ignored. For turbulent reacting flows the density is not constant; hence, another non-linear unclosed correlation arises which is ρ′u′i. In order to eliminate this term, a density-wighted averaging which is called Favre averaging is defined: ũi ≡ ρui ρ (2.9) Using the Favre version of the Reynold decomposition, i.e. ui = ũi + u′′i , one can derive the Favre-averaged Navier-Stokes equations for compressible flows: ρ ( ∂ũi ∂t + ũj ∂ũi ∂xj ) = − ∂p ∂xi + ∂ ∂xj (τ̃ij − ρũ′′i u′′j ) (2.10) where ρũ′′i u ′′ j is the compressible Reynolds stress tensor and τ̃ij is the vis- cous stress term and is usually ignored compared to the Reynolds stress. 6 2.1. Turbulence The challenge in all RANS turbulence models is to provide a modelling clo- sure for the Reynolds stress tensor based on the known mean variables. The most comprehensive RANS turbulence model is the Reynolds Stress Model (RSM). In this model, a transport equation is solved for the Reynolds stress term to directly obtain this tensor. This transport equation has several other unclosed terms which are not easy to close. Also, this model is com- putationally expensive and is not always numerically stable. The most popular RANS turbulence models are two-equation models which rely on the eddy viscosity concept. In this approach, the effect of tur- bulent eddies, which is the momentum transfer, is modelled with a turbulent viscosity. Hence, the Reynolds stress term can be written as: ρ̄ũ′′i u ′′ j = −µt ( ∂ũi ∂xj + ∂ũj ∂xi − 2 3 ∂ũl ∂xl δij ) + 2 3 ρ̄k̃δij (2.11) Several different models such as k− ε and k− ω are proposed to obtain the turbulent viscosity µt. In the k − ε turbulence model this quantity is given as: µt = Cµρ̄ k̃2 ε̃ (2.12) where Cµ is a model constant. Transport equations are solved for k̃ and ε̃. These transport equations in Standard k−ε turbulence model are as follows: ∂ ∂t (ρ̄k̃) + ∂ ∂xj (ρ̄ũj k̃) = ∂ ∂xj ( µt σk ∂k̃ ∂xj ) − ρ̄ũ′′i u′′j ∂ũi ∂xj − ρ̄ε̃ (2.13) ∂ ∂t (ρ̄ε̃) + ∂ ∂xj (ρ̄ũj ε̃) = ∂ ∂xj ( µt σε ∂ε̃ ∂xj ) − Cε1 ε̃ k̃ ρ̄ũ′′i u ′′ j ∂ũi ∂xj − C2 ρ̄ ε̃2 k̃ (2.14) where σk, σε, Cε1 and Cε2 are also model constants. Although RANS models are not as accurate as LES and DNS, they are not computationally expen- sive. Hence, RANS turbulence models - especially two-equations models like k − ε – are still the method of choice for practical purposes. 2.1.2 Large-Eddy Simulation Large-eddy simulation (LES) is a powerful turbulence model which provides a much more accurate and realistic representation of turbulence compared to RANS models. In LES, the spatially filtered governing equations of the flow field are solved and the effect of sub-filter scales (SFS) is modelled. In other words, the large scales of turbulence are captured and the small scales are 7 2.1. Turbulence modelled. Hence, modelling error in LES is substantially lower than RANS turbulence models in which all scales of turbulence are modelled. Defining G(~x) as a spatial low-pass filter function at location ~x, the filtered velocity field is: 〈ui(~x, t)〉 = ∫ V ui(~x′, t)G(~x− ~x′)d~x′ (2.15) where 〈ui〉 is the filtered velocity component in the i-th direction. An ex- ample for G(~x) is the Gaussian filter: G(~x) = ( 6 pi∆2 )1/2 exp ( −6~x · ~x ∆2 ) (2.16) where ∆ is the filter scale. In LES, the filtered Navier-Stokes equation is solved. Unlike the averaging operator in RANS, which commutes with the spatial derivatives in the Navier-Stokes equations, the filtering operator does not necessarily have this property: 〈∂ui ∂xj 〉 6= ∂〈ui〉 ∂xj (2.17) for a spatially-invariant filter it is easy to show that the spatial derivative commutes with the filtering operator. The problem is when the filter scale is not homogeneous. Ghosal and Moin [51] have shown that it is possible to define filter in a way that commutes with the derivatives in space at the second order accuracy. Higher-order commuting filters have also been proposed [155, 156], but there is no general agreement on the significance of these filters [50]. Implicit filtering is the most widely used approach in which the governing equations of the flow are not explicitly filtered according to equation 2.15. In this approach, the numerical scheme and the grid resolution act as a filtering operator. This approach is preferred in complex geometries with anisotropic grid where using an explicit filter is not easy and involves commutation errors. In such cases, the characteristic LES filter size is the local cell size which in three dimensions can be defined as: ∆ ≡ (∆x1∆x2∆x3)1/3 (2.18) where ∆xi is the grid spacing in i-th direction. Filtering the Navier-Stokes equations for an incompressible flow field gives the following equation: ρ ( ∂〈ui〉 ∂t + ∂〈ui〉〈uj〉 ∂xj ) = −∂〈p〉 ∂xi + ∂ ∂xj ( µ ∂〈ui〉 ∂xj − ρ〈uiuj〉+ ρ〈ui〉〈uj〉 ) (2.19) 8 2.1. Turbulence where σij ≡ ρ〈uiuj〉 − ρ〈ui〉〈uj〉 is the the sub-filter scale turbulent stress tensor. Several different LES sub-filter scale models have been proposed for this term [136]. Unlike RANS where the viscous stress tensor is small compared to the Reynolds stress, the viscous stress in LES µ∂〈ui〉∂xj is not negligible compared to the sub-filter scale turbulent stress tensor. Similar to RANS models for turbulent reacting flows, a density-weighted or Favre-filtered velocity field can also be defined: ũi ≡ 〈ρui〉〈ρ〉 (2.20) using Favre-averaging, the Navier-Stokes equation for compressible flows reads: 〈ρ〉 ( ∂ũi ∂t + ũj ∂ũi ∂xj ) = −∂〈p〉 ∂xi + ∂ ∂xj (τ̃ij − 〈ρuiuj〉+ 〈ρ〉ũiũj) (2.21) where σij ≡ 〈ρuiuj〉 − 〈ρ〉ũiũj is the sub-filter scale turbulent stress tensor and τ̃ij is the viscous stress tensor, which is defined as: τ̃ij ≡ −〈µ〉 ( ∂ũi ∂xj + ∂ũj ∂xi − 2 3 ∂ũl ∂xl δij ) (2.22) The eddy-viscosity hypothesis (also known as Boussinesq hypothesis) is also used in LES for modelling the sub-filter scale turbulent stress tensor. This approach for a compressible flow is as follows: σij − 13σllδij = −2〈ρ〉νt ( 〈Sij〉 − 13δij〈Sll〉 ) (2.23) where 〈Sij〉 is the filtered rate of strain tensor: 〈Sij〉 = 12 ( ∂〈ui〉 ∂xj + ∂〈uj〉 ∂xi ) (2.24) The first LES sub-filter scale model was proposed by Smagorinsky [141] based on the above eddy-viscosity approach and the following equation for the turbulent viscosity: νt = (Cs∆)2 √ 2〈Sij〉〈Sij〉 (2.25) where ∆ is the filter scale and Cs is the Smagorinsky coefficient. This coefficient is between 0.1 to 0.23 and can be calculated dynamically [49]. 9 2.2. Chemistry Another model for the sub-filter scale turbulent viscosity is the k-equation model. In this approach, one transport equation is solved for sub-grid scale turbulent kinetic energy k∆ = 1/2(ũiuj − ũiũj) [168]: ∂ ∂t (〈ρ〉k∆) + ∂ ∂xj (〈ρ〉ũjk∆) = ∂ ∂xj [ 〈ρ〉(ν + νt Prk ) ∂k∆ ∂xj ] (2.26) −σij ∂ũi ∂xj − 〈ρ〉Cεk 3/2 ∆ ∆ where Prk and Cε are modelling constants and νt is calculated using the sub-grid scale turbulent kinetic energy: νt = Cνk 1/2 ∆ ∆ (2.27) where Cν is a modelling constant. The k-equation sub-grid scale LES model depends on three modelling constants. Schmidt and Schumann proposed the values of Cε = 0.845, Cν = 0.0856 and Prk = 0.25 [138]. These parameters can also be computed using a dynamic model [82]. The Smagorinsky and k-equation models are categorized as functional approaches for LES. In functional approaches, the effect of sub-filter scales on the resolved scales is modelled. Another group of models are the struc- tural models which directly calculates the sub-filter scales [48, 136]. In these models, either the sub-filter scale turbulent stress tensor σij or the sub-filter scale velocity fluctuations u′i ≡ ui − 〈ui〉 are calculated. Several different models are available for this purpose: models based on the scale similar- ity hypothesis, models that use a formal series expansion and models that directly solve transport equations for σij [136]. Another modeling approach in LES is the Monotonically Implicit Large- Eddy Simulation (MILES) [11, 44]. In this approach, the sub-filter scale model is not explicitly implemented in the governing equations of turbulent flows. Hence, MILES relies on the numerical dissipation of the discretization scheme to represent the effect of sub-filter scales. 2.2 Chemistry The most simple representation of chemistry is an infinitely fast chemistry model [79]. In this model, the reactants are assumed to convert to products instantaneously. Hence; the local composition is determined by turbulent transport. A better model is to use chemical equilibrium which is based on the minimization of the Gibbs free energy of the local mixture. In this 10 2.2. Chemistry model, the equilibrium concentration of radicals such as CO and NO can be found. However, this model is also assuming a fast combustion process. For slow reactions finite rate chemistry happens where chemical kinetics and turbulent mixing are fully coupled. Finite rate chemistry was first described as a global one-step reaction from reactants to the final products. However, in reality many different species are present in combustion and they relate to each other through elementary reactions. GRI-Mech 3.0 [142] is a sample detailed chemistry mechanism for combustion of methane, which is a simple hydrocarbon. This mechanism contains 325 elementary reactions between 53 different species. Some of these species last for a very short time compared to others that have a longer time scale. Hence, chemistry is a multi-scale phenomena which is shown in Figure 2.1. This figure also shows that the chemistry times are over a larger range compared to the turbulent flow time scales and the resulting governing equations are stiff. chemical kinetics time scales 10−9 s 10 s10−7 s 10−5 s 10−3 s 10−1 s turbulence time scales Figure 2.1: A sample of range of time scales in a turbulent reacting flow. As the computing facilities have become more powerful these days, use of detailed chemistry in laminar flames is increasing. However, these are mainly limited to laminar flames especially in one-dimensional and two-dimensional simulations. Three-dimensional CFD modeling of turbulent reacting flows still require reduced chemistry models, especially in practical application. 2.2.1 Reduced Chemistry Chemistry was traditionally reduced using Quasi-Steady State Assumption (QSSA) for species and Partial Equilibrium Assumption (PEA) for reac- tions [26, 114]. The outcome of these methods is a skeletal mechanism or a 11 2.2. Chemistry reduced mechanism with several global steps that involve fewer species and fewer reactions compared to the detailed chemistry. However, these reduced mechanisms may still remain stiff [98]. Moreover, these methods are not efficient and their accuracy is not known prior to application [56]. More novel techniques for chemistry reduction were proposed such as Computational Singular Perturbations (CSP) [90], direct relation graph [92], Rate-Controlled Constrained Equilibrium (RCCE) [78] and Intrinsic Low- Dimensional Manifold (ILDM) [99]. CSP and ILDM are based on dynamical systems-based approaches for time scale analysis. In these methods, the main objective is to find fast and slow time scales and eliminate the slow time scales so that the system lays down on a lower-dimensional manifold in the composition space. This is shown schematically in Figure 2.2. Figure 2.2: Schematic representation of low-dimensional manifold technique for chemistry reduction. Low-Dimensional Manifolds are efficient methods for chemistry reduc- tion, especially for simulation purposes [99]. The initial work of Maas and Pope [99] was based on finding an intrinsic lower-dimensional manifold 12 2.2. Chemistry (ILDM) using a Jacobian analysis of the governing equations for scalars in the absence of transport. The lower-dimensional manifold can be obtained by several other methods [52, 68, 129, 133, 154]. Trajectory-Generated Lower-Dimensional Manifold (TGLDM) [68, 129] is one of these methods in which a prototype flame is used to build a lower-dimensional manifold. In premixed combustion, a one-dimensional laminar premixed flame can be used to generate trajectories. This method is called Flame Prolongation of ILDM (FPI) [52] or the Flame-Generated Manifold (FGM) [154]. The lower-dimensional manifold is tabulated and used in the multi- dimensional reacting flow simulations. Both the chemical reaction source terms and species mass fractions are tabulated as a function of a few con- trolling variables for which transport equations is usually solved in the real simulation. The storage and retrieval of these data can be performed using advanced numerical algorithms such as artificial neural networks [29] and in situ adaptive tabulation (ISAT) [127]. 2.2.2 Laminar Flames with Detailed Chemistry The interaction between chemistry and momentum transport categorizes the reacting flows into three different types: non-premixed, premixed and par- tially premixed flames. In non-premixed flames, fuel and oxidizer come from two different streams as shown in Figure 2.3a. Here, the chemical reaction happens where the two steams are mixing and the mixture is flammable. Hence, heat release is mainly controlled by hydrodynamic mixing. In pre- mixed flames, fuel and oxidizer are mixed before combustion and the chem- ical reaction happens at a self-sustained propagating zone as shown in Fig- ure 2.3b. Partially-premixed flames are in between these two cases. They can have a premixed propagating front, or diffusion characteristics. In some circumstance, both phenomena exists as a triple flame which was first ob- served experimentally by Phillips [118]. Premixed flames have a more complicated nature compared to non- premixed flames. In premixed flames, there is a strong interaction between molecular transport and chemical kinetics. This interaction can be first studied in laminar flames by neglecting the complexity due to turbulence. Also, the structure of laminar flames is important in turbulent reacting flows with a low turbulence intensity where the local structure of the flame is still laminar. All premixed flames have three different zones: preheat zone, inner layer and oxidation layer. In the preheat zone temperature increases due to heat conduction from the inner layer. When the temperature is high enough the 13 2.2. Chemistry (a) Non-premixed (b) Premixed Figure 2.3: Schematic of two different types of reacting flows: non-premixed and premixed flames. chemical reactions happen at the inner layer. In this layer intermediate radicals are present and oxidize the fuel to form CO and H2 species. The main heat of combustion is released in this layer. The final products are then formed in the post flame oxidation layer. These three different zones are shown in Figure 2.4 for a steady one-dimensional laminar premixed flame. The governing equations for a steady one-dimensional laminar premixed flame, which is the simplest laminar premixed flame, are as follows [123]: ρu = ρuSL = constant (2.28a) ρu dYk dx = −dJk dx + ω̇k (2.28b) ρcpu dT dx = d dx (λ dT dx )− K∑ k=1 cpkJk dT dx − K∑ k=1 ω̇khk (2.28c) where the subscript k refers to the k-th species, ρ is the density of the mixture, ρu is the density of the unburnt mixture, SL is the laminar flame speed, u is velocity, Yk is the mass fraction, ω̇k is the chemical production rate, T is temperature, cp is the heat capacity of the mixture at constant pressure, λ is the thermal conductivity, hk is the total enthalpy and Jk is the molecular diffusion flux of the k-th species. In case of reacting flows with air as an oxidizer, the mixture is highly dilute; hence, Jk can be approximated as follows [89]: Jk = −ρDk dYk dx (2.29) 14 2.2. Chemistry 0 1 2 3 4 5 6 7 8 9 x 10−3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x(m) T/T equilibrium YCH 4 * 10 YH * 10,000 Inner Layer Oxidation LayerPre−heat Layer Figure 2.4: Chemical structure of one-dimensional unstrained laminar pre- mixed flame for a lean methane-air mixture. These equations can be solved with detailed chemistry using software tools such as Cantera [55]. The thermal thickness of a one-dimensional premixed flame can be defined as: `F ≡ Tb − TudT dx |max (2.30) This quantity can be represented graphically as shown in Figure 2.5. The thickness of the inner layer of a premixed flame `δ is almost one order of magnitude less than the thermal thickness. Another important parameter in laminar premixed flames is the lam- inar flame speed. In case of a one-dimensional flame, this speed SL de- pends on molecular diffusion and chemical kinetic properties of the mixture. Mauss and Peters demonstrated the importance of using detailed chemical kinetic mechanism in laminar flame speed calculations [103]. They showed that ignoring C2 species and related elementary reactions results in under- 15 2.2. Chemistry x T Tb Tu ℓF Figure 2.5: Schematic representation of laminar flame thickness definition. estimation of the laminar flame speed in rich mixtures. In one-dimensional laminar premixed flames, the cold mixture velocity upstream of the flame front relative to the flame is uniform and constant. In two and three-dimensional flames, the upstream velocity field is not uniform with respect to the flame; hence, the flame surface is hydrodynamically stretched. The flame stretch can be defined as the normalized rate of change of the flame area [123]: κ ≡ 1 A dA dt (2.31) this stretch term can be decomposed into a tangential strain aT and a cur- vature term: κ = aT + Sd∇ · ~n (2.32) where ~n is the flame front normal vector and∇·~n is the local flame curvature. In this equation, Sd is the flame displacement speed which differs from the one-dimensional unstrained laminar flame speed SL. Theoretical studies suggest that the Sd has a linear relationship with stretch κ for small values of stretch: Sd = SL − Laκ (2.33) where La is the Markstein length. This length indicates the effect of stretch on the flame speed which is usually determined by using experimental data [88]. Simple examples of stretched laminar premixed flames are steady stagnation point flame and unsteady spherical flame. Figure 2.6 shows the stagnation 16 2.3. Turbulence-Chemistry Interactions point premixed flame where reactants and hot products are set against each each other. Changing the inlet speeds of these two streams results in stretch- ing the flame by changing the velocity in the flame tangential plane. In an unsteady spherical flame, the flame front is a sphere with increasing radius. In this flame, the reactants are outside of the sphere and the combustion products are inside. The stretch in this flame is due to flame front curvature. 0 0 x y Products u v Fuel + Air Stagnation Plane Flame Figure 2.6: Schematic of a stagnation point laminar premixed flame. 2.3 Turbulence-Chemistry Interactions The interaction between turbulence and chemistry is complex phenomenon due to the multi-scale and non-linear nature of both turbulence and chem- istry. In this section, first the physics of this interaction is explained and then the modelling approaches are briefly discussed. 2.3.1 Regimes in Premixed Turbulent Combustion Damköhler was the first to characterize the interaction between different scales of turbulence and chemistry in terms of turbulent combustion regimes [33]. He categorized this interaction into two different regimes: large-scale and 17 2.3. Turbulence-Chemistry Interactions small-scale turbulence. In large-scale turbulence regime, the interaction be- tween turbulence and chemistry is limited to moving the locally laminar flames while small-scale turbulence can locally modify the flame so that the turbulent diffusion is dominant compared to molecular diffusion. The first regime is also called the flamelet regime, which is especially important in tur- bulent premixed combustion. For low turbulence intensity, the flame front is perturbed and wrinkled by turbulent eddies. As the turbulence intensity in- creases, the premixed flame becomes corrugated as shown in Figure 2.7. For higher turbulence intensities where the small scales are of the same size of the reaction zone, the validity of Damköhler’s idea has been questioned [39, 117]. Peters has argued that the small scale turbulent fluctuations cause heat loss and temperature fluctuations which lead to extinction of laminar flamelets [117]. Hence, the only effect of small-scale turbulence on the local structure of laminar flamelets is extinction. His argument is based on the experimen- tal results of Mansour et al. [100]. On the other hand, other experimental and Direct Numerical Simulation results [4, 40, 58] questioned the validity of the flamelet assumption for high turbulence intensities. This is an open question in the turbulent combustion field. Figure 2.7: Interaction between turbulence and premixed flame. A more comprehensive phenomenological explanation of the interaction between different scales of turbulence and chemistry can be expressed in terms of different regime diagrams. The turbulent premixed regime dia- 18 2.3. Turbulence-Chemistry Interactions gram was proposed by several researchers including Bray [16], Borghi [10], Williams [165], Peters [115], Abdel-Gayed and Bradley [2] and Pitsch and Duchamp de Lageneste [122]. The most popular regime diagram comes from the work of Borghi [10] which is called the Borghi diagram. In this regime diagram, turbulence in- tensity over the laminar flame speed u′/SL is plotted against the ratio of the integral length scale over the laminar flame thickness `/`F on a logarithmic scale. This diagram is shown in Figure 2.8. 10−1 100 101 102 103 104 10−1 100 101 102 103 ℓ/ℓF u ′ / S L Laminar Flames Wrinkled flamelets Corrugated flamelets Ka = 1 Well−stirred Reactor Da = 1 Re = 1 Distributed Reaction Zone Figure 2.8: Different regimes in turbulent premixed combustion [10]. Assuming equal diffusivity D for all species, the laminar Schmidt number can be defined as Sc ≡ ν/D, where ν is the kinematic viscosity of the mixture. The laminar flame thickness is of order of `F ∼ D/SL where SL is the laminar flame speed. Also, the flame time scale can be taken as tF = D/S2L. A turbulent Reynolds number is defined based on the turbulence intensity u′ and turbulent integral length scale `: Re ≡ u ′` ν (2.34) 19 2.3. Turbulence-Chemistry Interactions Assuming unity Schmidt number, one can re-write the Reynolds number in the following way: Re = u′` D = ( u′ SL )( ` `F ) (2.35) For turbulent Reynolds number of less than one, the flame is laminar as shown in Figure 2.8. One of the important non-dimensional parameters in turbulent combustion is the Damköhler number which is defined as the large scale turbulent time scale over the large scale chemistry time scale. Here, the large scale turbulent time scale can be defined as tT = `/u′ and the flame time scale corresponds to a large time scale in laminar premixed flame. Therefore, the Danköhler number is: Da ≡ tT tF = ( u′ SL )−1( ` `F ) (2.36) For a large Damköhler number, chemical effects are not disturbed by the turbulent fluctuations and the flame is locally laminar. As the turbulent intensity increases the flamelets are wrinkled and then corrugated until the smallest turbulent scale, the Kolmogorov scale, is of the order of the flame time scale. This regime can be further characterized by another important non-dimensional number called Karlovitz number, which is defined as the ratio of the chemistry time scale to the Kolmogorov time scale: Ka ≡ tF tη = `2F η2 (2.37) one can easily show the following relationship between Damköhler, Karlovitz and the Reynolds number: Re = Da2Ka2 (2.38) Using this equation and equations 2.35 and 2.36, one can write the following equation for the Karlovitz number in terms of u′/SL and `/`F : Ka = ( u′ SL )3/2( ` `F )−1/2 (2.39) Based on the definition of the Damköhler and Karlovitz numbers, the four different regimes shown in the Borghi diagram (Figure 2.8) can be explained. An important regime in this diagram is the distributed reaction zone regime for Ka > 1 and Da < 1. Peters [117] has shown that distributed burning 20 2.3. Turbulence-Chemistry Interactions does not happen in this regime. The reason is that the small scale Kol- mogorov eddies, which are of the order of the thermal thickness of the flame `F , can penetrate into the preheat zone of the locally laminar premixed flame, but they vanish and cannot perturb the inner layer to destroy the structure of the underlying laminar flame. He defined a second Karlovitz number based on the inner layer thickness: Kaδ ≡ ` 2 δ η2 (2.40) He modified the Borghi diagram and proposed a different one as shown in Figure 2.9. The main difference between this diagram and the original Borghi’s diagram is thin reaction zone regime for Ka > 1 and Kaδ < 1 where the inner zone of the laminar flamelets are not disturbed and the small eddies can only penetrate the preheat layer. For Kaδ > 1, the dilemma of distributed burning or broken reaction exists which was explained before. 10−1 100 101 102 103 104 10−1 100 101 102 103 ℓ/ℓF u ′ / S L Laminar Flames Wrinkled flamelets Corrugated flamelets Re = 1 Ka = 1 Distributed burning/ Broken Reaction Zone Kaδ = 1 Da = 1 Thin Reaction Zone Figure 2.9: Modified turbulent premixed combustion regime diagram [116]. 21 2.3. Turbulence-Chemistry Interactions 2.3.2 Turbulent Combustion Modelling in RANS RANS turbulent combustion models provide closure for the mean chemical reaction source term in the transport equations for species mass fractions. The transport equation for instantaneous species mass fractions reads: ∂ ∂t (ρYk) + ∂ ∂xj (ρujYk) = ∂ ∂xj ( Dk ∂Yk ∂xj ) + ω̇k (2.41) where ω̇k is the chemical reaction source term which is a highly non-linear function of the scalar field. Therefore, mean reaction source terms cannot be correctly computed by evaluating them from averaged scalars: ω̇k(T, Yk, ρ) 6= ω̇k(T̄ , Y k, ρ̄) (2.42) Although researchers have been investing significant effort in this field, there is still no reliable method for predicting mean heat release rates in premixed turbulent combustion [18]. In order to characterize the extent to which combustion has been completed in turbulent premixed flames, most models are implemented using either a progress variable approach or a G-equation approach [123]. The first approach is based on a progress variable which is 0 in reactants and 1 in products. An example of this variable can be a product species mass fraction normalized by the equilibrium mass fraction. Hence, the transport equation for instantaneous progress variable c is similar to equation 2.41: ∂ ∂t (ρc) + ∂ ∂xj (ρujc) = ∂ ∂xj ( Dc ∂c ∂xj ) + ω̇c (2.43) the time-averaged version of this equation is solved in reacting flow simu- lations which involves the unclosed source term ω̇c. Several turbulent com- bustion models were proposed for closure of this term, and are discussed in the following sections. Bray Moss Libby (BML) model The Bray Moss Libby (BML) [19] model is a simple algebraic turbulent combustion model. The basic assumption in the BML combustion model is that the flame is infinitely thin. Under the BML assumptions, the progress variable is either 0 or 1; thus, the probability density function of c would be a combination of two Dirac delta functions: P (c∗) = αδ(c∗) + βδ(c∗). By using this PDF, the following relation for mean reaction is obtained: ω̇c = ρ̄̃c 2cm − 1 (2.44) 22 2.3. Turbulence-Chemistry Interactions where cm is a number of order one and ̃c is the scalar dissipation rate. This is a key parameter in turbulent combustion modelling and are discussed in the next chapter. Equation 2.44 is similar to the one obtained using the eddy break-up model proposed by Spalding [144], albeit more meaningful. Flame Surface Density (FSD) model Another widely-used model for ω̇c is the Flame Surface Density (FSD) model. Here, the mean reaction rates is assumed to be proportional to the flame surface area per unit volume multiplied by the laminar flame speed: ω̇c = ρuSLI0Σ̄ (2.45) where ρu is the fresh gas density, I0 is the stretch factor and Σ̄ is the average flame surface density inside a control volume. Several algebraic models and transport equations have been proposed to calculate this parameter, as summarized by Poinsot and Veynante [123]. Presumed Conditional Moment (PCM) model Presumed Conditional Moment (PCM) model [15, 52, 154] is another turbu- lent combustion model in which the local conditional structure of the flame is presumed to be that of a laminar flame. This model is then combined with a presumed PDF approach to get the mean reaction rates from the presumed conditional ones: ω̇c = ∫ 1 0 ω̇(c∗)P (c∗)dc∗ (2.46) where ω̇(c∗) comes from the solution of one-dimensional laminar canonical flame. In this model, the detailed chemistry is easily coupled with the multi- dimensional turbulent reacting flow simulation. G-equation model The G-equation approach [80] is based on calculating a G field by using a level-set approach. This parameter ranges from −∞ to +∞ and the flame is located at G = 0. The transport equations for G reads: ∂G ∂t + ũi ∂G ∂xi = ST |∇G| (2.47) where ST is the turbulent flame speed for which several empirical correlations have been proposed [39]. However, this is not a well-defined parameter 23 2.3. Turbulence-Chemistry Interactions and there are considerable uncertainties in the experimentally measured turbulent flame speed [123]. Transported PDF model All the above models are valid in the flamelet regime. However, the trans- ported PDF model and Conditional Moment Closure (CMC) are more gen- eral. In the transported PDF model, a transport equation for the joint probability density function of velocity and composition is solved. In a more simplified way, only the composition PDF transport equation is derived and solved. In both cases, an integral-differential equation with more than 4 independent variables is obtained which makes the regular grid-based CFD models impracticable for solving this equation. Statistical methods such as Lagrangian Monte Carlo particle method are the numerical methods of choice. The transported PDF model provides a perfect closure for chemical source terms. However, the diffusion terms are unclosed in the transport equation for the joint PDF [125] which is the main modelling challenge in this method. This challenge is of vital importance in a turbulent premixed flame, where a strict coupling between turbulent mixing and chemistry takes place. The transported PDF model [125] is a more computationally expen- sive approach compared to other turbulent combustion models. A recent comprehensive review of the progress of this method can be found in [65]. This method was first applied to premixed flames by Anand and Pope [3] where a transport equation for the joint probability density function of ve- locity and and a progress variable is solved. An idealized one-dimensional turbulent premixed flame was solved and the final turbulent flame speed re- sults were compared against the experimental data. This method has been used in several simulations of turbulent premixed flames in laboratory-scale burners [22, 97, 108, 146, 160] and in a practical combustor [73]. Conditional Moment Closure (CMC) Conditional moment closure was first formulated for premixed flames by Bilger [7]. CMC gives a less constrictive closure for the chemical source terms compared to flamelet models, in that it provides a statistical closure that does not rely on the flame being “thin”. Moreover, it requires less computational time than the transported PDF model. In first-order CMC, it is assumed that the conditional mean of reaction rates can be evaluated by using the conditionally averaged scalars: ω̇c(T, Yk, ρ)|c∗ ≈ ω̇c(T |c∗, Yk|c∗, ρ|c∗) (2.48) 24 2.3. Turbulence-Chemistry Interactions and then the unconditional mean reaction rate can be easily computed by knowing the PDF of progress variable: ω̇c = ∫ 1 0 ω̇c|c∗ P (c∗)dc∗ (2.49) This model is widely-used and well-developed in non-premixed turbulent combustion [83]. Swaminathan & Bilger [147] have shown – in an a priori sense using DNS data – that equation 2.48 can be used to predict mean chemical source terms in premixed flames. However, in order to use equa- tion 2.48, one needs to first calculate both the conditional scalar field and the PDF of reaction progress variable. A presumed PDF is used for the re- action progress variable and the conditional scalar field is computed by solv- ing transport equations. These transport equations involve several unclosed terms which is not easy to provide accurate closure in premixed combustion. Only two attempts to date have been made to use CMC in actual tur- bulent reacting flow simulations. In the first of those, Smith [143] could not obtain a converged solution. In the second, Martin et al. [102] used CMC to simulate a dump combustor, however their results are limited to a qualitative comparison with experimental data. The application of CMC to turbulent premixed combustion is still challenging and remains an area of research. Multiple Mapping Conditioning (MMC) Multiple Mapping Conditioning (MMC) is another model for turbulence- chemistry interactions. This model is a combination of CMC and the trans- ported PDF model proposed by Klimenko and Pope [84]. MMC does not have the limitations of the transported PDF model and provides a perfect closure for chemical reaction source terms. This model has been applied to simulate non-premixed flames [30, 161]; however, it has not yet been extended to premixed combustion. Conditional Source-term Estimation (CSE) Conditional Source-term Estimation (CSE) model is very similar to CMC in terms of the basic closure remedy for chemical reaction source term (equa- tion 2.48). However, the conditionally averaged scalar field is calculated by solving an integral equation using an inverse method. This method was first proposed by Bushe and Steiner [20] for LES of non-premixed flames. It is also used in RANS simulations of turbulent non-premixed flames [68, 105]. 25 2.3. Turbulence-Chemistry Interactions Jin et al. [76] showed that CSE can also be used in premixed flames by doing an a priori test with DNS data. Using this method in an actual CFD code and comparing the results with experiment is the subject of this research. 2.3.3 Turbulent Combustion Modelling in LES In LES of turbulent reacting flows – especially at high Reynolds and Damköhler numbers – chemical reactions happen at unresolved scales. As a result, large eddy simulation of turbulent reacting flows is still a modelling challenge. Nevertheless, LES is still preferred to RANS, particularly in the presence of large scale motions [128]. A variety of combustion models have been proposed for LES of turbulent premixed flames. Apart from the artificially thickened flame model [31] and linear eddy model [25], other LES combus- tion models have generally been adopted from similar RANS models [121]. Presumed Conditional Moments (PCM) [38, 45], Flame Surface Density (FSD) model [9, 63] and the G-equation model [69, 71, 81, 85] have been successfully modified and used in LES. The transported PDF model is also extended from RANS to LES by introducing the idea of Filtered Density Function (FDF) [53, 126]. This method is successfully used in LES of turbulent non-premixed flames in a laboratory-scale burner [140]. The transported PDF model has also been used in LES of premixed combustion [167]. Conditional Source-term Estimation (CSE) model was first proposed as a CMC turbulent combustion model for LES [20]. This model is used in LES of non-premixed laboratory-scale burner [145, 163]. The original CMC model in which transport equations are solved for the conditionally-averaged scalar field is also used as a sub-filter scale turbulent combustion model for LES [47, 109, 110, 151]. All these attempts are limited to turbulent non- premixed flames. CMC has not been used in turbulent premixed flames as yet. 26 Chapter 3 PCM-FPI Model in RANS Simulation of Turbulent Premixed Flames Presumed Conditional Moments (PCM) approach is a flamelet turbulent combustion model which assumes that the local structure of the turbulent flame as a function of a progress variable is that of a laminar flame. This local structure can be the Flame Prolongation of ILDM (FPI) chemistry model which was described in section 2.2.1. Using a presumed probability density function (PDF) approach, one can couple the flamelet chemistry with the turbulent flame. This method is a common approach in non-premixed flamelet mod- elling [113]. There, the definition of a conserved scalar and a coordinate transformation from physical space to conserved scalar space results in flamelet equations. The shape of the PDF of this conserved scalar – mix- ture fraction – is often approximated with reasonable accuracy by a β- function [96]. This distribution is a function of the mean and variance of the mixture fraction for which transport equations can be solved. In turbulent premixed flames, the governing equations cannot be trans- formed to a conserved scalar space. Instead, the structure of a laminar un- strained premixed flame can be used as a flamelet as proposed by Bradley et al. [14, 15], Van Oijen et al. [72, 154] and Gicquel et al. [52]. This method can be used with Reynolds-Averaged Navier-Stokes (RANS) turbu- lence models [42, 134] and with Large Eddy Simulation (LES) turbulence model and Direct Numerical Simulation (DNS) [38, 45]. Several researchers including Bradley et al. [13], De Goey et al. [34] and Kolla and Swami- nathan [86] have suggested models to account for the effect of fluid dynamic stretch on flamelets. In this chapter, FPI chemistry along with two different functional forms for the presumed PDF are applied to simulate the Bunsen burner experi- ments of Chen et al. [28]. This experiment has been simulated by using both 27 3.1. Theory a RANS turbulence model [5, 67, 97, 131, 139] and LES [35, 85]. In the fol- lowing sections, first the theory and the governing equations are discussed, then a brief review of the CFD approach are described and the results of the simulation are presented. 3.1 Theory 3.1.1 Presumed PDF The FPI representation of chemistry is obtained by solving transport equa- tions for conservation of momentum, energy and species mass fraction in a steady one-dimensional laminar unstrained premixed flame. The FPI chem- istry can then be tabulated for different equivalence ratios as a function of a reaction progress variable. The effect of radiation heat transfer can also be considered by adding a source term to the energy equation [41] and an additional dimension to the table. For an adiabatic premixed flame at a specific mixture fraction, the chemical reaction rates are only functions of the progress variable. Knowing the PDF of progress variable, the averaged reaction source term in the transport equation of the k-th species can be calculated by a simple integration: ω̇k = ∫ 1 0 ω̇k|c∗P (c∗)dc∗ = ∫ 1 0 ω̇FPIk (c ∗)P (c∗)dc∗ (3.1) where ω̇FPIk (c ∗) is the chemical source-term from the one-dimensional FPI calculation expressed as a function of the conditioning variable c∗, related to the progress variable c. Several different functional forms for P (c∗) have been used in combina- tion with FPI chemistry. They all are formed using mean and variance of progress variable. For high variances, there are either reactants or products at every point of the reacting field; in this extreme, the PDF becomes two delta functions [134]. This regime corresponds to a large Damköhler num- ber, defined as the ratio of a turbulent mixing time scale to a chemical time scale. For extremely low Damköhler numbers, the reaction zone is thick and there is a very low value of variance. In such cases, the PDF takes a narrow Gaussian shape. Since the β-PDF can match these two extremes, this form has found some favor in premixed flame calculations [42]. Bray et al. [18] proposed a laminar flamelet PDF defined as: P (c∗;x, t) = Aδ(c∗) +Bf(c∗) + Cδ(1− c∗) (3.2) 28 3.1. Theory where δ(c∗) is the Dirac delta function and f(c∗) is calculated using the solution of an unstrained laminar flame: f(c∗) = 1 |∇c∗| . (3.3) The constants A, B and C are calculated from the following equations:∫ 1 0 P (c∗)dc∗ = A+B ∫ 1− f(c∗)dc∗ + C = 1 (3.4a) ∫ 1 0 c∗P (c∗)dc∗ = B ∫ 1− c∗f(c∗)dc∗ + C = c̄ (3.4b)∫ 1 0 c∗ 2 P (c∗)dc∗ = B ∫ 1− c∗ 2 f(c∗)dc∗ + C = c′2 + c̄2 (3.4c) where is a very small number. Alternatively, a Favre PDF can be deter- mined by c̃ and c̃′′2. Favre PDF is defined as: P̃ (c∗;x, t) = ρ|c∗ ρ̄ P (c∗;x, t) (3.5) where ρ̄ is the mean density and ρ|c∗ is the conditional density which is calculated from the solution of an unstrained laminar premixed flame. In either case, the constants A, B and C should all be positive; hence the following constraints exist for mean and variance: c̄(c̄− 1) + c′2 I2 − I1 ≥ 0 (3.6a) c̄I2 − (c̄2 + c′2)I1 I2 − I1 ≥ 0 (3.6b) I2 − I1 − c′2I0 − c̄(c̄− 1)I0 − c̄I2 + (c′2 + c̄2)I1 I2 − I1 ≥ 0 (3.6c) where, I0 = ∫ 1− f(c ∗)dc∗, I1 = ∫ 1− c ∗f(c∗)dc∗ and I2 = ∫ 1− c ∗2f(c∗)dc∗. As a result, this PDF can be generated only for a limited range of mean and variance. As an example, this range is shown in figure 3.1 in a shaded area for stoichiometric methane-air chemistry. This figure clearly shows that this PDF is only applicable to high Damköhler number flames as stated by Bray et al. [18]. The experimental data of Chen et al. [28] are also depicted in figure 3.1. These data show that for this flame mean and variance 29 3.1. Theory of reaction progress variable does not fall into the applicable range of the above PDF. This range can also be depicted in a Favre-averaged graph (c̃′′2 vs. c̃); however, since the mean experimental data come with Reynolds- averaged quantities, this map is shown using the Reynolds-averaged mean and variance. Domingo et al. [38] proposed a combination of this PDF with a β-PDF for low variances in LES, the so called FSD-PDF sub-grid scale closure model. Figure 3.1 shows that using this model in RANS would result in using the β-PDF in the whole flame brush. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c̄ c′ 2 Figure 3.1: Realizable values of mean and variance of reaction progress variable. The shaded area shows the applicable range of the PDF of equation 3.2 for a stoichiometric methane-air flame. Symbols denote the experimental data for Flame F3 of the Bunsen burner of Chen et al. [28] at axial positions x/D = 8.5 (+), x/D = 6.5 (♦), x/D = 4.5 () and x/D = 2.5 (◦). In order to overcome this limitation, a modified version of the laminar flamelet PDF has been proposed by Jin et al. [76] which is called the modified laminar flamelet PDF (MLF-PDF). In this modification, four possible cases are considered: 1. The laminar flamelet PDF (equation(3.2)) 2. The laminar flamelet PDF with only one delta function at c∗ = 0 and 30 3.1. Theory then the same laminar inner distribution until c∗ = c2 < 1 3. The laminar flamelet PDF starting from c∗ = c1 > 0 and ending with a delta function at c∗ = 1 4. The inner distribution function form c1 > 0 to c2 < 1 Figure 3.2: Four possible cases of the modified laminar flamelet PDF These cases are depicted schematically in figure 3.2. In each case, the values of only three of A, B, C, c1 and c2 are needed, and these values are, again, chosen to match mean and variance of progress variable and to meet the constraint that ∫ 1 0 P (c ∗)dc∗ = 1. For a given mixture of fuel and oxidizer, the MLF-PDF can be calculated a priori for a range of mean and variance – from 0 to 1 for the mean and from 0 to 0.25 for the variance. The procedure for calculating the PDF is as follows: 31 3.1. Theory 1. All possible numerical distributions for P (c∗;x, t) are calculated by changing the three constants in each of the four possible cases. For example, in the second case of the PDF, A and c2 vary between 0 and 1. Then, B is set so that ∫ 1 0 P (c ∗)dc∗ = A+B ∫ c2 0 f(c ∗)dc∗ = 1. 2. The mean density is calculated using ρ̄ = ∫ 1 0 ρ FPI(c∗)P (c∗)dc∗. Thus, the Favre-PDF, c̃ and c̃′′2 are computed. 3. The calculated numerical distributions for P (c∗;x, t) and P̃ (c∗;x, t) are tabulated as functions of c̃ and c̃′′2. In order to optimize this look-up table, the domain is meshed using Delaunay triangulation. Now, given the mean and variance of reaction progress variable, the probability density function of the reaction progress variable can be interpolated using the calculated PDFs. Two samples of the c̃ - c̃′′2 domain are shown in figures 3.3 and 3.4 for methane-air mixtures at two different equivalence ratios. These figures clearly show that for a very low value of variance, only case 4 occurs. For high variance, depending on the mean reaction progress variable, a PDF generated using either case 1, 2 or 3 can occur. If the mean value of progress variable is near 0, a PDF of case 2 is more likely, and if the mean is near 1, a PDF of case 3 is more likely. The MLF-PDF was tested in the a priori sense against DNS data of a freely propagating planar premixed flame [76] and found to be dramati- cally superior to the β-PDF in predicting the shape of the PDF because it incorporates the effects of the chemistry on that shape. In this work, the ability of this PDF to provide the statistics of the reaction progress variable in highly turbulent Bunsen flames is compared to the widely used β-PDF. The purpose of this work is primarily to show that the new PDF can be used to obtain meaningful predictions of a turbulent premixed flame that are at least competitive with those obtained using existing methods. 3.1.2 Transport Equations for c̃ and c̃′′2 For a range of mean and variance of reaction progress variable and a partic- ular chemistry, the MLF-PDF can be formed and the mean reaction rates integrated and stored a priori as an FPI-PDF table. In order to use this table in a reacting flow CFD code, transport equations for mean reaction progress variable and its variance must be solved. The instantaneous trans- 32 3.1. Theory 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c̃ c̃” 2 Case 4 Case 2 Case 3 Case 1 Figure 3.3: Range of the applicability of each of the four possible cases of the modified laminar flamelet PDF for a lean mixture of methane and air at equivalence ratio of 0.6. port equation for reaction progress variable is: ∂ρc ∂t + ∂ ∂xi (ρuic) = ∂ ∂xi ( ρD ∂c ∂xi ) + ω̇c (3.7) where D is the molecular diffusivity, ui is velocity and ω̇c is the chemical source-term of the progress variable. If we average the above equation, the transport equation for the Favre-averaged field c̃ ≡ ρc/ρ is: ∂ρ̄c̃ ∂t + ∂ ∂xi (ρ̄ũic̃) = − ∂ ∂xi ( ρ̄ũ′′i c′′ ) + ∂ ∂xi ( ρ̄Dc ∂c̃ ∂xi ) + ω̇c (3.8) where c′′ ≡ c− c̃. The above equation has two unclosed terms: ω̇c and ũ′′i c′′. The latter of these is the turbulent scalar flux, which is usually modelled by using a gradient assumption for a passive scalar: ũ′′i c′′ = − νt Sct ∂c̃ ∂xi (3.9) 33 3.1. Theory 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 c̃ c̃” 2 Case 4 Case 2 Case 3 Case 1 Figure 3.4: Range of the applicability of each of the four possible cases of the modified laminar flamelet PDF for a stoichiometric methane-air mixture. where, νt is the turbulent dynamic viscosity and Sct is the turbulent Schmidt number. For a reactive scalar, the turbulent scalar flux and the mean progress variable vector are not necessarily in the same direction, i.e., ũ′′i c′′ ∂ec ∂xi can be either positive or negative. When it is positive, the so-called counter- gradient diffusion phenomenon occurs. This phenomenon is due to the inter- mittency of the flame front and takes place if the heat release effect is more dominant than turbulence. Increasing the turbulence makes the flame front wrinkled and corrugated in which case this correlation becomes negative and gradient diffusion occurs. The Bray number is a non-dimensional number which can determine whether gradient or counter-gradient diffusion happens. This number is less than one for a gradient diffusion and bigger than one for counter-gradient diffusion. For a constant pressure premixed flame, this number is [158]: B = τbSL 2αu′ (3.10) where, SL is laminar flame speed, u′ is the turbulence intensity and α is an empirical constant of order unity. An appropriate value for α can be 34 3.1. Theory Figure 3.5: Gradient vs. counter-gradient diffusion in turbulent premixed flames estimated as a function of the ratio of turbulent length scale and laminar flame thickness [158]. τb is the heat release index defined as τb ≡ (Tb − Tu)/Tu where, Tb and Tu are the temperature of burned and unburned gases respectively. When the Bray number is very small, a gradient assumption can be used to model the turbulent scalar flux. When the Bray number fluctuates around unity, both gradient and counter-gradient diffusion occurs inside the flame brush. For such flames, using a gradient model is not appropriate. Several algebraic closures have been proposed for turbulent scalar flux in both RANS [158, 171] and LES [152] that can switch sign between gradient and counter-gradient regimes. A better approach for modelling this unclosed term is second moment models in which a transport equation is solved for turbulent scalar flux [37]. In our work, a gradient assumption is used be- cause the Bray number is smaller than one everywhere inside the turbulent premixed flames of Chen et al. The transport equation for c̃′′2, with some simplifications, is [123]: ∂ρ̄c̃′′2 ∂t + ∂ ∂xi ( ρ̄ũic̃′′2 ) = ∂ ∂xi ( ρ̄ νT Sc1 ∂c̃′′2 ∂xi ) − 2ρ̄ νT Sc2 ∂c̃ ∂xi ∂c̃ ∂xi (3.11) −2ρDc∂c ′′ ∂xi ∂c′′ ∂xi + 2c′′ω̇c The last two terms in the above equation are unclosed. The latter of these 35 3.1. Theory is the correlation of reaction rate and fluctuations in progress variable. This term can be easily computed knowing the PDF: c′′ω̇c = (c− c̃) ω̇c = ∫ 1 0 (c∗ − c̃) ω̇c(c∗)P (c∗)dc∗ (3.12) The other unclosed term, ρ̄̃c ≡ ρDc ∂c′′∂xi ∂c ′′ ∂xi , is approximately equal to the Favre-averaged scalar dissipation rate which is an indication of local mixing. For a passive scalar, a good approximation for the scalar dissipation rate can be achieved by a linear relaxation: ρDc∂c ′′ ∂xi ∂c′′ ∂xi = CDρ̄ ε̃ k̃ c̃′′2 (3.13) where k̃ is turbulent kinetic energy, ε̃ is rate of dissipation of the turbulent ki- netic energy and CD is the ratio of turbulence time scale (k̃/ε̃) to scalar time scale (c̃′′2/̃c). For a reactive scalar, on the other hand, the effect of turbu- lence on micro-mixing of the scalar field is more complicated. Hence, the best approach for modelling scalar dissipation rate is deriving and solving a bal- ance equation. Several attempts can be found in the literature for deriving a transport equation for scalar dissipation rate [23, 24, 101, 107, 148, 149]. For a large Damköhler number flame, production and dissipation terms of such transport equations are dominant. Hence, an algebraic model can be derived for scalar dissipation rate. The first algebraic model for scalar dissipation rate in premixed flames was proposed by Mantel and Borghi [101]: ̃c = ( 1 + 2CεcSL 3 √ k̃ )( CD ε̃ k̃ ) c̃′′2, Cεc = 0.1, CD = 0.21 (3.14) where SL is the laminar flame speed. This model was modified by Swami- nathan and Bray [148]: ̃c = ( 1 + 2CεcSL 3 √ k̃ )( CDc SL `F + CD ε̃ k̃ ) c̃′′2, CDc = 0.24. (3.15) where `F is the laminar flame thickness. Kolla et al. have also proposed an algebraic model for scalar dissipation rate in premixed flames [87]: ̃c = 1 β′ ( 2K∗c SL `F + [C3 − τbC4DaL] ε̃ k̃ ) c̃′′2 (3.16) 36 3.1. Theory The model constants are given in [87]. Another algebraic model was pro- posed by Vervisch et al. which is based on the concept of flame surface density [157]: ρ̃c = 1 2 (2cm − 1) ρuSLΞ|∇c| c̃ ′′2 c̃(1− c̃) (3.17) where ρu is the density of the unburned gases, Ξ is the flame wrinkling factor and cm is defined as: cm = ∫ 1 0 c ∗ρω̇cP (c∗)dc∗∫ 1 0 ρω̇cP (c ∗)dc∗ (3.18) The value of cm can be assumed to be close to its BML value cm ≈ 3/4 [157]. The flame wrinkling factor is defined as Ξ = |∇c|/|∇c|. Here, |∇c| is the generalized flame surface density, Σ, which stays close to the flame surface density, Σ(c∗), under a flamelet assumption [159]. Several transport equa- tions have been proposed for flame surface density [123]. For low Damköhler number premixed flames, equation 3.13 is still valid based on the arguments of O’Young and Bilger [112] and Chen and Bil- ger [27]. They assert that the time scale ratio, CD, approaches a constant value in this regime for high Damköhler number premixed flames. Neverthe- less, there is not any universal value for the time scale ratio. Swaminathan and Bray [148] showed that a value of CD = 0.25 makes the results of equa- tion 3.13 close to their DNS data of a planar turbulent premixed flame in the flamelet regime. Vervisch et al. [157] computed the time scale ratio based on their DNS data of a turbulent premixed V-flame. They showed that CD is more than 2 in the whole flame brush at one specific axial location. Also, they showed that this value has quite significant fluctuations. The experi- mental results of Chen and Bilger [27] show that this time scale ratio seems to converge to a value between 2 and 3 for low values of Damköhler number. They conducted their experiments on a Bunsen burner using two different hydrocarbon fuels. Lindstedt and Vaos [97] conducted a parametric study on the effect of changing CD in their simulation of the premixed Bunsen flames of Chen et al [28]. Their simulation – based on transported PDF turbulent combustion model – shows that one of the flames, which has the highest Reynolds number, extinguishes for CD = 0.5. They conclude that the time scale ratio should be bigger than one to obtain a stable solution. In this work several models for scalar dissipation rate are attempted in simulation of the highly strained premixed flames of Chen et al.[28]. The results of applying different models are discussed. 37 3.2. Summary of the Governing Equations 3.2 Summary of the Governing Equations A summary of the governing equations for the RANS solver with the im- plemented models is given in this section. The equation for conservation of total mass is: ∂ρ̄ ∂t + ∂ ∂xj (ρ̄ũj) = 0 (3.19) The conservation of momentum reads: ∂ ∂t (ρ̄ũi) + ∂ ∂xj (ρ̄ũiũj) = − ∂p̄ ∂xi − ∂ ∂xj ( ρ̄ũ′′i u ′′ j ) (3.20) where ũ′′i u ′′ j is the Reynolds stress tensor and is modelled using the eddy viscosity concept: ρ̄ũ′′i u ′′ j = −µt ( ∂ũi ∂xj + ∂ũj ∂xi − 2 3 ∂ũk ∂xk δij ) + 2 3 ρ̄k̃δij (3.21a) µt = Cµρ̄ k̃2 ε̃ (3.21b) where transport equations are solved for k̃ and ε̃. These transport equations in Standard k − ε turbulence model are as follows: ∂ ∂t (ρ̄k̃) + ∂ ∂xj (ρ̄ũj k̃) = ∂ ∂xj ( µt σk ∂k̃ ∂xj ) − ρ̄ũ′′i u′′j ∂ũi ∂xj − ρ̄ε̃ (3.22) ∂ ∂t (ρ̄ε̃) + ∂ ∂xj (ρ̄ũj ε̃) = ∂ ∂xj ( µt σε ∂ε̃ ∂xj ) − Cε1 ε̃ k̃ ρ̄ũ′′i u ′′ j ∂ũi ∂xj − C2 ρ̄ ε̃2 k̃ (3.23) Energy equation is stated through the conservation of total enthalpy: ∂ ∂t (ρ̄h̃) + ∂ ∂xj (ρ̄ũj h̃) = ∂ ∂xj ( µt Pr ∂h̃ ∂xj ) (3.24) Transport equation for species mass fraction is as follows: ∂ ∂t (ρ̄Ỹk) + ∂ ∂xj (ρ̄ũj Ỹk) = ∂ ∂xj ( µt Sct ∂Ỹk ∂xj ) + ω̇k (3.25a) ω̇k = ∫ 1 0 ω̇k FPI(c∗)P (c∗)dc∗ (3.25b) 38 3.3. Validation Case In addition to the species transport equations, a transport equation is solved for the combustion progress variable in a premixed flame to better represent the state of combustion. This transport equation is the normalized version of the transport equation for CO2 mass fraction: ∂ ∂t (ρ̄c̃) + ∂ ∂xj (ρ̄ũj c̃) = ∂ ∂xj ( µt Sct ∂c̃ ∂xj ) + ω̇c (3.26) PCM-FPI combustion model requires solving a transport equation for vari- ance of progress variable: ∂ρ̄c̃′′2 ∂t + ∂ ∂xj ( ρ̄ũj c̃′′2 ) = ∂ ∂xj ( ρ̄ νt Sc1 ∂c̃′′2 ∂xj ) + 2ρ̄ νt Sc2 ∂c̃ ∂xj ∂c̃ ∂xj − 2ρ̄̃c + 2c′′ω̇c (3.27a) c′′ω̇c = ∫ 1 0 (c∗ − c̃) ω̇FPIc (c∗)P (c∗)dc∗ (3.27b) In order to close the scalar dissipation rate using equation 3.17, a trans- port equation for Flame Surface Density needs to be solved. The following transport equation is an example based on a Coherent Flame Model [123]. ∂Σ ∂t + ∂ ∂xj ( ũjΣ ) = ∂ ∂xj ( νt ScΣ ∂Σ ∂xj ) + α0 ε̃ k̃ Σ− β0SL Σ 2 1− c̃ (3.28) The modelling constants of these equations are given in Table 3.1. α0 and β0 are set based on the recommendation of reference [131]. Table 3.1: Model coefficients Cµ Cε1 Cε2 σk σε Pr Sct Sc1 Sc2 ScΣ α0 β0 0.09 1.6 1.92 1.0 1.3 1.0 0.7 0.7 0.7 0.7 1.7 1.0 3.3 Validation Case The premixed Bunsen flames of Chen et al.[28] have been simulated in this work. This burner is shown schematically in figure 1 of reference [28]. In this burner, stoichiometric methane-air mixture enters through a central jet at three different Reynolds numbers. The central jet is surrounded with hot pilot products. The diameter of the central jet, D, is 12 mm. The pilot 39 3.4. Numerical Simulation Code Table 3.2: Characteristics of the turbulent Bunsen burner of Chen et al.[28] Flame ReD u0 [m/s] k0 [m2/s2] Da Ka u′/SL τb Bb F1 52,000 65 12.70 1.2 11.0 9.9 6.46 0.41 F2 40,000 50 10.80 1.5 7.3 9.1 6.46 0.44 F3 24,000 30 3.82 2.5 3.4 5.4 6.46 0.75 b Bray number is calculated according to equation 3.10 using α = 0.8. stream comes from a perforated plate with 1175 holes of 1 mm each. The diameter of the pilot annulus is 68 mm. The exit cold flow from each pilot hole is 84 cm/s. Thus, the cold pilot has an average velocity of 0.22 m/s and the hot pilot a velocity of 1.5 m/s, taking into account the effect of density reduction due to heat release. The pilot is also surrounded by an outer cold air flow of 0.22 m/s. A summary of the characteristics of the three flames are shown in Ta- ble 3.2. The values of the Damköhler and Karlovitz numbers show that the flames are all in the thin reaction zone regime. The values of the Bray numbers are all less than one. These Bray numbers have been calculated using the characteristic turbulent velocity fluctuations, u′, at the nozzle exit. Therefore, the reaction zone – which is in the shear layer of the central jet and the pilot flame – should have a substantially higher u′, hence, a lower value for the Bray number. Thus, a gradient diffusion model for the scalar flux is appropriate. 3.4 Numerical Simulation Code The transport equations for c̃, c̃′′2, k̃ and ε̃ are solved along with the Reynolds averaged Navier-Stokes equations. Energy conservation is applied in solving a transport equation for the average total enthalpy assuming an adiabatic flame. Transport equations for the average mass fraction of major species, i.e. CH4, O2, H2O, CO2, H2 and CO, are solved using gradient assumption for turbulent scalar flux and FPI-PDF technique for modelling mean reaction rates. The major species mass fractions can also be retrieved using the FPI manifold rather than solving a transport equation for them. However, looking up the reaction rates and solving transport equations for major species was easier to implement in the existing reacting flow code, es- pecially for the current test case which will be discussed in the next section. FlameMaster [120] was used to obtain the solution of the one-dimensional 40 3.5. Computational Domain and Boundary Conditions laminar unstrained premixed flame for the FPI table with a detailed GRI- Mech 2.11 mechanism [12] for methane-air chemistry. A finite volume low- Mach number pressure-based approach with a PISO algorithm for pressure correction is used for solving the Navier-Stokes equations using C++ ob- jects available in the OpenFOAM package [111]. Flux calculations are done using limited extrapolation with TVD and NVD [75] schemes. In each it- eration, all the transport equations are solved sequentially and implicitly. The system of equations are solved using the Bi-Conjugate Gradient method with preconditioning, with the exception of the Poisson equation for pres- sure which is solved using a geometric algebraic multi-grid method for better convergence. 3.5 Computational Domain and Boundary Conditions The radial profiles of the mean velocity and mean turbulent kinetic energy at the inlet of the central jet have been measured in the experiments and these are taken as the inlet boundary conditions. The dissipation rate in the central jet is set using this assumption [96]: ε̃ = C 3/4 µ k̃3/2 lm (3.29) where Cµ = 0.09 and lm is the Prandtl mixing length. This length scale is lm = κy for a pipe with κ = 0.41; y = D − r is the distance from the wall. The thin wall of the inlet pipe has been modelled as a slip wall. The computational domain is axi-symmetric with outlet boundary conditions for the outer boundaries (figure 3.6). Coarse and fine grids have been used to ensure grid independence of the solution. The boundary condition for reaction progress variable is set to zero in the central jet and one in the hot pilot inlet stream. There is not any physical meaning for reaction progress variable in the air stream surrounding the pilot stream since it contains no fuel. If it is set to zero, then the model would erroneously predict a reaction taking place in the shear layer between the pilot stream and the surrounding air stream. To solve this problem, the reaction progress variable is set to one in this air stream. Since, the reaction rates are zero for c̃ = 1 there is no consumption or production rates for chemical species in case of looking up the reaction rates and solving transport equations for major species. However, if the mass fraction retrieval was the method of choice, then the code would assume that the surrounding 41 3.5. Computational Domain and Boundary Conditions Figure 3.6: Schematic of the computational domain and boundary condi- tions Figure 3.7: Computational grid in OpenFOAM air is also a product mixture (since c̃ = 1). Therefore, the composition of this air stream changes as the code iterates and the code diverges. There are two possible solutions for this problem in order to be able to look up species mass fractions rather than solving a transport equation for them: 1. Limit the computational domain to the main stream and the hot pi- lot stream. In this case, the fluid dynamic effect of surrounding air is missing. Moreover, the outlet boundary gets close to the computa- tional domain – causing non-physical perturbations and convergence problems. 2. Solve a transport equation for Favre-average mixture fraction z̃. This 42 3.6. Results parameter would be a passive scalar taking a value of zero in the air stream and unity in the pilot and the main central stream. Using this parameter the mass fraction of every species can be calculated easily: Ỹk = z̃Y FPIk + (1− z̃)Y airk . The second solution was applied to the code and the results were similar to the approach of solving transport equations for species mass fractions. 3.6 Results 3.6.1 Effect of Different Scalar Dissipation Rate Models As mentioned in the theory section, the best approach for modelling scalar dissipation rate is to solve a transport equation for it. The transport equa- tion given by Swaminathan and Bray [148] was added to the solver. However, no converged solution was obtained, despite several attempts at stabilizing the solution; even when some of the source terms were changed according to suggestions given by Chakraborty et al. [23] a converged solution could not be achieved. This may be because of the fact that turbulence is not well-modelled with a simple two-equation turbulence model. At the next level, algebraic models were applied. The simple algebraic model of equation 3.13 was applied to the simulation of flame F3. As men- tioned above, there is not any universal value for the time scale ratio CD. Thus, the effect of varying this parameter on the results was investigated. Simulations were performed with three different values for the time scale ratio: CD =0.25, 0.5 and 1.0. The axial distribution of the progress variable at the center line of the flame is shown in figure 3.8 for these three different values of CD. This figure shows that the flame speed is directly related to this parameter. As the time scale ratio decreases, the flame speed decreases; hence, the flame length increases. Since the solution is very sensitive to the value of CD and there is no reliable method for selecting this constant, equation 3.13 is not considered a proper closure for scalar dissipation rate. In an attempt to improve the predictions, the model given by Vervisch et al. (equation 3.17) was implemented. The results of using this model are shown in figure 3.9. The model given by Swaminathan and Bray [148] (equation 3.15) was also applied to the simulation of flame F3. The results are shown in fig- ure 3.10. The algebraic model of Kolla et al. [87] (equation 3.16) was also implemented, but the final flame is substantially shorter than is indicated by the experiments. This might be due to using the k − ε model which 43 3.6. Results 0 0.25 0.5 0.75 1 c̃ 0 5 10 15 20 25 30 35 0 0.25 0.5 0.75 x/D c̃ C D = 0.25 C D = 0.5 C D = 1.0 C D = 0.25 C D = 0.5 C D = 1.0 Modified Laminar Flamelet PDF β - PDF Figure 3.8: Effect of changing the time scale ratio on axial distribution of reaction progress variable for flame F3. results in a poor prediction of turbulent kinetic energy in this flame. A comparison between figures 3.9 and 3.10 shows that the algebraic model of equation 3.17 gives better results compared to the prediction obtained us- ing equation 3.15, especially for species mass fractions. However, in both cases the modified laminar flamelet PDF performs better than the β-PDF. Therefore, it appears that, regardless of the scalar dissipation rate model used, the predictions obtained using the new PDF are generally superior to those obtained with the widely-used β-PDF, although neither PDF approach leads to particularly compelling agreement with the experiments. For the discussion that follows, we use equation 3.17 for modelling scalar dissipation rate. 3.6.2 Turbulence Model In premixed flames, velocity fluctuations in the flame brush are not only due to turbulence, but also to the effect of flame front intermittency which causes additional fluctuations [159]. Assuming a very thin flame, the BML theory gives the following expression for velocity fluctuations in the flame 44 3.6. Results 0 1 2 0 0.5 1 1.5 r/D Ũ / U 0 0 1 2 0 2 4 6 r/D k̃ / k 0 0 0.5 1 0 0.5 1 1.5 r/D T̃ / T b 0 0.5 1 0 0.05 0.1 0.15 r/D Ỹ H 2 O 0 0.5 1 0 0.02 0.04 r/D Ỹ C H 4 0 0.5 1 0 0.02 0.04 0.06 0.08 r/D Ỹ C O Figure 3.9: Radial profile of several quantities calculated by using the al- gebraic model of equation 3.17 for scalar dissipation rate. The results are for flame F3 at x/D = 8.5. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. brush [17]: ũ′′2i = c̃(1− c̃)(ūib − ūiu)2 + (1− c̃)u′2iu + c̃u′2ib (3.30) where ūib and ūiu are the conditionally averaged velocities in burned and unburned gases in the i−direction. The first term is the effect of intermit- tency which is not accounted for in the k − ε turbulence model. However, experimental results show that the velocity difference ūib − ūiu decreases as the turbulence intensity increases [17]; therefore, the discrepancy between velocity fluctuations predicted by a k − ε turbulence model and the actual velocity fluctuations decreases as the turbulence intensity increases. For a laminar unstrained premixed flame of stoichiometric methane-air mixture, this velocity difference is almost 2.4 m/s. Thus, the effect of intermittency on the mean turbulent kinetic energy is of order of unity. Hence, using such a model for simulating the highly turbulent premixed flames of Chen et al.[28] is reasonable. The standard k − ε does not have enough dissipation and overpredicts the turbulent kinetic energy in axi-symmetric round jets [124]. Thus, an- other important issue is choosing a proper correction for the standard k− ε 45 3.6. Results 0 1 2 0 0.5 1 1.5 r/D Ũ / U 0 0 1 2 0 2 4 6 r/D k̃ / k 0 0 0.5 1 0 0.5 1 1.5 r/D T̃ / T b 0 0.5 1 0 0.05 0.1 0.15 r/D Ỹ H 2 O 0 0.5 1 0 0.01 0.02 0.05 0.04 0.05 r/D Ỹ C H 4 0 0.5 1 0 0.02 0.04 0.06 0.08 r/D Ỹ C O Figure 3.10: Radial profile of several quantities calculated by using the algebraic model of equation 3.15 for scalar dissipation rate. The results are for flame F3 at x/D = 8.5. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. turbulence model to be able to resolve the round-jet/plane-jet anomaly in this axi-symmetric geometry. One way to address this is to simply increase Cε1 to 1.6. The other possibility, proposed by Pope [124], is to add an extra term to the ε̃ transport equation. This term is vortex-stretching invariant with a constant of Cε3 = 0.79 chosen to produce a correct spreading rate for a simple cold round jet. Figure 3.11 compares the performance of these two approaches in the cold (ie. non-reacting) and hot flow results of flame F3. The correction of Pope gives a good prediction of both the mean axial velocity and turbulent kinetic energy in the cold flow case. However, it has a substantial under-prediction in the reacting case. Herrmann [67] interpreted this under-prediction as a result of not having a mechanism in k − ε to account for the effect of intermittency in velocity fluctuations. However, as mentioned above, this value is of order of unity. Vervisch et al. [157] used Cε3 = 0.5 to overcome this limitation. This approach has no theoretical basis; thus, in this work the first approach (using a value of Cε1 = 1.6) is chosen. Figure 3.11 shows that this final choice still under-predicts the turbulent kinetic energy in the hot flow case. This is likely due to neglecting the effect of flame-generated 46 3.6. Results turbulence. 0 0.5 1 1.5 2 0 0.5 1 1.5 r/D Ū / U 0 Cold 0 0.5 1 1.5 2 0 2 4 6 8 10 r/D k̄ / k 0 Cold 0 0.5 1 1.5 2 0 0.5 1 1.5 r/D Ũ / U 0 Hot 0 0.5 1 1.5 2 0 2 4 6 8 10 r/D k̃ / k 0 Hot Figure 3.11: Radial profile of mean axial velocity and mean turbulent kinetic energy for cold and reacting flow of Flame F3 at x/D = 8.5. Symbols denote experimental data [28], “-·-” is Standard k− ε, “—” is Standard k− ε with Pope’s correction and “- -” is Standard k − ε with Cε1 = 1.6 The main source of the flame-generated turbulence is the pressure di- latation term in the transport equation for the Favre-averaged turbulent kinetic energy [170]. In an attempt to improve the turbulent kinetic energy results, the model of Zhang and Rutland [170] for this term was adopted. In this model, the pressure dilatation term is closed by using the flame surface density function: p′ ∂u′′i ∂xi = ρuSL (∆u)2 2 c̃Σ (3.31) where ρu is the density of the unburned gases, SL is the unstrained laminar flame speed and ∆u is the change of the normal component of velocity across the flame. The results of using this correction are shown in figure 3.12. This 47 3.6. Results figure shows that using this model leads to a substantial over-prediction of turbulent kinetic energy especially at lower axial locations. Thus, for the remainder of the work to be described here, this correction is not added to the turbulent kinetic energy equation and a simple standard k−ε turbulence model with Cε1 = 1.6 is used. It is worth noting that, with this model for turbulence, there is still considerable discrepancy between the predicted value of turbulent kinetic energy and the experiments. However, the primary aim of this work is to show that the new proposed model for the PDF of reaction progress variable can produce usable predictions of premixed flames and to compare this model with the widely used β-PDF. Therefore, as long as this model is kept consistent between simulations with different PDFs, the effect of using different PDFs can be shown. 0 0.5 1 1.5 2 x/D=4.5 0 0.5 1 1.5 r/D x/D=6.5 0 0.5 1 1.5 0 5 10 x/D=8.5 r/D k̃ / k 0 Figure 3.12: The effect of using the model of Zhang and Rutland [170] in the Standard k− ε turbulence model with Cε1 = 1.6 in reacting flow simulation of flame F3 at three axial locations. Symbols denote experimental data [28]. 3.6.3 Discussion Figure 3.13 shows the temperature contour results in flame F3 for the two different PDFs. Figures 3.14 and 3.15 show the production rate and the variance of the reaction progress variable. These figures show that the β- PDF is giving a shorter flame compared to the MLF-PDF. Detailed comparison between the results of using the β-PDF and the modified laminar flamelet PDF are shown in figures 3.16-3.23. Both PDFs lead to over-predictions of the mean axial velocity as shown in figure 3.16, but the modified laminar flamelet PDF gives better predictions. Figure 3.17 shows the results of using these two PDFs in predicting the mean turbulent kinetic energy. The radial profiles of this quantity in the case of flame F3 have two maxima in the experimental data at lower axial locations. One peak is on the unburned side of the flame brush and the other one on the burned side. The inner peak vanishes at downstream 48 3.6. Results (a) β PDF (b) Modified laminar flamelet PDF Figure 3.13: Temperature contours in flame F3. locations. Also for higher Reynolds number flames, there is only one peak in turbulent kinetic energy profiles. The simulations predict these peaks and the modified laminar flamelet PDF gives a relatively better estimation of the radial location of them. The results using both PDFs have under- predictions for turbulent kinetic energy, likely due to using a simple two- equation turbulence model. The temperature field calculations are shown in figure 3.18. This figure shows that the calculated mean temperature field is higher than the actual value. This may be due to neglecting radiation heat transfer and the effect of cold air entrainment in partial premixing in the chemistry model. Using the β-PDF causes a higher over-prediction, especially farther downstream while the modified laminar flamelet PDF has a reasonable prediction. Tak- ing the effect of radiation heat transfer into account, as suggested by Fiorina [41], would likely improve the results. Another possible source of the discrep- ancy between the numerical and experimental results of temperature profiles could be the uncertainty in the inlet conditions, especially the pilot tempera- ture. Also to be considered is the modelling error for scalar dissipation rate. 49 3.6. Results (a) β PDF (b) Modified laminar flamelet PDF Figure 3.14: Production rate of reaction progress variable in flame F3. As discussed before, this parameter has a substantial effect on the scalar field. This is due to the direct dependence of the turbulent flame speed on the scalar dissipation rate value. The turbulent flame speed determines the flame length and shape. The radial profiles of temperature suggest that the flame speed was slightly over-predicted; hence, the flame length was under- predicted. Under-prediction of the flame length causes early heat release at lower axial locations. This can explain the over-heating effect, especially at downstream locations. Using an improved model for scalar dissipation rate ought to improve the results. It should be stressed that the accuracy of the models for scalar dissipation rate strongly depend on the accuracy of the underlying turbulence model. RANS turbulence models have known shortcomings in predicting the turbulent field in reacting flows, especially premixed flames. Therefore, using an LES turbulence model is likely to lead to considerable improvements in the results. Radial profiles of mean CH4 mass fraction are shown in figure 3.19. This figure shows that the modified laminar flamelet PDF gives good predictions while the β-PDF results in under-estimation of this quantity. The O2 mass 50 3.6. Results (a) β PDF (b) Modified laminar flamelet PDF Figure 3.15: Contours of variance of reaction progress variable in flame F3. fraction has also a similar trend as shown in figure 3.20. This fact also implies that the mean flame front shape is not predicted very well by using the β-PDF. Figure 3.22 compares the performance of the two methods in predicting mean H2O mass fraction. While the modified laminar flamelet PDF results indicate a good prediction, using the β-PDF results in an over-prediction, especially for the locations farther downstream. Both methods deviate from the experimental data for mean CO mass fraction, as shown in figure 3.23. This figure shows that the β-PDF has a substantial over-prediction of the mean CO mass fraction while the modified laminar flamelet PDF performs better. This is in part due to the over- prediction of temperature field. Moreover, a good estimation of slow species like CO and NO requires a good model for turbulence-chemistry interaction and the simple flamelet model implemented here may be limited in this regard. 51 3.6. Results 0 0.5 1 1.5 0 0.5 1 r/D Ũ / U 0 F1 0 0.5 1 1.5 2 0 0.5 1 Ũ / U 0 F1 0 0.5 1 1.5 2 0 0.5 1 Ũ / U 0 F1 0 0.5 1 1.5 2 0 0.5 1 1.5 F1 Ũ / U 0 0 0.5 1 1.5 r/D F2 0 0.5 1 1.5 2 F2 0 0.5 1 1.5 2 F2 0 0.5 1 1.5 2 F2 0 0.5 1 1.5 2 x /D =2 .5 r/D F3 0 0.5 1 1.5 2 x /D =4 .5F3 0 0.5 1 1.5 2 x /D =6 .5F3 0 0.5 1 1.5 2 F3 x /D =8 .5 Figure 3.16: Radial profile of mean axial velocity. Symbols denote experi- mental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 52 3.6. Results 0 0.5 1 1.5 0 4 r/D k̃ / k 0 F1 0 0.5 1 1.5 0 4 8 k̃ / k 0 F1 0 0.5 1 1.5 4 8 k̃ / k 0 F1 0 0.5 1 1.5 0 4 8 12 F1 k̃ / k 0 0 0.5 1 1.5 r/D F2 0 0.5 1 1.5 F2 0 0.5 1 1.5 F2 0 0.5 1 1.5 F2 0 0.5 1 1.5 2 x /D =2 .5 r/D F3 0 0.5 1 1.5 2 x /D =4 .5F3 0 0.5 1 1.5 2 x /D =6 .5F3 0 0.5 1 1.5 2 F3 x /D =8 .5 Figure 3.17: Radial profile of mean turbulent kinetic energy. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 53 3.6. Results 0 0.25 0.5 0.75 1 0 0.5 1 r/D T̃ / T b F1 0 0.25 0.5 0.75 1 0 0.5 1 T̃ / T b F1 0 0.25 0.5 0.75 1 0 0.5 1 T̃ / T b F1 0 0.25 0.5 0.75 1 0 0.5 1 1.5 T̃ / T b F1 0 0.25 0.5 0.75 1 r/D F2 0 0.25 0.5 0.75 1 F2 0 0.25 0.5 0.75 1 F2 0 0.25 0.5 0.75 1 F2 0 0.25 0.5 0.75 1 r/D x /D =2 .5F3 0 0.25 0.5 0.75 1 x /D =4 .5F3 0 0.25 0.5 0.75 1 x /D =6 .5F3 0 0.25 0.5 0.75 1 F3 x /D =8 .5 Figure 3.18: Radial profile of reduced temperature. Tb is the adiabatic flame temperature, Tb = 2248. Symbols denote experimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 54 3.6. Results 0 0.5 1 1.5 0 0.03 r/D Ỹ C H 4 0 0.5 1 1.5 0 0.03 Ỹ C H 4 0 0.5 1 1.5 0 0.03 Ỹ C H 4 0 0.5 1 1.5 0 0.03 0.06 Ỹ C H 4 F1 0 0.5 1 1.5 r/D 0 0.5 1 1.5 0 0.5 1 1.5 0 0.5 1 1.5 F2 0 0.5 1 1.5 r/D x /D =2 .5 0 0.5 1 1.5 x /D =4 .5 0 0.5 1 1.5 x /D =6 .5 0 0.5 1 1.5 F3 x /D =8 .5 Figure 3.19: Radial profile of methane mass fraction. Symbols denote ex- perimental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 55 3.6. Results 0 0.5 1 1.5 0 0.15 r/D Ỹ O 2 0 0.5 1 1.5 0 0.15Ỹ O 2 0 0.5 1 1.5 0 0.15Ỹ O 2 0 0.5 1 1.5 0 0.15 0.3 Ỹ O 2 F1 0 0.5 1 1.5 r/D 0 0.5 1 1.5 0 0.5 1 1.5 0 0.5 1 1.5 F2 0 0.5 1 1.5 r/D x /D =2 .5 0 0.5 1 1.5 x /D =4 .5 0 0.5 1 1.5 x /D =6 .5 0 0.5 1 1.5 F3 x /D =8 .5 Figure 3.20: Radial profile of O2 mass fraction. Symbols denote experimen- tal data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 56 3.6. Results 0 0.5 1 1.5 0 0.1 r/D Ỹ C O 2 0 0.5 1 1.5 0 0.1 Ỹ C O 2 0 0.5 1 1.5 0 0.1 Ỹ C O 2 0 0.5 1 1.5 0 0.1 0.2 Ỹ C O 2 F1 0 0.5 1 1.5 r/D 0 0.5 1 1.5 0 0.5 1 1.5 0 0.5 1 1.5 F2 0 0.5 1 1.5 r/D x /D =2 .5 0 0.5 1 1.5 x /D =4 .5 0 0.5 1 1.5 x /D =6 .5 0 0.5 1 1.5 F3 x /D =8 .5 Figure 3.21: Radial profile of CO2 mass fraction. Symbols denote experi- mental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 57 3.6. Results 0 0.5 1 1.5 0 0.085 r/D Ỹ H 2 O 0 0.5 1 1.5 0 0.085 Ỹ H 2 O 0 0.5 1 1.5 0 0.085 Ỹ H 2 O 0 0.5 1 1.5 0 0.085 0.17 Ỹ H 2 O F1 0 0.5 1 1.5 r/D 0 0.5 1 1.5 0 0.5 1 1.5 0 0.5 1 1.5 F2 0 0.5 1 1.5 r/D x /D =2 .5 0 0.5 1 1.5 x /D =4 .5 0 0.5 1 1.5 x /D =6 .5 0 0.5 1 1.5 F3 x /D =8 .5 Figure 3.22: Radial profile of water mass fraction. Symbols denote exper- imental data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 58 3.6. Results 0 0.5 1 1.5 0 0.03 r/D Ỹ C O 0 0.5 1 1.5 0 0.03 Ỹ C O 0 0.5 1 1.5 0 0.03 Ỹ C O 0 0.5 1 1.5 0 0.03 0.06 Ỹ C O F1 0 0.5 1 1.5 r/D 0 0.5 1 1.5 0 0.5 1 1.5 0 0.5 1 1.5 F2 0 0.5 1 1.5 r/D x /D =2 .5 0 0.5 1 1.5 x /D =4 .5 0 0.5 1 1.5 x /D =6 .5 0 0.5 1 1.5 F3 x /D =8 .5 Figure 3.23: Radial profile of CO mass fraction. Symbols denote experimen- tal data [28], “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 59 3.7. Summary 3.7 Summary In this chapter the PCM-FPI turbulent combustion model is used to simulate three premixed Bunsen flames at three different Reynolds numbers. These simulations are done using two different models for the presumed PDF of the reaction progress variable; a β-PDF and a MLF-PDF. The results of the MLF-PDF are better compared to the β-PDF. However, since these results are obtained using a simple RANS turbulence model, the superiority of the the new PDF over the β-PDF is not completely ascertained by this work only. There is a possibility that the new PDF is less sensitive to the errors in the turbulent kinetic energy field compared to the β-PDF. However, as Bray et al. proved in an analytical work [18] and Jin et al. showed by using the DNS data [76], the PDF which is formed based on the chemistry is superior to the β-PDF in that it includes the effects of different chemical reaction rates on the shape of the PDF. This work can be regarded as a complement to the previous results. 60 Chapter 4 CSE in RANS Simulation of Turbulent Premixed Flames 4.1 Introduction In a laminar flamelet model, a turbulent flame is modelled as an ensem- ble of laminar flames. The only effect of turbulence on chemistry that can be considered is the effect of local turbulent stretch. This is only truly valid in the “flamelet regime” in which chemistry time scales are smaller than the turbulence time scales. Conditional Moment Closure (CMC) [7] turbulent combustion model is not limited to the flamelet regime as de- scribed in section 2.3.2. Bushe and Steiner [20] proposed a new CMC-type technique based on solving an inverse problem called Conditional Source- term Estimation (CSE). This model is shown to work very well for diffusion flames [68, 145, 163]. Jin et al. [76] have shown that the CSE version of CMC works in pre- mixed flames in an a priori analysis using DNS data. In this chapter, the extension of their CSE variant of CMC to turbulent premixed flame simu- lations is studied. Numerical challenges of applying CSE are investigated. A Trajectory Generated Low-Dimensional Manifold (TGLDM) model has been used for chemistry reduction. The final CSE-TGLDM combustion model has been used to simulate a turbulent premixed Bunsen flame. 4.2 Theory 4.2.1 Conditional Source-term Estimation Chemical reaction source terms are highly non-linear functions of the scalar field. As a result; one cannot obtain mean chemical reaction source terms from the mean scalar field: ω̇k 6= ω̇k(T , Yk, ρ) (4.1) 61 4.2. Theory Conditional Source-term Estimation uses the CMC hypothesis for closing chemical source terms. In first-order CMC hypothesis, the conditional- averaged chemical source terms are closed by conditional-averaged scalars to suppress the fluctuations and obtain a better estimate of the mean reaction rates [6]: ω̇k|c∗ ≈ ω̇k(T |c∗, Yk|c∗, ρ|c∗) (4.2) where c∗ is a conditioning variable, T |c∗ is the conditional temperature, Yk|c∗ is the conditional mass fraction and ρ|c∗ is the conditional density. ω̇k(T, Yk, ρ) is given by the chemistry model. The unconditional mean chem- ical source term is then obtained by integrating the conditional value over the PDF of the conditioning variable: ω̇k = ∫ 1 0 ω̇k|c∗P (c∗)dc∗ (4.3) where P (c∗) is the probability density function of the conditioning variable. A shape for this function is presumed which is formed by knowing the mean and variance of the conditioning variable. A function composed of two delta functions at c∗ = 0 and c∗ = 1 [134], the β-function and a modified laminar flamelet PDF [76] are among the options. It was shown in chapter 3 that a modified laminar flamelet PDF is a good model for the statistics of the conditioning variable in the RANS simulation of premixed flames. The conditioning variable – which is zero in reactants and becomes one in products – can be the normalized mass fraction of one or a combination of product species, or normalized temperature. In this work the normalized CO2 mass fraction is chosen. The conventional CMC approach is based on solving transport equations for conditional mass fractions. Such transport equations contain several unclosed terms. Closure of these terms are a modelling challenge especially in premixed flames. In CSE, these values are obtained by inverting the following integral equation [20]: Ỹk(~x) = ∫ 1 0 Yk|c∗(~x, c∗)P̃ (c∗; ~x)dc∗ (4.4) where ~x is the spatial coordinate, Ỹk ≡ ρYk/ρ is the density-weighted (Favre- averaged) mass fraction of k-th species and P̃ (c∗; ~x) is the density-weighted PDF of the conditioning variable. The basic assumption of CSE is that the conditional averages are spatially invariant within a known ensemble of 62 4.2. Theory points. This presumption is also commonly taken in CMC combustion model to reduce the computational cost of the simulation. Using this assumption, the above relation can be rewritten in the following form: Ỹk(~xj) = ∫ 1 0 Yk|c∗(c∗)P̃ (c∗; ~xj)dc∗ j ∈ D (4.5) where ~xj is the spatial coordinate of j-th point in the ensemble D. In equation (4.5) the values of Ỹk(~xj) are known – transport equations are solved for them – and the values of Yk|c∗ are unknown. This equation must be solved at every time step in the reacting flow solver to provide the conditional averages and then equation (4.2) is used to close the chemical reaction source terms. 4.2.2 Ill-posed Inverse Problem Equation (4.5) is Fredholm integral equation of the first kind with P̃ (c∗; ~xj) as the kernel. The deconvolution of this integral equation is not a well-posed problem. Singular Value Decomposition (SVD) is a powerful tool to identify the nature of the ill-posedness of the Fredholm integral equation of the first kind. Equation (4.5) can be written in discrete form using m bins for the progress variable: ~b = A~α, Aji = ∫ c2 c1 P̃ (c∗; ~xj)dc∗ (4.6) where bj = Ỹk(~xj) and αi = Yk|c∗i is the conditional average in the i-th bin. The integration to obtain Aji is done over the i-th bin with c1 and c2 as the lower and upper limit. Using the modified laminar flamelet PDF, this integration can be done a priori and tabulated before the simulation. In order to have a better numerical accuracy the Cumulative Distribution Function (CDF) which is more monotone than the PDF is employed instead of the PDF to evaluate the above integral. The singular value decomposition of an n×m matrix A (n > m) is: A = U σ1 σ2 . . . σm 0 . . . 0 VT (4.7) 63 4.2. Theory where σi are the singular values arranged in decreasing order. ~Ui and ~Vi are column vectors of U and V which are orthonormal, i.e. (~UTi · ~Uj) = (~V Ti · ~Vj) = δij . The solution of equation 4.6 can be obtained using the following back-substitution: ~α = m∑ i=1 ~UTi ·~b σi ~Vi (4.8) If matrix A has zero singular values, then the above equation gives di- vision by zero. In this case, A is rank-deficient and has a nontrivial null space. As a result, an infinite set of solutions exists that solve ~b = A~α. According to Hadamard’s definition of a well-posed problem[59], this is type (ii) ill-posedness where the solution is not unique. This situation is likely to happen for an infinitely fast chemistry or a flame in the flamelet regime. This is due to the shape of the probability density function of the reaction progress variable which has two dominant delta functions in these combus- tion regimes. Therefore, some part of the the integral domain in equation 4.6 is effectively zero and A is very sparse. In these circumstances, the flamelet solution is sufficient and the integral equation is not required to be inverted. On the other hand, if matrix A has nonzero but small singular values, equa- tion 4.8 shows that these singular values can amplify any artifact in ~b and result in a non-physical solution. These artifacts are virtually always present due to modelling and discretization errors. According to Hadamard, this is type (iii) ill-posedness where the solution is unstable to small perturbations in the system. There are several regularization approaches to stabilize type (iii) ill- posedness. Tikhonov proposed a regularization approach based on an a priori knowledge of the solution [150]. One way to implement this approach is to solve the following least-squares problem to obtain the solution of equa- tion 4.6: ~α ≈ arg min ~α∗ ∥∥∥∥[AλI ] ~α∗ − [ ~b λ~α0 ]∥∥∥∥2 2 (4.9) where || · ||2 denotes the L2-norm of a vector and I is the identity matrix. In this equation, ~α0 is an a priori knowledge of the solution and λ is the reg- ularization parameter. In zeroth-order (standard) Tikhonov ~α0 is assumed to be zero which means that a small solution is required. In the first-order Tikhonov, a smooth solution is promoted: ~α ≈ arg min ~α∗ ∥∥∥∥[AλL ] ~α∗ − [ ~b 0 ]∥∥∥∥2 2 (4.10) 64 4.2. Theory where L is a discrete approximation to a derivative operator. Bushe and Steiner [20] used the first-order Tikhonov with a second-order derivative operator for L in a CSE simulation of a non-premixed flame. They used λ = Tr(ATA)/Tr(I) for the regularization parameter where Tr is the trace of the matrix. Grout et al. [57] performed their CSE simulation of a non- premixed flame using equation 4.9 for Tikhonov regularization approach. They conducted an unsteady RANS simulation of an igniting jet and they used the previous time-step solution for ~α0. Jin et al. [76] proposed to use the solution of an unstrained one-dimensional laminar premixed flame for ~α0 in a premixed version of CSE. In this way, the combustion model can switch to the flamelet solution in case of having type (ii) ill-posedness. This will render a robust combustion model to handle both flamelet and non-flamelet combustion regimes. Another widely-used method of regularizing integral equations is Trun- cated Singular Value Decomposition (TSVD) approach [60]. This method simply truncates the expansion of equation 4.8: ~α ≈ k∑ i=1 ~UTi ·~b σi ~Vi (4.11) where k < m is the regularization parameter. This method eliminates the effect of small singular values that make the solution unstable. Thus, in- creasing and decreasing k results in less or more filtering of the solution. Equation 4.11 can be expressed in terms of filter functions: ~α ≈ m∑ i=1 fi ~UTi ·~b σi ~Vi (4.12) where filter functions for the TSVD method is: f(σi) = { 1 if i ≤ k 0 if i > k (4.13) using the filter functions, the standard Tikhonov regularization method (equation 4.9 with ~α0=0) can also be expressed in terms of equation 4.12 with the following filter functions: f(σi) = σ2i σ2i + λ2 (4.14) These two filter functions are shown in figure 4.1. This figure clearly shows the difference between Tikhonov and TSVD regularization methods. While 65 4.2. Theory both methods filter out small singular components, TSVD has a sharp cut- off threshold and Tikhonov smoothly cancels the effect of small singular values. 10−2 10−1 100 101 102 0 0.2 0.4 0.6 0.8 1 σ f (σ ) λ = 1 Tikhonov TSVD Figure 4.1: Filter functions for Tikhonov and TSVD regularization methods 4.2.3 Parameter-Choice Methods There is an optimum value for the regularization parameters in both Tikhonov and TSVD methods. A variety of parameter selection methods exist which can be used to rigorously obtain the regularization parameter [62]. L-curve method and Generalized Cross-Validation (GCV) are among the best which have the property to ’squeeze out’ as much information as possible from the kernel and the ~b vector. Also, these methods avoids using too much reg- ularization especially in Tikhonov approach. In the L-curve approach, the norm of the solution vs. the residual norm is plotted in a logarithmic scale. This graph looks like an L shape and the corner of this L shape correspond to the optimum regularization parameter [61]. An operational definition of the corner of the L-curve can be define as 66 4.2. Theory the maximum curvature of the following curve: (ξ(θ), ζ(θ)) = ( log||A~α−~b||22, log||~α||22 ) , θ ∈ {λ, k} (4.15) the curvature is defined as: κ(θ) = ξ′ζ ′′ − ξ′′ζ ′ [(ξ′)2 + (ζ ′)2]3/2 (4.16) In Generalized Cross-Validation (GCV) method [54], the optimum regular- ization parameter is the minimizer of the following function: Υ(λ) = ||A~α−~b||22 (m−∑mi=1 fi)2 (4.17) where fi are the filter functions (equation 4.13 and 4.14) which depend on the singular values. Therefore, GCV method requires performing a singular value decomposition and solving an optimization problem whereas the L- curve method only requires solving an optimization problem. Hence, GCV is more computationally expensive compared to the L-curve method. 4.2.4 Nonlinearity The discussion in section 4.2.2 for the ill-posedness is valid for a linear Fredholm integral equation with a well-defined kernel. Here the kernel of the inversion is approximated with a functional form which is updated from one time step to another time step. Jin et al. [76] have shown that the modified laminar flamelet PDF gives a good approximation for the PDF and results in a good estimation of the conditional moments. However, this is only done in an a priori sense using the DNS data which issues a linear inverse problem. The same scenario happens in every individual time step in the reacting flow solver. Nevertheless, the outcome of the inversion changes the scalar field and this change results in a different kernel for the next time step. Consequently, a non-linear inverse problem is encountered through several time steps. The stability of a non-linear inverse problem is inherently a more complicated mathematical subject. In this study, numerical experiment is used to address this issue. 4.2.5 Chemistry Reduction Trajectory Generated Low-Dimensional Manifold is one way to obtain the lower-dimensional manifold which was first proposed by Pope & Maas [129]. 67 4.2. Theory Van Oijen and de Goey [154] used this idea to build a chemistry database for premixed combustion based on the idea of Flame-Generated Manifold (FGM). Huang & Bushe [68] developed TGLDM for non-premixed flames. Basically, the TGLDM method results in an invariant low-dimensional man- ifold for chemistry in composition space [133]. This means that if one point is initially on the low-dimensional manifold, it stays on the manifold due to having generated the lower-dimensional manifold based on the detailed chemical kinetic trajectories. The composition of different species and reac- tion rates in every trajectory at different states are stored and then tabulated as a function of a few progress variables. Huang & Bushe [68] used the autoignition trajectories to build the TGLDM database. The species mass fractions in an autoignition trajec- tory are controlled by the following equation: ρ ∂Yk ∂t = ω̇k (4.18) The reaction rates in an autoignition manifold are different from those in a premixed flame because of the effect of mass transport. Also, these tra- jectories could not be created for initial temperatures below the autoignition temperature by simply integrating the above equation; special treatments are required [162]. Van Oijen and de Goey [154] solved the governing equations of a one- dimensional laminar unstrained premixed flame (FGM method) to generate trajectories. The mass fractions transport equation in this steady flame is: ρuSL ∂Yk ∂x − ∂ ∂x ( ρDk ∂Yk ∂x ) = ω̇k (4.19) where ρu is the density of unburned gases, SL is the laminar flame speed and Dk is the molecular diffusion coefficient. This method takes the effect of mass transport into account and results in a chemistry database suitable for premixed flame simulations. No special treatment is required for the low-temperature cases. Different trajectories in this method are obtained from different inlet boundary conditions for the one-dimensional premixed flame at x = 0. These inlet conditions are calculated by either doing an element balance for a mixture of some major species or by calculating the constrained equilibrium composition. The TGLDM trajectories for a stoi- chiometric methane-air mixture in the CO2–H2O and CO2–CO domains are shown in Figure 4.2. Also, the mass fraction of H2 and the reaction rate of CH4 are shown in Figures 4.3 and 4.4 as a function of CO2 and H2O mass fractions. 68 4.2. Theory 0 0.05 0.1 0.15 0 0.05 0.1 YCO 2 Y H 2O 0.05 0.1 0.15 0 0.02 0.04 0.06 0.08 0.1 YCO 2 Y CO Figure 4.2: Trajectory Generated Low-Dimensional Manifold. Different tra- jectories are generated using the solution of an unstrained laminar premixed flame. The TGLDM database is usually parametrized as a function of two input parameters. CO2 mass fraction is chosen for the first entry and either H2O or CO for the second entry. This results in saving a considerable computational time in combination with the CSE combustion model. As mentioned in the previous section, the combustion progress variable in the turbulent flame is chosen to be the normalized CO2 mass fraction, i.e c∗ = YCO2/YCO2eq . Thus, the conditional average of the CO2 entry to the TGLDM table, i.e YCO2 |c∗ is a known linear function. As a result, the only unknown parameter required for closure of chemical reaction source is YH2O|c∗ or YCO|c∗. 4.2.6 Governing Equations A reacting flow field is obtained by solving several transport equations for velocity field and scalar field. Navier-Stokes equations conserve momentum. Total mass is conserved throughout the domain by applying the continuity equation. The concentration of every species can be directly retrieved from the TGLDM manifold or can be tracked by solving a transport equation with the chemical source terms obtained form the manifold. At least two transport equations for species, on which the chemistry parametrization is based on, should be solved. The conservation of energy is usually applied by solving a transport equation for temperature or enthalpy [123]. Here, these equations are solved with a RANS approach to model the effect of turbulent fluctuations in the flow field. A simple two-equation k− ε turbulence model 69 4.2. Theory 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.5 1 1.5 x 10−3 Y H 2 OYCO 2 Y H 2 Equilibrium Figure 4.3: H2 mass fraction for different trajectories. is used. The transport equation for a density-weighted average species mass fraction is: ∂ ∂t ( ρ̄Ỹk ) + ∂ ∂xi ( ρ̄ũiỸk ) = − ∂ ∂xi ( ρ̄ũ′′i Y ′′ k ) + ∂ ∂xi ( ρ̄Dk ∂Ỹk ∂xi ) + ω̇k, (4.20) where Y ′′k ≡ Yk − Ỹk. The above equation has two unclosed terms: ω̇k and ũ′′i Y ′′ k . The latter of these is the turbulent scalar flux, modelled by using a simple gradient transport assumption. As mentioned in chapter 3, this assumption is valid when the Bray number is less than one. For Bray number bigger than one, the counter-gradient diffusion happens in which the gradient diffusion assumption is no longer valid [158]. The CSE-TGLDM combustion model is used to model the mean chemical reaction source terms as discussed in section 4.2.1. As mentioned above, this combustion model relies on a presumed shape for the probability density function of the conditioning variable. This requires knowing both mean and variance of the conditioning variable. Two additional transport equations 70 4.2. Theory 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 −800 −600 −400 −200 0 Y H 2 O YCO 2 ω̇ C H 4 Equilibrium Figure 4.4: CH4 reaction rates for different trajectories. are solved for mean and variance of c. Similar to chapter 3, the transport equation for c̃”2 is as follows: ∂ ∂t ( ρ̄c̃′′2 ) + ∂ ∂xi ( ρ̄ũic̃′′2 ) = ∂ ∂xi ( ρ̄ νT Sc1 ∂c̃′′2 ∂xi ) +2ρ̄ νT Sc2 ∂c̃ ∂xi ∂c̃ ∂xi − 2ρDc∂c ′′ ∂xi ∂c′′ ∂xi + 2c′′ω̇c (4.21) where νT is turbulent viscosity. All Schmidt numbers are set to 0.7. Two terms in the above equation are unclosed. The correlation between the source term and the fluctuations of the conditioning variable c′′ω̇c can be modelled by knowing the PDF: c′′ω̇c = (c− c̃)ω̇c = ∫ 1 0 (c∗ − c̃)ω̇c|c∗P (c∗)dc∗. (4.22) The Favre-averaged scalar dissipation rate ρD ∂c′′∂xi ∂c′′ ∂xi is modelled by an al- 71 4.3. Implementation gebraic equation proposed by Vervisch et al. [157]: ρD ∂c′′ ∂xi ∂c′′ ∂xi ≈ 1 2 (2cm − 1) ρuSLΣ c̃ ′′2 c̃(1− c̃) , (4.23) where SL is the laminar flame speed and ρu is the density of unburned gases. A value of 0.75 is assumed for cm [157]. To obtain the generalized flame surface density, Σ, a transport equation is solved based on the Coherent Flame Model (CFM)[123]. This transport equation is as follows: ∂Σ ∂t + ∂ ∂xi ( ũiΣ ) = ∂ ∂xi ( νt ScΣ ∂Σ ∂xi ) + α0 ε̃ k̃ Σ− β0SL Σ 2 1− c̃ , (4.24) where α0 and β0 are model constants, k̃ is the mean turbulent kinetic energy and ε̃ is the mean turbulent dissipation rate. Figure 4.5 summarizes the structure of the CSE-TGLDM combustion model as discussed above. As shown in this flowchart, the CFD code pro- vides the mean and variance of the progress variable to the PDF sub-model and the unconditional average of the H2O mass fraction (or any other species like CO). The PDF sub-model constructs the kernel of the inverse problem (A) for the inversion process in the CSE module. CSE calculates the condi- tional averages from the unconditional average and the kernel. Here, since the conditioning variable is the CO2 mass fraction, the conditional average of CO2 mass fraction is known. Otherwise, the CFD code ought to sup- ply the unconditional average of the CO2 mass fraction to the CSE module and the CSE code should also calculate the conditional average of carbon monoxide. This is the case in non-premixed combustion where the condi- tioning variable is mixture fraction [68, 163]. YH2O|c∗ and YCO2 |c∗ are then supplied to the TGLDM chemistry database and the conditional reaction rates are retrieved. The PDF module calculates the unconditional mean reaction rates and supply that to the CFD code. 4.3 Implementation 4.3.1 Validation Case The Bunsen burner of Chen et al.[28] is selected for validation of the CSE- TGLDM combustion model. As mentioned in chapter 3 section 3.3, in this burner, the stoichiometric mixture of methane and air enters as a central cold jet which has a diameter of 12mm. This central jet is surrounded by 72 4.3. Implementation Figure 4.5: Interaction between different physical sub-models in a premixed CSE-TGLDM reacting flow simulation code. hot pilot products of a methane-air flame at the same equivalence ratio for stabilization. Only flame F3 of this burner is simulated. The Bray number is 0.75 [137]; therefore, the counter-gradient diffusion phenomena does not happen in this flame. 4.3.2 Simulation Code The Navier-Stokes equations are solved with a finite volume low-Mach num- ber pressure-based method using OpenFOAM [111]. A Standard k − ε tur- bulence model is used with Cε1 = 1.6 to account for extra dissipation due to vortex-stretching effects in the axisymmetric round jet. All transport equa- tions are solved using an implicit time integration approach. Second-order TVD and NVD schemes available in OpenFOAM are used for flux calcula- tions. The inlet condition for axial velocity and turbulent kinetic energy is set based on the profiles measured in experiment. The simulation domain is two-dimensional and axisymmetric. The chemistry database is generated by one-dimensional unstrained lam- inar premixed flame calculations using Cantera [55]. A stoichiometric mix- ture of methane and air is considered in order to simulate the Bunsen 73 4.4. Results and Discussion burner of Chen et al. [28]. The GRI-MECH 3.0 [142] is used as a detailed chemical kinetic mechanism for methane-air combustion along with a multi- component diffusion model [36]. The inlet condition for different trajectories is calculated using an element balance on five major species: CH4,O2, CO2, H2O and CO. Each TGLDM trajectory is parametrized by a piecewise cubic polynomial as a function of YCO2 which is the first entry of the chemistry database. For a given CO2 mass fraction, other species mass fractions and reaction rates are obtained from all trajectories. This set of data is then parametrized by another cubic polynomial as a function of the second entry to the chemistry database which is either H2O or CO mass fractions. Thus, given the first and the second entries, all species mass fractions and chemi- cal reaction source terms can be retrieved from the manifold. As mentioned before, since the conditioning variable is the normalized CO2 mass fraction, the conditional average of the CO2 mass fraction is known. Therefore, the chemistry database is only required to be tabulated at distinct values of CO2 to reduce computational time. 4.3.3 Grid Resolution Study Three different grids with 8120, 4060 and 2030 computational cells are used to evaluate the numerical uncertainties and ensure grid independency. Fig- ure 4.6 shows the effect of using these grids on the temperature profiles at four different axial locations. This figure shows that as the grid is refined from 2030 cells to 4060 cells, the numerical results become closer to the ex- perimental profiles. Further refining from 4060 cells to 8120 does not have any significant effect. Hence, the results of the grid with 8120 cells are con- sidered converged. During the rest of this chapter the results of the finest grid are shown. 4.4 Results and Discussion It is assumed that the conditional averages are homogeneous throughout the computational domain and a single profile for the conditional average of H2O or CO is calculated. The computational domain is initialized with the solution of a flamelet solver to be able to get a converged solution more quickly. The singular values of the kernel of the deconvolution is used to diagnose the nature of ill-posedness. This is done at several time steps until convergence including the beginning of the simulation using the initial conditions. Figure 4.7 shows the singular values for both the initial solution and the converged solution. These singular values are all nonzero; hence, 74 4.4. Results and Discussion 0 0.5 1 0 0.5 1 r/D T̃ / T b x/D=2.5 0 0.5 1 0 0.5 1 r/D T̃ / T b x/D=4.5 0 0.5 1 0 0.5 1 r/D x/D=6.5 T̃ / T b 0 0.5 1 0 0.5 1 r/D x/D=8.5 T̃ / T b Experiment 8120 cells 4060 cells 2030 cells Figure 4.6: Grid convergence test matrix A does not have nontrivial null space and the solution ought to be unique. However, the singular values range over several orders of magnitude and the small singular values – which are of order 10−2 and 10−3 – makes the solution unstable and result in ill-posedness type (iii). Therefore, the problem need to be regularized. This criteria is checked at several different time steps to make sure that the nature of the problem does not change. 4.4.1 Regularization and Inversion Results Two different methods have been used for regularization: Tikhonov and TSVD. The initial guess for the solution in the Tikhonov case is first con- sidered to be zero, i.e. ~α0 = 0 and then the solution of a one-dimensional laminar unstrained premixed flame is used. The resulting least-square prob- lem is solved using a non-negative least-square method [93]. As mentioned in the theory section, λ is the regularization parameter in the Tikhonov 75 4.4. Results and Discussion 0 5 10 15 20 25 30 35 40 45 50 10−4 10−3 10−2 10−1 100 101 102 Index Si ng ul ar V al ue s Initial Field Converged Field Figure 4.7: Singular Values of the kernel of the inverse problem for the initial and converged field. method and k in the TSVD method. The L-curve approach is used to find the regularization parameters. Figures 4.8 and 4.9 shows two L-curves for the two regularization meth- ods. These curves are the result of changing the regularization parameters in order to find the optimum value for them. The L-curve for the TSVD method is discrete because k is an integer varying from 0 to 50. Increasing k decreases the filtering of the solution and decreasing k causes less number of singular values to construct the solution; hence, increases the filtering effect in this regularization method. On the other hand, the L-curve for Tikhonov method is continuous because λ can be any number bigger than zero. In this method, increasing λ increases the filtering effect. In Figure 4.8 the flamelet solution is used for ~α0. This figure shows that as λ increases the residual does not increase beyond a certain limit. This infers that the a priori knowledge of the solution is already a good estimate of the solu- tion. This is likely due to the fact that flame F3 is in the lower margin of the thin reaction zone regime and the presumed conditional moments are 76 4.4. Results and Discussion 10−1 100 10−1 100 101 ||A~α – ~b||2 ||~α || 2 Figure 4.8: L-curves for regularization parameters. The solid line represents the Tikhonov regularization using the laminar flamelet solution for ~α0 and the symbols are for TSVD. good estimates of the conditional averages. Figure 4.10 depicts the decon- volution results. This figure shows that different regularization approaches give similar results, specially the Tikhonov approach with two different ini- tial guesses. This figure shows that the solution is almost independent of the a priori estimation. These solutions are based on the optimal value of the regularization parameters. The solution of both regularization meth- ods are not initially smooth; an explicit smoothing filter is applied on the conditional averages. This could be done implicitly by combining both the first-order and second-order Tikhonov regularization and using a first-order or a second-order derivative matrix for L in equation (4.10). Figure 4.10 shows that the estimated conditional moments are close to the laminar flamelet solution in this case. However, there is a noticeable difference between the unstrained laminar solution and the CSE results. It appears that this difference is due to the effect of turbulence on the solution. Turbulence stretches a premixed flame and can eventually penetrate into it 77 4.4. Results and Discussion 10−1 100 101 10−1 100 101 ||A~α – ~b||2 ||~α || 2 Figure 4.9: L-curves for regularization parameters. The solid line represents the Tikhonov regularization without any guess for the solution, i.e. ~α0 = 0 and the symbols are for TSVD. as the turbulence intensity increases. This figure shows that turbulence ef- fects make the conditional average of water mass fraction to be lower than an unstrained flamelet. Unfortunately, there are not any available experi- mental data for validation of this effect in the current case. However, the canonical system of laminar premixed stagnation flame can give an insight into the effect of turbulence induced stretch on a premixed flame. There are several configurations for a premixed stagnation flame [91]. Hawkes and Chen [64] have shown that the reactant-to-product (RtP) configuration reproduces the DNS results well. A schematic of this flame is shown in Fig- ure 2.6. Cantera [55] can be used to obtain different strained flamelets in RtP configuration with stoichiometric mixture of methane-air. Figure 4.11 shows the velocity field in this canonical flame. Figure 4.12 shows the effect of increasing the strain rate on the conditional average of H2O mass fraction. This figure demonstrates that as the strain rate increases, the conditional average of H2O mass fraction becomes lower than the unstrained flamelet 78 4.4. Results and Discussion 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 c∗ Y H 2 O |c∗ Figure 4.10: Conditional average of water mass fraction. Symbols denote the unstrained laminar flamelet solution, “—” is the result of using Tikhonov regularization with the flamelet solution as an initial guess, “-·-” is obtained from using the Tikhonov regularization without any initial guess and “- -” depicts the result of TSVD approach. values; a trend that is in agreement with CSE predictions. It should be emphasized, however, that this effect in CSE might not necessarily be due to turbulent strain. It could be due to some other effect such as unsteady or curvature effects or even be a fortunate artifact due to using a simple turbulence model. Unfortunately, conditional averages were not measured in the experiments. 4.4.2 Simulation vs. Experiment Simulation is run using both H2O and CO mass fractions for TGLDM parametrization in combination with CO2 mass fraction. They show simi- lar results; hence, the rest of this section is based on a TGLDM chemistry database which uses CO2 and H2O mass fractions as the entries. 79 4.4. Results and Discussion −15 −10 −5 0 5 10 −4 −3 −2 −1 0 1 2 3 4 x/ℓf u (m / s) Stagnation Point Expansion in Flame Figure 4.11: Axial velocity in a strained laminar premixed flame. The radial profiles of different quantities were measured and averaged at several axial locations. The simulation results are compared with these pro- files as shown in figures 4.13-4.21. Figure 4.13 shows the comparison between the experimental results and numerical predictions for mean axial velocity at four axial locations. The mean velocity field is captured well. The pre- diction of the CSE combustion model for this quantity is mostly better than the PCM-FPI flamelet model results presented in chapter 3. Radial profiles of turbulent kinetic energy are shown in figure 4.14. This figure shows that the turbulent kinetic energy profiles have two distinct peaks at lower axial locations. These two peaks are clearly captured with the CSE model at correct locations as opposed to the flamelet model which is not accurately predicting these two peaks as shown in figure 3.17. For farther downstream locations, the CSE model is giving an almost perfect prediction at x/D = 6.5 while the flamelet model is not capturing the correct magnitude and radial location for the turbulent kinetic energy peak at this axial location. All in all, the results of the CSE combustion model are generally better than the flamelet model for axial velocity and turbulent kinetic energy. On the other hand, these results are obtained using a simple two-equation RANS turbu- 80 4.4. Results and Discussion 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 c∗ Y H 2 O |c∗ Increased strain Figure 4.12: Effect of increasing aerodynamic strain on Conditional average of H2O mass fraction in a laminar premixed flame of stoichiometric methane- air mixture in RtP configuration. lence model. There is a possibility that the CSE combustion model has less sensitivity to the errors in the turbulence model, or there exists favourable error contributions to the superiority of the CSE combustion model. In this way, the advantage of the CSE combustion model over the flamelet mod- elling is not completely established and further investigations using a better turbulence model is recommended. The radial profile of mean temperature is shown in Figure 4.15. There exists more deviations from the experimental results at lower axial locations. This is likely due to the uncertainty in the inlet condition of the tempera- ture which is not measured in experiments. Also, the effect of using a RANS model and neglecting the radiation heat transfer may in part contribute to these discrepancies. The CSE results are better compared to the flamelet results in figure 3.18 in the previous chapter especially at downstream loca- tions. Radial distributions of four major species mass fractions are shown in 81 4.4. Results and Discussion 0 0.5 1 1.5 0.5 1 1.5 r/D Ũ / U 0 x/D=2.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=4.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=6.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=8.5 Figure 4.13: Radial profiles of mean axial velocity at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. U0 is the mean axial velocity at the inlet. figures 4.16-4.19. These major species are obtained by solving transport equations for them. This figure shows that the major species have been captured with reasonable accuracy with CSE turbulent combustion model. These results also infer that the flame brush and the turbulent flame speed are captured very well. Similar trends were observed in chapter 3 with the PCM-FPI flamelet model using the MLF-PDF. Prediction of two radical species mass fractions are shown in figure 4.20 and figure 4.21. These radical species have been directly retrieved from the chemistry lower dimensional manifold. There is a notable difference between the numerical and experimental results for the radical species mass fractions, most likely due to the inherent limitations of the RANS model. 82 4.4. Results and Discussion 0 0.5 1 1.5 0 3 6 r/D k̃ / k 0 x/D=2.5 0 0.5 1 1.5 0 3 6 r/D k̃ / k 0 x/D=4.5 0 0.5 1 1.5 0 3 6 r/D k̃ / k 0 x/D=6.5 0 0.5 1 1.5 0 3 6 r/D k̃ / k 0 x/D=8.5 Figure 4.14: Radial profiles of mean turbulent kinetic energy at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experi- mental data and solid lines are the simulation results using CSE combustion model. k0 is the mean turbulent kinetic energy at the inlet. Correct estimation of radical species in the current simulation relies mainly on the variance field. There are fundamental limitations in RANS to model the interaction between turbulence and scalar transport in a premixed flame. Using a better turbulence model, like LES, would likely improve the results in this regard. The simulation results do not show a substantial superiority of the CSE combustion model over a PCM flamelet model. Moreover, estimating the conditional averages is more computationally expensive than simply assum- ing a flamelet distribution for them. This might raise questions about ex- ploiting this model in practical applications. However, it is worth mentioning that the current simulation case is very close to the flamelet regime; thus, 83 4.4. Results and Discussion 0 0.25 0.5 0.75 0 0.5 1 r/D T̃ / T b x/D=2.5 0 0.25 0.5 0.75 0 0.5 1 r/D T̃ / T b x/D=4.5 0 0.25 0.5 0.75 0 0.5 1 r/D x/D=6.5 T̃ / T b 0 0.25 0.5 0.75 1 0 0.5 1 r/D x/D=8.5 T̃ / T b Figure 4.15: Radial profiles of mean temperature at four different axial locations for turbulent Bunsen flame F3 [28], Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. Tb is the temperature of the hot pilot products at the inlet. presuming a flamelet like shape for conditional averages is not far from real- ity. On the other hand, CSE does not rely on the limiting assumption of the flamelet models. Therefore, using this model far outside the flamelet regime has the potential to show a better representation of turbulence-chemistry in- teractions. This work is at least a successful detailed extension of using the first-order CMC hypothesis in premixed combustion. CMC and its CSE rep- resentations are well-established combustion models in non-premixed flames and this work is a start for using this model in premixed combustion mode. 84 4.5. Summary 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=2.5 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=4.5 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=6.5 0 0.5 1 r/D Ỹ C H 4 x/D=8.5 Figure 4.16: Radial profiles of mean CH4 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 4.5 Summary In this chapter the CSE turbulent combustion model along with the TGLDM chemistry model is used to simulate a Bunsen flame. This is done with a simple RANS turbulence model. The results of the CSE model are mainly better compared to the flamelet model results in the chapter 3 within the limited RANS context. The main outcome of this chapter is to prove that the CSE turbulent combustion model is stable and converges to a physically meaningful results. This is an important conclusion due to the inherent non- linearity of the underlying inverse problem in the CSE turbulent combustion model. 85 4.5. Summary 0 0.5 1 0.15 0.3 r/D Ỹ O 2 x/D=2.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=4.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=6.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=8.5 Figure 4.17: Radial profiles of mean O2 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 86 4.5. Summary 0 0.5 1 0.085 0.17 r/D Ỹ C O 2 x/D=2.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=4.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=6.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=8.5 Figure 4.18: Radial profiles of mean CO2 mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 87 4.5. Summary 0 0.5 1 0.08 0.16 r/D Ỹ H 2 O x/D=2.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=4.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=6.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=8.5 Figure 4.19: Radial profiles of mean H2O mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 88 4.5. Summary 0 0.5 1 1.5 0.02 0.04 r/D Ỹ O H × 10 x/D=2.5 0 0.5 1 1.5 r/D Ỹ O H × 10 x/D=4.5 0 0.5 1 1.5 0 0.02 0.04 r/D Ỹ O H × 10 x/D=6.5 0 0.5 1 1.5 0 0.02 0.04 r/D Ỹ O H × 10 x/D=8.5 Figure 4.20: Radial profiles of mean OH mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 89 4.5. Summary 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=2.5 0 0.5 1 0 0.015 0.03 r/D Ỹ C O x/D=4.5 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=6.5 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=8.5 Figure 4.21: Radial profiles of mean CO mass fraction at four different axial locations for turbulent Bunsen flame F3 [28]. Symbols denote experimental data and solid lines are the simulation results using CSE combustion model. 90 Chapter 5 PCM-FPI in LES Simulation of Turbulent Premixed Flames 5.1 Introduction Flamelet models are easier to implement and computationally less expen- sive compared to transported PDF and CMC combustion models. This is of vital importance for large-eddy simulations which are computationally- demanding and require high performance parallel computing. The accuracy of flamelet models is satisfactory in the flamelet regime – an important practical regime, particularly for premixed flames. Hence, different variants of flamelet models have been used in commercial and academic simulation tools. In this chapter the Presumed Conditional Moments (PCM) turbu- lent combustion model is used in combination with a FPI chemistry. The model requires a presumed PDF for reaction progress variable. In chapter 3 it was shown that the Modified Laminar Flamelet PDF (MLF-PDF) is a better presumed PDF model compared to the widely-used β-PDF. In this chapter, the MLF-PDF and the β-PDF are used in large-eddy simulation of a turbulent premixed Bunsen burner and the results are compared with experiments. 5.2 Formulation In large-eddy simulation, the three-dimensional unsteady large scale features of the flow field are captured and the small scales are filtered. Defining G(~x) as a spatially-invariant low-pass filter function, the resolved portion of every quantity φ can be expressed as: 〈φ(~x, t)〉 = ∫ V φ(~x′, t)G(~x− ~x′)d~x′ (5.1) 91 5.2. Formulation where 〈φ〉 is the filtered quantity. A density-weighting or Favre-filtering can also be defined as φ̃ = 〈ρφ〉/〈ρ〉. The transport equations for the large scale quantities can be obtained by applying the above filtering operation to the governing equations of a reacting flow. The filtered conservation equation for mass fraction of species k assuming constant diffusivity is as follows: ∂ ∂t ( 〈ρ〉Ỹk ) + ∂ ∂xi ( 〈ρ〉ũiỸk ) = ∂ ∂xi ( 〈ρ〉Dk ∂Ỹk ∂xi ) − (5.2) ∂ ∂xi ( 〈ρuiYk〉 − 〈ρ〉ũiỸk ) + 〈ω̇k〉 where, 〈ρ〉 is the filtered mixture density, Ỹk is the mass fraction of species k, t is time, xi is the spatial coordinate in i-direction, ũi is the filtered velocity in i-direction, Dk is the molecular diffusivity, (〈ρuiYk〉 − 〈ρ〉ũiỸk) is the unresolved turbulent scalar flux and 〈ω̇k〉 is the filtered chemical reaction source term. The filtered chemical reaction source term is the result of interactions between chemical reactions and the sub-grid scale turbulence fluctuations. This interaction in the PCM-FPI flamelet model is reflected through a model for the filtered probability density function (also known as FDF [53]) of a reaction progress variable. The filtered PDF of reaction progress variable c(~x, t) is defined as [32]: P (c∗; ~x, t) ≡ ∫ V δ [ c∗ − c(~x′, t)]G(~x− ~x′)d~x′ (5.3) where δ is the Dirac delta function. The conditional filtered value of every quantity φ(~x, t) can also be defined knowing the PDF: 〈φ(~x, t)|c∗〉 ≡ ∫ V φ(~x ′, t)δ [ c∗ − c(~x′, t)]G(~x− ~x′)d~x′ P (c∗; ~x, t) (5.4) If the above equation is expressed for the chemical reaction source term and is integrated in the progress variable space, the following expression is obtained: 〈ω̇k(~x, t)〉 = ∫ 1 0 〈ω̇k(~x, t)|c∗〉 P (c∗; ~x, t)dc∗ (5.5) In the PCM-FPI combustion model the conditional averages are assumed to be the laminar flamelet values. These values come from the FPI chem- istry model. In the FPI chemistry model a steady one-dimensional laminar unstrained premixed flame is calculated using detailed chemistry. Chemical reaction source terms and species mass fractions are tabulated as a function 92 5.2. Formulation of a progress variable c∗ = Yc/Y eq c . Any linear combination of species mass fraction can be used for Yc so far as it changes monotonically from reactants to products. In chapter 3 and 4, Yc = YCO2 was chosen. However, this progress variable does not change monotonically for rich mixtures. A better candidate is Yc = YCO2 +YCO, which is chosen in this study [41]. The PCM- FPI combustion model requires a functional form for the the FDF which is formed knowing the filtered progress variable c̃ and sub-grid scale variance of the reaction progress variable c̃v ≡ c̃c− c̃c̃ at each point. This parameter can be normalized with maximum possible variance and is called segregation factor Sc = c̃v/(c̃(1− c̃)). Based on these assumptions the filtered chemical reaction source term at each point in space and time (~x, t) is: 〈ω̇k〉 ≈ ∫ 1 0 ωk(c∗)FPI P (c∗; c̃, c̃v)dc∗ (5.6) the above integration can be done a priori and 〈ω̇k〉 can be stored in a table as a function of c̃ and c̃v. The Favre-filtered species mass fractions can also be tabulated using Equation 5.6 with the Favre FDF P̃ . In this way there is no need to solve a transport equation for every species mass fractions; the pre-computed table can be used instead. The β-function has been used as a PDF model in large-eddy simulation of turbulent reacting flows with PCM-FPI combustion model [42, 66]. How- ever, as mentioned in Chapter 3 this PDF is ad hoc and used primarily be- cause it can recover the limits of extremely high and extremely low variance. While it has been used very successfully in non-premied combustion, in pre- mixed combustion, it over-predicts the progress variable chemical reaction source term and does not always give satisfactory results in LES [43, 95]. As mentioned in Chapter 3 the modified laminar flamelet PDF is intrinsi- cally a better presumed PDF compared to the β-PDF. This PDF presumes four different shapes depending on the value of filtered progress variable and sub-filter scale fluctuations of this parameter. This is shown in Figure 5.1. MLF-PDF is not ad hoc and is affected by changes to the chemical kinetic mechanism and how these affect the shape of the laminar premixed flame [76]. 5.2.1 Filtered Laminar Flame Speed Fiorina et al. [43] have studied the ability of the β-PDF to recover the filtered laminar flame speed. This is of particular importance when all the turbulence scales are captured in the grid scale and there is no sub-grid scale wrinkling. In this case, the filtered flame speed S∆ is equal to the laminar 93 5.2. Formulation Figure 5.1: Four possible shapes for the Modified Laminar Flamelet PDF for lean methane-air premixed flame. Sc is the segregation factor Sc = c̃v/(c̃(1− c̃)) flame speed SL. In a canonical case of a steady one-dimensional laminar premixed flame the following equation is valid: ρuSL dc dx = d dx ( ρD dc dx ) + ω̇c(x). (5.7) where ρu is the density of unburned gases, ρ is the mixture density, c is the progress variable and x is the spatial coordinate which is changed from −∞ to ∞. Equation.5.7 can be integrated over the entire spatial coordinate to find the laminar flame speed: ρuSL = ∫ +∞ −∞ ω̇c(x)dx (5.8) If the laminar flame is filtered with the LES filter, one should recover the laminar flame speed (figure 5.2). This can be shown by filtering equation. 5.7 with a spatial filter of size ∆ and integrated to obtain the filtered flame speed 94 5.2. Formulation S∆ [43]: ρuS∆ ≡ ∫ +∞ −∞ 〈ω̇c(x)〉 dx = ρuSL. (5.9) where 〈ω̇c〉 is the filtered reaction rate which is obtained from the com- Figure 5.2: Schematic of laminar vs. filtered flame speed. bustion model. In PCM-FPI combustion model the filtered reaction rate is tabulated as a function of filtered reaction progress variable and sub-filter fluctuations of the reaction progress variable and depends on the form of the PDF. In case of using the β-PDF, the filtered flame speed is not equal to the the laminar flame speed. This can be shown by filtering the solution of equation. 5.7 with a one-dimensional Gassian filter [43]: G(x) = ( 6 pi∆2 )1/2 exp ( −6x 2 ∆2 ) (5.10) where ∆ is the filter size. The filtered progress variable and the sub-filter scale variance can be calculated for different filter sizes ∆. Using these two values, the filtered reaction rate is retrieved from the PCM-FPI table and equation. 5.9 is used to calculate the filtered flame speed S∆. Using both the the β-PDF and the MLF-PDF the filtered flame speed is calcu- lated and shown in Fig. 5.3. This figure shows that the β-PDF does not satisfy Eq. (5.9) while the MLF-PDF gives the filtered flame speed with a relatively good accuracy. These results suggests that the MLF-PDF ought to perform better in LES calculations of premixed turbulent flames in the flamelet regime. 95 5.2. Formulation 100 101 102 0 0.5 1 1.5 2 2.5 3 ∆/ℓF S ∆ /S L β−PDF MLF−PDF Figure 5.3: Filtered laminar flame speed for different filter sizes using β-PDF and MLF-PDF. `F is the laminar flame thickness. 5.2.2 LES SGS Closures In the PCM-FPI combustion model, two transport equations for c̃ and c̃v are solved in addition to the continuity, momentum and energy equations. The transport equation for c̃ is similar to Eq. (5.3). In this work a k-equation SGS model is used where one transport equation is solved for sub-grid scale turbulent kinetic energy k∆ = 0.5(ũiuj − ũiũj) [168]: ∂ ∂t (〈ρ〉k∆) + ∂ ∂xj (〈ρ〉ũjk∆) = ∂ ∂xj [ 〈ρ〉(ν + νt Prk ) ∂k∆ ∂xj ] (5.11) −σij ∂ũi ∂xj − 〈ρ〉Cεk 3/2 ∆ ∆ where ν is the laminar dynamic viscosity of the mixture, νt is the turbulent dynamic viscosity, Prk is a turbulent Prandtle number, ∆ is the filter scale, Cε is the a modelling constant and σij is the sub-grid scale shear stress which 96 5.2. Formulation is obtained through the eddy-viscosity hypothesis: σij = −〈ρ〉νt ( ∂ũi ∂xj + ∂ũj ∂xi − 2 3 ∂ũl ∂xl δij ) + 2 3 〈ρ〉k∆δij (5.12) here νt is calculated using the sub-grid scale turbulent kinetic energy in the k-equation model: νt = Cνk 1/2 ∆ ∆ (5.13) where Cν is a modelling constant. The k-equation sub-grid scale LES model depends on three modelling constants. In this work, the values of Cε = 0.845, Cν = 0.0856 and Prk = 0.25 is chosen [138]. These parameters can also be computed using a dynamic model [82]. The unresolved turbulent scalar flux in Eq. (5.3) is closed using a gradient assumption: σc = 〈ρuiYk〉 − 〈ρ〉ũiỸk ≈ −〈ρ〉 νtSct ∂Ỹk ∂xi (5.14) where Sct is a turbulent Schmidt number and νt is the SGS turbulent vis- cosity. Equation (5.14) is not valid when counter-gradient diffusion occurs. However, σc is small in LES compared to RANS, because in LES the large scale portion of the turbulent scalar flux is captured. Also, as the grid resolution increases, the uncertainty in modelling this term decreases in LES [123]. Moreover, counter-gradient diffusion does not happen when the Bray number is less than one. The transport equation for c̃v after using the gradient diffusion assump- tion for unresolved scalar fluxes with a single Schmidt number is as fol- lows [38]: ∂ ∂t (〈ρ〉c̃v) + ∂ ∂xi (〈ρ〉ũic̃v) = ∂ ∂xi [ 〈ρ〉(D + νt Sct ) ∂c̃v ∂xi ] (5.15) +2〈ρ〉 νt Sct ∂c̃ ∂xi ∂c̃ ∂xi − 2〈ρc〉+ 2 (〈ω̇cc〉 − 〈ω̇c〉c̃) where 〈ω̇c〉 is the filtered chemical reaction source term in the progress vari- able transport equation; 〈ω̇cc〉 is also unclosed and is computed by integrat- ing over the filtered PDF and stored; 〈ρc〉 is the SGS scalar dissipation rate and can be modelled using a linear relaxation: 〈ρc〉 = 〈ρ〉 νtSct c̃v ∆2 (5.16) where ∆ is the filter scale. The turbulent Schmidt number is set to Sct = 1.0 [63]. 97 5.3. Implementation 5.3 Implementation The chemistry library was generated by solving a one-dimensional unstrained laminar premixed flame calculated with Cantera [55]. Continuity, momen- tum and energy equations are solved along with a transport equation for the SGS turbulent kinetic energy. Two additional transport equations are solved for the filtered and SGS fluctuation of reaction progress variable. Species mass fractions are tracked by solving transport equations. These transport equations are solved using a compressible density-based approach available in the CFFC code [46]. In this code, the temporal derivatives are solved using a second-order Runge-Kutta scheme and spatial derivatives are solved via a second-order finite volume approach. The inviscid fluxes are calculated using limited linear reconstruction and a Reimann-solver-based flux calculation approach. Viscous fluxes are computed via a hybrid average gradient-diamond path. The details of the computational approach can be found in [46]. A multi-block hexahedral mesh is used and all the compu- tations are done in parallel using a domain decomposition approach with MPI. The premixed Bunsen flame of Chen et al. [28] is simulated in this work. In this burner the stoichiometric mixture of methane and air enters as a central cold jet at different Reynolds numbers. This central jet is surrounded by hot pilot products of a methane-air flame at the same equivalence ratio for stabilization. This geometry is modelled with approximately 1.6 million hexahedral computational cells (figure 5.4) and is run on 128 processors. In this study flame F3 of this burner is simulated which has a mean inlet velocity of 30 m/s. For this flame the Bray number is less than one as mentioned in chapter 3; hence, the gradient diffusion assumption should be valid. 5.3.1 Turbulent Inlet Boundary Condition In order to obtain a realistic turbulent field in LES and DNS, the turbu- lent inlet boundary condition needs to be specified with accurate turbulent fluctuations. These fluctuations mainly depend on the upstream conditions. Therefore, the most accurate method for generating inlet turbulent bound- ary condition is to run a precursor simulation of the upstream flow field. This is the so called ”recycling method”. This method is applicable to the cases where the upstream condition is known and can be easily simulated. Oth- erwise, synthetic turbulence is generated and used as the inflow condition. This synthetic turbulence can be homogeneous isotropic turbulence based 98 5.3. Implementation Figure 5.4: Velocity and vorticity contours in flame F3. on a presumed turbulence spectrum [135] or inhomogeneous turbulence gen- erated based on the statistical information of the upstream condition from experiment such as first and second moments of the velocity field [74, 94]. The central jet in the Bunsen burner of Chen et al. [28] is a circular pipe. Hence, a precursor large-eddy simulation of a periodic pipe is perfomed to generate the turbulent inlet velocity filed. The pipe has a length-to-diameter ratio of 20 and is periodic in the axial direction. The no-slip conditin is set on walls and the turbulent boundary layer is captured. In order to correctly simulate the turbulent coherent structures in large-eddy simulation of a pipe flow, the first grid point should be at y+ < 2 [119]. Here, y+ is the distance from the wall normalized by viscosity ν and friction velocity uτ = U √ f/8, where U is the mean velocity and f is the Darcy friction factor. Also, the grid spacing should be of order ∆z+ = 50− 150 and r∆θ+ = 15− 40 [119]. The final mesh is multi-block and structured with approximately 1.1 million cells. OpenFOAM [111] CFD solver is used for large-eddy simulation of this periodic pipe. The computational domain is initialized with parabolic lam- inar velocity profile. Since the Reynolds number is high enough, numerical perturbations causes instability and transition to a fully-developed turbu- lent flow. Figure 5.7 shows the L2-norm of the residual of velocity field in the domain as a function of time. This figure shows the transition from an initial laminar field to a stationary turbulent regime with significant velocity fluctuations. The contours of the magnitude of velocity and sub-filter scale turbulent kinetic energy in the pipe is shown in figure 5.6. The lateral velocity vectors 99 5.4. Results and Discussion 0 0.05 0.1 0.15 0.2 0.25 10−6 10−5 10−4 10−3 10−2 10−1 100 Time(s) ||U n + 1 − U n || Turbulent Laminar Figure 5.5: The temporal variation of the L2-norm of the residual of the velocity field in LES of a periodic pipe. at one cross-section of the pipe is shown in figure 5.7. This figure shows that the turbulent boundary layer is captured within the computational grid. A more qualitative comparison between the profile of the averaged veloc- ity field and the ”Law of the wall” is shown in figure 5.8. This figure shows that the LES results are very similar to reality in both the viscous sublayer and the turbulent outer layer. The velocity field at the outlet surface is stored and recycled at the inlet boundary condition of the Bunsen burner during the actual simulations. 5.4 Results and Discussion Simulation is run with the initialization of a cylinder of reactants inside the simulation domain otherwise filled with products. The simulation is run until the total heat release in the domain achieves a statistically stationary condition. The statistics are collected thereafter over several large eddy turn-over time. Figure 5.9 shows the instantaneous contours of density and 100 5.4. Results and Discussion Figure 5.6: Contours of mean velocity and sub-grid scale turbulent kinetic energy in the LES simulation of a periodic pipe figure 5.10 shows the instantaneous contours of velocity and vorticity field. This figure shows the vorticity generation at the shear layer between the central jet and the pilot stream. This is due to the large difference between their velocity. The instantaneous contours of two radical species are shown in figure 5.11. This figure shows the effect of the large scale turbulent structures on formation and distribution of radicals at different parts of the flame. This phenomena cannot be captured with RANS turbulence models. The radial profiles of different quantities were measured in the experi- ments at four axial locations. Figure 5.12 shows the comparison between the numerical and experimental results of the axial velocity profiles. The devi- ations at downstream locations are likely due to using a simple k-equation sub-filter scale LES model. This figure shows that the β-PDF performs 101 5.4. Results and Discussion Figure 5.7: Crosswise velocity vectors in LES of a periodic pipe. slightly better compared to the MLF-PDF. Figure 5.13 shows the compar- ison between the experimental results and numerical predictions for mean temperature. Predictions are satisfactory except at the lowest axial loca- tions x/D = 2.5. This is likely due to uncertainty in the inlet boundary condition for the pilot stream. This figure shows that the β-PDF over- predicts the temperature due to over-estimation of the reaction rates and local over-heating. Radial profiles of major species mass fractions are shown in figures 5.14- 5.17. These figures show that the major species have been captured with reasonable accuracy with both the β-PDF and the MLF-PDF. Figures 5.18 shows the radial profiles of OH mass fraction, a radical species that is hard to predict. This figure shows that both PDFs under-estimate the OH mass fractions; however, the results of using the MLF-PDF are closer to the ex- perimental results. Both methods have a reasonably good prediction of the CO mass fraction as shown in figure 5.19. This figure shows that the MLF- 102 5.5. Summary 100 101 102 103 0 5 10 15 20 25 y+ U + LES U+ = y+ U+ = 5.0 + 2.5 ln y+ Figure 5.8: Mean velocity profiles of LES simulation of a periodic pipe. PDF provides a better prediction of the trend in the experimental results at farther downstream locations. At these locations the CO mass fraction increases, has a peak at r/D = 0.5 and then decreases. The MLF-PDF cap- tures this trend while in the β-PDF results the CO mass fractions remain constant up to r/D = 1.0 before decreasing. 5.5 Summary In this chapter, the PCM-FPI turbulent combustion model is used in large- eddy simulation of a Bunsen flame. Two different presumed PDF models are used in this work: β-PDF and MLF-PDF. The results of MLF-PDF in prediction of species mass fractions are slightly better compared to the β-PDF results. While the β-PDF is ad hoc, the MLF-PDF is based on the chemistry and has a much better prediction of the filtered flame speed as shown in section 5.2. Therefore, the MLF-PDF is a better candidate for the 103 5.5. Summary statistical representation of the reaction progress variable in LES. 104 5.5. Summary Figure 5.9: Instantaneous contour of density in kg/m3 for flame F3. 105 5.5. Summary Velocity (m/s) −0.02 0 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 5 10 15 20 25 30 35 Vorticity (1/s) −0.02 0 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 2000 4000 6000 8000 10000 12000 Figure 5.10: Instantaneous velocity and vorticity contours in flame F3. 106 5.5. Summary ỸCO −0.02 0 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 0.005 0.01 0.015 ỸH2 −0.02 0 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 1 2 3 4 x 10−4 Figure 5.11: Instantaneous contours of two radical species in flame F3. 107 5.5. Summary 0 0.5 1 1.5 0.5 1 1.5 r/D Ũ / U 0 x/D=2.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=4.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=6.5 0 0.5 1 1.5 0 0.5 1 1.5 r/D Ũ / U 0 x/D=8.5 Figure 5.12: Radial profiles of mean axial velocity at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. U0 is the mean axial velocity at the inlet. 108 5.5. Summary 0 0.25 0.5 0.75 0 0.5 1 r/D T̃ / T b x/D=2.5 0 0.25 0.5 0.75 0 0.5 1 r/D T̃ / T b x/D=4.5 0 0.25 0.5 0.75 0 0.5 1 r/D x/D=6.5 T̃ / T b 0 0.25 0.5 0.75 1 0 0.5 1 r/D x/D=8.5 T̃ / T b Figure 5.13: Radial profiles of mean temperature at four different axial locations for turbulent bunsen flame F3 [28]. Tb is the adiabatic flame tem- perature, Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 109 5.5. Summary 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=2.5 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=4.5 0 0.5 1 0 0.03 0.06 r/D Ỹ C H 4 x/D=6.5 0 0.5 1 r/D Ỹ C H 4 x/D=8.5 Figure 5.14: Radial profiles of mean CH4 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 110 5.5. Summary 0 0.5 1 0.15 0.3 r/D Ỹ O 2 x/D=2.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=4.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=6.5 0 0.5 1 0 0.15 0.3 r/D Ỹ O 2 x/D=8.5 Figure 5.15: Radial profiles of mean O2 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 111 5.5. Summary 0 0.5 1 0.085 0.17 r/D Ỹ C O 2 x/D=2.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=4.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=6.5 0 0.5 1 0 0.085 0.17 r/D Ỹ C O 2 x/D=8.5 Figure 5.16: Radial profiles of mean CO2 mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 112 5.5. Summary 0 0.5 1 0.08 0.16 r/D Ỹ H 2 O x/D=2.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=4.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=6.5 0 0.5 1 0 0.08 0.16 r/D Ỹ H 2 O x/D=8.5 Figure 5.17: Radial profiles of mean H2O mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 113 5.5. Summary 0 0.5 1 1.5 0.02 0.04 r/D Ỹ O H × 10 x/D=2.5 0 0.5 1 1.5 r/D Ỹ O H × 10 x/D=4.5 0 0.5 1 1.5 0 0.02 0.04 r/D Ỹ O H × 10 x/D=6.5 0 0.5 1 1.5 0 0.02 0.04 r/D Ỹ O H × 10 x/D=8.5 Figure 5.18: Radial profiles of mean OH mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 114 5.5. Summary 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=2.5 0 0.5 1 0 0.015 0.03 r/D Ỹ C O x/D=4.5 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=6.5 0 0.5 1 1.5 0 0.015 0.03 r/D Ỹ C O x/D=8.5 Figure 5.19: Radial profiles of mean CO mass fraction at four different axial locations for turbulent bunsen flame F3 [28]. Symbols denote experimental data, “- -” is the result of using the β-PDF and “—” is the result of using the modified laminar flamelet PDF. 115 Chapter 6 CSE in LES Simulation of Turbulent Premixed Flames Application of the conditional source-term estimation turbulent combus- tion model to large-eddy simulation is more challenging compared to RANS turbulence models. This is mainly due to parallel processing and high- performance computing requirements. Implementation of CSE in a parallel code is more difficult compared to the PCM-FPI flamelet model. The reason is that flamelet models are local i.e. the reaction rates at one computational cell only depends on the properties of that cell, whereas in CSE an ensemble of computational cells contribute to the reaction rates of that one cell. In this chapter, CSE is implemented in the LES context. The detailed chemical kinetics is reduced using a TGLDM model similar to chapter 4. The results section consists of using the CSE-TGLDM model of chapter 4 in two different Bunsen burners and the main outcome of these simulations is stability and convergence to meaningful solutions. 6.1 Theory As described in section 4.2.5, the outcome of the TGLDM chemistry reduc- tion method is a function of two control variables in a premixed flame: ω̇k(T, ρ, Yk) ≈ ω̇k(c1, c2) (6.1) where in methane-air combustion c1 = YCO2 and c2 = YH2O. Using the probability density function of c1, one can write the following equation for the filtered reaction rates: 〈ω̇k〉 = ∫ 1 0 〈ω̇k|c∗1〉P (c∗1)dc∗1 (6.2) The conditionally-filtered reaction rates in the above equation can be ap- proximated using the first-order Conditional Moment Closure (CMC) hy- pothesis [7] with TGLDM chemistry: 〈ω̇k|c∗1〉 ≈ ω̇k(c∗1, 〈c2|c∗1〉) (6.3) 116 6.1. Theory where 〈c2|c∗1〉 is the conditionally-filtered value of c2 conditioned on c1. A transport equation is solved for 〈c2|c∗1〉 in the conventional CMC ap- proach [83]. However, as mentioned in chapter 4 such transport equation has closure issues, especially in premixed combustion where there is a strong coupling between transport and chemistry. Conditional Source-term Esti- mation (CSE) [20] is another approach to obtain the conditional scalar field as described in chapter 4. This is done by solving the following integral equation: c̃2(~x) = ∫ 1 0 〈c2|c∗1〉P̃ (c∗1; ~x)dc∗1 (6.4) In this equation, c̃2 is the Favre-filtered value of c2 for which a transport equation is solved; P̃ (c∗1; ~x) is obtained from a presumed PDF model. As- suming that 〈c2|c∗1〉 is homogeneous within an ensemble of points, the above integral equation can be discretized and inverted to get 〈c2|c∗1〉. Similar to the RANS analysis in chapter 4, equation 6.4 is discretized using m bins in the c∗1 space and is written for n points in the following matrix form: ~b = A~α (6.5) where bj = c̃2(~xj) for j = 1 : n and αi for i = 1 : m is the discrete form of 〈c2|c∗1〉; A is called the design matrix and its elements are calculated from integration of the PDF in each bin for every point in the ensemble. This inverse problem is ill-posed; hence, the Tikhonov regularization technique is used: ~α = arg min ~α∗ ( ||A~α∗ −~b||22 + λ2||~α∗ − ~α0||22 ) (6.6) where || ||2 is the L2-norm of a vector, λ is the regularization parameter and ~α0 is an a priori knowledge of the solution. As mentioned in chapter 4, the solution of a one-dimensional laminar premixed flame can be used for ~α0. In this way, CSE is capable of returning the flamelet results in the flamelet regime. In this work, two different jet flames are simulated. In an axi-symmetric round jet the variation of the conditional averages in the tangential direction is ignored. Hence, different CSE ensembles can be formed with computa- tional cells at the similar axial locations and the first-order CMC hypothesis is used. The L-curve method can be used to find an optimum value for the regularization parameter [61]. 117 6.2. Implementation 6.1.1 Stability The CSE inverse problem is linear in each time step. Nevertheless, the solution of the inverse problem in one time step affects both the design matrix and the right hand side vector ~b in the next time step. Therefore, the inverse problem is non-linear over consecutive time steps. The stability of a non-linear inverse problem cannot be determined easily and a numerical experiment needs to be done to demonstrate stability. 6.2 Implementation Two different cases were attempted in this work. The first case is flame F3 of the Bunsen burner of Chen et al. [28]. The second case is the Bunsen burner of Yuen and Gulder [169] which is referred to as the Gulder burner. Both cases have premixed methane-air coming from a central jet surrounded by hot pilot product. The central jet in the Chen burner is a stoichiometric mixture with a mean velocity of 30 m/s. The pilot has a lower velocity of 1.5 m/s which causes shear-generated turbulence in addition to the incoming turbulence which has similar statistics to a pipe flow. The Gulder burner has a lean CH4-air mixture at φ = 0.7 and the pilot has a high velocity. Hence, the main source of turbulence is from the incoming central jet which is grid-generated turbulence. This burner has a much higher turbulence intensity compared to the first case. The TGLDM chemistry library was generated by solving one-dimensional unstrained laminar premixed flames using Cantera [55]. The library is tabu- lated as a function of two control variables: c1 = YCO2 and c2 = YH2O. The incoming turbulence in the Chen burner comes from a separate precursor simulation of a periodic pipe. In the Gulder burner, homogeneous isotropic turbulence is used for this purpose. In this work, the first-order CMC hypothesis is used. Tikhonov reg- ularization is used to solve the underlying inverse problem. The L-curve approach is used to optimally find the regularization parameter every 50 time steps. In order to use CSE with the TGLDM chemistry model in the multi-dimensional CFD code, seven transport equations are solved: con- tinuity, momentum and energy equations, the transport equation for the sub-filter scale turbulent kinetic energy and two transport equations for fil- tered control variables c̃1 and c̃2. Also another transport equation needs to be solved for sub-grid fluctuations of c1. These transport equations are solved using a compressible density-based approach available in CFFC [46]. The details of the computational approach can be found in [46]. 118 6.3. Results 6.3 Results 6.3.1 Test Case –I Based on LES results, the location of all the points in the turbulent premixed flame can be shown on a turbulent premixed regime diagram [122]. This is shown in figure 6.1 for flame F3 of the Chen burner. According to this figure, the Chen burner is partly in the flamelet regime and patly in the thin reaction zone regime. The singular values of the design matrix at several axial locations are shown in figure 6.2. This matrix is rank-deficient at all axial locations. Therefore, the CSE turbulent combustion model returns the flamelet conditonal averages as shown in figure 6.3. This flame is primarily in the flamelet regime (as seen in figure 6.1) and this behaviour is expected in that regime. Also, according to figure 6.1 the filter size ∆ is between 2 to 10 times larger than the laminar flame thickness. Using a finer grid might result in deviation of the conditional averages away from the flamelet ones. Figure 6.1: Location of the Chen burner on LES turbulent premixed regime diagram. 119 6.3. Results 0 10 20 30 40 50 10−8 10−6 10−4 10−2 100 102 104 Index Si ng ul ar V al ue s z/D=0 z/D=1 z/D=2 z/D=5 z/D=8 Figure 6.2: Singular values of the design matrix A for Chen burner at several different axial locations. Since all of the conditional averages calculated using CSE return the flamelet conditional average, CSE’s predictions of this flame are identical to those obtained using the PCM-FPI model, as presented in chapter 5. As such, thy will not be presented here. However, it is worth noting that these CSE simulations did prove that a stable solution can be obtained in LES with CSE. Also, these results show that for a flame in the flamelet regime, the CSE model correctly returns the flamelet model results. 6.3.2 Test Case – II The second test case is the Gulder burner which is mainly in the thin reaction zone regime as shown in figure 6.4. This is due to the high turbulence intensity coming through the inlet of this burner. A non-trivial fraction of the flame is in the Broken Reaction Zone regime, which suggests that a flamelet model would not be appropriate for this flame. Figure 6.5 shows the singular values of the design matrix at different axial locations. This figure shows that, unlike the Chen burner, the design matrix in this burner is not rank-deficient at lower axial locations. Therefore, the design matrix contains enough information to reconstruct the conditional 120 6.3. Results 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 c∗1 c 2 |c∗ 1 Figure 6.3: The conditional averges in Chen burner at different axial loca- tions. Symbolds denote the unstrained laminar flamelet results. averages independent of the flamelet solution and the effect of the a priori solution is only regularization. As a result, the conditional averages deviate from the flamelet solution at the lower axial locations, as shown in figure 6.6. Further downstream these conditional averages converge to the flamelet so- lution. This is likely due to the fact that the turbulence in this burner is generated using a grid; this type of turbulence decays exponentially down- stream of the grid. The actual conditional averages were not measured in the experiment. However, this type of behaviour, where the conditional av- erages tend to a straight line for high turbulence intensities was observed in a recent DNS of premixed flames in thin reaction zone and distributed reaction zone regimes [4]. Furthermore, similar to the discussion in section 4.4, these results can be justified based on the laminar strained flamelets in a counter flow configuration as shown in figure 6.7. The instantaneous contours of velocity and vorticity of the Gulder burner are shown in figure 6.8. This figure shows that the turbulence intensity of this flame is very high. Also, the vorticity contours indicate that the main source of turbulence is the incoming turbulence from the central jet. Figure 6.9 shows the filtered temperature and methane mass fraction at one 121 6.4. Summary Figure 6.4: Location of the Gulder burner on the premixed LES regime diagram. snapshot of the flame. The LES results for temperature were spatially averaged at different angles for each axial location to obtain the radial profiles for comparison with experimental results. The experimental radial profiles are the ensemble average of 300 snapshot of the flame [169]. These results are shown in Fig. 6.10. 6.4 Summary In this chapter, the CSE-TGLDM turbulent combustion model is used in large-eddy simulation of two different turbulent premixed flames. One flame is primarily in the flamelet regime and the other one in the thin reaction zone regime and partly in the broken reaction zone regime. The results show that CSE returns the flamelet conditional averages in the flamelet regime. Outside of the flamelet regime, the final conditional averages are different from the flamelet one in a way that is consistent with the reported values in literature. The main outcome of these simulations is that CSE in premixed combustion is stable and converges to meaningful results. As 122 6.4. Summary 0 10 20 30 40 50 10−8 10−6 10−4 10−2 100 102 Index Si ng ul ar V al ue s x/D=0 x/D=1 x/D=2 x/D=5 x/D=8 Figure 6.5: Singular values of the design matrix A for Gulder burner at several different axial locations. mentioned in section 6.1.1, the stability of the CSE inverse problem cannot be mathematically demonstrated easily due to the non-linear nature of the problem. Hence, CSE was implemented in a simulation code to perform numerical experiments and the stability of the method was substantiated. 123 6.4. Summary 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 c∗1 c 2 |c∗ 1 Increased z/D Figure 6.6: Filtered control variable c2 conditioned on c1 for the Gulder burner. 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 c1 c 2 Increased strain Figure 6.7: c2 vs. c1 in the counter-flow premixed flame configuration. 124 6.4. Summary Figure 6.8: Instantaneous velocity and vorticity contours in Gulder burner. 125 6.4. Summary Figure 6.9: Instantaneous temperature and methane mass fraction contours in Gulder burner. 126 6.4. Summary −1 −0.5 0 0.5 1 0 0.5 1 r/D c̃ T −1 −0.5 0 0.5 1 0 0.5 1 r/D c̃ T z/D=7 z/D=5 Figure 6.10: Radial profiles of the normalized temperature c̃T = (T̃ − Tu)/(Tb − Tu) at two axial locations. Solid line is the averaged LES simula- tion results and symbols denote experimental results [169]. 127 Chapter 7 Conclusions and Future Work 7.1 Conclusions In this chapter, the conclusions of the main four chapters of this thesis are first explained and then the final conclusion is stated. 7.1.1 Chapter 3 In this chapter first the PCM flamelet model with an FPI representation of chemistry was employed to simulate a highly-strained turbulent premixed flame. A progress variable approach for flamelet modelling using a presumed shape for the probability density function of reaction progress variable was used to calculate the mean reaction rate. The effects of using the β-function for the PDF and Modified Laminar Flamelet PDF (MLF-PDF) were com- pared. There are some noticeable deviations between the experimental data and the numerical simulations using both the β-PDF and the modified lam- inar flamelet PDF. These deviations appear mostly to be due to using a simple RANS model for turbulent fluctuations. However, the results of this work indicate that the new PDF is a viable option for the statistical de- scription of the reaction progress variable in a turbulent premixed flame. This new PDF is based on the chemistry; thus, it is an intrinsically better choice than the ad hoc β-PDF. That the new PDF gives results that are superior to those obtained using the β-PDF with only a modest additional computational cost (in a pre-processing operation only) suggests that it is preferable to the available alternatives. 7.1.2 Chapter 4 In this chapter the Conditional Source-term Estimation in combination with a TGLDM model for chemistry was applied to RANS simulation of the same turbulent premixed flame. This combustion model successfully converged to 128 7.1. Conclusions a physical solution. Hence, it was proved by numerical experiment that, at least for the flame studied here, the underlying nonlinear inverse problem in CSE has a unique and stable solution. The results show that this model can predict the mean velocity and temperature field with reasonable accuracy. Major species mass fractions are also captured well, while the predictions for OH and CO radicals show some discrepancy from the experimental val- ues. This was attributed to modelling errors introduced by using a RANS turbulence model. The results of the CSE combustion model were close to the flamelet results in the current test case. It is expected to see the advan- tages of the CSE combustion model over the flamelet models become more pronounced under conditions that depart from the flamelet regime. 7.1.3 Chapter 5 In this chapter the PCM-FPI turbulent combustion model was used in large- eddy simulation of the same turbulent premixed flame. The results of using the Modified Laminar Flamelet PDF (MLF-PDF) were compared with the widely-used β-PDF. The results with the MLF-PDF were marginally better compared to those from the β-PDF. Also, it was demonstrated that predic- tion of the filtered laminar flame speed using the MLF-PDF is superior to that obtained using the β-PDF especially when the filter scale is bigger than the laminar flame thickness. Clearly, the MLF-PDF should be the preferred option for presumed PDF modelling of turbulent premixed flames. 7.1.4 Chapter 6 In this chapter the Conditional Source-term Estimation turbulent combus- tion model was used in large-eddy simulation of turbulent premixed flames. This model was coupled with a TGLDM model for chemistry. The first-order Conditional Moment Closure (CMC) hypothesis was used in CSE to close the chemical reaction source terms. The conditional moments were calcu- lated by solving an integral equation. This is a non-linear inverse problem, where the solution of one time step affects the inversion in the next time step. In this work, it was shown that this algorithm is stable and converges to meaningful results both in low and high turbulence intensity regimes. 7.1.5 Summary of Accomplishments The main outcomes of this thesis can be summarized as follows: 129 7.2. Future Work • The modified laminar flamelet PDF model was first applied to simula- tion of turbulent premixed flames in this work. This new PDF model represent an improvement over the β-PDF. Further improvement using the linear eddy model (LEM) might lead to better predictions. • It was first shown in this work that the the CSE turbulent combustion model for premixed combustion appears to be stable and converges to physically meaningful results. • CSE seems to provide results consistent with our expectations. It defaults to flamelet-like behaviour in the flamelet regime. Outside of the flamelet regime, it provides deviations from flamelet behaviour that are consistent with other results in the literature. 7.2 Future Work The simulations in this work are limited to a few turbulent premixed Bun- sen flames. It is suggested to use the CSE model in simulation of other laboratory-scale burners such as bluff-body stabilized burners and swirl burners. The turbulent coherent structures in these burners differ from those in a Bunsen burner. Hence, they can provide different testing con- ditions for the models proposed in this work. Premixed swirl burners are similar to the actual gas turbine combustion chambers and have the poten- tial to reach a very high turbulence intensity. This can provide conditions that are far from the flamelet regime where the flamelet assumption is no longer valid. Hence, simulation of such burners will be a useful test before modelling actual gas turbine combustors which is the ultimate purpose of these models. The current work only uses methane as the fuel. Simulation of burners with different fuels such as propane and biofuels is recommended. Premixed propane flames are more strain sensitive and biofuels can have a slower chemistry time scales. Hence, a stronger coupling exists between turbulence and chemistry which is a good testing condition for CSE. The LES sub-filter scale model used is a simple k-equation model. It is suggested to use a dynamic sub-filter scale LES model to improve the results. The Bayesian inverse problem which is explained in Appendix A has not been implemented yet. Using this approach along with the second- order CMC hypothesis (as explained in Appendix B) might improve the current version of CSE. This improvement depends on using an efficient numerical algorithm. 130 7.2. Future Work Finally, after extending CSE from non-premixed to premixed combus- tion, the next step is to use this model in partially-premixed flames. In order to simulate these flames, the conditional averages will likely have to be doubly conditioned on progress variable and mixture fraction. This will make challenges in solving the inverse problem to find the conditional scalar field. 131 References [1] A Special Report of Working Group III of the Intergovernmental Panel on Climate Change, Montreal, Canada. Carbon Dioxide Capture and Storage Summary for Policymakers, 2005. [2] R. G. Abdel-Gayed and D. Bradley. Combustion regimes and the straining of turbulent premixed flames. Combustion and Flame, 76:213–218, 1989. [3] M. S. Anand and S. B. Pope. Calculations of premixed turbulent flames by PDF methods. Combustion and Flame, 67:127–142, 1987. [4] A. J. Aspden, M. S. Day, and J. B. Bell. Turbulence-flame interac- tions in lean premixed hydrogen: transition to the distributed burning regime. Journal of Fluid mechanics, 680:287–320, 2011. [5] M. Bidi, M.R.H. Nobari, and M.S. Avval. A numerical investigation of turbulent premixed methane-air combustion in a cylindrical chamber. Combustion Science and Technology, 179:1841–1865, 2007. [6] R. W. Bilger. Conditional moment closure for turbulent reacting flow. Physics of Fluids, 5:436–444, 1993. [7] R. W. Bilger. Conditional moment closure modelling and advanced laser measurements. In Takeno T, editor. Turbulence and molec- ular processes in combustion (Sixth Toyota Conference), Amster- dam:Elsevier, pages 267–285, 1993. [8] R. W. Bilger. The role of combustion technology in the 21st century. In T. Echekki and E. Mastorakos, editors, Turbulent Combustion Mod- eling, pages 3–18. Springer, 2012. [9] M. Boger, D. Veynante, H. Boughanem, and A. Trouve. Direct numer- ical simulation analysis of flame surface density concept for large-eddy simulation of turbulent premixed combustion. In 27th Symposium (In- ternational) on Combustion/The combustion Institute, pages 917–925, 1998. 132 References [10] R. Borghi. On the structure of turbulent premixed flames. In C. Bruno and C. Casci, editors, Recent Advances in Aeronautical Science, pages 117–138. Pergamon, London, 1985. [11] J. P. Boris, F. F. Grinstein, E. S. Oran, and R. S. Kolbe. New insights into large eddy simulation. Fluid Dynamics Research, 10:199–228, 1992. [12] C.T. Bowman, R.K. Hanson, W.C. Gardiner, V. Lissianski, M. Fren- klach, M. Goldenberg, G.P. Smith, D.R. Crosley, and D.M. Golden. An optimized detailed chemical reaction mechanism for methane combus- tion and NO formation and reburning. Technical Report GRI-97/0020, Gas Research Institute, Chicago, IL, USA, 1997. [13] D. Bradley, P. H. Gaskell, and X. J. Gu. Application of a Reynolds stress, stretched flamelet, mathematical model to computations of tur- bulent burning velocities and comparison with experiments. Combus- tion and Flame, 96:221–248, 1994. [14] D. Bradley, P. H. Gaskell, X. J. Gu, M. Lawes, and M. J. Scott. Premixed turbulent flame stability and no formation in a lean-burn swirl burner. Combustion and Flame, 115:515–538, 1998. [15] D. Bradley, L. K. Kwa, A. K. C. Lau, M. Missaghi, and S. B. Chin. Laminar flamelet modeling of recirculating premixed methane and propane-air combustion. Combustion and Flame, 71:109–122, 1988. [16] K. N. C. Bray. Turbulent flows with premixed reactants. In Turbulent Reacting Flows, pages 115–183. Topics in Applied Physics, Springer, 1980. [17] K. N. C. Bray. Turbulent transport in flames. Proc. R. Soc. Lond. A, 451:231–256, 1995. [18] K. N. C. Bray, M. Champion, P. A. Libby, and N. Swaminathan. Fi- nite rate chemistry and presumed PDF models for premixed turbulent combustion. Combustion and Flame, 146:665–673, 2006. [19] K. N. C. Bray, P. A. Libby, and J. B. Moss. Unified modelling ap- proach for premixed turbulent combustion - Part I: General formula- tion. Combustion and Flame, 61:87–102, 1985. 133 References [20] W. K. Bushe and H. Steiner. Conditional moment closure for large- eddy simulation of non-premixed turbulent reacting flows. Physics of Fluids, 11:1896–1906, 1999. [21] J. Cannon and B. Shivamoggi. Mathematical and Physical Theory of Turbulence. Chapman & Hall/CRC, 2006. [22] S. M. Cannon, B. S. Brewster, and L. D. Smooth. PDF modeling of lean premixed combustion using in situ tabulated chemistry. Combus- tion and Flame, 119:233–252, 1999. [23] N. Chakraborty, J.W. Rogerson, and N. Swaminathan. A priori as- sessment of closures for scalar dissipation rate transport in turbulent premixed flames using direct numerical simulation. Physics of Fluids, 20:045106, 2008. [24] N. Chakraborty and N. Swaminathan. Influence of the damköhler number on turbulence-scalar interaction in premixed flames. ii. model development. Physics of Fluids, 19:045104, 2007. [25] V. K. Chakravarthy and S. Menon. Large-eddy simulation of turbu- lent premixed flames in the flamelet regime. Combustion Science and Technology, 162:175–222, 2001. [26] Y. C. Chen. A general procedure for constructing reduced reaction mechanisms with given independent relations. Combustion Science and Technology, 57:89–94, 1988. [27] Y. C. Chen and R. W. Bilger. Experimental investigation of three- dimensional flame front structure in premixed turbulent combustion-I: hydrocarbon/air bunsen flames. Combustion and Flame, 131:400–435, 2002. [28] Y. C. Chen, N. Peters, G. A. Schneemann, N. Wruck, U. Renz, and M. S. Mansour. The detailed flame structure of highly stretched turbu- lent premixed methane-air flames. Combustion and Flame, 107:223– 244, 1996. [29] F. C. Christo, A. R. Masri, and E. M. Nobet. An integrated PDF/neural network approach for simulating turbulent reacting sys- tems. Proceedings of the Combustion Institute, 26:43–48, 1996. 134 References [30] M. J. Cleary, A. Y. Klimenko, J. Janicka, and M. Pfitzner. A sparse- lagrangian multiple mapping conditioning model for turbulent diffu- sion flames. Proceedings of the Combustion Institute, 32:1499–1507, 2009. [31] F. Colin, D. Veynante, and T. Poinsot. A thickened flame model for large-eddy simulation of premixed turbulent combustion. Physics of Fluids, 12:1843–1863, 2000. [32] P. J. Colucci, F. A. Jaberi, P. Givi, and S. B. Pope. Filtered density function for large-eddy simulation of turbulent reacting flows. Physics of Fluids, 10:499–515, 1998. [33] G. Damköhler. Der einfluss der turbulenz auf die flammengeschwindi- genkeit in gasgemischen. Z. Electrochem, 46:601–652, 1940. [34] L. P. H. de Goey and J. H. M. ten Thije Boonkkamp. A flamelet description of premixed laminar flames and the relation with flame stretch. Combustion and Flame, 119:253–271, 1999. [35] L. Duchamp de Lageneste and H. Pitsch. A level-set approach to large-eddy simulation of premixed turbulent combustion. Center for Turbulence Research Annual Research Briefs, pages 105–116, 2000. [36] G. Dixon-Lewis. Flame structure and flame reaction kinetics, II: Transport phenomena in multi-component systems. Proc. R. Soc. A, 307:111–135, 1968. [37] P. Domingo and K.N.C. Bray. Laminar flamelet expressions for pres- sure fluctuation terms in second moment models of premixed turbulent combustion. Combustion and Flame, 121:555–574, 2000. [38] P. Domingo, Luc Vervisch, S. Payet, and R. Hauguel. DNS of a pre- mixed turbulent V-flame and LES of a ducted flame using a FSD-PDF subgrid scale closure with FPI-tabulated chemistry. Combustion and Flame, 143:566–586, 2005. [39] J. F. Driscoll. Turbulent premixed combustion: Flamelet structure and its effect on turbulent burning velocities. Progress in Energy and Combustion Science, 34:91–134, 2008. [40] M. J. Dunn, A. R. Masri, R. W. Bilger, R. S. Barlow, and G. H. Wang. The compositional structure of highly turbulent piloted pre- 135 References mixed flames issuing into a hot coflow. Proceedings of the Combustion Institute, 32:1779–1786, 2009. [41] B. Fiorina, R. Baron, O. Gicquel, D. Thevenin, S. Carpentier, and N. Darabiha. Modelling non-adiabatic partially premixed flames us- ing flame-prolongation of ILDM. Combustion Theory and Modelling, 7:449–470, 2003. [42] B. Fiorina, O. Gicquel, L. Vervisch, S. Carpentier, and N. Dara- biha. Premixed turbulent combustion modeling using tabulated de- tailed chemistry and PDF. Proceedings of the Combustion Institute, 30:867–874, 2005. [43] B. Fiorina, R. Vicquelin, P. Auzillon, N. Darabiha, O. Gicquel, and D. Veynante. A filtered tabulated chemistry model for LES of pre- mixed combustion. Combustion and Flame, 157:465–475, 2010. [44] C. Fureby and F. F. Gristein. Monotonically integrated large-eddy simulation of free shear flows. AIAA Journal, 37:544–556, 1999. [45] J. Galpin, A. Naudin, L. Vervisch, C. Angelberger, O. Colin, and P. Domingo. Large-eddy simulation of a fuel-lean premixed turbulent swril burner. Combustion and Flame, 155:247–266, 2008. [46] X. Gao and C. P. T. Groth. A parallel solution - adaptive method for three-dimensional turbulent non-premixed combusting flows. Journal of Computational Physics, 229:3250–3275, 2010. [47] A. Garmory and E. Mastorakos. Capturing localized extinction in San- dia flame F with LES-CMC. Proceedings of the Combustion Institute, 33:1673–1680, 2011. [48] E. Garnier, N. Adams, and P. Sagaut. Large-Eddy Simulation for compressible flows. Springer, 2009. [49] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic sub-grid scale eddy viscosit model. Physics of Fluids, 3:1760–1765, 1991. [50] S. Ghosal. Analysis and control of errors in the numerical simulation of turbulence. In D. Drikakis and B. J. Geurts, editors, Turbulent Flow Compuations, pages 101–140. Kluwer Academic Publishers, 2004. 136 References [51] S. Ghosal and P. Moin. The basic equations fot the large-eddy simula- tion of turbulent flows in complex geometry. Journal of Computational Physics, 118:24–37, 1995. [52] O. Gicquel, N. Darabiha, and D. Thévenin. Laminar premixed hy- drogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion. Proceedings of the Combustion In- stitute, 28:1901–1908, 2000. [53] P. Givi. Model free simulations of turbulent reactive flows. Progress in Energy and Combustion Science, 15:1–107, 1989. [54] G. H. Golub, M. Heath, and G. Wahba. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21:215–223, 1979. [55] D. G. Goodwin. Cantera: Object-oriented software for reacting flows. http://www.cantera.org. [56] D. A. Goussis and U. Maas. Model reduction for combustion chem- istry. In T. Echekki and E. Mastorakos, editors, Turbulent Combustion Modelling, pages 193–220. Springer, 2011. [57] R. W. Grout, W. K. Bushe, and C. Blair. Predicting the ignition delay of turbulent methane jets using conditional source-term estimation. Combustion Theory and Modelling, 11:1009–1028, 2007. [58] A. Gruber, R. Sankaran, E. R. Hawkes, and J. H. Chen. Turbulent flame-wall interaction: a direct numerical study. Journal of Fluid Mechanics, 658:5–32, 2010. [59] J. Hadamard. Lectures on Couchy’s Problem in Linear Partial Differ- ential Equations. Yale University Press, 1923. [60] P. C. Hansen. The truncated SVD as a method for regularization. BIT Numerical Mathematics, 27:534–553, 1990. [61] P. C. Hansen. Numerical tools for analysis and solution of Fredholm integral equations of the first kind. Inverse Problems, 8:849–872, 1992. [62] P. C. Hansen. Rank-Deficient and Discrete Ill-Posed Problems. SIAM Press, 1998. 137 References [63] E. R. Hawkes and R. S. Cant. A flame surface density approach to large-eddy simulation of premixed turbulent combustion. Proceedings of the Combustion Institiute, 28:51–58, 2000. [64] E. R. Hawkes and J. H. Chen. Comparison of direct numerical simula- tion of lean premixed methane-air flames with strained laminar flame calculations. Combustion and Flame, 144:112–125, 2006. [65] D. C. Haworth. Progress in probability density function methods for turbulent reacting flows. Progress in Energy and Combustion Science, 36:168–259, 2010. [66] F. E. Hernandez-Perez, F. T. C. Yuen, C. P. T. Groth, and O. L. Gulder. LES of a laboratory-scale turbulent premixed bunsen flame using FSD, PCM-FPI and thickened flame models. Proceedings of the Combustion Institiute, 33:1365–1371, 2011. [67] M. Herrmann. Numerical simulation of turbulent bunsen flames with a level set flamelet model. Combustion and Flame, 145:357–375, 2006. [68] J. Huang and W. K. Bushe. Simulation of an igniting methane jet using conditional source-term estimation with a trajectory generated low-dimensional manifold. Combustion Theory and Modeling, 11:977– 1008, 2007. [69] Y. Huang, H. Sung, S. Hsieh, and V. Yang. Large-eddy simulation of combustion dynamics of lean-premixed swirl-stabilized combustor. Journal of Propulsion and Power, 19:782–794, 2003. [70] Jérôme Idier. Bayesian Approach to Inverse Problems. ISTE Ltd and John Wiley & Sons Inc, 2008. [71] H. G. Im, T. S. Lund, and J. H. Ferziger. Large-eddy simulation of turbulent front propagation with dynamic subgrid models. Physics of Fluids, 9:3826–3833, 1997. [72] F. A. Lammers J. A. van Oijen and L. P. H. de Goey. Modeling of com- plex premixed burner systems by using flamelet generated manifolds. Combustion and Flame, 127:2124–2134, 2001. [73] S. James, M. S. Anand, M. K. Razdan, and S. B. Pope. In situ detailed chemistry calculations in combustor flow analyses. Journal of Engineering for Gas Turbines and Power - Transactions of the ASME, 123:747–756, 2001. 138 References [74] N. Jarrin, B. Benhamadouche, D. Laurence, and R. Prosser. A synthetic-eddy-method for generating inflow conditions for large-eddy simulation. International Journal of Heat and Fluid Flow, 27:585–593, 2006. [75] H. Jasak, H.G. Weller, and A.D. Gosman. High resolution NVD dif- ferencing schemes for arbitrary unstructured meshes. International Journal for Numerical Methods in Fluids, 31:431–449, 1999. [76] B. Jin, R. Grout, and W. K. Bushe. Conditional source-term estima- tion as a method for chemical closure in premixed turbulent reacting flow. Flow, Turbulence and Combustion, available online, 2008. [77] J. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems. Springer, 2004. [78] J. Keck and D. Gillespie. Rate-controlled partial-equilibrium method for treating reacting gas mixtures. Combustion and Flame, 17:237– 241, 1971. [79] Robert J. Kee, Micheal E. Coltrin, and Peter Glarborg. Chemically Reacting Flow: Theory and Practice. John Wiley & Sons Inc, 2003. [80] A. R. Kerstein, W. Ashurt, and F. A. Williams. Field equation for interface propagation in an unsteady homogeneous flow field. Physical Review A, 37:2728–2731, 1988. [81] W. W. Kim and S. Menon. Numerical modeling of turbulent premixed flames in the thin-reaction-zones regime. Combustion Science and Technology, 160:119–150, 2000. [82] W. W. Kim, S. Menon, and H. C. Mongia. Large-eddy simulation of a gas turbine combustor flow. Combustion Science and Technology, 143:25–62, 1999. [83] A. Y. Klimenko and R. W. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25:595–687, 1999. [84] A. Y. Klimenko and S. B. Pope. The modeling of turbulent reac- tive flows based on multiple mapping conditioning. Physics of Fluids, 15:1907–1925, 2003. 139 References [85] E. Knudsen and H. Pitsch. A dynamic model for the turbulent burning velocity for large-eddy simulation of premixed combustion. Combus- tion and Flame, 154:740–760, 2008. [86] H. Kolla, , and N. Swaminathan. Strained flamelets for turbulent premixed flames, I: formulation and planar flame results. Combustion and Flame, 157:943–954, 2010. [87] H. Kolla, J.W. Rogerson, N. Chakraborty, and N. Swaminathan. Scalar dissipation modelling and its validation. Combustion Science and Technology, 181:518–535, 2009. [88] Kenneth K. Kuo. Principles of Combustion. John Wiley & Sons Inc, 2005. [89] Kenneth K. Kuo and Ragini Acharya. Fundamentals of Turbulent and Multiphase Combustion. John Wiley & Sons Inc, 2012. [90] S. H. Lam and D. A. Goussis. Understanding complex chemical ki- netics with computational singular perturbations. Proceedings of the Combustion Institute, 22:931–941, 1988. [91] C. K. Law. Combustion Physics. Cambridge University Press, 2006. [92] C. K. Law. On the applicability of direct relation graph to the reduc- tion of reaction mechanisms. Combustion and Flame, 146:472–483, 2006. [93] C. L. Lawson and R. J. Hanson. Solving least square problems. Prentice-Hall, 1974. [94] H Le, P Moin, and J Kim. Direct numerical simulation of turbulent flow over a backward-facing step. Journal of fluid mechanics, 330:349– 374, 1997. [95] G. Lecocq, S. Richard, O. Colin, and L. Vervisch. Hybrid presumed pdf and flame surface density appoaches for large-eddy simulation of premixed turbulent combustion. part 1: Formalism and simulation of a quasi-steady burner. Combustion and Flame, 158:1201–1214, 2011. [96] P.A. Libby and F.A. Williams. Turbulent Reacting Flows. Academic Press, London, 1994. 140 References [97] R. P. Lindstedt and E. M. Vaos. Transport PDF modeling of high- Reynolds-number premixed turbulent flames. Combustion and Flame, 145:495–511, 2006. [98] T. F. Lu and C. K. Law. Toward accommodating realistic chemistry in large-scale computations. Progress in Energy and Combustion Science, 35:192–215, 2009. [99] U. Maas and S. B. Pope. Simplifying chemical kinetics: Intrinsic low- dimensional manifolds in composition space. Combustion and Flame, 88:239–264, 1992. [100] M. S. Mansour, N. Peters, and Y. C. Chen. Investigation of scalar mixing in the thin reaction zones regime using a simultaneous CH- LIF/Rayleigh laser technique. In 27th Symposium (International) on Combustion/The combustion Institute, pages 767–773, 1998. [101] T. Mantel and R. Borghi. A new model of premixed wrinkled flame propagation based on scalar dissipation equation. Combustion and Flame, 96:443–457, 1994. [102] S. M. Martin, J. C. Kramlich, G. Kosaly, and J. J. Riley. The premixed conditional moment closure method applied to idealized lean premixed gas turbine combutors. Journal of Engineering for Gas Turbines and Power - Tansactions of ASME, 125:895–900, 2003. [103] F. Mauss and N. Peters. Reduced kinetic mechanisms for premixed methane-air flames. In N. Peters and B. Rogg, editors, Lecture Notes in Physics: Reduced Kinetics Mechanisms for Application in Combustion Systems, pages 58–75. Springer-Verlag, 1993. [104] S. McAllister, J. Y. Chen, and A. C. Fernandez-Pello. Fundamentals of Combustion Processes. Springer, 2011. [105] G. P. Mctaggart-Cowan, N. Wu, B. Jin, S. N. Rogak, M. H. Davy, and W. K. Bushe. Effects of fuel compostion on high-pressure non- premixed natural gas combustion. Combustion Science and Technol- ogy, 181:397–416, 2009. [106] P. Moin and K. Mahesh. Direct numerical simulation: A tool in turbu- lence research. Annual Review of Fluid Mechanics, 30:539–578, 1998. 141 References [107] A. Mura and R. Borghi. Towards an extended scalar dissipation equation for turbulent premixed combustion. Combustion and Flame, 133:193–196, 2003. [108] J. R. Nanduri, D. R. Parsons, S. L. Yilmaz, I. B. Celik, and P. A. Strakey. Assessment of RANS-based turbulent combustion models for prediction of emissions from lean premixed combustion of methane. Combustion Science and Technology, 182:794–821, 2010. [109] S. Navarro-Martinez and A. Kronenburg. LES-CMC simulation of a turbulent bluff-body flame. Proceedings of the Combustion Institute, 31:1721–1728, 2007. [110] S. Navarro-Martinez, A. Kronenburg, and F. Di Mare. Conditional moment closure for large eddy simulations. Flow Turbulence and Com- bustion, 75:245–247, 2005. [111] OpenFOAM. http://www.opencfd.co.uk/openfoam/, 2007. [112] F. O’Young and R. W. Bilger. Scalar gradient and related quantities in turbulent premixed flames. Combustion and Flame, 109:682–700, 1997. [113] N. Peters. Laminar diffusion flamelet models in non-premixed turbu- lent combustion. Progress in Energy and Combustion Science, 10:319– 339, 1984. [114] N. Peters. Numerical and asymptotic analysis of systematically re- duced reaction schemes for hydrocarbon flames. In R. Glowinski, B. Larrouturou, and R. Temam, editors, Numerical Simulation of Combustion Phenomena, pages 90–109. Springer-Verlag, 1985. [115] N. Peters. Laminar flamelet concepts in turbulent combustion. In 21st Symposium (International) on Combustion/The combustion Institute, pages 1231–1250, 1986. [116] N. Peters. The turbulent burning velocity for large-scale and small- scale turbulence. Journal of Fluid mechanics, 384:107–132, 1999. [117] N. Peters. Turbulent Combustion. Cambridge University Press, 2000. [118] H. Phillips. Flame in a buoyant methane layer. In 10th Symposium (International) on Combustion/The combustion Institute, pages 1277– 1250, 1965. 142 References [119] U. Piomelli and J. R. Chasnov. Large-eddy simulations: Theory and applications. In M. Hallback, D. S. Henningson, A. V. Johansson, and P. H. Alfredsson, editors, Turbulence and Transition Modelling, pages 269–336. Kluwer Academic Publishers, 1996. [120] H. Pitsch. FlameMaster, A C++ computer program for 0D combus- tion and 1D laminar flame calculations. [121] H. Pitsch. Large-eddy simulation of turbulent combustion. Annual Review of Fluid Mechanics, 38:453–482, 2006. [122] H. Pitsch and L. Duchamp De Lageneste. Large-eddy simulation of premixed turbulent combustion using a level-set approach. Proceedings of the Combustion Institute, 29:2001–2008, 2002. [123] Thierry Poinsot and Denis Veynanate. Theoretical and numerical com- bustion. R.T. Edwards, 2nd edition, 2005. [124] S. B. Pope. An explanation of the turbulent round-jet/plane-jet anomaly. AIAA Journal, 16:279–281, 1978. [125] S. B. Pope. PDF methods for turbulent reactive flows. Progress in Energy and Combustion Science, 11:119–192, 1985. [126] S. B. Pope. Computations of turbulent combustion: Progress and challenges. In 23th Symposium (International) on Combustion/The Combustion Institute, pages 591–612, 1991. [127] S. B. Pope. Computationally efficient implementation of combustion chemistry using in situ adaption tabulation. Combustion Theory and Modelling, 1:41–63, 1997. [128] S. B. Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New Journal of Physics, 6, 2004. [129] S. B. Pope and U. Maas. Simplifying chemical kinetics: Intrinsic low- dimensional manifolds. Technical Report 11, FDA, 1993. [130] Stephen B. Pope. Turbulent flows. Cambridge University Press, 2000. [131] R.O.S. Prasad and J.P. Gore. An evaluation of flame surface density models for turbulent premixed jet flames. Combustion and Flame, 116:1–14, 1999. 143 References [132] Published by Ministry of Justice, http://laws-lois.justice.gc.ca. On- Road Vehicle and Engine Emissions Regulations, 2012. [133] Z. Ren, S. B. Pope, A. Vladimirsky, and J. M. Guckenheimer. The invariant constrained equilibrium edge preimage curve method for the dimensional reduction of chemical kinetics. The Journal of Chemical Physics, 124:114111, 2006. [134] G. Ribert, M. Champion, O. Gicquel, N. Darabiha, and D. Veynante. Modelling non-adiabatic turbulent premixed reactive flows including tabulated chemistry. Combustion and Flame, 141:271–280, 2005. [135] R. S. Rogallo. Numerical experiments in homogeneous turbulence. Technical Report 81315, NASA Technical Memorandum, 1981. [136] Pierre Sagaut. Large-Eddy Simulation for Incompressible flows. Springer, 2006. [137] M. M. Salehi and W. K. Bushe. Presumed PDF modeling for RANS simulation of turbulent premixed flames. Combustion Theory and Modelling, 14:381–403, 2010. [138] H. Schmidt and U. Schumann. Coherent structure of the convective boundary layer derived from large-eddy simulations. Journal of Fluid Mechanics, 200:511–562, 1989. [139] E. Schneider, A. Sadiki, and J. Janicka. Modeling and 3D-simulation of the kinetic effects in post-flame region of turbulent premixed flames based on the g-equation approach. Flow, Turbulence and Combustion, 75:191–216, 2005. [140] M. R. H. Sheikhi, T. G. Drozda, P. Givi, F. A. Jaberi, and S. B. Pope. Large-eddy simulation of a turbulent nonpremixed piloted methane jet flame (Sandia flame D). Proceedings of the Combustion Institute, 30:549–556, 2005. [141] J. Smagorinsky. General circulation experiments with the primitive equations. Monthly Weather Review, 91:99–164, 1963. [142] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eite- neer, M. Goldenberg, C. T. Bowman, R. K. Hanson, S. Song, W. C. Gardiner, V. V. Jr. Lissianski, and Z. Qin. GRI-MECH 3.0. Technical report, http://www.me.berkeley.edu/gri mech. 144 References [143] N. S. A. Smith. Development of the conditional moment closure method for modeling turbulent combustion. PhD thesis, University of Sydney, Australia, 1994. [144] D. B. Spalding. Mixing and chemical reaction in steady confined tur- bulent flames. Proceedings of the Combustion Institute, 13:649–657, 1971. [145] H. Steiner and W. K. Bushe. Large-eddy simulation of a turbulent reacting jet with conditional source-term estimation. Physics of Fluids, 13:754–769, 2001. [146] M. Stöllinger and S. Heinz. PDF modeling and simulation of premixed turbulent combustion. Monte Carlo Methods Appl, 14:343–377, 2008. [147] N. Swaminathan and R. W. Bilger. Analysis of conditional moment closure for turbulent premixed flame. Combustion Theory and Mod- elling, 5:241–260, 2001. [148] N. Swaminathan and K.N.C Bray. Effect of dilatation on scalar dissi- pation in turbulent premixed flames. Combustion and Flame, 143:549– 565, 2005. [149] N. Swaminathan and R.W. Grout. Interaction of turbulence and scalar fields in premixed flames. Physics of Fluids, 18:045102, 2006. [150] A. Tikhonov and V Arsenin. Solution of incorrectly formulated prob- lems and the regularization method. Soviet Mathematics Doklady, 4:1035–1038, 1963. [151] A. Triantafyllidis, E. Mastorakos, and RLGM Eggels. Large-eddy sim- ulations of forced ignition of a non-premixed bluff-body methane flame with conditional moment closure. Combustion and Flame, 156:2328– 2345, 2009. [152] S. Tullis and R.S. Cant. Scalar transport modelling in large-eddy sim- ulation of turbulent premixed flames. Proceedings of the Combustion Institute, 29:2097–2104, 2002. [153] UNFCC: United Nations Framework Convention on Climate Change, www.unfccc.int. 1992. [154] J. A. van Oijen and L. P. H. de Goey. Modelling of premixed laminar flames using flamelet-generated manifolds. Combustion Science and Technology, 161:113–137, 2000. 145 References [155] O. V. Vasilyev. A general class of commutative filters for LES in complex geometries. Journal of Computational Physics, 146:82–104, 1998. [156] H. Van Der Ven. A family of large-eddy simulation filters with nonuni- form filter widths. Physics of Fluids, 7:1171–1172, 1995. [157] L. Vervisch, R. Hauguel, P. Domingo, and M. Rullaud. Three facets of turbulent combustion modelling: DNS of premixed v-flame, LES of lifted non-premixed flame and RANS of jet-flame. Journal of Turbu- lence, 5:1–36, 2004. [158] D. Veynante, A. Trouve, K. N. C. Bray, and T. Mantel. Gradient and counter-gradient scalar transport in turbulent premixed flames. Journal of Fluid Mechanics, 332:263–293, 1997. [159] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress in Energy and Combustion Science, 28:193–266, 2002. [160] W. Vicente, M. Salinas, Esteban Barrios, and C. Dopazo. PDF mod- eling of co and no formation in lean premixed methane flames. Com- bustion Science and Technology, 176:585–601, 2004. [161] K. Vogiatzakim, A. Kronenburg, M. J. Cleary, and J. H. Kent. Multi- ple mapping conditioning of turbulent jet diffusion flames. Proceedings of the Combustion Institute, 32:1679–1685, 2009. [162] M. Wang, A. Frisque, J. Huang, and W. K. Bushe. Trajectory gen- rated low-dimensional manifolds generated using the stochastic parti- cle method. Combustion Theory and Modeling, 12:249–267, 2008. [163] M. Wang, J. Huang, and W. K. Bushe. Simulation of a turbulent non-premixed flame using conditional source-term estimation with a trajectory generated low-dimensional manifold. Proceedings of The Combustion Institute, 31:1701–1709, 2007. [164] J. Warnatz, U. Maas, and R.W. Dibble. Combustion. Springer, 3 edition, 2000. [165] F. A. Williams. Turbulent combustion. In J. D. Buckmaster, editor, Mathematics of Combustion, pages 97–132. SIAM, 1985. [166] World Energy Outlook, International Energy Agency, Paris, France. 2008. 146 References [167] S. L. Yilmaz, M. B. Nik, P. Givi, and P. A. Strakey. Scalar filtered density function for large-eddy simulation of a bunsen burner. Journal of Propulsion and Power, 26:84–93, 2010. [168] A. Yoshizawa and K. Horiuti. A statistically-derived subgrid scale kinetic energy model for the large-eddy simulation of turbulent flows. Journal of the Physical Society of Japan, 54:2834–2839, 1985. [169] F. T. C. Yuen and O. L. Gulder. Premixed turbulent flame structure investigation by Rayleigh scattering in the thin reaction zone regime. Proceedings of the Combustion Institue, 32:1747–1754, 2009. [170] S. Zhang and C.J. Rutland. Premixed flame effects on turbulence and pressure-related terms. Combustion and Flame, 102:447–461, 1995. [171] V. L. Zimont and F. Biagioli. Gradient, counter-gradient transport and their transition in turbulent premixed flames. Combustion Theory and Modelling, 6:79–101, 2002. 147 Appendix A: Bayesian Inverse Problem An inverse problem can be formulated using deterministic optimization ap- proaches or through statistical Bayesian methods [70]. The solution of a deterministic approach is one point in the parametric solution space. How- ever, a Bayesian interface formulates the solution as an ensemble of random points in solution space with different probabilities. Therefore, a Bayesian interface can give a more comprehensive description of inverse problems. The basis of a Bayesian interface is Bayes’ theorem: p(~α|~b) = p( ~b|~α)p(~α) p(~b) (7.1) where ~α is the unknown in discrete form of the inverse problem i.e., A~α = ~b. In this equation, p(~α|~b) is the probability of the solution ~α being correct given the observation ~b which is called posterior probability; p(~b|~α) is the probability of the observation ~b actually happening assuming that the solu- tion is known. p(~b|~α) is also known as likelihood; p(~α) is called prior which contains the initial knowledge of the statistics of the solution; p(~b) acts as a normlizing term and can be calculated knowing the prior and the likelihood: p(~b) = ∫ p(~b|~α)p(~α)d~α. The information about the solution of an inverse problem comes from the linear operation of A~α = ~b. In reality, this linear system is contaminated with error. Assuming a Gaussian-distributed error, the likelihood can be stated as: p(~b|~α) = (2piσ2)−m/2exp ( − 1 2σ2 ||A~α−~b||22 ) (7.2) where || ||2 denotes the L2-norm of a vector. Based on this equation, the solution that minimizes the least-square error maximizes the likelihood. However, minimizing the least-square error is ill-posed and depending on the nature of ill-posedness, the least-square solution is either not unique or unstable. Therefore, the information regarding the likelihood needs to be 148 Appendix A: Bayesian Inverse Problem augmented with the prior knowledge of the statistics of the solution. As discussed in chapter 4, a laminar unstrained premixed flame can provide a priori knowledge of the solution i.e., ~α0. This information can be used to define the following prior: p(~α) = (2piγ2)−m/2exp ( − 1 2γ2 ||~α− ~α0||22 ) (7.3) This form of prior is also called a conjugate prior because it is mathemat- ically similar to the likelihood function; hence, it simplifies the analytical calculations. Combining equation 7.2 with 7.3 using the Bayes’ rule gives the posterior: p(~α|~b) ∝ exp ( −||A~α−~b||22 − λ2||~α− ~α0||22 ) (7.4) where λ2 ≡ σ2/γ2. The regularization parameters can also be selected in a rigorous probabilistic approach. The posterior distribution is the basis for statistical inference of the solution of the inverse problem. Several different estimators can be used to obtain the solution from the posterior distribu- tion [77]. The minimum variance Bayes estimator and the maximum a posterior (MAP) estimator are often used. The minimum variance Bayes estimator is the first moment of the solu- tion based on the posterior: ~̂α = ∫ ~αp(~α|~b)d~α (7.5) this equation involves the integration over a multi-dimensional function which is not usually analytically possible. It is common to use Monte Carlo techniques for sampling from the high-dimensional solution space to calcu- late the integration. The MAP solution is the most probable solution which is evaluated through an optimization problem: ~̂α = arg max ~α p(~α|~b) (7.6) In general, this optimization problem can be solved using a gradient de- scent method. For the simple posterior of equation 7.4, the MAP solution simplifies to the following mimization problem: ~̂α = arg min ~α ( ||A~α−~b||22 + λ2||~α− ~α0||22 ) (7.7) 149 Appendix A: Bayesian Inverse Problem this is equivalent to the Tikhonov solution of equation 4.9. The advantage of the Bayesian framework is that one can estimate the uncertainty in the soluion of the inverse problem by using the posterior: δ~̂α = ∫ (~α− ~̂α)2p(~α|~b)d~α (7.8) Now the solution of the inverse problem can be stated as ~̂α± δ~̂α. 150 Appendix B: Second-order CMC The TGLDM chemistry for premixed combustion is a function of at least two progress variables: ω̇k(T, ρ, Yk) ≈ ω̇k(c1, c2) (7.9) where these two progress variables are c1 = YCO2 and c2 = YH2O for methane-air combustion. As mentioned in chapter 6, the first-order CMC hypothesis for this chemistry reads: 〈ω̇k|c∗1〉 ≈ ω̇k(c∗1, 〈c2|c∗1〉) (7.10) This approximation can be further improved by taking the effect of higher conditional moments into account. This can be done using a Taylor series expansion of the reaction source terms [83]: ω̇k(c∗1, c ∗ 2) = ω̇k(c ∗ 1, c2|c∗1) + c′′2 ∂ω̇k ∂c2 ∣∣∣∣∣ c∗2=〈c2|c∗1〉 + 1 2 c′′22 ∂2ω̇k ∂2c2 ∣∣∣∣∣ c∗2=〈c2|c∗1〉 + TE (7.11) where c′′2 ≡ c∗2 − 〈c2|c∗1〉. If the above equation is filtered conditioned on c∗1, the following equation is obtained: 〈ω̇k|c∗1〉 ≈ ω̇k(c∗1, 〈c2|c∗1〉) + 1 2 〈c′′22 |c∗1〉 ∂2ω̇k ∂2c2 ∣∣∣∣∣ c∗2=〈c2|c∗1〉 (7.12) where 〈c2|c∗1〉 is the conditionally-filtered value of c2 conditioned on c1 and 〈c′′22 |c∗1〉 is the conditional variance of c2. Conditional Source-term Estima- tion model can be used to estimate 〈c2|c∗1〉 as explained in chapters 4 and 6. In CSE, an integral equation is solved assuming that 〈c2|c∗1〉 is constant in an ensemble. This integral equation in discrete form reduces to the following linear system: ~b = A~α (7.13) 151 Appendix B: Second-order CMC where ~α is the discrete form of 〈c2|c∗1〉. As explained in Appendix A, a Bayesian interface can be used to not only find the solution vector ~α, but also the uncertainty in the vector δ~α using equation 7.8. This uncertainty can be used to approximate 〈c′′22 |c∗1〉 in equation 7.12. In this way, the homogenity assumption of CSE can be relaxed. In other words, when the conditional average is not homogeneous in an ensemble, assuming a constant conditional average and using the first-order CMC hypothesis will not be a good approximation. A more accurate approximation can be obtained using second-moment closure. This is of particular importance in complex geometries where there is unlikely to be a direction of homogenity. 152
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical simulation of turbulent premixed flames with...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical simulation of turbulent premixed flames with conditional source-term estimation Salehi, Mohammad Mahdi 2012
pdf
Page Metadata
Item Metadata
Title | Numerical simulation of turbulent premixed flames with conditional source-term estimation |
Creator |
Salehi, Mohammad Mahdi |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Conditional Source-term Estimation (CSE) is a closure model for turbulence-chemistry interactions. This model is based on the conditional moment closure hypothesis for the chemical reaction source terms. The conditional scalar field is estimated by solving an integral equation using inverse methods. CSE was originally developed for - and has been used extensively in - non-premixed combustion. This work is the first application of this combustion model to predictive simulations of turbulent premixed flames. The underlying inverse problem is diagnosed with rigorous mathematical tools. CSE is coupled with a Trajectory Generated Low-Dimensional Manifold (TGLDM) model for chemistry. The CSE-TGLDM combustion model is used with both Reynolds-Averaged Navier-Stokes (RANS) and Large-Eddy Simulation (LES) turbulence models to simulate two different turbulent premixed flames. Also in this work, the Presumed Conditional Moment (PCM) turbulent combustion model is employed. This is a simple flamelet model which is used with the Flame Prolongation of ILDM (FPI) chemistry reduction technique. The PCM-FPI approach requires a presumption for the shape of the probability density function of reaction progress variable. Two shapes have been examined: the widely used beta-function and the Modified Laminar Flamelet PDF (MLF-PDF). This model is used in both RANS and large-eddy simulation of a turbulent premixed Bunsen burner. Radial distributions of the calculated temperature field, axial velocity and chemical species mass fraction have been compared with experimental data. This comparison shows that using the MLF-PDF leads to predictions that are similar, and often superior to those obtained using the beta-PDF. Given that the new PDF is based on the actual chemistry - as opposed to the ad hoc nature of the beta-PDF - these results suggest that it is a better choice for the statistical description of the reaction progress variable. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-07-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution 3.0 Unported |
DOI | 10.14288/1.0072895 |
URI | http://hdl.handle.net/2429/42775 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_salehi_mohammad.pdf [ 4.9MB ]
- Metadata
- JSON: 24-1.0072895.json
- JSON-LD: 24-1.0072895-ld.json
- RDF/XML (Pretty): 24-1.0072895-rdf.xml
- RDF/JSON: 24-1.0072895-rdf.json
- Turtle: 24-1.0072895-turtle.txt
- N-Triples: 24-1.0072895-rdf-ntriples.txt
- Original Record: 24-1.0072895-source.json
- Full Text
- 24-1.0072895-fulltext.txt
- Citation
- 24-1.0072895.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072895/manifest