Abstract ii Abstract 'Implementation of Conditional Source-term Estimation for the Prediction of Methane Ignition' The Laminar Flamelet Decomposition approach to Conditional Source-term Estimation is a recent development in combustion modeling that approximates the conditional averages of scalars as linear combinations of basis functions obtained from a library of unsteady solutions to the classical flamelet equations. Then an integral equation for a single scalar is inverted and solved for the weighting function of the basis functions. This weighting function is then used to obtain closure for the chemical source terms in the transport equations. Many parameters affect the solution of the Conditional Source-term Estimation method including the solution technique for inverting the integral equation; the number of solutions in theflameletlibrary; and the scalar used for the integral equation. The first part of this research investigated these parameters to optimize the solution for the Conditional Source-term Estimation. This is accomplished in the a priori sense by comparison to Direct Numerical Simulations of a simple turbulent flame. The second part of this work implemented Conditional Source-term Estimation into the commercially available multi-dimensional combustion software, KIVA-II by replacing the existing ignition mechanism. The model was used to predict the ignition of methane for different initial temperatures and compared to experimental results. The predicted ignition delay times were shorter than the experimental results but exhibited similar trends. Contents iii Contents Abstract ii Contents iii List of Tables vi List of Figures Nomenclature Acknowledgements 1 Introduction v i i xi xiv 1 1.1 Introduction 1 1.2 Objective 2 1.3 Thesis Outline 2 2 Background 4 2.1 Introduction 4 2.2 Combustion 4 2.3 Conservation Equations 5 2.3.1 Conservation of Mass 5 2.3.2 Conservation of Momentum 6 2.3.3 Species Balance 2.3.4 Conservation of Energy 7 2.3.5 Chemical Kinetics 8 2.4 Turbulence 2.4.1 Averaging Concepts 7 9 10 Contents 2.4.2 Turbulence Simulation 2.5 Turbulent Combustion Models . . 3 iv 12 '. 14 2.5.1 Mixture Fraction 15 2.5.2 Mixture Fraction Variance 16 2.5.3 Scalar Dissipation Rate 17 2.5.4 Probability Density Function 18 2.5.5 Fast Chemistry 20 2.5.6 Laminar Diffusion Flamelets 21 2.5.7 Conditional Moment Closure 26 2.5.8 Conditional Source-term Estimation 28 2.6 Summary 34 A priori Testing 35 3.1 Introduction 35 3.2 Formulation 35 3.3 Methodology 38 3.4 Linear Regularization Parameters 40 3.4.1 Results of Linear Regularization Parameters 3.5 Optimum Number of Flamelets 3.5.1 Results of Optimum Number of Flamelets 3.6 Optimum Scalar for Decomposition 3.6.1 Results of Optimum Scalar for Decomposition 3.7 Conclusions and Recommendations . . . 41 45 46 51 53 60 4 Implementation 62 4.1 Introduction 62 4.2 KIVA-II - Multi-Dimensional Model for Reactive Flows . . . 62 4.2.1 Jet Simulations with KIVA 4.3 Implementation of Conditional Source-term Estimation . . . 64 65 4.3.1 Numerical Model 65 4.3.2 Unsteady Flamelet Libraries 67 4.3.3 Shock Tube Simulations of Methane Ignition 69 Contents 4.4 Results and Discussion 72 4.4.1 Shock Tube Simulations of Methane Ignition 4.4.2 Temperature Dependence of Methane Ignition 76 4.4.3 Location of Ignition 81 Conclusions and Recommendations 84 Conclusions and Recommendations 87 4.5 5 v . . . . . 5.1 A priori Testing of Conditional Source-term Estimation . . . 5.2 Implementation of Conditional Source-term Estimation 5.3 Recommendations . . . 72 87 88 90 Bibliography 91 A A d d i t i o n a l Converged Statistics 96 B 97 A d d i t i o n a l Scalar Combinations Comparisions C V a l i d a t i o n o f M i x t u r e F r a c t i o n a n d Variance 99 D G r i d Convergence 105 E T i m e E v o l u t i o n of Temperature Contours 107 F Temperature, M i x t u r e Fraction and Variance 110 G Reaction R a t e and Temperature Contours 113 List of Tables vi List of Tables 3.1 Comparison of global reaction rates for different smoothing matrices 45 3.2 Comparison of global reaction rates for libraries of various sizes 50 3.3 Comparison of global reaction rates for converged statistics of reaction rate 1 3.4 Comparison of global reaction rates for different scalars . . . 51 56 3.5 Comparison of global reaction rates for combinations of scalars 59 4.1 Ignition characteristics based on the ignition criterion of 2000 K 4.2 Peak reaction rate characteristics 81 84 A . l Comparison of global reaction rates for converged statistics of reaction rate 2 96 A. 2 Comparison of global reaction rates for converged statistics of reaction rate 3 96 B. l Additional Comparison of global reaction rate for combinations of scalars 98 List of Figures vii List of Figures 2.1 Energy spectrum 13 2.2 Fast chemistry relations between mass fractions and mixture fraction 2.3 Borghi diagram 20 25 3.1 Weighting coefficient vectors for each H matrix of the linear regularization parameters 42 3.2 Reaction rate 1 predictions for different smoothing matrices in the linear regularization 43 3.3 Reaction rate 2 predictions for different smoothing matrices in the linear regularization 43 3.4 Reaction rate 3 predictions for different smoothing matrices in the linear regularization 44 3.5 Surface plot of temperature flamelet libraries shown in order of increasing scalar dissipation 46 3.6 Weighting coefficient vectors for the six flamelet libraries . . . 47 3.7 Reaction rate 1 predictions for different flamelet libraries . . . 48 3.8 Reaction rate 2 predictions for different flamelet libraries . . . 49 3.9 Reaction rate 3 predictions for different flamelet libraries . . . 50 3.10 Weighting coefficient vectors for the six scalars 3.11 Reaction rate 1 predictions for different scalars 53 54 3.12 Reaction rate 2 predictions for different scalars 55 3.13 Reaction rate 3 predictions for different scalars 55 3.14 Reaction rate 1 predictions for combinations of scalars . . . . 57 3.15 Reaction rate 2 predictions for combinations of scalars . . . . 58 List of Figures viii 3.16 Reaction rate 3 predictions for combinations of scalars . . . . 58 4.1 Unsteady round jet 64 4.2 Two dimensional mesh 66 4.3 Comparison of full and ignitingflameletlibraries 70 4.4 Temperature contours for the initial temperature of 1200 K using the fullflameletlibrary 73 4.5 Comparison between predicted and experimental ignition delay for fullflameletlibraries 74 4.6 Fullflameletlibrary weighting coefficient vector for ignition plane 75 4.7 Comparison of source terms for igniting and extinguishing flamelets 75 4.8 Temperature contours at different times during methane ignition 76 4.9 Ignitingflameletlibrary weighting coefficient vector for ignition plane 77 4.10 Comparison between predicted and experimental ignition delay for ignitingflameletlibraries 78 4.11 Light emission images from shock tube experiments for 1475 K 79 4.12 Temperature contours for 1475 K 79 4.13 Reaction rate contours for 1475 K 80 4.14 Temperature contours overlaid with mixture fraction and variance isopleths for 1200 K 82 4.15 Reaction rate and temperature contours at the onset of ignition 83 4.16 Reaction rates as functions of mixture fraction 83 C. 1 Predicted contours of mixture fraction for an unsteady methane jet 100 C.2 Comparison between predicted and experimental mixture fraction 101 List of Figures ix C.3 Predicted contours of mixture fraction variance for an unsteady methane jet 102 C. 4 Comparison between predicted and experimental mixture fraction variance 103 D. l Contours of fuel mass fraction for different grid spacings . . . 105 D. 2 Centreline fuel mass fractions for different grid spacings . . . 106 E. l Temperature contours at different times during methane ignition, 1000 K 107 E.2 Temperature contours at different times during methane ignition, 1100 K 108 E.3 Temperature contours at different times during methane ignition, 1300 K 108 E. 4 Temperature contours at different times during methane ignition, 1400 K 109 F. 1 Temperature contours overlaid with mixture fraction and variance contours for 1000 K 110 F.2 Temperature contours overlaid with mixture fraction and variance contours for 1100 K Ill F.3 Temperature contours overlaid with mixture fraction and variance contours for 1300 K Ill F. 4 Temperature contours overlaid with mixture fraction and variance contours for 1400 K 112 G. l Reaction rate and temperature contours at the onset of ignition for 1000 K 113 G.2 Reaction rates as functions of mixture fraction for 1000 K . . 114 G.3 Reaction rate and temperature contours at the onset of ignition for 1100 K . 114 G.4 Reaction rates as functions of mixture fractionfor1100 K . . 115 List of Figures x G.5 Reaction rate and temperature contours at the onset of ignition for 1300 K 115 G.6 Reaction rates as functions of mixture fraction for 1300 K . . 116 G.7 Reaction rate and temperature contours at the onset of ignition for 1400 K 116 G.8 Reaction rates as functions of mixture fraction for 1400 K . . 117 Nomenclature xi Nomenclature a CSE weighting coefficient vector m A n x m matrix in matrix equation B Smoothing matrix in linear regularization Bj Pre-exponential factor in Arrhenius Equation C Right hand side of matrix equation C\, C2, Cp, C , C<= ei 2 Model constants C Constant pressure specific heat (5^^) D Species molecular diffusivity i~) Da Damkohler number E Activation energy (kJ) g Gravitational forces (N) Gk Stress tensor in turbulent kinetic energy equations h Specific enthalpy (kJ/kg) H Smoothing characterization for linear regularization k Turbulent kinetic energy (fcJ) p a Jj Enthalpy diffusivity ( ^ ) Jj Species diffusivity (^) kb Backward reaction rate constant kf Forward reaction rate constant 10 Eddy length scale (m) 11 Laminar flame thickness (m) L Characteristic length (m) Le Lewis number 1 Nomenclature m Number of flamelets in C S E basis function rtij Forward reaction order Mi Chemical symbol representing i n Number of scalar realizations in the C S E integral equation rij Backward reaction order N Number of chemical species in kinetic reaction p Pressure P Probability density function Pr Prandtl number R Universal gas constant Re Reynolds number ih species (kPa) (j^) Number of reaction steps in kinetic reaction T i m e (s) T Temperature (JK) Tf Flamelet library for C S E basis function u Velocity ) v Velocity ) v" Turbulence intensity VL Laminar burning velocity ( ™ ) w Velocity ( a ) x Direction x Unknown vector in matrix equation Yk Species mass fraction for the k Z Mixture fraction m Z" th 2 Mixture fraction variance species xii Nomenclature Greek Letters Temperature exponent in Arrhenius equation a Represents the Beta probability density function Small increment of change 8 Kronecker delta € Dissipation of turbulent kinetic energy (^) 7 Weighting function for smoothing in linear regularization A Thermal diffusivity Viscosity ( ^ ) IH Eddy viscosity v' Stoichiometric coefficient for reactant species v" Stoichiometric coefficient for product species $m CSE basis function 9 Density ( ^ ) T Time scale (s) T Chemical time scale (s) ij Viscous stress tensor T Turbulent time scale (s) e Arbitrary scalar c T t Chemical reaction rate X Scalar dissipation (i) c Conditional value of mixture fraction xiii Acknowledgements xiv Acknowledgements The author would like to take this opportunity to acknowledge the people who played an important role during the realization of this project and the writing of this thesis. First, I would like to gratefully acknowledge the funding support from the Natural Sciences and Engineering Council of Canada and Westport Innovations Inc. I wish to express my appreciation to my supervisors, Dr. K. Bushe and Dr. G. Li for their invaluable assistance, guidance and encouragement throughout this research. Of course, I need to thank a long list of my office mates, J. Lozada, C. Larson, R. Grout, H. Jones, S.H. Park, and A. Frisque for their contributions and discussions of all things combustion. Most importantly, I would like to thank my family and friends, particularly my loving wife, Jennifer. Without her continued support and encouragement the completion of this thesis would not have been possible. It is to all my family and friends that I dedicate my thesis. 1 Chapter 1 Introduction 1.1 Introduction Diesel engines have long been the traditional power source for both the transportation and power generation industries but are also among the largest contributors to environmental pollution (Taylor [1]). Air pollution has been identified as a contributor to the depletion of the ozone layer and global warming through increased carbon dioxide emissions. Another very noticeable result of emissions from the combustion process is photochemical smog caused by solar radiation on air polluted with oxides of nitrogen and hydrocarbons. These emissions result in an increased danger to human health including a higher risk of cancer and pulmonary diseases (Sher [2]). One method to reduce the emissions from diesel engines is to burn cleaner burning fuels like natural gas or methane. Westport Innovations Inc. has developed a strategy called High Pressure Direct Injection (HPDI) that allows the diesel engine to operate on natural gas while maintaining comparable power output and efficiency as conventionally fuelled diesel engines. Natural gas is injected into the combustion chamber and ignited with a small amount of pilot diesel fuel. As experimental development costs are high, it is desirable to develop computer models to predict the ignition of methane to facilitate the design and development of the HPDI system. In the past quarter century, advances in computer technology have sparked the development of computer models to assist the combustion engineer. As such, models to account for fluid dynamics, turbulent flows, and chemical Chapter 1. Introduction 2 reactions have been developed for Computational Fluid Dynamics (CFD) codes and have become a valuable tool for design engineers (Bilger [3]) providing insight and understanding of the combustion process (Veynante and Vervisch [4]). Several models have been developed for turbulent combustion modeling including Laminar Diffusion Flamelet models (Peters [5]), Probability Density Function (PDF) methods (Pope [6]), Conditional Moment Closure (CMC) (Bilger and Klimenko [7]), and more recently Conditional Source-term Estimation (CSE) (Bushe and Steiner [8]). 1.2 Objective The objective of this research is to further explore the concept of Conditional Source-term Estimation. First the parameters affecting the solution of CSE were explored in an attempt to optimize the solution and improve the closure of the chemical source terms. Once the optimization process was complete, the CSE method was implemented into a commercial combustion modeling package and used to predict the ignition of methane under diesel engine-relevant conditions. 1.3 Thesis Outline In this chapter, the motivation for this research has been presented and discussed. The next chapter begins with a review of combustion theory including the governing equations and an introduction to turbulence modeling followed by a discussion of three different combustion models; Laminar Diffusion Flamelets, Conditional Moment Closure, and Conditional Source-term Estimation. The third chapter presents the results of further a priori testing of the Conditional Source-term Estimation method including optimizing the various parameters that affect the prediction of the chemical source terms. The implementation of the CSE method in a commercial software package Chapter 1. Introduction 3 and its use to predict the ignition of methane is described in the fourth chapter. The final chapter summarizes the findings of the preceding chapters and also proposes additional research to further the understanding of the Conditional Source-term Estimation method and its applicability to turbulent combustion modeling. 4 Chapter 2 Background: An Introduction to Turbulent Combustion Modeling 2.1 Introduction This chapter begins with a description of combustion and a review of combustion simulation principles including the governing equations. A more detailed look at several key combustion model developments including: Laminar Diffusion Flamelet models, Conditional Moment Closure and Conditional Source-term Estimation follows. This review forms the foundation for this research by defining the key terms and outlining the CSE methodology. It is not meant as an exhaustive review of combustion simulation research but as an overview for the context of this work. 2.2 Combustion Combustion is a chemical reaction that occurs when a fuel is oxidized to release energy. Typical fuels consist of hydrogen and carbon denoted by the form C H m n while the most common oxidizer is air with nitrogen behaving as an inert species. In the case of complete combustion, the fuel and oxidizer react together to produce water and carbon dioxide. However, complete combustion rarely occurs. In reality, the fuel and oxidizer are accompanied by additional components like sulphur and nitrogen that also react during Chapter 2. Background 5 the combustion process to form dangerous pollutants such as oxides of nitrogen (NO ), and oxides of sulphur (SO ), among others. x x There are two distinct modes of combustion in IC engines: premixed and non-premixed combustion (Heywood [9]). The first occurs typically in spark ignited (SI) engines where the fuel and oxidizer are thoroughly mixed before entering the combustion chamber and are ignited by an ignition source. The second mode, which occurs when the fuel and oxidizer are injected separately into the combustion chamber where they mix and burn, is commonly found in combustion chambers for gas turbines and diesel engines. This type of combustion is referred to as mixing limited or non-premixed combustion as the reaction rate depends on the mixing of the reactants. The objective of this research is to model methane ignition for high pressure direct injection (HPDI) in diesel engines and so the focus will be primarily on non-premixed combustion. 2.3 Conservation Equations Conservation equations are balance equations for conserved quantities, such as energy, and can be applied to a system of fixed identity with everything external to the system termed the surroundings and separated from the system by a boundary. Any gain in a property by the system must be offset by an equal loss of property to the surroundings (Qengel and Boles [10]). These same fundamental equations are used to describe turbulent reacting flows and can be modified to include chemical reactions, species balance and variable-property effects. 2.3.1 Conservation of Mass The conservation of mass states that mass can neither be formed nor destroyed (Warnatz et al.[ll]). Chapter 2. Background dp dpuj dt dxj 6 (2.1) 0 The first term represents the time, rate of change of mass while the second term represents the mass flux density of the system. 2.3.2 Conservation of Momentum The conservation of momentum states that momentum in a control volume is conserved (Warnatz et al. [11]). In the case of a viscousflowfield, this yields the Navier-Stokes equation: d(puj) dt d (pujUj) dp_ + dxj = QTJJ dxi dxj + P9 (2.2) The first term represents the time rate of change of momentum, the second term is the convective momentumfluxand the third term is the momentum flux due to pressure. The viscous stress tensor (Eq. 2.3) is the first term on the right hand side and the last term is gravitational effects or momentum due to body forces. (2.3) where p is the fluid viscosity and <5y is the Kronecker delta: i 1 for i = j 0 for i ^ j (2.4) Chapter 2. Background 2.3.3 7 Species Balance Additional equations can be derived for the species balance: dY dt P k | d{pujY ) dx, k = djf dx where Y represents the mass fraction of the k species and uj is the chem- th k k ical source term or the mass rate of change of the species due to chemical reaction and accounts for the production or consumption of the species. T h e species diffusivity is generally defined using Fick's Law: here D is the species molecular diffusivity. A simplification is to assume equal species diffusivities which leads to the unity Lewis number assumption which is generally valid for hydrocarbons (Warnatz et al. [11]): L e - p c ^ ^ - (2 7) where C is the constant pressure specific heat capacity and A is the thermal p diffusivity. T h e Lewis number is the ratio between the thermal and mass diffusivities. 2.3.4 Conservation of Energy T h e conservation of energy, shown in terms of the specific enthalpy, is also known as the First Law of Thermodynamics (Warnatz et al. [11]). dph dp d(pujh) _ 9 ( h \ T h e first term represents the time rate of change of specific enthalpy, the Chapter 2. Background 8 second is the time rate of change of pressure followed by the convective term. The enthalpy diffusive and viscous flux terms are on the right hand side and finally the heat generation due to radiation. Enthalpy diffusion is described by the Fourier Law: where Pr is the Prandtl number, or the ratio of momentum and thermal diffusivities: Pr=^- 2.3.5 (2.10) Chemical Kinetics To understand the chemical source term (oJk), consider a general system in- volving NR reactant species and Np product species and R reaction steps s that has the following form (Libby and Williams [12]): N Np R £ v' Mi i=l kj 4jMi, ^ j = 1,2,..., R , s (2.11) i=l where Mj is the chemical symbol for species i, v\j and v"j are the stoichiometric coefficients for species i as reactant and product respectively. The reaction would have a corresponding forward reaction rate constant (fc/) and backward reaction rate constant (kb). Both rate constants take the form of an Arrhenius equation with a pre-exponential factor (Bj) k = BjT e~~Br aj fj (2.12) where E j is the activation energy, T is the temperature and R is the a Chapter 2. Background 9 universal gas constant. The Arrhenius relationship can be developed from experiments using empirically determined values. The rate of change in species concentration (c) for species i in reaction r is given by the expression(Libby and Williams [12]) : (^)= fWi-»ri)h s'° ( - ) k ^ C ' 2 13 s=l where c is the concentrations of the S different species s. a From these reaction steps, an expression for the chemical source term can be derived (Eq. 2.14) by summation over the rate equations (Eq. 2.13) where rrij and rij are the overall orders of the forward and backward reactions (Libby and Williams [12]): ** = £ K - O 3=1 U ^p mj k=i \ (i n K • ii fc=i / (2.14) Clearly, the chemical reaction rate is a highly non-linear function of temperature and species mass fractions. 2.4 Turbulence Turbulence is the fluctuating and agitated motion of a flowing fluid and is a property of the flow, not of thefluid.It can be characterized by continuous fluctuations in the velocity about its mean (White [13]). In a mixture, these fluctuations in the velocity can lead tofluctuationsin scalars like density, temperature and composition. The Reynolds number (Eq. 2.15) is the ratio of dynamic to viscous forces and described by the velocity (u), a characteristic length (L) and thefluidviscosity (p,) and is used to determine whether Chapter 2. Background 10 or not aflowis turbulent. Jtes^ A (2.15) 1 An important dimensionless number used in turbulent combustion is the Damkohler number (Eq. 2.16) which is a comparison of turbulent (r ) and t chemical (T ) time scales (Veynante and Vervisch [4]). It describes the time c needed for chemical change compared to the time needed for change induced byfluidmotion. Da = - 2.4.1 (2.16) Averaging Concepts In a turbulentfieldwithfluctuatingscalars, it is often convenient to determine mean, or average, values for the scalars. Time averaging is defined as: — @ r 1 -fJ t+T @ t ^ d T - <- ) 2 17 The ensemble average is defined as: where N is the number of realizations. Afluctuatingscalar can then be written with its mean and fluctuating term (0') as: B(x,t) = e(x,t) + Q'(x,t) (2.19) Chapter 2. Background 11 If the transport equations are averaged using this concept (Eq. 2.19), many terms appear for which there is no closure. For example, if the density and velocity are broken into their respective means and fluctuations, a new velocity correlation term (p'u ) appears in the conservation of mass equation 1 that cannot be expressed from the known variables. Similarly, more unclosed terms appear for the conservation of momentum and conservation of energy equations. One method of reducing the number of unclosed terms from the transport equations is to use density-weighted, or Favre averaging: 0= P (2.20) The scalar can now be written as: e(x,t) = Q(x,t) + e"(x,t) (2.21) where 0 is the density weighted average and 0" is the density weighted fluctuation. While this may simplify the equations, several unclosed terms remain that tend to diffuse momentum and therefore act to increase the effective viscosity. As a result, the molecular viscosity in the transport equations is enhanced with an additional term called the eddy viscosity (fit)Another useful averaging concept is the conditional average which means only the realizations that meet a given criteria or condition are averaged. For non-premixed combustion, the conditional scalar of choice is usually the mixture fraction (27), discussed later. The conditional average (Qi) of a scalar Oj is expressed as: Chapter 2. Background Qi (C; x, t) = {Qi (x, t) \Z (x, t) = 0 12 (2.22) and is the average of 0 given that the mixture fraction has a value of C 2.4.2 Turbulence Simulation There are three strategies for simulating turbulence in computational fluid dynamics. 1. Direct Numerical Simulation (DNS) requires solving the transport equations while resolving all terms in the equations for all length and time scales including the Kolmogorov length and time scales without making any modeling approximations. 2. Large Eddy Simulations (LES) resolves the integral scales while using models for viscous dissipation and combustion in the unresolved scales. 3. Reynolds Averaged Navier-Stokes (RANS) models solve the ensemble or time averaged forms of the transport equations and model fluctuating dissipation and combustion. DNS is the most accurate solution as it resolves all length and time scales without resorting to modeling assumptions. However, it is also the most computationally expensive. At present, DNS is only able to determine solutions for low Reynolds Number flows that are of little practical interest. This is not to say that DNS does not provide useful results, in fact, DNS has the capability to extract or isolate specific phenomena that can be used to test and validate models for other turbulence modeling strategies. The basis for LES is to use the Kolmogorov hypothesis to model the eddy viscosity. The goal is to compute the largest structures in the flow while modeling the effects of the smaller scales. The large structures are Chapter 2. Background 13 dependent on the system geometry while the small structures tend to exhibit more universal behaviour. Examining the energy spectrum (Figure 2.1), LES resolves the scales to the left of the vertical line (integral scales) then models a turbulent viscosity that dissipates the equivalent amount of energy to the right of the vertical line (Kolmogorov scales). This results in a more efficient turbulence model as the smaller scales are simulated rather than averaged as in RANS. As well, LES provides valuable information on the resolved scales of motion. Figure 2.1: Energy spectrum RANS is the most common turbulence model found in commercial software packages. Either a time or ensemble average is calculated and then closure is required for any unclosed terms, usually the velocity fluctuation correlations (u^Uj). The principal effect of these unclosed terms on the fluid flow is to diffuse momentum which leads to the development of the eddy viscosity. One of the most common methods for modeling the eddy viscosity is to use the k — e model where transport equations are written for the turbulent kinetic energy (k) (Peters [14]): Chapter 2. Background 14 + G -pe k (2.23) and its dissipation (?) (Peters [14]): d(pe) at _d_ + dx where Gk is a function of the stress tensor: G k = + The eddy viscosity can be modeled using (Launder and Spalding. [15]): IH = In this case, the constants C , C M Cy. = 0.09 2.5 £ l C>7 • , and C C = 1.44 ei e 2 (2- are empirically determined. C = 1.92 t2 (2.27) Turbulent Combustion Models Turbulent combustion is a complex problem requiring the solution of the transport equations in a highly coupled system. Fluid dynamics describe the turbulent mixing and transport phenomena such as heat transfer and diffusion while the chemistry describes the production and consumption of species. Multi-phase systems and radiation can also add to the complexity Chapter 2. Background 15 of the problem. The objective of turbulent combustion modeling is to propose closures for the unknown terms that arise in the averaged transport equations. 2.5.1 Mixture Fraction Mixture fraction (Z) is an important scalar in non-premixed combustion which describes the state of mixedness of two or more fluids and is a measure of the stoichiometry. The mixture fraction is a conserved scalar that represents the fraction of gas that originated from the fuel stream and is conveniently defined to be unity in pure fuel and zero in pure oxidizer. A transport equation can be written for the Favre averaged mixture fraction (Peters [14]): (2.28) with the velocity correlation: (2.29) There is no chemical source term as the mixture fraction is a conserved scalar. Assuming equal diffusivities, the mixture fraction has a linear relationship with the species mass fractions. Many combustion models solve transport equations for the species mass fractions but not the mixture fraction. As such, an inert element such as argon (already present in air), can be included in the oxidizer stream to act as a tracer for the transport of the mixture fraction. Using the linear relationship between mixture fraction and the mass fraction of argon, the equation for mixture fraction becomes: Chapter 2. Background 16 Z=1 - & Yiro where Y aTi (2.30) is the massfractionof argon at time i and Y aro is the initial mass fraction of argon. 2.5.2 Mixture Fraction Variance In addition to the mean mixturefraction,a transport equation can be derived for the mixturefractionvariance (Z" ) (Peters [14]). 2 (2.31) recalling that the velocity correlation is: u " ' =- A ^ Z a (2.32) X is the scalar dissipation rate which dissipates fluctuations in scalars much the same way as viscosity dissipates fluctuations in the velocityfield.The mixturefractionvariance characterizes non-homogeneities and intermittencies in the mixing field and is a measure of the degree of mixing at unresolved scales. The mixturefractionand its variance provide a description of the turbulent mixing. Using empirically derived constants (Peters [5]), the transport equation for the mixturefractionvariance can be re-written: 17 Chapter 2. Background d(pZ" ) 2 dt d + dxj La*) _d_ dxj [pUj V2,eff OXj + cr pt ( dZ dZ o~2,t dxj dxj J - PX (2.33) with k 2 0"2,e// = 0.9 <T ,t = 1-0 2 Ci = 2.7 p t (2.34) and Cp has been previously denned (Eq. 2.27). The model: PX = pC lz^, x C = 2.0 x (2.35) has been suggested for the last term on the right hand side of Eq. 2.33 (Peters [14]). 2.5.3 Scalar Dissipation Rate The influence of theflowfieldon the mixture fraction and its variance is introduced through the scalar dissipation rate (Peters [5]). The scalar dissipation rate acts as an external parameter imposed on the flame structure and directly measures the decaying speed of fluctuations via turbulent micro-mixing. The scalar dissipation rate is a measure of the inverse of a characteristic diffusive time and as the scalar dissipation rate increases, the heat and mass transfer are enhanced (Veynante and Vervisch Chapter 2. Background 18 [4])In a simulation, the scalar dissipation (2.37) X = C7i can be estimated from a first order approximation of the gradient in the mixture fraction (Chapra and Canale [16]). 2.5.4 Probability Density Function The Probability Density Function (PDF) quantifies the probability of finding a scalar 8 between the values of 0 + 50 and 0 — (50 at a given location and time (xj,t). The PDF can provide a probabilistic description of the internal structure of a turbulent flame by focusing on its statistical properties (Veynante and Vervisch [4]). Allflamescan be parameterized by the fuel mass fraction Yf and the mixture fraction Z and completely described by the joint PDF of mixture fraction and fuel mass fraction: P(Y ,Z; ,t) f (2.38) Xj The conditional PDF is introduced to better study the fuel mass fraction at a given mixture fraction: P{Y , Z;x t) = Pi (Y \Z;x t)P f it f Jt 2 (Z;x ,t) 3 (2.39) In this case P<i{Z;xj,t) describes the statistics of the mixing fuel and oxidizer and P\ (Yf\Z;Xj,t) is linked to the internal structure of the flame (Veynante and Vervisch [4]). Chapter 2. Background 19 From the conditional PDF, it is easy to define the conditional average, that is the average of a quantity given a specific condition (e.g. average all the fuel mass fractions in the flowfieldwhere the mixture fraction is a specified value). (Y \Z{xj,t) = Q= pYflMYffcxjrfdYf f (2.40) A presumed shape, typically the /3-function, is often chosen to represent the PDF of mixture fraction (P (Z;Xj,t)){Cook and Riley [17]). The flexibility 2 of the /3-function exhibits the appropriate limiting behaviour observed in experiments from unmixed reactants (delta functions) to well-mixed reactants (Gaussian distribution) (Cook and Riley [17]). Once the values for the mixture fraction and its variance are known, it is then straightforward to capture the mixing of the fuel and oxidizer. For this reason, the /3-function parameterized by the mixture fraction and its variance is chosen as the presumed form of the PDF. The /3-PDF (P (Z;xj,t)) is given by: 2 7=r(o)r(6) p (Z;x ,t) = 1z - (i-z) - , a 2 l b 1 j ( 2 ' 4 1 ) The normalization parameter (7) is taken to be one (Libby and Williams [18]) and the parameters a (Eq. 2.42) and b (Eq. 2.43) are obtained from the mean (Z) and variance (Z" ) of the mixture fraction. 2 ^ ( F S M ) - b= I - a 1 ) (2 - 42) (2.43) The PDF can be tabulated in a library or calculated directly from a nu- Chapter 2. Background 20 merical algorithm (Press et al. [19]). In this implementation of the CSE method, the /3-function is calculated from the numerical algorithm for accuracy considerations and results in a trade-off in computation time for the sake of accuracy. 2.5.5 Fast Chemistry The fast chemistry assumption is the simplest approach that can be used to illustrate general concepts for combustion modeling. This method assumes that the chemical time scales are considerably smaller than the mixing time scales resulting in local chemical equilibrium, which leads to the maxim that "mixed is burned" (Warnatz et al. [11]). Often, the entire process is approximated by a single step reaction between fuel and oxidizer to produce products. Relating the species mass fractions and temperature to the mixture fraction leads to the development of linear relationships (Figure 2.2). Figure 2.2: Fast chemistry relations between mass fractions and mixture fraction Chapter 2. Background 20 merical algorithm (Press et al. [19]). In this implementation of the C S E method, the /3-function is calculated from the numerical algorithm for accuracy considerations and results i n a trade-off i n computation time for the sake of accuracy. 2.5.5 Fast Chemistry The fast chemistry assumption is the simplest approach that can be used to illustrate general concepts for combustion modeling. T h i s method assumes that the chemical time scales are considerably smaller than the mixing time scales resulting i n local chemical equilibrium, which leads to the m a x i m that "mixed is burned" (Warnatz et al. [11]). Often, the entire process is approximated by a single step reaction between fuel and oxidizer to produce products. Relating the species mass fractions and temperature to the mixture fraction leads to the development of linear relationships (Figure 2.2). Figure 2.2: Fast chemistry relations between mass fractions and mixture fraction Chapter 2. Background 22 semble of thin, one-dimensional reaction zones in a locally laminar mixing region (calledflamelets).Working with the conserved scalar concept, the transport equations can be re-written in terms of the mixture fraction and scalar dissipation rate which provides the coupling between the flow field and the chemistry that is missing in the Fast Chemistry model. Flamelet models can be viewed as an improvement over the Fast Chemistry model as thefiniterate chemistry effects are included through the additional parameter: the scalar dissipation rate (x)Assuming that the reaction zone is small compared to the small scales of turbulence and restricting the analysis to a small region around the flame, a coordinate transformation from physical space into mixture fraction space is performed. With the assumption of unity Lewis number, the classical flamelet equations can be developed (Peters [5]) for species mass fraction: dYi pd Yi 2 (2.46) and temperature: (2.47) Steady Flamelets Assuming that the scalar dissipation rate varies slowly enough such that the time dependence can be neglected (Peters [14]), the steady equation for species mass fraction becomes: P&Yi + 2dZ 2 uJi = 0. (2.48) If this condition holds, then a library of solutions to the steady flamelet Chapter 2. Background 23 equations can be developed with the solutions listed in order of increasing scalar dissipation rate. This library of solutions can then be used as a reference table for calculating species mass fractions and temperature and becomes a direct replacement for the equilibrium chemistry model in turbulent combustion simulations. However, this model does not capture the unsteady transitions associated with ignition and extinction. Unsteady Flamelets While it may be convenient to neglect the unsteady term in the classicalflameletequations, mostflamesof practical interest are neither onedimensional nor steady. Most turbulentflamesinclude unsteady effects such as ignition and quenching and require the inclusion of the unsteady terms (Pitsch [20]). Two of the more commonflameletmodels that consider the unsteady phenomena are the Lagrangian Flamelet Model (Pitsch et al. [21]) and the Eulerian Particle Flamelet Model (Coelho et al. [22]). Both of these models assume that theflameletis convected by the turbulent flow, and integrate the classical flamelet equations (Eqs. 2.46, 2.47) in time and account for the time variation of the scalar dissipation rate (x) to capture theflamelethistory. Often, the solution to the unsteadyflameletequations is done interactively with the calculation of theflowfield. In the Lagrangian Flamelet model, the time derivative (J^)fromthe flamelet equations is evaluated at constant mixture fraction. The flamelet coordinate system will move with the velocity of a point on the stoichiometric surface. Thus, the position of theflameletis uniquely related to time which corresponds to a Lagrangian treatment (Pitsch [20]). Closure of the flamelet equations is achieved by replacing the scalar dissipation rate and velocity with the instantaneous local conditionally filtered values. Consequently, the temperature and species mass fraction obtained through this method are also the conditional average values. Chapter 2. Background 24 The Eulerian Particle Flamelet model introduces fluid particles that are transported with the fluid. Each particle represents a uniqueflameletwith its own scalar dissipation rate history. Initially the particles are uniformly distributed through the domain and transport equations for the particles are solved along with theflowcalculations. At each time step, the scalar dissipation rate is calculated for each particle and assumed to have varied linearly from the previous time step to the present one. Knowing the variation of the scalar dissipation rate, the unsteadyflameletequations can then be integrated for the time increment (Coelho et al. [22]). Applications of Steady and Unsteady Flamelets Three regions of turbulent non-premixed combustion can be identified according to the Borghi diagram (Figure 2.3) (Warnatz et al. [11], Lentini [24] and Heywood [9]). Plotted is the turbulence intensity (v") normalized with the laminar burning velocity (VL) against the eddy length scale (lo) normalized with the laminar flame thickness ( I I ) - The three distinct regimes are (Veynante and Vervisch [4]): 1. The Flamelet Regime occurs when the chemical time scales are smaller than the Kolmogorov time scales. This regime is characterized by very thinflameswhose internal structure is not affected by turbulent motion. 2. The Perturbed Flamelet Regime is the second regime and occurs when the chemical time scales are larger than the Kolmogorov time scales but smaller than the largest turbulent time scales. In this regime, the flame interacts with the turbulent motion but maintains many of the flamelet features. 3. The last combustion regime is the Thickened Flame regime and occurs when the chemical times scales are larger than the largest turbulent Chapter 2. Background 25 time scales. This regime behaves as a well stirred reactor. log(lo/l ) L Figure 2.3: Borghi diagram Comparisons between predictionsfromsteady laminarflameletmodels with experimental results have shown good correlations for jet diffusion flames (Buriko et al. [25]) and can account for finite rate chemistry in turbulentflames(Lentini [24]). When the strain rate is constant or varies slowly compared to the chemistry time scales, the steadyflameletmodel provides a good approximation to the real solution for turbulent diffusion flames (Cuenot et al. [26]). However, a considerable amount of the prediction error in the steady laminarflameletmodel may be due to unsteady effects that have been neglected. In an effort to improve the predictions for emissions, and unsteady phenomena (ignition and extinction), many researchers have begun using unsteadyflameletmodels. The unsteady flamelet models have been tested in the a priori sense against DNS data and compared to experimental results. It was found that the unsteadyflameletmodel exhibited better performance and smaller errors compared to the steadyflameletmodel in DNSflowsex- Chapter 2. Background 26 hibiting intermittency (Bourlioux et al. [27]). The unsteadyflameletmodel was also found to yield better predictions near extinction in counterflow diffusionflames(Cuenot et al. [26]). The predictions for the species mass fractions were improved (Coelho et al. [22] and Pitsch [23]) by using the unsteadyflameletmodel rather than the steady flamelet model to simulate the experimental results from the Sandia piloted methane/air jet (Barlow et al. [28]). Both the steady and unsteadyflameletmodels provide a convenient decoupling of the chemistry andflowfield in turbulent reactingflows.When in the combustion regime where the assumptions of thin reaction zone and reactants mixing in a laminar fashion hold, the laminar diffusion flamelet model provides good approximations to both DNS and experimental data. However, the predictions from these models become less reliable in the combustion regimes where these assumptions are no longer valid. 2.5.7 Conditional Moment Closure Another prominent turbulent combustion model was independently developed by Bilger ([29]) and Klimenko ([30]) and has become known as Conditional Moment Closure (CMC). The CMC method involves deriving and closing transport equations of conditional averages of scalars based on the mixture fraction. The underlying assumptions of the CMC method is that fluctuations in the scalar quantities of interest can be associated with fluctuation of the mixture fraction. Thus, in a non-premixed combustion problem where fuel and oxidizer are mixing, the values of the concentrations and temperature within the mixingfieldare dependent on the local, instantaneous value of the mixture fraction. While the general formulation for CMC and Flamelet Models may have similarities, the underlying physical assumptions are considerably different. CMC incorporates diffusion as turbulent micromixing and decouples it from chemistry while the flamelet assumption uses a Chapter 2. Background 27 coupled representation of diffusion and reaction (Veynante and Vervisch [4]). The assumption of a thin reaction zone is not necessary for the development of the CMC equations leading to their validity across all combustion regimes. Following the derivation of the CMC equations, the conditional chemical reaction rate is closed by assuming that the fluctuations around the conditional scalars are small and then the chemical reaction rate is approximated using the conditionally averaged scalars (Bilger and Klimenko [7]): ^{P^YK)^ u! (HC, rJC, YKK) • (2-49) Applications o f Conditional M o m e n t Closure Although a relatively new development in turbulent combustion modeling, the CMC method has been applied to a wide range of practical problems with comparison to both DNS and experimental data. DNS investigations of the validity of the CMC method found that the overall predictions of the species mass fractions were good but the model was sensitive to the modeling of the scalar dissipation rate (Mell et al. [31]). This sensitivity was observed to increase for faster chemistry rates as the model approached the fast chemistry limit. Numerical studies have also been undertaken to determine the applicability of the CMC model for predicting lifted turbulent flames (Devaud et al. [32]). Comparisons made with lifted turbulentflameexperiments have found that the lift-off height can be reasonably predicted with this approach. Predictions of species mass fractions were shown to have good qualitative agreement with experimental data. In all, the CMC method was found to give satisfactory results for simulations of lifted turbulent flames. Chapter 2. Background 28 An important issue in the design of IC engines is the prediction of emissions. To that end, many researchers have investigated the ability of CMC to predict the formation of pollutants. In the investigation of hydrogen jet flames, it was determined that the CMC method gave useful predictions of the formation of Nitrous Oxide (NO) (Barlow et al. [33]). As well, it was found that both major and minor species and flame temperature could be reasonably predicted compared to experimental data using this approach (Fairweather et al. [34]). The CMC approach has been tested in practical applications and has demonstrated its ability to provide a convenient means for incorporating the kinetic effects into the turbulent flow calculations. This method adds another dimension, the conditioning variable (mixture fraction) to the solution of the problem resulting in an increase in computation time. However, this approach provides reasonable closure for the chemical source terms by hypothesizing that the conditional average of the chemical source terms is a function of the conditional average of the scalars. 2.5.8 Conditional Source-term Estimation CSE is a recent development in the area of turbulent combustion simulation (Bushe and Steiner [8], [35], [36]) and uses the same closure hypothesis for the chemical source terms as the CMC method. The CSE method presumes a shape for the PDF of the mixture fraction and then inverts integral equations for the conditional average of species mass fractions and temperature. These conditional averages are then substituted into the reaction rate expressions to obtain the conditional reaction rates thereby obtaining closure for the chemical source terms. In CMC, the transport equations for the conditional averages are solved explicitly which involves adding a new independent variable (mixture frac- Chapter 2. Background 29 tion) to the system of equations and leads to a high computation cost. CSE takes advantage of spatial homogeneity of the conditional averages observed in a mixing layer; the conditional averages of scalars vary in the direction of the mean flow but have a weak sensitivity in the direction normal to the fuel and oxidizer interface (Bushe and Steiner [8]). Thus, integral equations can be written and solved for the scalars of interest and closure for the chemical source terms is achieved with adequate accuracy and a reduction in computation time. Formulation The chemical source terms in the species balance equation (Eq. 2.5) are highly non-linear functions of the local temperature, density and species mass fractions. Consequently, using the average values of these scalars yields a poor approximation of the average reaction rates e.g.: (2.50) Instead, CSE uses the same closure for the chemical reaction rates as the CMC method (Eq. 2.49) (repeated here) (p, T, F ) |C « c^p|C, T|C, F/dC) K Using the conditional averages of the scalars gives reasonable approximations for conditional average of the reaction rates compared to DNS data (Bilger and Klimenko [7]). In both the Reynolds Averaged Navier-Stokes and Large Eddy Simulation paradigms, the transport equations are solved for the average values of the scalars. Integral equations can be written for the scalars (species mass fractions, temperature, etc.): (2.51) Chapter 2. Background 30 where O (x, t\Q represents the conditional average of the scalar. The PDF of the mixture fraction, represented by P(C.;x,t), describes the turbulent mixing of fuel and oxidizer, and is often approximated using the presumed shape of the ^-function, P (C; Z (x,t), Z" 2 (Cook and Riley [17]). Both the mixture fraction and its variance can be obtained by solving their respective transport equations. In a simulation, the average value of a scalar is known at many discrete points, therefore the conditional average is the only unknown in the integral equation. Assuming that the gradient over the mixture fraction in the conditional average is small (Bushe et al. [37]), the integral equation can be inverted and solved for the conditional average of the scalar using a technique such as linear regularization (Press et al. [19]). Integral equations can be written and solved for the conditional average of each scalar (temperature, species mass fraction and density) then substituted into reaction rate expressions to provide an approximation to the conditional average of the reaction rates based on the CMC closure hypothesis. The unconditional average of the reaction rates can then be obtained by integrating: idJ\CP(C;Z(x,t),z^(x,t)) * == f *(x,t) = p Jo >- (2.52) p\C thus obtaining closure to the chemical source terms for turbulent combustion modeling (Bushe and Steiner [8]). An alternative approach is to use a decomposition where the conditional average is approximated as a linear combination of basis functions. These basis functions should be chosen such that they exhibit the sharp gradients that are expected in the conditional averages of species mass fractions. Using this method, the conditional average of the temperature would be approximated by: Chapter 2. Background 31 (2.53) m where $ (£) is thera*"basis function and a would be the weighting coefTO m ficient vector for that basis function. The integral equation for temperature could then be re-written using this basis function decomposition: m which could be inverted and solved for the weighting coefficient vector (a ). m Rather than approximating the conditional average at each discrete mixture fraction point, one would solve for the weighting coefficient vector to be multiplied by the basis functions. An obvious choice for the basis functions is a library of steady solutions to the flamelet equations. However, these solutions fail to capture the underlying physics of local ignition or extinction. A better choice for the basis function would be a library of solutions to the unsteady flamelet equations (Pitsch [20]) as these solutions would represent physically realizable conditions for the flame and include ignition and extinction phenomena. The conditional average of each scalar can be approximated using the basis function decomposition and then substituted into the reaction rate expressions to close the chemical source terms using the CMC closure hypothesis (Bushe and Steiner [36]). The use of a library of unsteady solutions to the flamelet equations has become known as the Laminar Flamelet Decomposition (LFD) approach to CSE. If LFD is used, a third approach can be developed as the conditional average of the chemical source terms can be expressed as a linear combination of the source terms from the flamelet library: (2.55) m Chapter 2. Background 32 Rather than solving for a weighting coefficient vector (a ) for each scalar, m one pertinent scalar, such as temperature, would be selected and the integral equation for that scalar would be inverted to obtain a single weighting coefficient vector. The single weighting coefficient vector would then be used to approximate the conditional average of the chemical source terms directly (Eq. 2.55). Additionally, the conditional averages of the species mass fraction would be given by: (2.56) m where Y j (£) is the species basis function from the unsteady flamelet lim brary. Then the unconditional averages of the species mass fractions are obtained without solving transport equations but by integrating (Bushe and Steiner [36]): m Consequently, the chemical source terms can be closed and predictions for pollutant emissions can be obtained. Applications of Conditional Source-term Estimation CSE has been tested in the a priori sense with DNS data and against experimental results in both the RANS and LES paradigms. The first CSE strategy is to invert integral equations for the conditional average of each scalar and then approximate the conditional average of the chemical source terms. This model gave good approximations for the reaction rates and when using a second conditioning variable, could predict ignition and extinction phenomena in piloted flames (Bushe and Steiner [8]). This second conditioning variable was the scalar dissipation of mixture fraction which Chapter 2. Background 33 is not an ideal variable as both ignition and extinction phenomena exhibit time-lag between peaks in scalar dissipation. The CSE method does not require assumptions of fast chemistry or steady-state behaviour for the conditional averages. It was found to provide an acceptable closure model for non-premixed piloted turbulent jet flames working with reduced chemistry (Bushe and Steiner [35]). The testing of this method has also identified the shortcoming that the smoothing matrix used in the linear regularization technique causes the conditional average to be overly diffusive. This overly diffusive behaviour has a negligible effect in simple chemical kinetics but can lead to unacceptable errors in the predictions of the chemical source terms in more complex chemistry. The Laminar Flamelet Decomposition strategy for CSE provides reasonable closures for the chemical source terms. This method was actually shown to improve the predictions over the original CSE method (Bushe and Steiner [36]). As well, using only one scalar to solve for the weighting coefficient vector was found to offer considerable computational savings as individual species balance equations would not be solved and only one integral equation is being inverted resulting in an improvement in computation time. Although still early in its development, CSE has been shown to give good predictions of species mass fractions, temperatures and chemical reaction rates when compared to DNS data. The model assumes that the gradients in the conditional averages are small and presumes the shape of the PDF of mixture fraction. Closure for the chemical source terms is based on the CMC closure hypothesis. Chapter 2. Background 2.6 34 Summary In this chapter, the modeling of turbulent combustion has been reviewed. In particular, the equations of motion, chemical kinetics and the role of turbulence have been introduced. Three turbulent combustion models - Laminar Flamelet model, Conditional Moment Closure, and Conditional Source-term Estimation - were presented and discussed. The Laminar Flamelet Decomposition approach for CSE offers an alternative model for closing the chemical source terms that takes advantage of the CMC closure and advances in flamelet models and achieves a reduction in computation time. As such, further development and testing of the Laminar Flamelet Decomposition approach to CSE is the primary topic of this thesis and presented in the following chapters. 35 Chapter 3 A priori Testing of Conditional Source-term Estimation 3.1 Introduction The purpose of this investigation is to explore the effects of the solution parameters on the CSE method and compare them to DNS data in the a priori sense. In particular, the goals are to: 1. Test additional linear regularization parameters such as higher order smoothing matrices. 2. Perform tests to determine the number of flamelets that should be used in the flamelet library to optimize the computation time, and accuracy of solution. 3. Determine which combination of species mass fraction and temperature yields the best approximation of the reaction rates. 3.2 Formulation The underlying concept of CSE is to approximate conditional averages of the relevant scalars by inverting integral equations for these scalars. In the Laminar Flamelet Decomposition, the conditional average is approximated as a linear combination of basis functions and the integral equations are solved Chapter 3. A priori Testing 36 for the weighting coefficient vector which is then used to obtain closure for the chemical source terms. The weighting coefficient vector is truncated to eliminate negative values then rescaled to return the original values of the scalar. This is necessary because negative values in the weighting function would mean the negative of the forward reaction rate which is not the same as the reverse reaction rate. The solution technique for inverting the integral equation is described below. The integral equation for temperature is: %= f P (Z)T\ZdZ Jo (3.1) 1 n where P (Z) is a short form of the PDF of mixture fraction approximated n using the /3-function (P (C; Z (x, t), Z (x,*))). 712 The integral equation can be re-written as a linear combination of basis functions: Tn~ = £p n (Z) (^2 amTfm (Z)^ dZ (3.2) where T\Z = £ a T m / m (Z) (3.3) m and a is the weighting coefficient vector while T / is the basis function m m composed of unsteady solutions to theflameletequations. Eq. 3.1 is a Fredholm equation of the first kind (Kythe and Puri [38]). Assuming a quadrature of the form: J f(x)dx^ £W f J k k (3.4) Chapter 3. A priori Testing 37 where Wk is a weighting function, Eq. 3.2 is analogous to the matrix equation: (3.5) Ax = C where C= (3.6) \T J n and x is the weighting coefficient vector (3.7) x = ( a i . . . ctm) and A is an n x m matrix that is not necessarily square: I ... ^ P dZ n n (3.8) Tf m The matrix equation (Eq. 3.5) can be solved using a technique known as linear regularization (Press et al. [19] and Kythe and Puri [38]) or nonnegative least squares (Lawson and Hanson [39]). As the matrix equation (Eq. 3.5) is ill-posed, the equation is modified for ease of solution to become: (A A T + 7#) x= AC T (3.9) where yH is introduced to provide smoothing of the weighting coefficient vector (x) based on a priori knowledge of its shape. The weighting of the Chapter 3. A priori Testing 38 smoothing is given by 7 and its recommended form is (Press et al. [19]): Tr(A A) T ^ - . W • (3 - 10) where Tr represents the trace of the matrix and H is given by: H = B B. T (3.11) B is the smoothing matrix obtained by setting the first or higher order derivatives of the weighting coefficient vector to zero depending on the amount of smoothing required through a priori knowledge of the weighting coefficient vector. Once the weighting coefficient vector is determined, the conditional averages of each scalar can be approximated and then used to close the chemical source terms. 3.3 Methodology A MATLAB ([40]) program was developed to optimize the parameters that affect the solution of the Laminar Flamelet Decomposition approach to the CSE method. The results of the MATLAB program are compared to DNS data of turbulent non-premixed combustion that uses a realistic reduced chemical kinetic mechanism for combustion of a nitrogen/methane mixture with a nitrogen/oxygen mixture (Bushe et al. [37]). The kinetic mechanism for this flame is: Fuel + Oxidizer —• Intermediate + Products Intermediate + Oxidizer —• 2Product JV + Oxidizer -* 2NO 2 (3.12) Chapter 3. A priori Testing 39 The DNS data is of a turbulent temporal mixing layer. Its domain is 120 2 x 240 and the conditional averages of the scalars are organized into planes in mixture fraction space 240 points by 50 evenly spaced values of mixture fraction. The Reynolds number in the DNS database was initially around 60 in the coldfluidbut is reduced to 20 in the flame due to higher viscosity. The MATLAB program reads the temperature, mixture fraction, and mixture fraction variance from the DNS data, and presumes the PDF of the mixture fraction takes the shape of a /3-function. The program inverts the integral equations and solves for the weighting coefficient vector then calculates the chemical source terms for each step in the kinetic mechanism. The weighting coefficient.vector is truncated and re-scaled to eliminate negative values as a result of the physical constraint of non-negative reaction rates (a negative forward reaction rate is not the same as the reverse reaction rate). The predicted reaction rates are then integrated and compared to the global reaction rates from the DNS data. The comparison is made across planes of constant x through the mixing layer which are referred to as physical space. The error is calculated as the absolute error between the predicted reaction rate and the DNS reaction rate and is expressed as a percentage. Flamelet libraries were built for these a priori tests from unsteady solutions to theflameletequations and were initialized with the same data as the DNS. The scalar dissipation rate was increased until extinction occurred and the species mass fractions returned to a mixing only state. The flamelets are selected so that the library will contain realizations of a broad scope including burning, extinguishing and mixing. The libraries contained 45 flamelets except those used to determine the optimum number of flamelets. No significance is attached to the actual time or scalar dissipation of the flamelets; the intent is to represent many realizable conditions in the flame. Chapter 3. A priori Testing 3.4 40 Linear Regularization Parameters To test the linear regularization parameters, the MATLAB program is modified to determine the effect of using different smoothing matrices (73) to add varying amounts of smoothing to the system. The following smoothing matrices (Eqs. 3.13, 3.14, 3.15) are evaluated by calculating the respective weighting coefficient vectors and corresponding reaction rates. / - l l 0 .. o \ 0 -1 1 .. 0 0 ... 0 -1 1 / \ / -1 2 1 o\ 0 0 - 1 2 1 B = (3.13) 0 (3.14) 2 0 - 1 2 1 / / -1 B = 3 -3 1 o\ 0 0 - 1 3 - 3 1 0 3 0 -1 3 -3 (3.15) 1/ Thefirstmatrix (Eq. 3.13) sets thefirstderivatives to be zero and forces the solution to be more smooth. The next two matrices (Eqs. 3.14, 3.15) set the second and third derivatives to be zero respectively and enforce less restrictive smoothing. Each smoothing matrix yields a different H matrix and therefore a different character of smoothing to the solution of the matrix equation. H\ = Bi (3.16) Chapter 3. A priori Testing 41 H = BlB 2 (3.17) 2 H = B%B 3 (3.18) 3 However, the matrix H (Eq. 3.11) can also be a function of the H's determined from the smoothing matrices. For example: Hi - 0.5#i + 0.25H + 0.25# 2 (3.19) 3 where Hi is from Eq. 3.16, H is from Eq. 3.17, and H is from Eq. 3.18. 2 3.4.1 3 Results of Linear Regularization Parameters The reaction rates for the three kinetic steps are calculated for four different H matrices corresponding to the three B smoothing matrices and the linear combination of matrices (Eq. 3.19) then compared to the DNS database. The simplest H matrix (Hi) provides the smoothest weighting coefficient vector followed by the linear combination of H's (dominated by Hi) then H and Hz (Figure 3.1). The weighting coefficient vectors from H and H3 2 2 have a considerable number of values truncated to eliminate negative values. Truncating the weighting coefficient vectors will likely adversely affect the prediction of the chemical source terms. Figure 3.2 shows a comparison of the actual conditional average of the chemical reaction rates obtained from the DNS database for thefirstkinetic reaction (Eq. 3.12) and those calculated using the different smoothing matrices in physical space. The matrix that provides the smoothest weighting Chapter 3. A priori Testing 42 Flamelet number Figure 3.1: Weighting coefficient vectors for each H matrix of the linear regularization parameters coefficient vector (H\) also provides the best approximations for predicting the DNS reaction rate. As well, the linear combination of i/'s yields reasonable approximations of the reaction rates but this is due to the strong dependence on Hi. The matrices Hi and if3 provide poor approximations of the reaction rate since both greatly under-predict the rate. Comparing the predicted reaction rates and the DNS data for the second reaction (Figure 3.3) shows that Hi over-estimates the reaction rate. Hi under-predicts the reaction rates while H3 and Hi provide reasonable approximations to the reaction rate. The comparison between the predicted reaction rates and the DNS database for the third kinetic reaction are shown in Figure 3.4. This time, Hi and Hz over-predict the reaction rates compared to Hi and H4 which provide reasonable approximations. Chapter 3. A priori Testing 43 Physical space Figure 3.2: Reaction rate 1 predictions for different smoothing matrices in the linear regularization Physical space Figure 3.3: Reaction rate 2 predictions for different smoothing matrices in the linear regularization Chapter 3. A priori Testing 44 18 DNS H1 — H2 H3 — H4 -2 0 2 4 6 8 10 12 14 16 Physical space Figure 3.4: Reaction rate 3 predictions for different smoothing matrices in the linear regularization Integrating the predicted reaction rates across the mixing layer gives the global reaction rates for comparison. Table 3.1 shows the percent error between the predicted and the actual global reaction rates. Overall, H\ provides a reasonable estimate for the first reaction rate but H4 provides better approximations for reaction rates 2 and 3. Optimizing the linear combination of the smoothing matrices has the potential to improve the predictions of the reaction rates. As the solution is allowed to become less smooth, the prediction of the reaction rates becomes less reasonable. Physically, the weighting coefficient vector should be a smooth function; there should be a smooth transition between physical conditions in the flame. It is therefore reasonable that the parameters that yield the smoothest weighting coefficient vector also provide the best predictions of the reaction rates. As values are truncated from the weighting coefficient vector, it becomes less smooth resulting in poor Chapter 3. A priori Testing Matrix Reaction 1 Reaction 2 (% e r r o r ) (% e r r o r ) (% e r r o r ) Hi 3.15 26.34 15.17 H 76.69 55.17 50.98 H 50.99 11.97 47.58 H 9.42 10.820 2.22 2 3 4 45 Reaction 3 Table 3.1: Comparison of global reaction rates for different smoothing matrices predictions of the reaction rates. 3.5 Optimum Number of Flamelets The number of flamelets in the library for the Laminar Flamelet Decomposition will affect the prediction of the chemical source terms. Too small a selection of flamelets will not adequately capture sufficient physically realizable conditions in the flame while moreflameletsthan necessary will adversely affect the computation time. As such, it is necessary to determine an optimum number offlameletsto use in the library while sufficiently describing the physical conditions in the flame. To determine the optimum number offlameletsnecessary to reasonably estimate the chemical source terms, six flamelet libraries were built consisting of 5, 10, 25, 50, 75 and 100 unsteady solutions to theflameletequations (Eqs. 2.46 and 2.47). The libraries were initialized to match the initial conditions of the DNS database and then the scalar dissipation rate was increased until extinction occurred and the system returned to the pure mixing condition. Each library containsflameletsbetween burning and extinction but limited by the number offlamelets.The MATLAB program was modified to test the number offlameletsand the corresponding predicted chemical reaction rates are compared to the DNS data. Chapter 3. A priori Testing 46 To test for convergence of statistics, equally spaced cells were removed from the DNS data for temperature, mixture fraction and variance then used in the M A T L A B program to predict the three reaction rates. The number of cells tested are 60, 80, 120 and 240 (the maximum number of cells in the DNS database). 3.5.1 Results of Optimum Number of Flamelets The six temperature libraries are shown in Figure 3.5 indicating significant discontinuities are present in the 5, 10 and 25 flamelet libraries, particularly near extinction; the higher flamelet numbers. (d) Mixture fraction (•) Mixture fraction (f) Mixture fraction Figure 3.5: Surface plots of temperature flamelet libraries shown in order of increasing scalar dissipation: (a) 5 flamelets; (b) 10 flamelets; (c) 25 flamelets; (d) 50 flamelets; (e) 75 flamelets; (f) 100 flamelets The resulting weighting coefficient vectors are shown in Figure 3.6. The Chapter 3. A priori Testing 47 lines are scaled such that the markers indicate comparable flamelet numbers used in each library. The libraries containing 5 and 10 flamelets may may not contain sufficient flamelets to adequately represent the conditions in the flame. As the number of flamelets increases, the magnitude of each coefficient decreases. Flamelet number Figure 3.6: Weighting coefficient vectors for the six flamelet libraries The predicted chemical source terms based on each flamelet library are calculated using the M A T L A B program and are compared to the source terms from the DNS database. Figure 3.7 shows the predicted chemical source terms for the first reaction in the kinetic mechanism. It is shown that the libraries with 10 and 25 flamelets over-predict the reaction rates while the libraries with 50, 75 and 100 flamelets yield reasonable predictions of the reaction rate. The library with only 5 flamelets under-predicts the reaction rate but is a reasonable approximation. Figure 3.8 shows the predictions from each flamelet library for the second step in the kinetic mechanism. All of the libraries considerably over-predict Chapter 3. A priori Testing 48 DNS — — 5 flamelets 10 flamelets 25 flamelets — 50 flamelets — 75 flamelets - 100 flamelets 0 2 4 6 8 10 12 14 Physical space 16 Figure 3.7: Reaction rate 1 predictions for different flamelet libraries the reaction rate. The library with 50 flamelets provides the best prediction when compared to the other libraries. The over-prediction of the chemical source terms is likely due to the fact that the rates are sensitive to the intermediate mass fraction and the presence of extinction phenomena (Bushe and Steiner [36]). All six flamelet libraries are shown to under-predict the formation of oxides of nitrogen in Figure 3.9. Again, the library with 5 flamelets provides a reasonable approximation to the reaction rate. There is a considerable difference again between the libraries with 50, 75 and 100 flamelets compared to the predictions from the libraries with 10 and 25 flamelets. The libraries with more flamelets tend to capture more of the physical conditions in the flame and thus provide better predictions of the chemical source terms. It has been previously shown that both the C M C and C S E methods under-predict the formation of oxides of nitrogen in the presence of extinction phenomena (Bushe and Steiner [36] and Bushe et al. [37]). Chapter 3. A priori Testing 49 Figure 3.8: Reaction rate 2 predictions for different flamelet libraries Comparing percent error in the global reaction rates (Table 3.2) shows that the library with 50 flamelets provides the best overall predictions for the three reaction rates. In this case, all the libraries significantly over- predicted the second reaction rate, but the libraries with fewer flamelets over-predicted the reaction rate by over three times the value of the library with 50 flamelets. No improvement in the reaction rate predictions is observed when the number of flamelets is increased from 50 to 75 while the predictions obtained using 100 flamelets are slightly worse. There can be orders of magnitude difference between reaction rates for burning and mixing flamelets. As the magnitude of the weighting coefficient vector decreases with increasing number of flamelets, it may be possible that the smaller magnitude of the weighting coefficient vector coupled with the higher reaction rates could cause the accuracy of the chemical source terms predictions to diminish. Chapter 3. A priori Testing 50 Physical space Figure 3.9: Reaction rate 3 predictions for different flamelet libraries N u m b e r of flamelets Reaction 1 Reaction 2 (% error) (% error) (% error) 5 10.45 894.7 30.66 10 20.62 893.3 46.91 25 23.28 600.4 48.13 50 0.29 260.3 36.77 75 1.45 283.6 39.26 100 6.13 359.2 37.53 Reaction 3 Table 3.2: Comparison of global reaction rates for libraries of various sizes Chapter 3. A priori Testing 51 Based on thesefindings,the optimum number offlameletsto include in a flamelet library is around 50 provided that the conditions captured by the library closely match those to be simulated. Convergence of statistics were tested by calculating the global reaction rates based on modifying the DNS data to contain 60, 80, 120 and 240 cells. Comparisons of the percent error between the predicted global reaction rate and thefirstDNS reaction rate for the converged statistics are presented in Table 3.3 while additional converged statistics can be found in Appendix A. These tables demonstrate that the error is minimized when using 50 flamelets, regardless of the size of the domain. Number of flamelets 60 cells 80 cells 120 cells 240 cells (% error) (% error) (% error) (% error) 5 10.46 10.46 10.45 10.45 10 20.49 20.49 20.54 20.62 25 24.22 24.00 23.78 23.28 50 0.26 0.27 0.28 0.29 75 1.37 1.39 1.41 1.45 100 5.96 6.00 6.04 6.13 Table 3.3: Comparison of global reaction rates for converged statistics of reaction rate 1 3.6 Optimum Scalar for Decomposition In the initial examination of Laminar Flamelet Decomposition for CSE, only one scalar (temperature) was chosen and its integral equation was inverted for the weighting coefficient vector (Bushe and Steiner [36]). To determine which scalar provides the best approximation for the reaction rates, the integral equation for each scalar is inverted and used to solve for a different weighting coefficient vector to estimate the reaction rates. The predicted Chapter 3. A priori Testing 52 chemical source terms are then compared to the DNS reaction rates. In addition, integral equations for combinations of scalars are also evaluated in an effort to improve the predictions of the reaction rates. To create combinations of scalars, integral equations are written for each scalar to be used in the combination. If two scalars are used in the combination, the matrix equation (Eq.3.5) is rewritten such that: (3.20) where T and $ represent two different scalars. The weighting coefficient vector is: X = (ai . . . (3.21) dm) and A is now a 2n x m matrix: ( ... ^ ( P dZ n n \ T/m A= (3.22) P dZ n n where T / and ^ff represent the corresponding flamelet libraries for the m scalars T and m The new matrix equation from the combination of scalars is solved using the linear regularization technique. Chapter 3. A priori Testing 3.6.1 53 Results of Optimum Scalar for Decomposition The scalars evaluated are the five species mass fractions (Fuel, Oxidizer, Intermediate, Product and Nitric Oxide (NO)) appearing in the kinetic mechanism and the temperature. The weighting coefficient vectors for the six scalars can be seen in Figure 3.10. The weighting coefficient vectors for the temperature, fuel and product scalars are nearly identical whereas the other scalars (oxidizer, intermediate and NO) are considerably different. 0.15, 1 1 1 1 , 1 —, 1 , — Temperature — Fuel — Oxidizer -0.151- -0.21 0 ' 5 1 10 1 15 ' 1 20 ' 25 1 30 1 35 40 45 Flamelet number Figure 3.10: Weighting coefficient vectors for the six scalars Figure 3.11 shows a comparison of the predicted reaction rates for the first reaction in the kinetic mechanism. Using the weighting coefficient vector derived from the oxidizer greatly over-predicts the reaction rates. As the weighting coefficient vectors from the temperature, fuel and product exhibit similar shapes, their reaction rate predictions are also similar. Additionally, these three scalars provide the best approximations to the reaction rate. The intermediate and NO scalars tended to reasonably capture the reaction rate but provide an under-prediction. Chapter 3. A priori Testing 2.5 X10'r 54 5 DNS — — Temperature Fuel Oxidizer - Product — Intermediate 2 - NO 0 2 4 6 8 10 12 14 16 Physical space Figure 3.11: Reaction rate 1 predictions for different scalars In the comparison between the predicted reaction rates and the DNS database for the second reaction, the oxidizer again provides a significant over-prediction (Figure 3.12). The three similar weighting coefficient vectors derived from the temperature, fuel and product scalars also over-estimate the reaction rate. The NO and intermediate scalar yield reasonable predictions of the reaction rate. Comparisons of the predictions for the third reaction can be seen in Figure 3.13. The temperature, fuel, product and intermediate scalars provide reasonable approximations for the formation of oxides of nitrogen. The oxidizer and the N O scalar significantly under-predict the reaction rate. The reaction rates are integrated over the mixing plane to determine the global reaction rates and the percent error between the DNS database and the predictions can be seen in Table 3.4. The temperature scalar provides the best prediction for the first reaction rate followed by the fuel and prod- Chapter 3. A priori Testing Physical space Figure 3.12: Reaction rate 2 predictions for different scalars Physical space Figure 3.13: Reaction rate 3 predictions for different scalars Chapter 3. A priori Testing 56 uct scalars. The fuel and NO scalars provide the best predictions for the second reaction and the product provides the best predictions for the third reaction. The temperature, fuel, product and intermediate scalars provide predictions that are within 30% of the global reaction rates. The oxidizer does not provide acceptable predictions for the reaction rates. The NO provides acceptable predictions for the first two reactions but provides a poor prediction for the third reaction. Scalar Reaction 1 Reaction 2 Reaction 3 (% e r r o r ) (% e r r o r ) (% e r r o r ) Temperature 3.15 26.35 15.17 Fuel 4.78 13.7 17.48 Oxidizer 187.7 300.5 34.27 Product 4.86 23.10 8.35 Intermediate 18.36 14.32 13.87 NO 15.16 9.13 42.27 Table 3.4: Comparison of global reaction rates for different scalars Using combinations of scalars to solve for a weighting coefficient vector and close the chemical source terms may improve the reaction rate predictions. Six combinations of scalars are tested and presented here while the results of additional combinations can be found in Appendix B. The six combinations of scalars include: temperature and fuel (TF); temperature and oxidizer (TOx); product and intermediate (PI); fuel and intermediate (FI); temperature, product and intermediate (TPI); and temperature, fuel and product (TFP). Comparing the reaction rates predicted by the combinations of scalars for the first reaction can be seen in Figure 3.14. Similar to the predictions using only the oxidizer scalar, the combination of temperature and oxidizer over-predict the reaction rate. The combinations of temperature and fuel, Chapter 3. A priori Testing 57 and temperature, fuel and product give the best overall predictions. The other combinations tend to under-predict the reaction rates. Physical space Figure 3.14: Reaction rate 1 predictions for combinations of scalars Figure 3.15 shows the predictions for the second reaction rate compared to the DNS database. Again, the combination of temperature and oxidizer over-predict the reaction rate. The remaining combinations provide fair predictions of the reaction rate. A comparison between the DNS reaction rate and the predicted rates from the combinations of scalars for the third reaction is shown in Figure 3.16. This time, the combination of temperature and oxidizer under-predict the reaction rate. Similarly, the combinations of temperature and fuel, and temperature, fuel and product also under-predict the reaction rate. The combination of fuel and intermediate provide a good approximation. The remaining combinations provide fair predictions. Chapter 3. A priori Testing Figure 3.15: Reaction rate 2 predictions for combinations of scalars Physical space Figure 3.16: Reaction rate 3 predictions for combinations of scalars 58 Chapter 3. A priori Testing 59 Table 3.5 shows the percent error between the predicted and DNS global reaction rates. Overall, the combination of temperature, fuel and product provides the best approximations for all three reactions closely followed by the combination of temperature and fuel. The combination of temperature and oxidizer yields the worst predictions for the reaction rates, particularly the second reaction. The other combinations provide reasonable approximations for the reaction rates while the combination of fuel and intermediate provide the best prediction for the formation of oxides of nitrogen. Combination Reaction 3 Reaction 1 Reaction 2 (% error) (% error) (% error) TF 1.76 18.33 15.51 TOx 65.52 107.1 63.39 PI 27.30 23.55 20.21 FI 27.45 22.59 3.51 TPI 23.04 18.35 13.12 TFP 0.53 19.18 13.72 Table 3.5: Comparison of global reaction rates for combinations of scalars While improvements are made in the predictions of the global reaction rates, the additional coding and computation time involved with using combinations of scalars may not justify their use compared to using a single scalar. The above results show that inverting the integral equation for the single temperature scalar and solving for the weighting coefficient vector provides predictions that are comparable to the predictions offered by using combinations of scalars. All of these results are obtained for a plane in the DNS database when all of the scalars are present to some degree; in a combustion simulation, it is possible to have only some of the scalars present (e.g. fuel is absent before the start of injection). In that case, using the fuel scalar to solve Chapter 3. A priori Testing 60 for the weighting coefficient vector would provide a meaningless solution, therefore, it is better to invert for a scalar that is always present such as the temperature. 3.7 Conclusions and Recommendations In this chapter the effects of various parameters affecting the predictions from the CSE method using Laminar Flamelet Decomposition are investigated and has led to the following observations: 1. The lowest order smoothing matrix provides good predictions of the reaction rates compared to the DNS data. Higher order smoothing matrices can be used to predict the weighting coefficient vectors but require more truncation which leads to weaker predictions for the reaction rates. Linear combinations of smoothing matrices show promise for improving the reaction rate predictions and more work should be performed in an attempt to optimize the smoothing matrix. 2. Approximately 50 flamelets are necessary to sufficiently describe the physical conditions in the flame and provide reasonable predictions of the reaction rates. Comparing convergence of statistics from different grid spacings demonstrated that 50 flamelets was sufficient regardless of grid spacing. Fewer than 50 flamelets failed to adequately capture the physically realizable conditions in the flame although occasionally good predictions are obtained. The magnitude of the coefficients decreases as the number of flamelets in the library increases which led to high reaction rate flamelets having a greater influence on the predictions. 3. Inverting each scalar and solving for a weighting coefficient vector to predict the reaction rates found that the temperature, fuel and product provide good approximations. The intermediate scalar provides reasonable approximations for the global reaction rates. The oxidizer Chapter 3. A priori Testing 61 scalar offers the worst predictions for closing the chemical source terms. The NO scalar gave reasonable predictions for thefirsttwo reactions in the kinetic mechanism but gave a poor prediction for the third reaction rate. 4. The predictions from the combination of scalars tend to be an average of the predictions from the individual scalars. Combinations involving the individual scalars that provided good predictions (i.e. temperature, fuel and product) gave better predictions of the reaction rates. From the above investigations, it is recommended that the simplest smoothing matrix be used in conjunction with linear regularization while truncating and re-scaling the weighting coefficient vector. Theflameletlibrary should contain approximately 50 unsteady solutions to the flamelet equations to adequately describe the physical conditions in theflame.Finally, inverting the integral equation for the single temperature scalar provides good predictions for obtaining closure to the chemical source terms; temperature is the most convenient scalar to use as other scalars may not always be present in turbulent combustion simulations. 62 Chapter 4 Implementation of Conditional Source-term Estimation 4.1 Introduction This chapter describes the implementation of the CSE method (Bushe and Steiner [8],[35],[36]) into a commercial combustion package for the prediction of methane ignition under engine-relevant conditions. A brief introduction to the KIVA-II (Amsden et al. [41]) combustion software package is made followed by a description of the implementation of the CSE method. Finally, simulation results of methane ignition are presented and discussed. 4.2 KIVA-II - Multi-Dimensional Model for Reactive Flows KIVA-II (Amsden et al.[41]) is a commercially available, multi-dimensional model for chemically reactive flows specifically developed for engine applications. It uses a general numerical solution procedure capable of handling turbulent or laminar, subsonic or supersonic, and single or multiphase flows. Even though written primarily for internal combustion engines, the initial and boundary conditions, and mesh generation logic can be modified for other applications. The governing equations are discretized in both space and time and are solved using a finite volume method called the Arbitrary Chapter 4. Implementation 63 Lagrangian-Eulerian (ALE) method (Pracht [42]). An automated mesh generation tool is included with KrVA-II that can create a variety of piston and head shapes for engine applications. Incorporated with this tool are boundary conditions for inflow, outflow, wall, symmetry and periodic boundary conditions. Additionally, the moving mesh capabilities of the program make simulating moving pistons and valves possible. KIVA-II is a RANS model that uses the k — e approach to turbulence modeling by solving additional transport equations for the turbulent kinetic energy (Eq. 2.23) and its dissipation (Eq. 2.24). The scaling constants in the k — e transport equations are determined from experiments and theoretical considerations and can be modified by the user. KIVA-II can account for an arbitrary number of species and chemical reactions while distinguishing between slow reactions assumed to proceed kinetically and fast reactions assumed to be in equilibrium. The rate expressions for the kinetic reactions are assumed to take the Arrhenius form. The reaction progress rates are calculated assuming that every participating species is inert or appears on only one side of the reaction. It is assumed that the fuel species does not participate in any of the equilibrium reactions. KIVA-II calculates the change in the progress variable based on the forward and backward reaction rate coefficients and then predicts the individual species density and specific internal energy. The source code is included with KIVA-II which allows the user to modify or replace the various sub-models. As such, it is relatively straight-forward to implement improved models, such as CSE, to close the chemical source terms in the turbulent transport equations. Chapter 4. Implementation 4.2.1 64 Jet Simulations with K I V A A submerged jet is one that is spreading through a medium at rest and takes the form shown in Figure 4.1 (Abramovich [43]). The jet is composed of an initial region where the velocity is uniform followed by a transition region then the fully developed region. The fully developed region extends beyond approximately 20 nozzle diameters from the nozzle and the jet has the property of self-similarity, meaning that the non-dimensional velocity profile is independent of the distance from the nozzle. In unsteady jets, the tip is often treated as a vortex ball in front of a steady-state jet. Figure 4.1: Unsteady round jet Recently it was shown that the KIVA-II code could give nearly gridindependent solutions of unsteady gaseous jets when the grid size in the nozzle area was equal to \ the nozzle diameter and good results were obtained when the cell size was about ^ the nozzle diameter (Ouellette and Hill [44]). Chapter 4. Implementation 4.3 65 Implementation of Conditional Source-term Estimation The Laminar Flamelet Decomposition approach for CSE offers an improvement in the closure of the chemical source terms compared to other modeling strategies (Bushe and Steiner [36]). From the a priori testing discussed in Chapter 3, the Laminar Flamelet Decomposition approach for the CSE method will be implemented into KIVA-II. The unsteady flamelet library is built with 55 flamelets, and the integral equation for temperature is inverted and used to solve for the weighting coefficient vector. The lowest order smoothing matrix is used with the linear regularization technique along with truncating and re-scaling the weighting coefficient vector. Methane ignition delay times from the implementation will be tested against non-premixed shock tube experiments (Iaconis [45]). The following sections describe the various steps for the implementation and testing of the model. 4.3.1 Numerical Model For the purpose of validation and application to diesel engine-relevant conditions, the model is used to simulate the injection of methane into a hot cylinder similar to the experimental shock tube (Iaconis [45]). From symmetry, the three-dimensional, unsteady jet is modeled as a two-dimensional axisymmetric jet with a domain that is 50 mm in the radial direction and 100 mm in the axial direction. The injection is modeled using a pressure boundary inlet condition on a zone 0.5 mm wide at the top, left edge of the domain. This inlet condition is equivalent to the nozzle radius of the fuel injector used in shock tube experiments at the University of British Columbia (Iaconis [45]). A simple, cylindrical, structured mesh is constructed for the simulations (Figure 4.2) with the jet situated in the top left corner and walls placed suf- Chapter 4. Implementation 66 ficiently far from the jet so as to not interfere with its propagation. The left wall is set as an axis of symmetry boundary condition. The mesh is stretched in both directions; away from the axis of symmetry and away from the top boundary wall. Variations of the mesh include using one, two and five cells across the inlet boundary with corresponding refinements to the remaining mesh. Inlet and boundary conditions are specified to correspond to the particular experiments used for validation. 1 jijif ri Ififf T f f r EE Birr I'll II \ 1 : r ' 0 1 2 3 4 5 X Figure 4.2: Two dimensional mesh The C S E model is used to replace the existing ignition mechanism in KIVA-II. The weighting coefficient vector is calculated from inverting the integral equation for temperature and then used to calculate the chemical source terms (using the temperature source term as the basis function). The new source terms from the C S E method are then used to calculate the species density and specific internal energy for KIVA-II. The weighting coefficient vector is calculated at each time step for each row in the domain Chapter 4. Implementation 67 but does not account for its evolution with time and space. Essentially, the weighting coefficient vector for each row is calculated independently of its previous and neighbouring values. KIVA-II does not solve transport equations for the mixture fraction or its variance, therefore solution of these equations (Eqs. 2.30 and 2.33) were included with the CSE implementation. Validation of the mixture fraction and mixture fraction variance is discussed in Appendix C. The PDF of mixture fraction was implemented using a numerical algorithm to describe the internal structure of the flame and turbulent mixing. 4.3.2 Unsteady Flamelet Libraries The Laminar Flamelet Decomposition approach to CSE uses a library of solutions to the unsteady flamelet equations (Eqs. 2.46 and 2.47), rewritten here: Using a library of unsteady solutions to the flamelet equations is an innovative approach; to date most researchers have solved the unsteady flamelet equations in parallel using input from the flow calculations to solve the unsteadyflameletequations. In CSE, the conditional averages of the scalars are assumed to be a linear combination of basis functions and the library of unsteady solutions provides these basis functions. The advantage to using unsteady solutions is that they provide physically realizable conditions for theflameincluding local ignition and extinction phenomena. Chapter 4. Implementation 68 The kinetic mechanism used to develop the flamelet library is based on shock tube studies of premixed methane ignition under engine-relevant conditions (Huang et al. [46]). The mechanism contains 38 species undergoing 192 reactions and achieved a good agreement with experimental ignition delay data including the reproduction of the reversed-'S'-shape characteristic of the data. Using this mechanism, the flamelet library contains 39 scalars; 38 species plus temperature. Two sets offlameletlibraries were built and their ability to predict the ignition of methane under engine-relevant conditions is tested. The first set of libraries is referred to as Pull Flamelet Libraries while the second is referred to as Igniting Flamelet Libraries. The Full Flamelet Libraries are composed of approximately equal numbers of mixing, igniting, burning and extinguishingflameletsand are meant to describe the entire history of a flame. The Igniting Flamelet Library differs from the Full Flamelet Library in that it is composed of predominantly ignitingflameletswith relatively few burning flamelets. The initial condition of pure mixing is used to initialize the flamelet equations for both libraries. The scalar dissipation is increased exponentially with time until extinction occurred and the species mass fractions return to an unreacted, mixing only state. The libraries consist of 55 flamelets placed in chronological order encompassing realizations from a broad scope including mixing, ignition, burning and extinction. There is no significance attached to the actual time or scalar dissipation and theflameletsare referred to only by theirflameletnumber; the libraries form the basis functions to be used in the decomposition. Each library contains two parts: the temperature and the source term component. The temperature component forms the basis function for inverting and solving the integral equation for the weighting coefficient vector. The Chapter 4. Implementation 69 source term component is developed from the source term in the flamelet temperature equation (Eq. 2.47) and forms the basis function to close the chemical source terms (Eq. 2.55). The temperature is used to characterize theflameletsin the library with the initial conditions for theflameletsmatched to the shock-tube studies of temperature sensitivity (Iaconis [45]). Each set of libraries consists of five flamelet libraries, one for each initial condition (1000 - 1400 K). Figure 4.3 presents the temperature and source-term libraries for both the Full Flamelet and Igniting Flamelet Libraries for the initial temperature of 1200 K. Considerable differences can be seen between the different libraries. For the temperature libraries, there are significantly more burning flamelets and fewer ignitingflameletsin the Full Flamelet Library compared to the Igniting Flamelet Library. Accordingly, there are moreflameletswith higher source-terms in the reaction rate portion of the Full Flamelet Library. 4.3.3 Shock Tube Simulations of Methane Ignition Fundamental combustion experiments have been performed for premixed (Huang [47]) and non-premixed (Iaconis [45]) systems using a shock tube to study the ignition of methane under diesel engine-relevant conditions. In the non-premixed studies, a Westport Innovation J41 injector is used to inject fuel into the shock tube after the shock has reflected from the end-plate. Ignition is detected using light emission signals observed through a 10 mm diameter quartz window located in the shock tube end-plate. More recent experiments have provided images through a 200 x 20 mm quartz window parallel to the jet axis. Light emission is again used to detect ignition but the images encompass a region from the tip of the nozzle to approximately 2000 nozzle diameters downstream (Sullivan and Huang [48]). 70 Chapter 4. Implementation 10 (a) 20 30 40 50 Mixture fraction 10 (c) 20 30 10 „ 40 Mixture fraction (b) 50 10 (d) 20 30 40 50 Mixture fraction 20 30 40 50 Mixture fraction Figure 4.3: Comparison of full and igniting flamelet libraries (a) Full flamelet temperature library (b) Igniting flamelet temperature library (c) Full flamelet temperature source-term library (d) Igniting flamelet temperature source-term library Chapter 4. Implementation 71 These non-premixed shock tube experiments are used to test and validate the CSE implementation in KIVA-II and its ability to predict ignition trends with respect to different initial temperatures. The initial conditions are matched to the non-premixed experimental conditions: the injection duration is set to 1.0 milliseconds (ms), the chamber pressure is set to 30 atmospheres (atm), the injection pressure ratio is 3:1 and the temperature is varied between 1000 and 1400 Kelvin (K). The simulations are limited to the onset of ignition and two different conditions are used to identify ignition in the simulations. The first condition is the temperature in a control volume reaches the ignition limit of 2000 K as proposed by Bi and Agrawal ([49]). The second ignition criterion uses a normalized consumed fuel mass fraction: Y = c F m g " Y f (4.1) where Yf is the average fuel mass fraction in the control volume and Y i m X is the average fuel mass fraction for pure mixing in the control volume and is given by: Y = mix I' P(C;x,t)Y<,dZ Jo (4.2) where Yo represents the pure mixing line for the fuel mass fraction as a function of mixture fraction. Ignition is said to have occurred when this normalized consumed fuel mass fraction reaches a value of one (the fuel is completely consumed) in a control volume. Chapter 4. Implementation 4.4 72 Results and Discussion T h e Laminar Flamelet Decomposition approach for C S E is implemented into the KrVA-II code and the various steps are validated independently with experimental data as presented in the following sections. T h e code was tested for grid convergence and is discussed in Appendix D . 4.4.1 Shock Tube Simulations of Methane Ignition T h e ignition delay time of methane is strongly affected by the initial temperature (Iaconis [45]). A s such, simulations are performed using the C S E implemented code in K I V A - I I to predict the ignition of methane under enginerelevant conditions for initial temperatures between 1000 and 1400 K . Full Flamelet Libraries T h e Full Flamelet Libraries are used to calculate methane ignition for different temperatures. Temperature contours at four different times for the initial temperature of 1200 K are shown in Figure 4.4. Based on the first ignition criteria (control volume temperature reaching 2000 K ) , the 1200 K initial temperature case begins to ignite near the tip of the jet 0.35 ms after the start of injection. A comparison between the predicted ignition delay times to experimental data (Iaconis [45]) is presented in Figure 4.5. T h e Full Flamelet Libraries greatly under-predict the ignition delay time although the same trend of decreasing ignition delay time with increasing temperature is exhibited. However, the differences between the predicted ignition delay times (0.05 ms) is small for the change in initial temperature (100 K ) and consequently may not be significant. T h e slope of the ignition delay curve is representative of the activation energy. Performing a linear regression analysis of the data, the slope of the experimental curve is about 3060 K while the slope of the Chapter 4. Implementation 73 Figure 4.4: Temperature contours for the initial temperature of 1200 K using the full flamelet library (a) 0.05 ms after SOI (b) 0.15 ms after SOI (c) 0.25 ms after SOI (d) 0.35 ms after SOI simulations with the full flamelet library is about 570 K . To investigate the small change in the ignition delay times for changes in temperature, the weighting coefficient vector for the plane where ignition occurs is examined (Figure 4.6). The weighting coefficient vector indicates that the solution is selecting predominantly burning and extinguishing flamelets from the library (flamelet numbers 30 - 40 see Figure 4.3) which tend to have higher reaction rates compared to the igniting flamelets. Comparing the chemical source terms for an igniting flamelet (flamelet 15) to an extinguishing flamelet (flamelet 45) (Figure 4.7) shows that the source term from the extinguishing flamelet is higher than that of the igniting flamelet. The CSE solution method determines the weighting coefficient vector based on using all the flamelets in the basis function. As a result, the higher reaction rates from the extinguishing flamelets are influencing the basis function resulting in the over-prediction of the chemical source terms and the correspondingly short ignition delay times. Developing new flamelet libraries Chapter 4. Implementation 74 1000/Temperature (1000/K) Figure 4.5: Comparison between predicted and experimental ignition delay for full flamelet libraries consisting of predominantly igniting flamelets may provide better predictions for the ignition delay time as they would not be as influenced by the higher chemical source terms of the extinguishing flamelets. Igniting Flamelet Libraries Temperatures at successive times for the Igniting Flamelet Libraries for the initial temperature of 1200 K are presented in Figure 4.8. A hot spot, or ignition kernel, first appears at the trailing edge of the tip vortex and is then transported downstream with the jet. Similar images for the other initial temperatures can be found in Appendix E . T h e ignition kernel appears in a well mixed, but slightly fuel rich region. For the higher initial temperatures, the ignition site appears sooner and is closer to the nozzle. T h i s effect is i n agreement with combustion and jet theory as the fuel should ignite sooner for higher initial temperatures and is therefore closer to the nozzle as the jet has not penetrated as far into the chamber. Chapter 4. Implementation 75 0.03 Flamelet number Figure 4.6: Full flamelet library weighting coefficient vector for ignition plane — Flamelet 15 — Flamelet 45 3 2.5 2 - Mixture fraction Figure 4.7: Comparison of source terms for igniting and flamelets extinguishing Chapter 4. Implementation 0.5 (c) 1 1.5 2 Radial distance (cm) 2.5 0.5 (d) 1 1.5 76 2 2.5 Radial distance (cm) Figure 4.8: Temperature contours at different times during methane ignition, (a) 0.10 ms after SOI (b) 0.30 ms after SOI (c) 0.50 ms after SOI (d) 0.70 ms after SOI Examining the weighting coefficient vector (Figure 4.9) for the ignition plane in the 1200 K case shows that the solution is selecting mostly igniting flamelets (see Figure 4.3). As a result, the solutions predicted using the Igniting Flamelet Libraries are less dominated by the high reaction rates exhibited by extinguishing flamelets and show better agreement with the experimental data as is discussed in greater detail in the following section. From here on the Igniting Flamelet Libraries are used to predict the ignition of methane. 4.4.2 Temperature Dependence of Methane Ignition Five simulations to determine the sensitivity of methane ignition delay time to initial temperature (1000, 1100, 1200, 1300, 1400 K) are studied. Two ignition criteria are considered as discussed in Section 4.3.3. Ignition criterion 1 is determined using the control volume temperature reaching 2000 K Chapter 4. 78 Implementation Experimental -•©- Criterion 1 ~o- Criterion 2 I i 07 0 75 i 0.8 i 0.85 i 0.9 i 0.95 i 1 I 1.05 1000/Temperature (1000/K) Figure 4.10: Comparison between predicted and experimental ignition delay for igniting flamelet libraries slope of 3060 K . The slopes of the predicted ignition delay curves using the Igniting Flamelet Library are closer to the experimental results compared to the predictions from the Full Flamelet Library. The most recent experiments have produced images where the perspective is perpendicular to the axis of the jet and use light emission to determine the ignition delay time. Images from these experiments are shown in Figure 4.11. Comparable images from simulation results showing temperature contours are presented in Figure 4.12 while Figure 4.13 shows the reaction rate contours. The experimental images are two-dimensional images of threedimensional phenomena; so far, it cannot be ascertained whether or not the experimental images indicate ignition at the tip of the vortex or along the side of the jet. Although no correlation has been made between the recorded light emission and temperature, several qualitative comparisons can be made. Chapter 4. Implementation 79 Figure 4.11: Light emission images from shock tube experiments for 1475 K (a) 0.95 ms after SOI (b) 1.18 ms after SOI (c) 1.30 ms after SOI (d) 1.42 ms after SOI (e) 1.66 ms after SOI Figure 4.12: Temperature contours for 1475 K: (a) 0.05 ms after SOI (b) 0.25 ms after SOI (c) 0.50 ms after SOI (d) 0.90 ms after SOI Chapter 4. Implementation (c) Radial distance (cm) (d) 80 Radial distance (cm) Figure 4.13: Reaction rate contours for 1475 K: (a) 0.05 ms after SOI (b) 0.25 ms after SOI (c) 0.50 ms after SOI (d) 0.90 ms after SOI Ignition in the experimental images is recognized by the white cloud appearing in image (c) that increases in size and propagates downstream in images (d) and (e). The location of the ignition site in both the ex- perimental and simulation results are downstream from the nozzle. In the experiment, the ignition site is approximately 4.7 mm downstream and in the simulation, the ignition site (based on the 2000 K criterion) is about 1.8 mm downstream. Experimentally, the ignition occurs 1.2 ms after the start of injection while in the simulation ignition occurs 0.20 ms after the start of injection. The experimental images and the simulation results are similar in that the ignition site is near the tip and then propagates downstream. Comparing the trends in the simulated reaction rates to the experimental data identifies similar trends. The peak reaction rate originates near the tip of the jet and is also transported downstream. The structure is similar in the experiments and the simulations but the comparison in magnitudes of ignition delay is only mediocre. These results show early, promising similarities but require further investigation to develop better correlations before more rigorous conclusions can be made. Chapter 4. Implementation 81 4.4.3 Location of Ignition The following analysis is based on the ignition criterion of 2000 K. Table 4.1 shows the ignition delay times (ri „) for each initial temperature with 9 the ignition site and the corresponding mixture fraction, mixture fraction variance and the downstream location (x). All of the mixture fraction values are slightly to the fuel-rich side of stoichiometry for methane/air combustion and the low variance indicates that these are well mixed regions. As the initial temperature increases, the ignition delay time decreases and therefore, the ignition site is progressively closer to the nozzle. Ti(K) x (mm) Z Z7h T ign 1000 3.31 0.070 0.0016 7.5 1100 2.57 0.072 0.0024 4.5 1200 2.09 0.082 0.0036 2.5 1300 1.81 0.083 0.0043 2.0 1400 0.50 0.083 0.0037 1.0 Mean - 0.078 0.0031 - (ms) Table 4.1: Ignition characteristics based on the ignition criterion of 2000 K Plotting contours of mean temperature at the onset of ignition overlaid with the mean mixture fraction and variance from Table 4.1 shows that the isopleths of the mixture fraction and variance intersect at the peak temperature (Figure 4.14). This trend of the mean mixture fraction crossing the mean mixture fraction variance contour repeats for all the simulations as seen in Appendix F. The mixture fraction variance is related to the mixture fraction gradient and the intersection of these isopleths indicates that the gradient in mixture fraction is reaching a maximum while the isopleth of mixture fraction variance is nearing an inflection point. Chapter 4. Implementation 0.5 1 1.5 2 82 2.5 Radial distance (cm) Figure 4.14: Temperature contours overlaid with mixture fraction and variance isopleths for 1200 K In Figure 4.15 the peak reaction rate and the peak temperature at the onset of ignition are compared for an initial temperature of 1200 K . Additional images for the other simulations are included in Appendix G . The peak reaction rate is found on the trailing edge of the vortex but closer to the tip of the jet compared to the peak temperature. The maximum reaction rate in the domain is determined by plotting the reaction rate as a function of the mixture fraction in planes perpendicular to the jet axis. Figure 4.16 shows the reaction rates for various planes in the jet for the initial temperature of 1200 K; additional reaction rate plots can be found in Appendix G . The reaction rate is found to peak near the stoichiometric mixture fraction on the trailing edge of the jet vortex. The peak reaction rate occurs predominantly at the stoichiometric mixture fraction and in a well-mixed region, similar to the peak temperature (Table 4.2). The lower mixture fraction variance coinciding with the peak reaction rate indicates that the location of peak reaction rate has experi- Chapter 4. Implementation 83 14 Figure 4.15: Reaction rate and temperature contours at the onset of ignition: (a) Reaction rate (b) Temperature Mixture fraction Figure 4.16: Reaction rates as functions of mixture fraction Chapter 4. Implementation 84 enced more mixing between the fuel and oxidizer. The location of the peak reaction rate is in good agreement with combustion theory as the reactants will be present in ideal concentrations at stoichiometry, allowing the reaction to proceed at its maximum rate. Ti (K) x (mm) Z ,n z Reaction Rate (kJ/kmol s) 1000 4.72 0.055 0.00053 7.94 xlO 13 1100 3.74 0.065 0.00080 1.57 xlO 14 1200 2.92 0.059 0.0010 2.96 xlO 14 1300 2.74 0.071 0.00046 7.50 xlO 14 1400 2.24 0.055 0.00015 1.34 xlO 15 Mean - 0.061 0.00059 - Table 4.2: Peak reaction rate characteristics 4.5 Conclusions and Recommendations Laminar Flamelet Decomposition for the CSE method was incorporated into KIVA-II and used to simulate the ignition of methane in shock tube experiments for different initial temperatures. 1. Flamelet libraries that describe ignition were developed based on a realistic kinetic mechanism for the combustion of methane and air. Two sets of libraries were developed: the first set is referred to as Full Flamelet Libraries consisting of nearly equal numbers of mixing, igniting, burning and extinguishingflameletswhile the second set is referred to as Igniting Flamelet Libraries where the majority of flamelets are taken between the onset of ignition and the onset of burning. 2. The Full Flamelet Libraries under-predict the ignition delay time for methane and predict a small difference between ignition delay times for different temperatures. The source terms from the extinguishing Chapter 4. Implementation 85 flamelets in the Pull Flamelet Libraries have an undesired influence on the solution; extinguishingflameletshave considerably higher reaction rates compared to ignitingflameletsthat led to the over-prediction of the chemical source terms. Igniting Flamelet Libraries provide better predictions of the ignition delay time as the chemical source terms are not as influenced by the higher reaction rate of the burning and extinguishing flamelets.. 3. The CSE with LFD method is used to predict the ignition of methane in air. The model reasonably predicts the trend of the ignition delay time decreasing with increasing initial temperature. Two different ignition criteria are used to detect ignition in the simulations and the slopes of the temperature-ignition delay time curves are similar compared to experimental results. The slopes of the simulation curves range between 2240 and 2350 K while the experimental results have a slope of 3060 K. Ignition occurs near the trailing edge of the jet vortex and proceeds downstream with the jet. 4. Using the ignition criterion of 2000 K, the isopleths of the mean mixture fraction and variance intersect at the peak temperature. The peak temperature coincides with a well-mixed, fuel-rich region in the domain. 5. The peak reaction rate occurs ahead of the peak temperature in the jet and is characterized by a slightly fuel rich mean mixture fraction (0.061) and a low mean mixture fraction variance (0.00059). The code shows reasonable agreement with several combustion experiment results including mixture fraction and mixture fraction variance. More physical details need to be included for a better correlation with the experimental light emissions and temperature or reaction rate from the simulations. The computation time of the code may be improved by using a library for the PDF of mixture fraction rather than the numerical algorithm. Selecting mostly ignitingflameletsprovides better predictions of the ignition Chapter 4. Implementation 86 delay times. This suggests that a more sophisticated method of selecting theflameletsused in the basis function needs to be developed to allow the flamelet library to describe the evolution from the initial pure mixing case to the igniting case and then to burning and should be investigated. This implementation of LFD for CSE solves for the weighting coefficient vector for each time step and does not account for its progression in both time and space; further investigation of using the weighting coefficient vector from the previous time step and the surrounding neighbours to calculate the new weighting coefficient vector may offer some improvements to the prediction of the ignition delay time. 87 Chapter 5 Conclusions and Recommendations This research has been conducted to determine the applicability of Conditional Source-term Estimation for predicting the ignition of methane under diesel engine-relevant conditions. First, an investigation to determine the effects of the various parameters affecting the solution of Conditional Sourceterm Estimation was conducted. Conditional Source-term Estimation was then implemented into a commercial combustion modeling package and used to predict the ignition of methane. The findings of this study are discussed below. 5.1 A "priori Testing of Conditional Source-term Estimation The investigated parameters include the linear regularization parameters, the number of flamelets used in the library and the optimum scalar to be used to solve for the weighting coefficient vector. 1. The lower order smoothing matrices provide better approximations for the chemical source terms. Higher order smoothing matrices require substantially more truncation, which leads to poor predictions of the chemical source terms. 2. Approximately 50 solutions to the unsteadyflameletequations, or flamelets, are sufficient to describe the physically realizable conditions Chapter 5. Conclusions and Recommendations 88 in the turbulentflameand provide reasonable predictions of the reaction rates. Fewer than 50flameletsinadequately describes the flame while libraries with more than 50flameletsshow that as the magnitude of the weighting coefficient vector decreases, the higher reaction rate flamelets have a greater influence on the predictions. 3. Using temperature only to calculate the weighting coefficient vector provides good approximations for the chemical source terms. The fuel and product scalars also provide good approximations to the reaction rates but are not always present in an engine simulation. The oxidizer provides the worst approximation to the chemical reaction rates while the intermediate species and the nitrous oxide provide only fair approximations. 4. Combining scalars to improve the prediction of the reaction rates offers minimal benefit over using one scalar that provides good approximations. The predictions provided by using combinations of scalars tend to be an average of the predictions from the individual scalars. 5.2 Implementation of Conditional Source-term Estimation Conditional Source-term Estimation was implemented into KIVA-II and used to successfully predict the ignition of methane. The method incorporated the earlier findings and used 55 unsteadyflameletsto form the basis functions then inverted and solved the temperature integral equation for the weighting coefficient vector. 1. Conditional Source-term Estimation replaces the existing kinetic ignition mechanism in KIVA-II and successfully predicts the ignition of methane and reasonably predicts the trend between ignition delay time and initial temperature. 2. Using the Full Flamelet Libraries to predict the ignition of methane Chapter 5. Conclusions and Recommendations 89 found that the ignition delay time of methane was greatly underpredicted compared to experimental results. Upon further investigation, the higher reaction rates of the extinguishing flamelets cause on over-prediction of the chemical source terms. Additionally, the small differences observed between predicted ignition delay times with increasing temperature may not be significant. The slope of the predicted ignition delay curve was 570 K compared to 3060 K for the experimental results. 3. The Igniting Flamelet Libraries provide better predictions for the ignition delay time when compared to experimental results. Consisting of primarily igniting flamelets, the predicted reaction rates are not as influenced by the higher reaction rates of the extinguishing flamelets. The slope of the ignition delay-temperature curve is representative of the activation energy and the predicted slope ranged between 2240 K and 2350 K while the experimental results exhibited a slope of 3060 K. 4. Using an ignition limit of 2000 K, an ignition kernel originates at the trailing edge of the jet vortex and is transported downstream with the jet. This kernel coincides with a well-mixed, fuel-rich region with a mean mixture fraction of 0.078 and a mean mixture fraction variance of 0.0031. As the temperature increases, the kernel appears sooner and closer to the fuel inlet nozzle. 5. The peak reaction rate in the igniting jet occurs physically ahead of the peak temperature and closer to the tip of the jet. It is characterized by a mean mixture fraction of 0.061 and a mean mixture fraction variance of 0.00059. 6. Several similarities between the predicted ignition of methane and experimental results are found. However, no direct correlation has been made between either the simulated results and the light emissions from the experiment. Chapter 5.3 5. Conclusions and Recommendations 90 Recommendations To further test and validate the Conditional Source-term Estimation method, the following work needs investigation: 1. Develop a correlation between the light emissions from the experimental results and the simulation results for temperature or reaction rate. Perhaps an intermediate species would also provide a better correlation. 2. The code should be extended to provide predictions of pollutant emissions and tested against experimental results. 3. A method that would allow the flamelet library to evolve from the pure mixing state to an igniting state without selecting burning or extinguishing flamelets may alleviate the problems encountered with the Full Flamelet Libraries. 4. The weighting coefficient vector is calculated independently in both space and time along each row in the computational domain. A method that would allow the weighting coefficient vector to be calculated using its neighbouring values and values form previous time-steps may improve the prediction of the ignition delay time. 5. It may also be beneficial to replace with basis functions in the Conditional Source-term Estimation method with a different library altogether. Current research is in the process of developing stochastic models incorporating the mixture fraction and strain that may better represent the realizable physical conditions in a flame compared to unsteady solutions to the flamelet equations. 91 Bibliography [1] G. Taylor, Trucks and air emissions, Environment Canada, 2001. [2] E. Sher, Handbook of air pollution from internal combustion engines: pollutant formation and control, Academic Press, 1998, Ch. Environmental aspects of air pollution. [3] R. Bilger, Future progress in turbulent combustion research, Prog. Energy Combust. Science 26 (2000) 367. [4] D. Veynante, L. Vervisch, Turbulent combustion modeling, Prog. Energy Combust. Science 28 (2002) 193. [5] N. Peters, Laminar diffusion flamelet models in non-premixed turbulent combustion, Prog. Energy Combust. Sci. 10 (1984) 319. [6] S. Pope, PDF methods for turbulent reactive flows, Prog, in Energy Combust. Sci. 11 (1985) 119. [7] A. Klimenko, R. Bilger, Conditional moment closure for turbulent combustion, Prog. Energy Combust. Sci. 25 (1999) 595. [8] W. Bushe, H. Steiner, Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows, Physics of Fluids 11 (7) (1999) 1896. [9] J. Heywood, Internal Combustion Engine Fundamentals, McGraw Hill, 1988. [10] Y. Cengel, M. Boles, Thermodynamics: an engineering approach, McGraw Hill, 1989. Bibliography 92 [11] J. Warnatz, U. Maas, R. Dibble, Combustion. Physical and Chemical Fundamentals, Modeling and Simulation, Experiments, Pollutant Formation, Springer, 2001. [12] P. Libby, F. Williams, Topics in applied physics: Turbulent reacting flows, Springer-Verlag, 1980, Ch. Fundamental aspects. [13] F. White, Fluid mechanics, McGraw Hill, 1986. [14] N. Peters, Turbulent combustion, Cambridge University Press, 2000. [15] B. Launder, D. Spalding, Mathematical models of turbulence, Academic Press, 1972. [16] S. Chapra, R. Canale, Numerical methods for engineers, 2nd Edition, McGraw Hill, 1988. [17] A. Cook, J. Riley, A subgrid model for equilibrium chemistry in turbulent flows, Physics of Fluids 6 (1994) 2868. [18] P. Libby, F. Williams, Turbulent reactingflows,Academic Press, 1980. [19] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in FORTRAN 77: The art of scientific computing, Cambridge University Press, 1992. [20] H. Pitsch, Unsteadyflameletmodeling of differential diffusion in turbulent jet diffusionflames,Combust. Flame 123 (2000) 358. [21] H. Pitsch, M. Chen, N. Peters, Unsteadyflameletmodeling of turbulent hydrogen-air diffusionflames,in: Proceedings of the Twenty-Seventh Symposium (International) on Combustion. The Combustion Institute, 1998. [22] P. Coelho, N. Peters, Unsteady modelling of a piloted methane/air jet flame based on the Eulerian particleflameletmodel, Combust. Flame 124 (2001) 444. Bibliography 93 [23] H. Pitsch, Extended flamelet model for les of non-premixed combustion, in: Annual Research Briefs. Center for Turbulence Research, NASA Ames/Stanford University, 1998. [24] D. Lentini, Assessment of the stretched laminar flamelet approach for nonpremixed turbulent combustion, Combustion Science and Technology 100 (1994) 95. [25] Y. Buriko, V. Kuznetsov, D. Volkov, S. Zaitsev, A test of a flamelet model for turbulent nonpremixed combustion, Combust. Flame 96 (1994) 104. [26] B. Cuenot, F. Egolfopoulos, T. Poinsot, An unsteady laminar flamelet model for non-premixed combustion, Combustion Theory and Modelling 4 (2000) 77. [27] A. Bourlioux, O. Volkov, Validation of unsteady flamelet models for non-premixed turbulent combustion with intermittency, in: Proceedings of the Eleventh Annual Conference of the CFD Society of Canada, 2003. [28] R. Barlow, Piloted CH4 /airflameD comparisons, in: Proceedings of the Third International Workshop on Measurement and Computation of Turbulent Nonpremixed Flames, 1998. [29] R. Bilger, Conditional moment closure for turbulent reactingflow,Phys. Fluids 5 (1993) 436. [30] A. Klimenko, Multicomponent diffusion of various admixtures in turbulentflow,Fluid Dynamics 25 (1990) 327. [31] W. Mell, V. Nilsen, G. Kosaly, J. Riley, Investigation of closure models for nonpremixed turbulent reacting flows, Physics of Fluids 6 (1994) 1331. Bibliography 94 [32] C. Devaud, K. Bray, Assessment of the applicability of conditional moment closure to a lifted trubulentflame:firstorder model, Combust. Flame 132 (2003) 102. [33] R. Barlow, N. Smith, J.-Y. Chen, R. Bilger, Nitric oxide formation in dilute hydrogen jetflames:Isolation of the effects of radiation and turbulence-chemistry submodels, Combust. Flame 117 (1999) 4. [34] M. Fairweather, R. Woolley, First order conditional moment closure modeling of turbulent, nonpremixed hydrogenflames,Combust. Flame 133 (2003) 393. [35] W. Bushe, H. Steiner, Large eddy simulation of a turbulent reacting jet with conditional source-term estimation, Physics of Fluids 13 (3) (2001) 754. [36] W. Bushe, H. Steiner, Laminarflameletdecomposition for conditional source-term estimation, Physics of Fluids 15 (6) (2003) 1564. [37] W. Bushe, R. Bilger, G. Ruetsch, Incorporating realistic chemistry into direct numerical simulations of turbulent non-premixed combustion, in: Annual Research Briefs. Center for Turbulence Research, NASA Ames/Stanford University, 1997. [38] P. K. Kythe, P. Puri, Computational Methods for Linear Integral Equations, Birkhauser, 2002. [39] C. Lawson, R. Hanson, Solving Least Squares Problems, Prentice-Hall Inc., 1974. [40] MATLAB v6.5, The Mathworks Inc., 2002. [41] A. Amsden, P. O'Rourke, T. Butler, KIVA - IT. A computer program for chemically reactiveflowswith sprays, Los Alamos National Laboratory . Bibliography 95 [42] W. Pracht, Calculating three-dimensional fluid flows at all speeds with an eulerian-lagrangian computing mesh, Journal of Computing Physics 17 (1975) 132. [43] G. Abramovich, The theory of turbulent jets, M.I.T. Press, 1963. [44] P. Ouellette, P. Hill, Turbulent transient gas injections, ASME Journal of Fluids Engineering 122 (2000) 743. [45] J. Iaconis, An investigation of methane autoignition behaviour under diesel engine-relevant conditions, MASc Thesis (UBC), 2003. [46] J. Huang, P. Hill, W. Bushe, S. Munshi, Shock-tube study of methane ignition under engine-relevant conditions: experiments and modeling, Combust. Flame . [47] J. Huang, Experimental shock tube study of ignition promotion for methane under engine relevant conditions, MASc Thesis (UBC), 2001. [48] G. Sullivan, J. Huang, Personal communication, University of British Columbia, 2003. [49] H. Bi, A. Agrawal, Study of autoignition of natural gas in diesel environments using computationalfluiddynamics with detailed chemical kinetics, Combustion and Flame 113 (1998) 289. [50] A. Birch, D. Brown, M. Dodson, J. Thomas, The turbulent concentration field of a methane jet, Journal of Fluid Mechanics 8 (3) (1978) 431. [51] R. Schefer, R. Dibble, Mixture fraction field in a turbulent nonreacting propane jet, AIAA Journal 39 (2001) 64. [52] R. Alvani, M. Fairweather, Ignition characteristics of turbulent jet flows, in: Proceedings of the 7th annual UK conference onfluidmixing, 2002. [53] FLUENT v6., FLUENT Inc., 2002. 96 Appendix A Additional Converged Statistics N u m b e r of flamelets 60 cells 80 cells 120 cells 240 cells (% error) (% error) (% error) (% error) 5 891.0 892.1 893.2 894.7 10 883.3 886.6 889.0 893.3 25 592.0 594.7 596.6 600.4 50 261.2 261.1 260.9 260.3 75 282.4 283.2 283.4 283.6 100 357.3 358.4 358.7 359.2 Table A . l : Comparison of global reaction rates for converged statistics of reaction rate 2 N u m b e r of flamelets 60 cells 80 cells 120 cells 240 cells (% error) (% error) (% error) (% error) 5 30.54 30.56 30.60 30.66 10 46.54 46.63 46.71 46.91 25 24.22 24.00 23.78 48.13 50 37.12 37.40 36.96 36.77 75 39.34 39.33 39.30 39.26 100 37.41 37.43 37.46 37.53 Table A.2: Comparison of global reaction rates for converged statistics of reaction rate 3 Appendix B Additional Scalar Combinations Comparisions Appendix B. Additional Scalar Combinations Comparisons Combination 98 Reaction 1 Reaction 2 Reaction 3 (% e r r o r ) (% e r r o r ) (% e r r o r ) Temperature and NO 14.66 66.42 37.37 Temperature and Product 4.11 27.52 14.38 Temperature and Intermediate 97.84 85.20 92.05 Fuel and Oxidizer 4.47 17.75 8.95 Fuel and Product 3.78 13.86 11.90 Fuel and NO 2.49 37.05 45.3 Oxidizer and Product 97.88 123.9 61.67 Oxidizer and Intermediate 94.51 84.56 99.92 Oxidizer and NO 54.21 145.1 92.18 Product and NO 18.19 76.23 31.81 Temperature, Fuel and Oxidizer 7.38 22.26 0.23 Temperature, Fuel and Intermediate 23.42 18.56 13.47 Temperature, Fuel and NO 13.96 57.28 53.19 Temperature, Product and NO 49.92 44.41 18.67 Temperature, Intermediate and NO 81.95 112.5 79.18 Fuel, Product and Intermediate 77.02 67.96 98.20 Fuel, Product and NO 45.31 41.35 27.73 Fuel, Intermediate and NO 51.96 78.11 80.95 Product, Intermediate and NO 15.68 21.60 36.14 Table B.l: Additional Comparison of global reaction rate for combinations of scalars. 99 Appendix C Validation of Mixture Fraction and Variance The PDF of mixture fraction is used to describe the internal structure of the flame and provide an understanding of the mixing process in turbulent jets. The mixture fraction and mixture fraction variance are need to calculate the PDF but KIVA-II does not solve transport equations for mixture fraction or mixture fraction variance, therefore solution of these equations needs to be included. To validate the solution of the transport equations for mixture fraction and its variance, a non-reacting jet was simulated and the predicted fields were compared to experimental results. Many researchers have made detailed measurements of turbulent concentration parameters in constant-density turbulent jets to investigate the turbulent mixing and flow structure (Birch et al. [50], Schefer and Dibble [51], Alvani and Fairweather [52]). The results of these studies are used here to validate the predictions of mixture fraction and mixture fraction variance from the modified KIVA-II code. FLUENT (FLUENT [53]) is another CFD package that uses RANS models for turbulence and is also used to predict the centreline mixture fraction and variance for comparison with the KIVAII predictions. The mixture fraction and mixture fraction variance curves are similar for jets of different gasses provided the density ratio (^) is between 0.14 to 5.11. Further research has shown that the asymptotic centerline value is independent of jet density ratio (Schefer and Dibble [51]). With the modified Appendix C. Validation of Mixture Fraction and Variance 100 KIVA-II code, non-reacting methane jets are simulated where the density ratio is approximately 0.44 and thus within the acceptable range. Mixture Fraction The mixture fraction is calculated as a linear combination of mass fractions using the additional inert element argon. This is accomplished by adding argon to the initial conditions describing the oxidizer in KIVA-II then solving Eq. 2.30 for each control volume cell. Predicted contours of the mixture fraction for a methane jet are shown in Figure C . l . The contours follow those of the jet with a strong core of nearly pure fuel diffusing into the surrounding air. 0.5 1 1.5 2 2.5 Radial distance (cm) Figure C . l : Predicted contours of mixture fraction for an unsteady methane jet Figure C.2 shows comparisons between the predicted mixture fraction and experimental results for the jet centreline and several radial locations downstream of the nozzle. The centreline mixture fraction calculated with the modified KIVA-II code closely matches both the experimental results Appendix C. Validation of Mixture Fraction and Variance 101 and the F L U E N T prediction. The downstream radial profiles also show very good agreement between the predictions and the experimental results (Schefer and Dibble [51]). Overall, the predictions of mixture fraction are in good agreement with experimental results. Figure C.2: Comparison between predicted and experimental mixture fraction (a) Centreline comparison (b) 15 nozzle diameters downstream (c) 30 nozzle diameters downstream (d) 50 nozzle diameters downstream Mixture Fraction Variance Figure C.3 shows contours of mixture fraction variance for an unsteady methane jet and exhibits the expected trends. The mixture fraction variance is low in both pure fuel and pure oxidizer and increases to a peak where the fuel and oxidizer are mixing. Comparisons between the predicted values and experimental results along the centreline and three downstream radial locations are shown in Figure Appendix C. Validation of Mixture Fraction and Variance 102 10.045 0.04 0.035 0 03 ! 0.025 0.02 0.015 0.01 0.005 0.5 1 1.5 2 2.5 Radial distance (cm) Figure C.3: Predicted contours of mixture fraction variance for an unsteady methane jet C.4. The KIVA-II model and F L U E N T over-predict the peak mixture fraction variance along the centreline. The simulated radial profile 15 nozzle diameters downstream agrees best with the experimental results even though it over-predicts the radial penetration. Further downstream (30 and 50 nozzle diameters) the simulation over-predicts the peak mixture fraction variance while under-predicting the radial penetration. The poor predictions of the mixture fraction variance are probably due to inaccuracies in the fc — e turbulence model. As these results show fair agreement with experiment and better predictions compared to the F L U E N T centreline predictions, the model is expected to provide adequate predictions for the mixture fraction variance. The mixture fraction predictions from the KIVA model are in good agreement with experimental results while the predictions of the mixture fraction variance are mediocre. However, the predictions of the mixture fraction and its variance are deemed adequate for the purpose of implementing the C S E Appendix C. Validation of Mixture Fraction and Variance 0.04 103 0.02 r 1 2 3 4 Radial distance (X/D) (d) (C) 1 2 3 4 Radial distance (X/D) Figure C.4: Comparison between predicted and experimental mixture fraction variance (a) Centreline comparison (b) 15 nozzle diameters downstream (c) 30 nozzle diameters downstream (d) 50 nozzle diameters downstream Appendix C. Validation of Mixture Fraction and Variance method to predict the ignition of methane. 104 105 Appendix D Grid Convergence The grid independence of the KIVA-II code with CSE solution was tested by comparing solutions with successively refined grids. Several meshes were built for the domain where the number of cells across the inlet boundary were varied. Thefirstmesh used a grid spacing equivalent to \ the nozzle diameter followed by grid spacings of \ and ^ the nozzle diameter. A correspondingly refined stretched mesh was used for the remainder of the domain. For comparison, the fuel mass fraction at a time of 0.40 ms after the start of injection is used. Figure D.2 shows the centreline values of the fuel mass fractions for the three different grid spacings. Although an improvement in the accuracy of the predictions can be achieved by reducing the grid spacing, the \ of the nozzle diameter grid spacing gives satisfactory accuracy and further improvement is relatively small. When the grid spacing is reduced to ^ the nozzle diameter, the computation time increased considerably. For example, the simulations to reach 0.40 milliseconds after the start of injection with a grid spacing of j took approximately 2.5 hours whereas the same simulation for a grid spacing of ^ took nearly 8 hours . Consequently, it was decided l to use a grid spacing of \ for the simulations as a trade-off in accuracy and for the sake of computation time. 1 computations were performed using a Pentium III 550 MHz processor 106 — gricMra D — grid=1/4 D — grid*1/10 D Axial distance (Z/D) Figure D . l : Centreline fuel mass fractions for different grid spacings Appendix E. Time Evolution of Temperature Contours 107 Appendix E Time Evolution of Temperature Contours 0.5 1 1.5 2 2.5 0.5 (a) Radial distance (cm) (b) (c) Radial distance (cm) (d) 1 1 5 2 2.5 Radial distance (cm) 0.5 1 1.5 2 2.5 Radial distance (cm) Figure E . l : Temperature contours at different times during methane ignition, 1000 K (a) 0.15 ms after SOI (b) 0.35 ms after SOI (c) 0.55 ms after SOI (d) 0.75 ms after SOI Appendix E. Time Evolution of Temperature Contours 108 Figure E.2: Temperature contours at different times during methane ignition, 1100 K: (a) 0.10 ms after SOI (b) 0.30 ms after SOI (c) 0.50 ms after SOI (d) 0.70 ms after SOI Figure E.3: Temperature contours at different times during methane ignition, 1300 K: (a) 0.10 ms after SOI (b) 0.15 ms after SOI (c) 0.20 ms after SOI (d) 0.25 ms after SOI Appendix E. Time Evolution of Temperature Contours 05 (c) 1 1.5 2 Radial distance (cm) 2.6 0.5 ( d ) 1 1.5 2 109 2.5 Radial distance (cm) Figure E.4: Temperature contours at different times during methane ignition, 1400 K: (a) 0.05 ms after SOI (b) 0.10 ms after SOI (c) 0.15 ms after SOI (d) 0.20 ms after SOI 110 Appendix F Temperature, Mixture Fraction and Mixture Fraction Variance Contours at Ignition 0-5 1 1.5 2 2.5 Radial distance (cm) Figure F . l : Temperature contours overlaid with mixture fraction and variance contours for 1000 K. Appendix F. Temperature, Mixture Fraction and Variance 0.5 1 1.5 2 111 25 Radial distance (cm) Figure F.2: Temperature contours overlaid with mixture fraction and variance contours for 1100 K . 0.5 1 1.5 2 2.5 Radial distance (cm) Figure F.3: Temperature contours overlaid with mixture fraction and variance contours for 1300 K . Appendix F. Temperature, Mixture Fraction and Variance Radial distance (cm) Figure F.4: Temperature contours overlaid with mixture fraction and ance contours for 1400 K . 113 Appendix G Reaction Rate and Temperature Contours 0.5 (a) 1 1.5 2 Radial distance (cm) 2.5 0.5 (b) 1 1.5 2 2.5 Radial distance (cm) Figure G.l: Reaction rate and temperature contours at the onset of ignition for 1000 K: (a) OReaction rate, (b) Temperature. Appendix G. Reaction Rate and Temperature Contours 114 Figure G.3: Reaction rate and temperature contours at the onset of ignition for 1100 K: (a) Reaction rate, (b) Temperature. Appendix G. Reaction Rate and Temperature Contours 115 Figure G.5: Reaction rate and temperature contours at the onset of ignition for 1300 K: (a) Reaction rate, (b) Temperature. Appendix G. Reaction Rate and Temperature Contours 116 Mixture fraction Figure G.6: Reaction rates as functions of mixture fraction for 1300 K Figure G.7: Reaction rate and temperature contours at the onset of ignition for 1400 K: (a) Reaction rate, (b) Temperature. Appendix G. Reaction Rate and Temperature Contours 117 Figure G.8: Reaction rates as functions of mixture fraction for 1400 K
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Implementation of conditional source-term estimation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Implementation of conditional source-term estimation for the prediction of methane ignition Blair, C. W. 2003
pdf
Page Metadata
Item Metadata
Title | Implementation of conditional source-term estimation for the prediction of methane ignition |
Creator |
Blair, C. W. |
Date Issued | 2003 |
Description | 'Implementation of Conditional Source-term Estimation for the Prediction of Methane Ignition' The Laminar Flamelet Decomposition approach to Conditional Source-term Estimation is a recent development in combustion modeling that approximates the conditional averages of scalars as linear combinations of basis functions obtained from a library of unsteady solutions to the classical flamelet equations. Then an integral equation for a single scalar is inverted and solved for the weighting function of the basis functions. This weighting function is then used to obtain closure for the chemical source terms in the transport equations. Many parameters affect the solution of the Conditional Source-term Estimation method including the solution technique for inverting the integral equation; the number of solutions in the flamelet library; and the scalar used for the integral equation. The first part of this research investigated these parameters to optimize the solution for the Conditional Source-term Estimation. This is accomplished in the a priori sense by comparison to Direct Numerical Simulations of a simple turbulent flame. The second part of this work implemented Conditional Source-term Estimation into the commercially available multi-dimensional combustion software, KIVA-II by replacing the existing ignition mechanism. The model was used to predict the ignition of methane for different initial temperatures and compared to experimental results. The predicted ignition delay times were shorter than the experimental results but exhibited similar trends. |
Extent | 13635562 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0081001 |
URI | http://hdl.handle.net/2429/15085 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-0079.pdf [ 13MB ]
- Metadata
- JSON: 831-1.0081001.json
- JSON-LD: 831-1.0081001-ld.json
- RDF/XML (Pretty): 831-1.0081001-rdf.xml
- RDF/JSON: 831-1.0081001-rdf.json
- Turtle: 831-1.0081001-turtle.txt
- N-Triples: 831-1.0081001-rdf-ntriples.txt
- Original Record: 831-1.0081001-source.json
- Full Text
- 831-1.0081001-fulltext.txt
- Citation
- 831-1.0081001.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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0081001/manifest