Combustion Modeling with Conditional Source-Term Estimation and Laminar Flamelet Decomposition by Ray W.S. Grout B . A . S c , The University of British Columbia, 2002 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F THE REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF APPLIED SCIENCE in The Faculty of Graduate Studies (Department of Mechanical Engineering) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A July 8, 2004 © Ray W.S. Grout, 2004 THE UNIVERSITY OF BRITISH C O L U M B I A FACULTY OF G R A D U A T E STUDIES Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Name of Author (please < ^ 3 / o Sr/2JZX=><J' / _( Date. (dd/mm/yyyy print) Title of Thesis: Degree: Year: Department of The University of British Columbia Vancouver, B C Canada grad.ubc.ca/forms/?formlD=THS page 1 of 1 last updated: 20-Jul-04 11 ABSTRACT Conditional source-term estimation w i t h laminar flamelet decomposition has been utilized to model the mean chemical source term i n a predictive R A N S simulation of two different problems. W i t h this model, integral equations relating the unconditional mean temperature and species fields to the conditional means are used to determine the flame structure as a combination of basis functions formed from solutions to the unsteady laminar flamelet equations. In the simulation of the Sandia F l a m e ' D ' , a well studied co-flowing piloted jet flame w i t h a steady average solution, a converged solution was obtained which captured the trends i n the temperature and major species, although the nitric oxide prediction overestimated the peak concentration by a substantial margin. Simulation of the autoignition process of non-premixed methane produced simulation predictions in excellent agreement w i t h the experimental data. Using a new, more stringent, criteria to define ignition t h a n earlier studies, the effect of ambient temperature on the ignition delay was captured, as was the expected physical behaviour prior to ignition. T h e errors i n the simulation of the co-flowing piloted jet can be attributed to a large degree to the lack of ability of the basis functions used to account for the effect of the pilot flame. T h a t the autoignition simulation was much more successful highlights the importance of including all relevant physics i n the library of basis functions. Earlier formulations of the model were extended by expanding the sample space for the basis functions to include a wide variety of solutions to the laminar flamelet equations; i n order to distinguish between the larger set, a new method of stabilizing the numerics was found to be necessary and the number of scalars used to determine the flame structure was increased. T h e formulation is easily extensible to libraries of basis functions which are capable of including effects such as pilot and edge flames not captured by the laminar flamelet solutions. iii CONTENTS Abstract ii Contents iii List of Tables v List of Figures vi Nomenclature viii Acknowledgements x Introduction xi 1 Background 1.1 Turbulent Combustion 1.2 Flamelet Models 1.3 Conditional Moment Closure 1.4 Conditional Source-Term Estimation References 2 Simulation of the Sandia Flame 'D' 2.1 2.2 2.3 2.4 Introduction P r o b l e m Description Background Formulation 2.4.1 Chemical Source Term 2.4.2 F l o w Solution 2.5 Results and Discussion 2.6 Conclusion References 1 2 3 12 17 25 30 31 32 32 35 35 36 37 41 44 CONTENTS iv 46 3 Non-Premixed Methane Auto-Ignition 3.1 3.2 Introduction 47 Formulation 48 3.2.1 C S E Combustion M o d e l 48 3.2.2 F l o w F i e l d 52 3.2.3 Basis Functions 53 3.3 Results and Discussion 58 3.3.1 Predictions of Ignition Delay 58 3.3.2 Relation to P r i o r Studies 60 3.3.3 Behaviour P r i o r to Ignition and Identification of Ignition 64 3.3.4 Effect of Turbulence on Ignition Delay 68 3.4 Conclusion 68 References 71 4 Discussion and Conclusion 73 References 82 83 A RANS Fields \ V LIST OF TABLES 3.1 3.2 A l l o c a t i o n of Variance for Ignition Delay Location and M i x t u r e Fraction at tj, 59 60 vi LIST OF FIGURES 1.1 1.2 1.3 2.1 2.2 2.3 2.4 2.5 3.1 3.2 3.3 3.4 3.5 3.6 3.7 Borghi Diagram Scatter P l o t of Temperature vs. M i x t u r e Fraction for Sandia F l a m e ' D \ x/D = 30 [1] Relation Between Temperature and M i x t u r e Fraction, L o g of Scalar Dissipation. Portion of Figures 1 and 2 i n [9] Selected Basis Function Solutions Predicted Temperature F i e l d C S E / L F D Validation - Comparison of Conditional Averages to Experimental Results (Solid lines predictions, symbols experimental measurements) C S E / L F D Validation - Comparison of A x i a l and R a d i a l P r o files to Experimental Results (Solid lines predictions, symbols experimental measurements) N i t r i c Oxide Prediction 13 17 34 38 39 42 43 Flamelet Solutions for Temperature and Enthalpy Source T e r m Flamelet Solutions for Species Mass Fraction - Yco and YCH OH Flamelet Solutions for Species Mass Fraction - Contour Plots Flamelet Temperature [K] Contours for Complete L i b r a r y . . . Effect of Ambient Temperature on Ignition Delay Effect of Ambient Temperature on Ignition Delay - C o m p a r i son P r i o r Works Evolution of Conditional Averages in M i x t u r e Fraction Space, 54 55 56 57 58 T 65 3 ox = U50K 3.8 Evolution of Peak Increase in Conditional Averages 3.9 Comparison of Results for Different Ignition C r i t e r i a 3.10 R a t i o Between Turbulent and L a m i n a r Ignition Delays 4.1 8 Flamelet Solutions for T and Yco . . . . 61 67 67 69 76 LIST OF A.l FIGURES R A N S Temperature and C H Creation Contours at td 84 Vlll NOMENCLATURE Symbols Notation Velocity Spatial Coordinate Temperature Species Mass Fraction M i x t u r e Fraction Sample space for mixture fraction Z" M i x t u r e Fraction Fluctuations (Variance) Le Lewis Number k T h e r m a l Conductivity D Diffusivity Cp Heat capacity Ka K a r l o v i t z Number Da Damkohler Number tk Kolmogorov time scale t Chemical time scale tj Integral time scale t f Non-premixed reference time scale lf Non-premixed reference length scale u f Non-premixed reference velocity scale v K i n e m a t i c Viscosity e Dissipation rate of turbulent kinetic energy p Density cjj Source term for species i X Scalar Dissipation R a t e a A r b i t r a r y Scalar a Coefficient for flamelet n ~a Vector of coefficients ~u ~x T Y Z £ 2 c re re Te n (ai,...,a ) N ipj k Integral of p(C') across discretization interval t ~v Vector M Matrix (T\Q Conditional average of quantity to left of |, conditioned o n conditioning variable taking the value to the right of | T Favre average T Average ||T|| Magnitude M* Conjugate X*(zT) Functional dependence of scalar dissipation on mixture fraction n* Continuous extension of n A Relative penalty i n minimization NOMENCLATURE Subscripts i n j k Species index Basis function index Computational cell (physical space) F i n i t e difference point (mixture fraction space) ix Abbreviations Pdf CMC CSE RANS LES DNS Probability density function Conditional Moment Closure Conditional Source-term E s t i m a t i o n Reynolds Averaged Navier-Stokes Large E d d y Simulation Direct Numerical Simulation X ACKNOWLEDGEMENTS I would like to thank the many people who supported and encouraged me throughout this project. F i r s t , I am grateful to my supervisor, D r . W . K e n d a l Bushe, who not only guided my research efforts, but shortly after we met convinced me that the best way to describe turbulent combustion simulation is " S i x pounds of cool i n a five pound sack". A l s o , thanks are due to my friends and family (especially my parents) who kept my morale high, and life i n the office interesting. A m o n g the many people who have helped me out along the way are C o l i n , C h a d , and Jorge, who were k i n d enough to proofread this thesis. I a m also indebted to C o l i n for fighting the good fight w i t h K I V A and sparing me the trials of adding transport equations to that particular code. A s well, I would perhaps not have undertaken this project had it not been for the encouragement of Dave M u m f o r d and D r . H i l l at Westport Innovations while I was an engineering undergrad. T h e financial assistance of Westport Innovations, the N a t u r a l Science and Engineering Research C o u n c i l of Canada, and the B C Advanced Systems Institute was gratefully appreciated. Finally, I would like to express my gratitude to A m a n d a , my loving wife, who has given me her unwavering support, caring, and encouragement throughout this project as always. It is to her that I dedicate this thesis. Ray Grout J u l y 2004 INTRODUCTION INTRODUCTION xii There is significant motivation to improve the performance of combustionbased energy conversion devices. A s we enter the twenty-first century, our society's demand for energy is extraordinary. One need only recall the effects of the recent power failures i n eastern N o r t h A m e r i c a , or remember the plight of California during the hot summer months to realize that satisfying the demand for energy is essential to maintain our way of life. Awareness of the impact of our energy consumption on the environment has led to efforts to reduce the harm caused by energy conversion devices, either through concern for harmful emissions resulting i n tightening E P A and European acceptable limits for automotive emissions, or broad base policy initiatives such as the K y o t o A c c o r d on greenhouse gas reduction. A t the same time, the demand for affordable energy i n developing economies has stretched the world's energy resources. O n all fronts, it is obvious that we must improve our primary means of energy conversion — controlled combustion — i n terms of b o t h efficiency and minimizing the release of harmful emissions. Realizing this objective necessitates accurate tools for modeling the combustion process. T h e design of devices which control the combustion process for greatest efficiency and least environmental impact can only be done w i t h a thorough understanding of the physics at work and an effective means of evaluating design alternatives at one's disposal. Experimentation, valuable for confirming that the final design performs as expected, is limited as a useful design tool. W h i l e it is possible to build a prototype of, for example, a proposed combustion environment's geometry, obtaining the precise temperature distribution throughout that environment during the combustion event is unfeasible experimentally. Numerical simulations, based on models developed and validated rigorously using fundamental experiments, hold the possibility to economically meet the need for detailed information. Development of adequate models involves incorporating the effects of a multitude of complex and coupled physical phenomena. In non-premixed gas phase turbulent combustion, the subject of this work, there are three m a i n areas of interest which are tightly coupled. T h e macro-mixing behaviour of the flow can be treated using conventional computational fluid d y n a m ics ( C F D ) techniques. Chemical kinetic mechanisms form an entire field of ongoing research; high quality detailed mechanisms are available only for some fuel/oxidizer combinations, tuned by comparison to a wide variety of experimental data. M i c r o - m i x i n g phenomena, and the resulting turbulencechemistry interactions, are not directly captured by a C F D solution for ensemble or spatial averages. Since direct numerical simulation ( D N S ) (which INTRODUCTION xiii does capture these effects) is beyond current computing limits for either realistic geometries or chemical kinetic mechanisms, models to account for these effects must be provided. T h e subject of this thesis is one such model, conditional source-term estimation ( C S E ) . After a general formulation of the model i n the latter part of the next chapter, two papers are presented which describe the application of the model w i t h little adjustment to two diverse problems using commercial C F D tools commonly employed by industry. T h e first, an extended version of that presented at the 2004 Canadian Section of the Combustion Institute Spring Technical M e e t i n g , addresses the simulation of a n atmospheric pressure co-flowing piloted methane jet. T h e second considers the auto-ignition of non-premixed methane under high pressure and moderate temperature such as found i n a diesel engine, and is scheduled for submission i n the near future. These two applications illustrate the flexibility of the model, a n d the wide range of conditions under which it has the potential to function effectively. A concluding discussion regarding the strengths and weaknesses of the model, as well as how it relates to the other common approaches t o combustion modeling, will follow. 1 Submitted to CI/CS 2004 Spring Technical Meeting as: Grout, R. and Bushe, W.K., Analysis of the Sandia Flame 'D' Using an Implementation of Conditional Source-term Estimation in a Commercial RANS Solver, CI/CS Spring Technical Meeting (2004). BACKGROUND CHAPTER 1.1 1. BACKGROUND 2 Turbulent Combustion In order to perform a useful computational analysis of a reacting flow, it is necessary to account for energy release and species production/consumption due to chemical reaction. Even i n the absence of transport phenomena, this is a daunting task due to the complexity of the reaction mechanisms necessary to describe even the simplest of combustion processes. T y p i c a l reaction mechanisms for the combustion of methane i n air may involve 40-50 species, and 200-300 elementary reactions [14, 15]. Solution of such a react i o n mechanism for a fixed composition requires solving a system of ordinary differential equations ( O D E s ) involving one equation for each species, plus temperature, which may evolve on timescales varying by several orders of magnitude. W h e n transport is included, a numerical solution to the NavierStokes and species evolution equations is physically accurate a n d theoretically obtainable using state of the art techniques for a fully resolved flow. However, if the flow is turbulent fully resolved simulations are precluded by computational limitations. T h e typical approach is to not fully resolve the flow, and provide a model for the unresolved effects. However, the methods typically used to account for the unresolved effects of turbulence i n a nonreacting flow are inadequate for turbulent combustion. Numerical simulation of practical turbulent flows solve the governing equations for average values, i n either an ensemble averaged ( R A N S ) or spatially filtered ( L E S ) sense, a n d provide a model for the mean effect of the unresolved (fluctuating or sub-grid scale) components based on the available mean values. W h i l e this approach is reasonably effective at approximating these effects i n the eddy viscosity, attempting to estimate chemical source terms directly from mean quantities is hopeless — i n a highly non-linear Arrhenius reaction rate expression, the fluctuations about the mean for temperature and species concentration w i l l dominate the result. Assuming that the unresolved variation for all dependent variables of i n terest is correlated, it is possible to relate the unresolved variation i n all scalars to a representative scalar. It is then necessary to know the probability density function of only that representative scalar. In the literature, this is done i n different ways to obtain two of the most popular modeling approaches for non-premixed turbulent combustion — L a m i n a r Flamelet m o d eling ( L F M ) and Conditional Moment Closure ( C M C ) . T h e representative scalar most appropriate for non-premixed combustion is termed the mixture fraction, Z. It is a conserved scalar i n the flow, often assumed to be trans- CHAPTER 1. BACKGROUND 3 ported w i t h unity Lewis number, and represents the degree of mixedness. Neither created nor destroyed, conserved scalars are free from the encumbrance of a source term. T h i s feature allows accurate transport equations to be written for the mean and higher moments of conserved scalars. Expressing the mean values of unconserved dependent scalars (temperature, species mass fractions) i n terms of the conserved scalar(s) allows the known higher moments of the conserved scalar to be used to improve the prediction of the mean effects of chemical reaction. T h e probability density function for the mixture fraction is commonly assumed to be a B e t a pdf. W i t h such a presumption, only the first two moments need to be estimated from the flow. L a m i n a r Flamelet models and Conditional Moment Closure, while sharing the use of a conserved scalar, differ dramatically i n their conceptual origins. A s always, there is a trade-off between the performance of the model i n terms of accuracy and economy. L a m i n a r Flamelet models are conceptually straightforward, simple to implement, and easily incorporate the effects of complex chemistry w i t h minor computational cost. However, the range of applicability is limited, and many of the subtleties of the combustion process are not accounted for. Conditional Moment Closure methods are much more robust i n their range of application, have a rigorously justifiable derivation, and are much more computationally expensive, particularly when realistic chemistry is employed. Also, the necessary modeling of many unclosed terms i n the resulting equations introduces modeling assumptions of questionable accuracy. A s a n alternative to these two presumed-pdf / conserved scalar models, a class of models known as ' p d f models exist, for which the evolution of the joint probability density function (pdf) of a l l scalars and velocities is solved [36]. In such a formulation, the chemical terms are closed exactly — it is the turbulent transport and molecular mixing which must be modeled. A s there is little from this model which is essential to the body of this thesis, it w i l l not be further discussed. T h e balance of this chapter w i l l first review flamelet modeling and C M C , then finish w i t h the development of the conditional source-term estimation model which is the central topic of this thesis. 1.2 Flamelet Models Pioneered for use i n non-premixed flows by Peters [29], flamelet models use a coordinate transform to express the reactive scalars (species concentration CHAPTER 1. BACKGROUND 4 and temperature) in terms of a conserved scalar. B y assuming a laminar structure, it is possible to close the transformed equations i n a form decoupled from the original (~x,t) coordinates. T h e chemistry problem can then be solved i n conserved scalar space only. T h i s results i n a highly cost-effective method, but unfortunately it is only strictly valid where it is justifiable to make the assumption of the laminar structure. S t i l l , flamelet modeling is conceptually attractive and serves to provide an introduction to the use of a conserved scalar and presumed pdf in combustion modeling. A s well, the flavour of conditional source-term estimation explored i n this work makes use of solutions to the flamelet equations. In a general turbulent reacting flow, the key scalars (velocity, temperature, and species mass fraction) are naturally functions of a spatial and temporal coordinate system; Tt = ll{-x,t), T = T(x,t), Yi = YiCx,t) (1.1) T h e mixture fraction (Z) is also a dependent variable naturally a function of (~x,t). A l t h o u g h an expression for the mixture fraction can be fabricated from the individual species concentrations by combining them so as to produce a conserved scalar and normalizing such that the resulting scalar takes a value of 1 i n the fuel stream and 0 in pure oxidizer, convention in the simulation literature now follows a more convenient definition. T h i s definition is simply that the mixture fraction is transported as a conserved scalar w i t h unity Lewis number [4, 44]; that is, it satisfies E q . 1.5 w i t h Co = 0 and D = -J^r. Boundary conditions are Z = 0 in the oxidizer stream and Z = 1 i n the fuel stream. A l t h o u g h some would argue that this definition lacks a clear physical interpretation, there is little point i n debating the issue. F i r s t , for the purposes of the flamelet models (and, as we shall see later, C M C ) , the physical interpretation is irrelevant beyond the existence of a strong correlation between the mixture fraction and the other dependent variables [34]. Second, although it will involve a computationally difficult deconvolution for a reacting flow, by virtue of the existence of the correlation between dependent variables it is conceptually possible to see that any alternate definition of the mixture fraction (for instance, one based on elemental conservation) must be expressible i n terms of this now widely accepted definition. T h e 'flamelet equations' are developed i n detail by Peters [29] for a globally unity Lewis number. A n extended formulation incorporating the effects of differential diffusion was evaluated by P i t s c h [32] for a test case involving a CHAPTER 1. BACKGROUND 5 turbulent jet flame where molecular diffusivity was of the same order of magnitude as turbulent diffusivity i n the near field. A s only modest improvement resulted from the substantial effort to address differential diffusion, the brief synopsis of the method presented here will consider only equal diffusivities. N o t i n g the correlation between the mixture fraction and the other scalars of interest, Peters observes that it should be possible to perform a Croccotransformation for the governing equations, which expresses dependent variables of interest i n terms of other dependent variables. If this is done, and the dependence on the original independent variables is eliminated entirely, it is then strictly valid to solve for the field Z(~x , t) and then map the other dependent variables into (~x, t) space from their solution in (Z, t) space, although some iteration may be necessary to capture the influence of the mapping on the Z(~x,t) solution, e.g., that due to density changes related to increased temperature. T o attempt to do this, Peters considers the coordinate system defined by Z (implicitly normal to the stoichiometric interface), Z2 = X2, Z3 = X3, T = t, and develops the following transformation rules: d_ d_ dZd_ = dt~ dr d 9 dxi d d k dZ + ' ] (1.3) 3 2 ~ dZ k { dxi dZ _ dx dtdZ + d dx dZ' { k " } He then applies them to the species transport equation (the diffusion term has been modeled using Fick's L a w i.e.: V • ji « — pDVYi) dY : _^ + pH • VYi - V • (pDVYi) = pui (1.5) resulting in the transformed species transport equation: T h e operator R is defined in [29], and involves derivatives w i t h respect to the higher dimensions (zT ,zT ). T h e first observation one may make is that the result, E q . 1.6, has not been to transform the species transport equation from a function of (~x,t) to one which is a function of (Z, t), but rather (Z, Z%, Z3, r ) . T h i s is obviously a mathematical necessity — we began w i t h a 4-dimensional coordinate system, so certainly 4 new orthogonal 2 3 CHAPTER 1. BACKGROUND 6 coordinates are necessary for a valid transformation. However, claiming that the variation along the coordinates normal to the gradient of mixture fraction is negligible, Peters discards the terms included i n R(Y{), and hence the dependence on Z , Z . In doing so, it is necessary to assume that the structure of the reaction zone is 1-dimensional; here we see the first use of the assumption that a turbulent flame can be likened to a laminar counterflow flame, which does have a 1-dimensional structure. Second, the appearance of the scalar dissipation, x = 27J (VZ • VZ), is noteworthy. T h i s implies that the attempt to eliminate the dependence on ~x was incomplete: the term's presence indicates coupling between the two coordinate systems. In order to solve the flamelet equations in mixture fraction space, a model for x must be provided. T h e typical model is to assume that x — XoX*(^), where x*(%) is some presumed functional form for the dependence on mixture fraction — the choice varies between formulations. E a r l y formulations utilize the form obtained by Peters ([29]) from an analytical solution for the spatial variation of mixture fraction in an infinite 1-D mixing layer, w i t h boundary conditions Z(y —»• —oo) = 0, Z(y —> oo) = 1. T h e resulting functional form is: 2 3 -2[erfc" (2Z)] 1 X*(Z) ^ = — e (1.7) -j. r ' 2 -2[erfc-l(2Z )] 0 2 V ' Alternatively, from an analytical solution to a symmetric one-dimensional m i x i n g layer where the fuel originates from the plane of symmetry, P i t s c h et al. [34] obtain: * • < * > - f £ i ™ Despite the observation by K i m et al. [19] that the latter form is inconsistent w i t h the transport equation for the pdf of Z , it now commonly used, e.g., in [32]. Other forms for this relationship are typically used i n C M C formulations, and w i l l be discussed later. Once such a functional form has been chosen, xo is the only remaining coupling to the flow field. Finally, it is important to note that the Eulerian time has been replaced by a Lagrangian 'flamelet time', r. A t this point, flamelet models spread out into several different directions w i t h varying success. Before considering some of the approaches, it is worthwhile to note the implication of the assumptions made thus far. In disregarding the variation along coordinates non-normal to the stoichiometric interface, there is an assumption that the local structure of the turbulent CHAPTER 1. BACKGROUND 7 flame resembles a laminar counterflow flame [29]. Physically, this requires that the turbulent structure of the flow field must not influence the internal flame structure. T h e only allowed effect is that the turbulence can strain the interface, increasing the local gradients and the rate of mixing. For this to be valid, the smallest scale of the turbulent structure, the Kolmogorov length, must be larger than the thickness of the reaction zone. A n alternative statement is that the chemical reactions must occur on a time scale shorter t h a n the Kolmogorov time scale. T h i s can be visualized on the B o r g h i diagram, where the ratio of the turbulence intensity to the laminar burning velocity is plotted against the ratio of large scale turbulent eddy size to the flame thickness on a logarithmic scale. Such a plot is given in Figure 1.1, adapted from Lentini [23]. A l t h o u g h laminar burning velocities are reasonable for premixed combustion, for non-premixed combustion Lentini proposes using revised definitions for the reference scales; specifically: Iref = \/vtc U ref = yj^ (1.9) Constant values of two dimensionless numbers, the K a r l o v i t z number (Ka) and the Damkohler number (Da) appear as straight lines on such a plot. T h e K a r l o v i t z number is defined as the ratio of chemical time scale of the laminar flame to the turbulent (Kolmogorov) time scale [46]: Ka=^ t k T h e Damkohler number is defined as the ratio between the (integral) time scale and the chemical time scale [46]: flow Da^. tc Using Lentini's revised definitions of the non-premixed length and time scales, tc i n the above definitions is replaced by trer = defined as i n E q . 1.9. In keeping w i t h the earlier criteria, strictly speaking the validity of the flamelet models require t < tk, or Ka < 1. T h i s region is the portion of Figure 1.1 below the line shown for Ka = 1. There is also a regime where the turbulence interacts w i t h the flame, but the flow does not completely overwhelm the flame. T h i s region, commonly known as the 'Thickened flame' or 'perturbed flamelet' regime falls in the region of Figure 1.1 between the c CHAPTER 1. BACKGROUND togfy Figure 1.1: Borghi Diagram 8 CHAPTER 1. BACKGROUND 9 lines Ka — 1 and Da — 1. For Da < 1, the flow field completely distorts the flame structure and the flamelet hypothesis is senseless. Lentini evaluated the chemical time scale and compared it to the flow time scale during flamelet computations for a turbulent burner. W i t h i n the flamelet regime, he obtained reasonable results. However, just outside the flamelet regime i n the perturbed flamelet regime Lentini noticed under-prediction i n minor species such as OH of more t h a n 50%. T h e implication is that while conceptually useful and convenient for providing engineering level solutions, basic flamelet models should be treated w i t h caution due to a limited regime of applicability. W i t h this in m i n d , let us now t u r n our attention to the various types of flamelet models. T h e most straightforward implementation of flamelet models is to neglect transient effects, such as extinction and re-ignition, and solve for steady solutions to E q . 1.6. Setting ^ = 0 and neglecting R(Yi) leads to a set of ordinary differential equations in mixture fraction space for the component species mass fractions and (including the corresponding equation obtained from a transformation of the energy equation) temperature, i n which the scalar dissipation (xo) appears as a parameter. Solutions to the system represent a balance between mixing (associated w i t h the given scalar dissipation) and chemical reaction. To use such a model, the typical approach is to solve the set of equations in mixture fraction space, and tabulate the resulting solutions in a 2-D lookup table as functions of Z , x- T o then recover the mean value of an arbitrary scalar, a, the problem is simply to evaluate: (1.10) A common assumption (e.g., [23, 29]) is to assume statistical independence of Z and xo- T h e n , the flow field is solved for sufficient moments of the mixture fraction and scalar dissipation to construct a pdf of these two quantities using a presumed form. Species concentrations and temperature are obtained from the lookup table using the pdf's to weight the individual flamelet solutions. A typical choice for the presumed form for the mixture fraction is the B e t a pdf (e.g., [24, 44]); although early (limited) use was made of a 'clipped-gaussian' pdf [25] it has lost favour. T h e B e t a pdf choice has been tested by W a l l et al. [45]; it has also been noted by L i e w et al. [24] (by reference to Jones, 1979) that the choice of the pdf for the conserved scalar is of relatively minor importance. To construct the B e t a pdf, transport equa- CHAPTER 1. BACKGROUND 10 tions are typically solved in physical space for the mean and variance of mixture fraction. T h e parameters of the pdf are then chosen to match the first two moments. For the scalar dissipation, it is well established that the correct functional form of the pdf is the log-normal pdf [23, 24]. A discussion of models available for the mean scalar dissipation rate are presented by Sanders and Gokalp [39]. T h e simplest algebraic model is based on an assumption of equality for the length and time scales between mechanical and scalar turbulence, and leads to the approximation: * k where c has a value of approximately 2. T h e second parameter for the lognormal pdf is often quite poorly estimated; based on observations made by Sreenivasan et al. (1977), Lentini [23] assigns the value a = 2, while L i e w et al. [24] develops an algebraic relationship containing the Kolmogorov and integral scales, as well as various constants (some of which depend on the flow geometry). In the commercial C F D code Fluent [13], the L a m i n a r Flamelet model implementation disregards the second parameter entirely and instead assumes a delta-pdf, based on the mean only. Despite the attractiveness of solving the chemistry problem before, and independent of, the flow problem, the use of steady flamelet modeling i n this way has some serious drawbacks. A s pointed out by Barths et al. [3], the history of flamelet solutions is relevant — it is as though the solutions have a 'memory' of their prior experiences, which, although not infinitely long, invalidates the steady state assumption for some cases. There has been some discussion regarding the appropriate criteria on which to assess the validity of the assumption of steadiness [30, 41, 42]. Regardless, recent developments i n the flamelet area have concentrated on unsteady models. Certainly, there are processes — notably ignition and extinction — for which an unsteady model is essential. Unsteady flamelet models have the feature of retaining the time derivative i n E q . 1.6 and attempting to relate the Lagrangian r to an Eulerian time/space coordinate system. T h i s is frequently done by introducing marker particles for each flamelet into the system and allowing them to be transported by an advection/diffusion equation. A separate type of marker particle is assigned for each independent flamelet history allowed, and a new transport equation added. A t each point in space, the portion of each flamelet solution present is taken as the fraction of the total number of marker particles present which are assigned to that flamelet. Developed x CHAPTER 1. BACKGROUND 11 jointly by Barths, Pitsch, and Peters [3] as the 'Representative Interactive Flamelet' ( R I F ) , a good description of what is now known as the ' E u l e r i a n Particle Flamelet M o d e l ' can be found i n [2]. In this model, the flow field is updated for flow and turbulence quantities, as well as scalar transport of mixture fraction, its variance, the scalar dissipation, and the transport equations for each marker particle. Following the flow time step, the scalar dissipation is extracted and passed to the flamelet code, which then advances the flamelet solutions in time. In this arrangement, the flamelet code can take much smaller time steps than the flow to resolve the chemical timescale. T h e species concentrations are then integrated over the pdf of mixture fraction (weighted by the particle fraction representing the corresponding flamelet) and used to determine the mixture enthalpy. T h i s enthalpy is then used to determine a new temperature to be fed back to the flow calculation. Transport equations for energy or individual concentrations are not solved i n physical space; aside from the temperature interaction, there is no feedback from the flamelet code to the solution i n physical space. These methods have been applied to problems such as the treatment of internal combustion engines [2] and m i l d combustion burners in furnaces [10] w i t h reasonable success at predicting temperature and major species fields. W h i l e the majority of the use of flamelet models is i n R A N S simulations, the model has also been applied to L E S computations. C o o k et al. [11] describe an application of steady flamelets i n L E S , using much the same approach as the earlier description; the various species concentrations are determined through a lookup table based on mixture fraction and mean scalar dissipation. A n implementation making use of an unsteady formulation is given by P i t s c h [33], P i t s c h and Steiner [35]. In the latter work, the formulation and solution procedure closely resembles that commonly employed for the modeled C M C equations which will be described in the next section. T h e flamelet formulation, however, leaves the specific omissions from the equations unclear. B r o a d assumptions regarding structure early in the derivation, as we have seen, limit the applicability of the method. Corrections introduced to arrive at equations similar to the C M C counterparts do produce improvement, but it is unclear whether any additional unclosed terms may have been overlooked. CHAPTER 1. BACKGROUND 1.3 12 Conditional Moment Closure Conditional Moment Closure ( C M C ) methods seek to reduce the magnitude of the fluctuations about the mean quantities by conditioning the means o n one or more closely correlated variables. Proposed independently by Bilger and K l i m e n k o and unified i n [22], the reward for reduced fluctuations about the conditional mean quantities is that simple first order closure for the chemical source term is often sufficient. T h e dramatically different derivations suggested by K l i m e n k o [21] (based on pdf methods) and Bilger [5] (applying a decomposition to conventional transport equations) producing the same basic equations lends substantial credibility to the method. Typically, the mixture fraction is used as a conditioning variable. A scatter plot of temperature against mixture fraction is shown i n Figure 1.2, obtained from experimental measurements of the Sandia F l a m e ' D ' [1] at 30 diameters downstream from the nozzle; this flame will be described i n detail when it is considered i n Chapter 2. A t the axial location shown, only moderate to low levels of extinction are expected. It is apparent that there is a significant correlation between temperature and the mixture fraction, and the variation about the conditional mean (shown i n black) when mixture fraction is used as a conditioning variable will be much smaller t h a n that about the unconditional mean. It is worthwhile to briefly consider the origin of the C M C equations, so that the modeling assumptions are clearly visible. T o derive the C M C equations using the decomposition approach, Bilger makes use of the decomposition: Y t) = Q(Z(1t, t), t, t) + Y"Cg, y). (1.11) Y" is the conditional fluctuation, and Q = (Y\Z = £); we will persist i n the use of Q i n this section only for brevity and consistency w i t h Bilger's notation. T h e quantity {Y\Z = C) should be interpreted as the conditional average of Y for the conditioning variable Z taking the value of £; £ forms the sample space for Z. To obtain the various terms appearing i n a general species transport equation, as i n E q . 1.5, E q . 1.11 is differentiated as CHAPTER 1. BACKGROUND 13 2200 Mixture Fraction Figure 1.2: Scatter P l o t of Temperature vs. M i x t u r e Fraction for Sandia F l a m e ' D \ x/D = 30 [1] CHAPTER 1. 14 BACKGROUND necessary: dY = dQ dt dQdZ dt dC dt V r = VQ + ^ V • (pDVY) dY^ dt V Z + VY" = V • (pDVQ) + | | v pD(VZf ^ + ' (1.13) • (pDVZ) (I-") DVZ.V^ d + K P + V • (pDVY") Substituting the above into E q . 1.5, using that Z also satisfies E q . 1.5 w i t h LJ = 0 to simplify, and taking the conditional average, Bilger arrives at the 'unclosed' C M C evolution equation for species concentration: <p|C> ?r + ' 0*10 (p\0 \ • V Q- i - s / * ~e <P|C> (xlO ^ Q xrib/ \/vi\>/ _Qr2 = (p\Q(u\0 +e Q +e Y (1.15) where: e Q = ( V • (pDVQ) + pDVZ • V ^ | C ) / dY" _^ \ {p— + pH • VY" - V • (DpVY")\(^ (1.16) (1.17) In [22], it is proposed that e^, representing diffusive effects, is small when Re is large, and can be neglected. The proposed closure for ey is: .. v• « P I C > ( u " Y " \ Q P ( 0 ) . . v• « P I C > P ( 0 ) , „ „ Where a gradient transport model has been used to close the term ( u "Y"\Q, similar to that used i n the development of the transport equation for turbulent kinetic energy, where the term u'u' is closed by the approximation u'u' « —D Vu. T h i s approximation is commonly used to close terms involving a correlation of fluctuations i n velocity and a scalar, e.g., i n [31]. T h e t CHAPTER 1. BACKGROUND 15 assumptions used for closure up to this point, along w i t h a linear approximation for the conditional average of velocity, have been broadly adopted, e.g.,in [12, 16, 28]. T h e final closure approximation is the functional form for the conditional scalar dissipation. K i m et al. [16, 18] and Devaud and B r a y [12] all make reference to an expression obtained by G i r i m a j i (1992) by consideration of the pdf transport equation. K i m and H u h [16] also make reference to the counterflow mixing model used i n flamelet calculations. In their work on spray autoignition, K i m and H u h [20] use what they refer to as an 'amplitude mapping closure' which produces a functional form resembling the counterflow mixing model; the same model is employed by Mastorakos and Wright [28], which also treats spray autoignition. T h e computational cost associated w i t h solving E q . 1.15 using the described approaches to obtain closure is significant. However, the conditional averages tend to vary slowly i n space. In fact, if the correlation between the dependent scalars and the conditioning variable were perfect, then there would be no need to retain a spatial dimension! T h i s has been done, e.g., by Swaminathan and Bilger [42], and is known as 0-dimensional C M C . In a homogeneous field, the conditional transport terms vanish, and the C M C species evolution equation reduces to a form reminiscent of the flamelet equations: ~i d + (xlO = WO (i-w) A l t h o u g h this form of the C M C equations bears resemblance to the flamelet equations, there is an important difference i n proceeding from the current state. W h i l e the assumptions made early in the derivation of the flamelet models disregard much of the physical behaviour, the rigorously justifiable derivation of the C M C equations, up to the point where closure assumptions are made, requires only the fundamental assumption that the variations i n the various scalars are well correlated w i t h the conditioning variable. T h i s means that the framework to reconsider or improve the closure for a given term is available. A l l of the necessary physics are present, and we are consciously aware of the assumptions made to neglect or model various terms. T h i s makes the model ideally suited to improvements not only i n modeling, but in basic understanding as well. W h i l e quite successful i n some circumstances, variations about the conditional averages in violation of the fundamental assumption can lead to errors in first moment closure for the conditional reaction rates. T h e reme- CHAPTER 1. BACKGROUND 16 dies for this vary; allowing for spatial variation of the conditional averages, i n 1 or more dimensions, can alleviate the difficulty. T h i s is successful as large variations about the conditional means are often the result of localized effects, such as extinction due to a small region of locally high strain. Extensions allowing for second-order closure have also been developed, such as described by K i m et al. [17], Mastorakos and Bilger [27]. However, this requires modeling of several unclosed terms i n the transport equation for the conditional fluctuations i n temperature; presently the models proposed are somewhat dubitable. Where significant local extinction is present, at scales unresolvable by any reasonable spatial discretization, the variation about the conditional average based on a single conditioning variable may be extreme; the left hand portion of Figure 1.3 is a scatter plot of temperature vs. mixture fraction from a direct numerical simulation of a flame w i t h substantial levels of local extinction performed by Sripakagorn given as Figure 1 i n [9]. In a case where the variation about the conditional average is significant, as i n Figure 1.3, an attractive alternative to estimating higher moments is to use a second conditioning variable. To be effective, it must be well correlated w i t h the scalar variations not correlated w i t h mixture fraction. W o r k using a scalar field defined such that its gradient is orthogonal to mixture fraction, the ideal behaviour for a second conditioning variable, has been undertaken w i t h promising results [26]. In a realistic combustion simulation, the choice of second conditioning variable is unclear. C h a et al. [9] attempted to use the scalar dissipation, after dividing by a presumed functional dependence on mixture fraction to produce a variable somewhat independent of mixture fraction (as shown by Bushe and Steiner [7]). However, as explained earlier by Bushe and Steiner [7], although this quantity has many of the attributes desirable in a second conditioning variable it w i l l ultimately prove unsuitable. A s the strain on a flame increases, the local gradients become steeper and the rate of mixing increases. For a m i l d increase i n mixing rates a corresponding increase in reaction rates, and a moderate change i n the species conditional averages is observed. W h e n the rate of mixing becomes too great, the temperature and intermediates are mixed away too quickly to be replaced by chemical reaction and the flame is extinguished. Once a flame has been extinguished, there is no longer a relation between the temperature and scalar dissipation — regardless of how the flame is strained, the temperature w i l l not increase until re-ignition occurs. Re-ignition involves either the rate of strain dropping below the level at which autoignition is likely, or an external event such as a passing edge flame [38]. Further complicating matters, there CHAPTER 1. BACKGROUND 17 is a temporal history effect. A n increase i n scalar dissipation to a level that would otherwise cause extinction can be prevented from doing so if the scalar dissipation is lowered again before the temperature and intermediate concentrations have dropped too far. T h e various scalar dissipation perturbation possibilities give rise to a plethora of possible scalar fields for the same scalar dissipation magnitude. T h e right side of Figure 1.3, which forms a portion of Figure 2 i n [9], illustrates the multiplicity of possible temperatures for a given scalar dissipation at constant mixture fraction. In the figure, which was extracted from the same field as the left portion of Figure 1.3, all of the points plotted have a constant mixture fraction, so an ideal conditioning variable would have a more or less unique relationship w i t h temperature. Obviously, this is not so for the scalar dissipation. Z 'Og Xst Figure 1.3: Relation Between Temperature and M i x t u r e Fraction, L o g of Scalar Dissipation. Portion of Figures 1 and 2 i n [9] T h e preceding discussion suggests that there are two factors to consider when evaluating reaction rates: the mixture fraction, and the elusive additional quantity that somehow represents the cumulative effects of the m i x i n g field on the flame throughout its recent history up to a n d including the present. In the next section the detailed derivation of an approach which has the potential to incorporate these two considerations w i l l be considered. 1.4 Conditional Source-Term Estimation C o n d i t i o n a l source-term estimation ( C S E ) ultimately uses the same first order closure based on conditional means to evaluate the chemical source terms CHAPTER 1. BACKGROUND 18 employed i n C M C . Rather than solving directly for the conditional averages, the flow is formally solved for the unconditional means i n physical space. T h e conditional means, which are assumed constant i n space over ensembles of cells, are related to the unconditional means and the pdf of the conditioning variables through integral equations which are then inverted for a discreet approximation to the conditional means. It is these conditional means which are then used to close the chemical source term; the mean reaction rates are then recovered by integration over the pdf of the conditioning variables. T h i s method was tested in an a priori sense against a D N S database [7] and was found to be effective when little extinction was present using only the mixture fraction as a conditioning variable. T h e addition of a modified form of the scalar dissipation as a second conditioning variable produced success when extinction was present. Later, the method was used to obtain closure for the chemistry in a L E S calculation of the Sandia F l a m e ' D ' using 2-step chemistry w i t h substantial success [40]. M o r e recently, Bushe and Steiner [8] proposed that the conditional means and source terms be decomposed into a linear combination of solutions to the L a m i n a r Flamelet equations. In this formulation, the integral equations for the mean temperature only are inverted to determine the appropriate combination of flamelet solutions; the conditional source terms associated w i t h the flamelets are then combined and integrated over the pdf of mixture fraction to close the enthalpy source term. T h i s method produced promising results in an a priori test i n [8]. Later, B l a i r [6] used this formulation in a K I V A - I I simulation of a t u r b u lent methane jet injected into hot air. However, B l a i r ' s implementation was unable to successfully predict the ignition delay when compared to experimental data. In the remainder of this section, the general formulation of C S E w i t h L a m i n a r Flamelet Decomposition used throughout this thesis w i l l be presented. Following Bushe and Steiner [7], first note that the unconditional average of any scalar quantity can be related to the conditional average of that quantity by employing the probability density function of the conditioning variable. Using a 'Favre density function' p ((; Z, Z" ^j 2 [22], the relationship then holds for Favre averages: (1.20) A s per [8], decomposition of the conditional average i n question into a sum- CHAPTER 1. BACKGROUND 19 mation of basis functions, e.g., N <T|C> « 5 > n <T|0 (1.21) n n=l can be expected to be valid for finite N provided the basis functions (T\Q capture characteristics of the true conditional average. Notably, it must be possible to reconstruct the potentially sharp gradients w i t h respect to mixture fraction. Substitution of E q . 1.21 into E q . 1.20 produces: n N (1.22) Jo .n=l A t this point, we can see that E q . 1.22 is a Fredholm E q u a t i o n of the first kind: g(t)= j K(t,s)f(s)ds, (1.23) Ja with s = ( a= 0 b= 1 (1.24) 9 = f (1.25) K=p(Q (1.26) N (1.27) / = 5> (T|C> . n n Choosing an ensemble of J cells over which a =i..JVj=i..j is constant, we can write E q . 1.22 numerous times, and form a system of equations resembling (18.4.5) in [37] (adapting to match notation): n N 9j (1.28) Jo n=l G i v e n J distinct 'observations' of Qj (possibly contaminated w i t h noise), each w i t h a different kernel p(Q which filters the underlying function for different aspects of its behaviour, Press et al. [37] suggest replacing the integral by a quadrature and solving the resulting system of equations in a least squares CHAPTER 1. 20 BACKGROUND sense to obtain an approximation of the underlying function at the quadrature points. A s the resulting system of equations is typically ill-posed (it could only be well posed if each of the N observations came conveniently supplied w i t h a response kernel that was a delta function at each of the quadrature points), the recommendation i n [37] is to add a priori information i n the form of smoothing to stabilize the solution. In our case, we do not actually wish to solve for the function / ( £ ) = ^2n=i ™ CHOn ^ quadrature points i n £-space, but rather for the basis function coefficients a . A s w i l l be discussed i n the following chapters, that a „ should have any elementary functional dependence on the position i n the library, n , is unphysical. Instead, a more physical a priori belief is that, if started from a physically realizable condition, the solution should evolve smoothly i n time. T u r n i n g now t o the algebra necessary to accomplish this, we revisit a rearranged version of E q . 1.22 a n d switch to vector/matrix notation: a a n (1.29) Jo Forming a i f - p o i n t quadrature i n £-space, / P ( C ) / ( C K « 5>(C*)/(CkKAfc =i J o fc K 1 /" 5 / E K Cfc+1 /(a) = £ > / ( & ) • (i-3o) p(C)dC fe=i fc=i A one-sided approximation (e.g., tbx = /^_ p(C')^C) is substituted for il>i,4>K- U s i n g E q . 1.30 to discretize E q . 1.29 a n d writing it for each of J computational cells i n an ensemble over which we assume ~a is spatially homogeneous i n the fashion of E q . 1.28, neglecting rij, we obtain: 1 ^j=l,fc=l • • • ^1,K f(T\O =l, =l k n a \ T j J \ ... jjjcj V (T\C) ,i K (T\C)K,N) (1.31) or, where A is defined by comparison to E q . 1.31, Aa. (1.32) CHAPTER 1. BACKGROUND 21 W e can also write an equation resembling E q . 1.20 for any scalar of i n terest. Repeating the above procedure for an arbitrary species mass fraction Yi, and defining: /^I0 f c = 1 , n = 1 ••• Bi = (Yi\O hN (1.33) we can also form the system of equations: Yi « Bi~a. (1.34) Following Bushe and Steiner [8], we suppose that it is possible to reconstruct the flame of interest from a linear combination of N basis functions, where the n basis function solution for (T\() , ( l i | C ) „ are corresponding scalar profiles from the same quasi-physical flame. In this situation, solution of either E q . 1.32 or E q . 1.34 for any i should yield the same ~a. However, we must recall the ill-posedness of these equations; not only is there a question of existence, there is a significant question regarding uniqueness. Later, i n Sections 2.3 and 3.2.3, it will be illustrated that using a typical library of basis functions created from solutions to the Laminar Flamelet equations there are many solutions where the part of the solutions described by (T\Q is nearly identical, yet the conditional averages for the various species mass fractions vary widely. To extract as much information as possible from the flow, it is proposed that the system t o be solved be expanded by solving, simultaneously, E q . 1.32 and E q . 1.34, using as many species mass fractions as necessary to differentiate between the basis functions i n the library. W e now have: th n n M~a ?s R, (1.35) with: ( A \ M = R = \ B J I (1.36) \YiJ Despite our hope that E q . 1.35 contains sufficient information t o select the most suitable combination of basis functions, it is still potentially illposed. We have doubts b o t h about the uniqueness of the solution and the CHAPTER 1. BACKGROUND 22 existence of a suitable solution. T h e ill-posedness is due in part to neglecting the noise rtj (caused by accumulated errors in the closure, poor predictions of the pdf of mixture fraction, and errors due to discretization E q . 1.30, etc.), and i n part due to insufficient filtering of the underlying function for different aspects of its behaviour by the kernel p(Q- There is also the possibility that it is impossible to reconstruct the conditions in the flow exactly due to the finite size of the library of basis functions. To address these issues, we t u r n here to the more general Tikhonov R e g ularization [43], rather than the Linear Regularization recommended i n [37] and used by Bushe and Steiner [8]. W h i l e b o t h of these regularization approaches include introducing a priori information into the solution, the general case allows us more freedom in choosing that information. T h i s allows us to more carefully introduce our prior belief about the solution behaviour and avoid contamination of the solution by ill-chosen and non-physical constraints. Tikhonov suggests reformulation of the problem as a minimization problem, where we simultaneously minimize the residual i n the problem we wish to solve as well as the deviation of the solution from an expected value: (1.37) T h i s reformulation as a minimization problem addresses the question of existence of a solution. However, it means that we must take care i n assembling E q . 1.35 to normalize b o t h sides of each equation, lest the m i n i m i z a tion procedure place greater emphasis on satisfying those equations involving scalars w i t h larger absolute magnitudes. T h e alternative formulation of the minimization which perhaps makes this more clear is: (1.38) Considered i n this form, it is obvious that we can adjust the influence of each scalar by adjusting the coefficients A .../+i. However, since we wish to satisfy the equations for all scalars w i t h the same weight, we can accomplish our objective simply by normalizing E q . 1.32 and E q . 1.34 by the m a x i m u m present so that all of the scalars fall in the range [0,1]. Returning to the original minimization problem E q . 1.37, the solution in a least squares sense 0 CHAPTER 1. BACKGROUND 23 is given by the solution to the linear system of equations: ( M * M + A J ) ~a = M*~R + A c * . (1.39) 0 W e now have an opportunity to provide an arbitrary initial solution of , and allow the minimization to select a new solution ~a which deviates from the original only to the extent that the current flow field, incorporated into M, justifies. W h i l e the smoothing used for stabilization i n earlier work by Bushe and Steiner [8] implies a functional dependence of a on n, the claim i n this work is that such functional dependence is unphysical. Insistence o n such dependence implies that if, for example, the order i n which the basis functions are stored or the spacing between the basis functions i n terms of flamelet time is altered, the physics captured will be altered. A s w i l l be discussed later, this is clearly unacceptable. In the two applications of this method discussed i n Chapters 2 and 3, the stabilizing solution ~a° will be chosen to be the solution for "0? from the previous time step or iteration as appropriate. 0 n T h e above formulation has a great deal i n common w i t h that proposed by Bushe and Steiner [7, 8] and used by B l a i r [6], Steiner and Bushe [40]. However, there are significant differences, which are believed to be responsible i n part for the relative success of the implementations presented i n this thesis. T h e two most notable which have appeared so far include: • Regularization Method. A l l of the previous formulations imposed smoothing on the quantity sought by the integral inversion. In [40], progress variables were solved for, which do tend to be smooth i n mixture fraction space. However, this necessitates inversion for a separate progress variable for each reaction i n the mechanism used, which is a prohibitively large set for realistic chemistry. In [6, 8], the smoothing was for the coefficient i n basis function index space; the drawbacks to this have already been touched on, and will be discussed further i n the following chapters. • Inversion Based on Multiple Species. W h i l e early work attempted to determine the combination of flamelet solutions that best describe the state of the system based on temperature alone, this formulation allows the determination to be made to ensure agreement w i t h the mean fields for minor species, which discriminate much more clearly between the flamelets. CHAPTER 1. BACKGROUND 24 A d d i t i o n a l differences between this and previous implementations are present; however, they are of a problem-specific nature; further discussion will be left until the final chapter. 25 REFERENCES REFERENCES [1] Proceedings of the Third International Workshop on Measurement and Computation of Turbulent Nonpremixed Flames, Boulder, Colorado, 1998. U R L http://www.ca.sandia.gov/TNF/3rdWorkshop/TNF3.html. Section 2 - Sandia / Darmstadt P i l o t e d CH /Air Flame D. 4 [2] H . B a r t h s , C . Hasse, et al. Simulation of combustion i n direct injection diesel engines using an eulerian particle flamelet model. In TwentyEighty Symposium (International) on Combustion, pages 1161-1168. T h e Combustion Institute, 2000. [3] H . B a r t h s , H . Pitsch, and N . Peters. 3d simulation of d i diesel combustion and pollutant formation using a two-component reference fuel. Oil & Gas Science and Technology: Revue de I'institut francais du petrole, 54:233-244, 1999. [4] R. Bilger. Turbulent jet diffusion flames. Progress in Energy and Combustion Science, l(2-3):87-109, 1976. [5] R. Bilger. C o n d i t i o n a l moment clsoure for turbulent reacting Physics of Fluids A - Fluid Dynamics, 5(2):436-444, 1993. flow. [6] C . B l a i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of B r i t i s h C o l u m b i a , 2003. [7] W . K . Bushe and H . Steiner. Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11 (7): 1896-1906, J u l y 1999. [8] W . K . Bushe and H . Steiner. L a m i n a r flamelet decomposition for conditional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [9] C M . C h a , G . Kosaly, and H . Pitsch. Modeling extinction and reignition i n turbulent nonpremixed combustion using a doubly-conditional moment closure approach. Physics of Fluids, 13(12):3824-3834, 2001. REFERENCES 26 [10] P. Coelho and N . Peters. Numerical simulation of a m i l d combustion burner. Combustion and Flame, 124:503-518, 2001. [11] A . W . Cook, J . J . Riley, and G . Kosaly. A laminar flamelet approach to subgrid-scale chemistry i n turbulent flows. Combustion and Flame, 109:332-341, 1997. [12] C . Devaud and K . Bray. Assessment of the applicability of conditional moment closure to a lifted turbulent flame: first order model. Combustion and Flame, 132:102-114, 2003. [13] Fluent 6.1. Fluent Inc., 2003. [14] D. M . G . Gregory P. S m i t h et al. Gri-mech 3.0. U R L h t t p : //www. me. b e r k e l e y . e d u / g r i - m e c h / . [15] J . Huang, P. H i l l , W . Bushe, and S. M u n s h i . Shock-tube study of methane ignition under engine-relevant conditions: experiments a n d modeling. Combustion and Flame, 136(l-2):25-42, 2004. [16] S. H . K i m and K . Y . H u h . Use of the conditional moment closure model to predict N O formation i n a turbulent C H 4 / H 2 flame over a bluff-body. Combustion and Flame, 130:94-111, 2002. [17] S. H . K i m , K . Y . H u h , and R. W . Bilger. Second-order conditional moment closure modeling of local extinction a n d reignition i n turbulent non-premixed hydrocarbon flames. Proceedings of the Combustion Institute, 29:2131-2137, 2002. [18] S. H . K i m , K . Y . H u h , and R. A . Fraser. Modeling autoignition of a turbulent methane jet by the conditional moment closure model. Proceedings of the Combustion Institute, 28:185-191, 2000. [19] S . - K . K i m , S . - M . K a n g , and Y . - M . K i m . Flamelet modeling for combustion processes a n d N O x formation i n the turbulent nonpremixed C O / H 2 / N 2 jet flames. Combustion Science and Technology, 168:47-83, 2001. [20] W . T . K i m and K . Y . H u h . Numerical simulation of spray autoignition by the first-order conditional moment closure model. Proceedings of the Combustion Institute, 29:569-576, 2002. REFERENCES 27 [21] A . Klimenko. Multicomponent diffusion of various admixtures i n turbulent flow. Fluid Dynamics, 25(3):327-334, 1990. [22] A . K l i m e n k o and R. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [23] D. Lentini. Assessment of the stretched laminar flamelet approach for nonpremixed turbulent combustion. Combustion Science and Technology, 100:95-122, 1994. [24] S. Liew, K . Bray, and J . Moss. A stretched laminar flamelet model of turbulent nonpremixed combustion. Combustion and Flame, 56:199213, 1984. [25] F . Lockwood and A . Naguib. T h e prediction of the fluctuations i n the properties of free, round-jet, turbulent, diffusion flames. Combustion and Flame, 24:109-124, 1975. [26] J . Lozada-Ramirez. Conditional moment closure for methane oxidation using two conditional variables and stochastic processes. Master's thesis, Univesity of B r i t i s h C o l u m b i a , 2004. [27] E . Mastorakos and R. Bilger. Second order conditional moment closure for the autoignition of turbulent flows. Physics of Fluids, 10(6), June 1998. Letter. [28] E . Mastorakos a n d Y . Wright. Simulations of turbulent spray autoignition w i t h elliptic conditional moment closure. In Proceeding of the European Combustion Meeting, 2003. [29] N . Peters. L a m i n a r diffusion flamelet models i n non-premixed turbulent combustion. Progress in Energy and Combustion Science, 10:319-339, 1984. [30] N . Peters. Comment and reply on the "Assessment of combustion and submodels for turbulent nonpremixed hydrocarbon flames". Combustion and Flame, 116:675-676, 1999. [31] N . Peters. Turbulent Combustion. Y o r k , 2000. Cambridge University Press, N e w REFERENCES 28 [32] H . Pitsch. Unsteady flamelet modeling of differential diffusion i n turbulent jet diffusion flames. Combustion and Flame, 123(3):358-374, 2000. [33] H . Pitsch. Improved pollutant predictions i n large-eddy simulations of turbulent non-premixed combustion by considering scalar dissipation rate fluctuations. Proceedings of the Combustion Institute, 29:1971-1978, 2003. [34] H . Pitsch, M . C h e n , and N . Peters. Unsteady flamelet modeling of turbulent hydrogen-air diffusion flames. In Twenty-Seventh Symposium (International) on Combustion, pages 1057-1064. T h e Combustion Institute, 1998. [35] H . P i t s c h and H . Steiner. Scalar mixing and dissipation rate i n largeeddy simulations of non-premixed turbulent combustion. Proceedings of the Combustion Institute, 28:41-49, 2000. [36] S. Pope. P d f methods for turbulent reactive flows. Progress in Energy and Combustion Science, 11:119-192, 1985. [37] W . Press, S. Teukolsky, W . Vetterling, a n d B . Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, U K , 1992. [38] G . Ruetsch, L. Vervisch, and A . L i n a n . Effects of heat release on triple flames. Physics of Fluids, 7(6): 1447-1454, June 1995. [39] J . Sanders and I. Gokalp. Scalar dissipation rate modelling i n variable density turbulent axisymmetric jets a n d diffusion flames. Physics of Fluids, 10(4):938-948, A p r i l 1998. [40] H . Steiner and W . Bushe. Large eddy simulation of a turbulent reacting jet w i t h conditional source-term estimation. Physics of Fluids, 13(3): 754-769, M a r c h 2001. [41] N . Swaminathan. Flamelet regime i n non-premixed combustion. Combustion and Flame, 129:217-219, 2002. [42] N . Swaminathan and R. Bilger. Assessment of combustion and submodels for turbulent nonpremixed hydrocarbon flames. Combustion and Flame, 116:519-545, 1999. REFERENCES 29 [43] A . Tikhonov and V . I. Arsenin. Solutions of ill-posed problems. Halsted Press, Washington, W i n s t o n , New York, 1977. Translation of M e t o d y resheniia nekorrektnykh zadach. [44] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress in Energy and Combustion Science, 28:193-266, 2002. [45] C . W a l l , B . J . Boersma, and P. M o i n . A n evaluation of the assumed beta probability density function subgrid-scale model for large eddy simulation of nonpremixed, turbulent combustion w i t h heat release. Physics of Fluids, 12(10):2522-2529, October 2000. [46] J . Warnatz, U . Mass, and R. Dibble. Combustion: physical and chemical fundmentals, modeling and simulation, experiments, pollutant formation. Springer-Verlag, New York, 2nd edition, 1999. 2. SIMULATION OF T H E SANDIA FLAME 'D' CHAPTER 2.1 2. SIMULATION OF THE SANDIA FLAME 'D' 31 Introduction Computer simulations of turbulent reacting flows seeking to solve for mean quantities i n either a filtered (Large E d d y Simulation ( L E S ) ) or Reynolds Average ( R A N S ) sense are hampered by the closure problem i n turbulent combustion, that is, determining the mean source-terms due to chemical reaction. T h e difficulty arises because unresolved fluctuations about mean quantities dominate the reaction rate expressions. Conditional moment closure ( C M C ) [54] seeks to address this issue by solving for conditional means where the conditioning variable accounts for much of the unresolved variation, leading to significantly smaller fluctuations about the mean quantities. T h i s allows for much improved closure of the source-term compared to evaluating the source-term w i t h unconditional means. However, this method requires adding another dimension to the solution for each conditioning variable leading to a substantial increase i n computational cost. Further, although applicable to arbitrarily complex chemistry i n principle, an additional transport equation must be solved for each species added. Conditional source-term estimation ( C S E ) [50, 58] a l lows use of the C M C hypothesis i n the L E S paradigm by solving for only unconditional means and then inverting integral equations to determine the conditional means necessary to close the source-terms. T h i s method can also be applied i n the R A N S paradigm resulting i n substantial computational savings compared to traditional C M C . However, using realistic chemistry still results i n a prohibitively large number of transport equations. L a m i n a r flamelet models [55] have a substantial advantage i n terms of computational efficiency i n that the chemistry is solved a priori and tabulated i n terms of representative scalars such as the mixture fraction used as a conditioning variable i n C M C . However, these models are, strictly speaking, valid only w i t h i n the flamelet regime and there is considerable concern regarding their applicability outside this regime. A s well, except i n limited circumstances, it is impossible to represent a flame of practical interest w i t h a single flamelet solution. A recently proposed extension to C S E , where the conditional averages are solved for i n terms of basis functions formed from varied solutions to the flamelet equations [51], has the potential to allow for the inclusion of arbitrarily complex chemistry while solving transport equations for only a small number of key scalars. T h i s proposal has been tested by comparison to results of a D N S simulation i n the a priori sense and the R A N S paradigm. CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 32 In this chapter, the method has been adapted for use i n a commercial R A N S code and used i n a predictive simulation. 2.2 Problem Description T h e arrangement considered i n this paper is a piloted methane-air jet flame. T h e geometry is that of the Sydney University piloted burner; it consists of a 7.2mm diameter fuel nozzle, surrounded by a 18.2mm diameter pilot burner. Outside the pilot burner is a uniform rising air column. T h i s geometry has been utilized to generate extensive experimental data for a variety of flames by Sandia N a t i o n a l Laboratories [47, 48]. T h e specific conditions used here are those of Sandia F l a m e ' D ' . In this flame, the fuel is 25% CH a n d 75% air, and the jet Reynolds number is 22,400. Under these conditions, there is only a small probability of local extinction occurring [47]. T h e pilot flame is a premixed C2H2 / H2 / CO2 / N2 / air flame used to stabilize the diffusion flame on the burner. 4 2.3 Background Conditional source-term estimation w i t h laminar flamelet decomposition [51] seeks to close the chemical source-term by inversion of integral equations for the conditional averages of key scalars to determine the combination of flamelet solutions necessary to construct an accurate representation of the flame. T h e conditional source-term is then obtained from the corresponding flamelet solutions, and integrated over the Favre density function (probability density function w i t h moments matching the Favre averaged mean and variance) of mixture fraction to obtain a mean source-term. These integral equations arise by representing the unconditional mean temperature i n terms of the conditional averages and Favre density function of mixture fraction: (2.1) G i v e n unique pairs of T and p(Z) for which the physical structure of the flame can be considered constant, a system of integral equations for the conditional average (T\Q can be constructed. T h e choice of the ensemble of cells over which the structure of the flame is constant i n a C F D computation is left to CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 33 the user; the specifics for this application will be discussed later. We then decompose the conditional averages into a summation of a finite number of basis functions: (T|C) « f > (T|C) , (2.2) B n=l where N need not be large if the basis functions T (£) are chosen to closely resemble the expected solution. Substituting this into Eq. 2.1 allows us to form a set of integral equations which may be solved for a : n n f>\_ (£zLi«n(TK) p (<;)<K\ n i ( 2 3 ) In searching for appropriate basis functions, an obvious choice is solutions to the unsteady laminar flamelet equations [55]: P——=u(Z,t) + —-—x{Z,t)-Qz2- (-) 2 4 These solutions produce, within the assumptions of the model, conditional averages for each scalar of interest and can be solved a priori using arbitrarily complex chemistry. If a variety of scalar dissipation histories are incorporated, the likelihood of obtaining solutions which resemble the physical flame structure is strong. The likelihood is so strong, in fact, that it is reasonable to suppose that it is possible to reconstruct the physical structure of the flame from some linear combination of basis functions obtained from a varied pool of solutions to the laminar flamelet equations. If this is in fact the case, the vector 7* could in principle be obtained by performing the inversion of the system Eq. 2.1 based on any scalar in the flow. In practice, it turns out that the ~a that results from inversion using various scalars is not necessarily unique, and the solution obtained is often highly dependent on the particular scalar used as well as the nature of the a priori information used. For example, in Figure 2.1, the temperature, carbon monoxide mass fraction, and enthalpy source term conditional averages are shown for three solutions from the library of basis functions described later. For these three solutions obtained at different times during the generation of the basis functions, the conditional averages of temperature are very similar, and might all fit the flow field equally well. In this case, it would fall to the a priori information to make the choice betuween the solutions. Obviously, as the source terms CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 34 differ dramatically, it is critical to choose the correct solution — this is why it is critical that the a priori information not be arbitrary, but rather provide realistic insight into the physical process. For the three solutions shown, the conditional averages of carbon monoxide are quite different; matching this scalar also to the flow field will dramatically reduce the reliance on the a priori information. Figure 2.1: Selected Basis Function Solutions To solve E q . 2.3, we note that the equation is a Fredholm equation of the first kind and replace the integral w i t h a numerical approximation [56]. W i t h this observation, it is straightforward to cast E q . 2.3 into a system of approximate linear equations of the form M~a « R, where the number of rows i n M need not equal the number of rows i n "a , and the problem may be ill-posed (implying possible problems related to existence, uniqueness, and continuity). T h e approximate nature of the equality adds further difficulties in that the data on the right hand side are contaminated w i t h noise due to the accumulated errors i n the scalar fields. In order to address the existence issue, the problem can be reformulated as a minimization problem, while the addition of a priori information encourages a unique solution. If we desired continuity of the solution, it could be encouraged through the a priori information; however, a discontinuous solution does not preclude obtaining a numerical solution as existence and uniqueness issues do. Solving E q . 2.3 for a , the mean chemical source-term can then be recovered i n the manner of E q . 2.1. Once a suitable choice for ~a has been made it is possible to obtain, using the same 7* i n the forward problem, the conditional averages for any arbitrary species. T h e conditional source terms for a l l species can then be n CHAPTER 2. SIMULATION OF THE SANDIA FLAME D' l 35 computed using first conditional moment closure. For efficiency, these source terms can also be determined i n advance and obtained from a table. 2.4 2.4.1 Formulation Chemical Source Term In this work, closure of the chemical source term was found by extending C S E w i t h laminar flamelet decomposition. Following the earlier method, we formed our system of integral equations relating the unconditional means to our chosen basis functions. However, we noted that i n the context of a realistic simulation covering a wide range of conditions, i n order to represent the physics present i n the library of flamelet solutions, the conditional average of temperature alone is insufficient to discriminate between flamelets. W e therefore augment our system of integral equations w i t h similar relations for the mass fraction of CO: ( f, \ (2.5) Sla (Y o\Q p{QdC YcOi V / /o«»<r|C>„p(C)dC \ n C n 1 ) For this case, where the configuration under consideration is the Sandia F l a m e ' D ' [47], we have chosen those cells which fall on planes normal to the centerline of the jet to form these ensembles and obtain two equations i n our system for each cell i n these ensembles. O u r hypothesis, supported by inspection of the library used and noting that these two scalars are sufficient to distinguish between all of the basis functions used, is that this pair is sufficient to determine ~a, and hence fully describe the structure of the system. W e employ Tikhonov Regularization [59], where the functional to be m i n imized takes the form: m i n { | | M ~ a - R\\ + A||c? - 7 ^ ° | | } (2.6) ( M * M + XI) 7 ^ = M*~R + A a ° . (2.7) 2 2 which has the least-squares solution: CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 36 In E q . 2.6, there is a balance between minimization of the error i n 'fitting' the unknown function to the data and incorporating a priori information about the function. T h e 'regularization parameter' A was chosen to a d d just enough a priori information to produce a well-behaved solution. W e explicitly supplied a guess — the solution from the previous time step — and allowed the solution to deviate from the expected value only when there was sufficient evidence to do so. T h i s choice is supported by the argument that the structure of the flame can only change at a finite rate; therefore, the flame structure at a given time step must depend on the structure at the previous step. In a flow configuration that is steady i n the mean sense such as this the choice is even more logical. A s well, beginning w i t h a physically realizable state, even i n the absence of any useful data for the inversion process a physical solution w i l l result. Other regularization possibilities (e.g., smoothing) lack this justification and may not lead to a physically realizable solution. For example, B l a i r [49] found a tendency to solutions incorporating a substantial number of negative coefficients using a smoothing operator for the a priori information, which the present choice does not encourage. 2.4.2 Flow Solution A R A N S simulation of the flame was performed using the k — e turbulence model w i t h standard model constants, except for the value of C = 1.8, i n preference to the standard value of C = 1.94 to obtain better agreement w i t h the experimentally measured mixing field near the nozzle exit. In a d dition to the Navier-Stokes and energy equations, scalar transport equations were solved for the mean mixture fraction, mixture fraction variance, and mass fractions of CO and NO. T h e 2-d axisymmetric mesh used encompassed 130 diameters i n the axial direction, and 90 diameters i n the radial direction. 120 cells i n the axial dimension, stretched w i t h a consecutive ratio of 1.04, and 78 cells i n the radial direction, stretched a variable consecutive ratio ranging from 1.05 — 1.09, were used. T h e size and stretching of the mesh was chosen to closely match that used by H i n z et al. [47], which also used a k — e model. A uniform velocity was used for the inlet boundary conditions, and a pressure outlet for the far field. T h e solution was obtained using the iterative solver i n Fluent [52]. e 2 £ 2 CHAPTER 2.5 2. SIMULATION OF THE SANDIA FLAME 'D' 37 Results and Discussion Solutions to the laminar flamelet equations with boundary conditions matching the fuel and oxidizer streams in the experiment were computed using a detailed methane-air reaction mechanism [53] at a variety of scalar dissipation rates. To build the library, the temperature and scalar fields were initialized with a mixing solution with boundary conditions at Z = 0 matching the coflow, at Z — 1 matching the fuel stream, and at Z — Z u matching the pilot stream. The mean scalar dissipation history was estimated by solving the non-reacting flow field and using the equal scales model mentioned in Section 1.2 to extract the mean scalar dissipation at the stoichiometric interface as a function of distance from the nozzle. Assuming that the stoichiometric interface was propagated downstream at the mean axial velocity of the stoichiometric interface, an approximate mean scalar dissipation history in Lagrangian time was obtained. This was then modeled using an approximate functional form for input into the flamelet evolution code. After the flamelet solutions reached steady state at the far field mean scalar dissipation, the scalar dissipation was increased linearly to generate solutions corresponding to higher scalar dissipation rates. The solution was sampled at regular intervals to build the library for use in the following R A N S simulation. However, when solutions corresponding to flamelets near extinction events were included in the library, the R A N S solution became unstable. ^ It is hypothesized that this approach can not capture strongly transient phenomena in a steady solution. For example, if extinction occurs, one would expect that re-ignition would occur a short time later, and some temporal oscillation would be observed before reaching a (Eulerian) steady solution; in an iterative solution method, this is impossible. Furthermore, solutions to the flamelet equations are incapable of capturing the appropriate re-ignition mechanism, viz, an edge flame [57]. This causes difficulty in including an extinguished flamelet solution, as there is no mechanism by which the flow can re-ignite at a later (Lagrangian) time. Instead, once extinction has occurred, it quickly propagates downstream until the entire domain is extinguished. Although this is a valid solution, it is obviously not what we hope for. To alleviate these issues, only 'burning' flamelets — those with moderate scalar dissipation rates close to the mean scalar dissipation in the flow — were included in the library. In the near field region, it is difficult to adequately describe the pilot premixed flame with decomposition, so the solution was constrained to a sinp CHAPTER 2. SIMULATION OF THE SANDIA FLAME TJ>' 38 gle arbitrary basis function near the nozzle. This basis function effectively took the place of the pilot in stabilizing the diffusion flame; although the flamelet selected did in fact stabilize the flame, the errors in the near field were substantial, and possibly influenced the downstream solution dramatically. Indeed, the solution was observed to be highly sensitive to the choice of the basis function enforced in the entrance region. Using the formulation for the flow and chemistry described above, the simulation was run until it achieved steady state in all scalars and velocity components, as well as a stationary solution for the vector ~a on all planes. A surface plot of the temperature field is shown in 2.2. A comparison of the predictions of the model to experimental data supplied in [47, 48] is given in Figure 2.3 in mixture fraction space, and in Figure 2.4 in physical space. Figure 2.2: Predicted Temperature Field The model does a reasonably good job of capturing trends in both physical and mixture fraction space despite the severe restriction of having only burning basis functions available for the decomposition. However, it is apparent that the predicted combustion occurs too close to the burner nozzle. While it is possible that this is due in part to the highly constrained solution in the near field, it is also quite likely a symptom of the lack of extinction and re-ignition in the library of basis functions. This would indicate that, to improve the simulation performance, an unsteady solution method and some means of circumventing the inability of the basic flamelet method to account for edge flame propagation is necessary. Despite the errors in the temperature and carbon monoxide fields, it is interesting to consider the predictions of the nitric oxide fields, to see if trends in pollutants can be predicted. As the temperature is over-predicted, it would be expected that the NO concentration would also be over-predicted. Figure 2.5, showing the predicted NO mass fraction along the centreline, shows that this is indeed the case. The two predictions shown correspond to the two CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 39 x/d=30 x/d=30 2500 2000 1500 1000 500 0.4 0.6 Mixture fraction 0.4 0.6 Mixture fraction x/d=45 x/d=45 2500 2000 1500 1000 0.2 0.4 0.6 Mixture fraction 0.8 Figure 2.3: C S E / L F D Validation - Comparison of Conditional Averages to Experimental Results (Solid lines predictions, symbols experimental measurements) different methods of obtaining a species field. The dashed line, referred to as 'Transport Equation' corresponds to the field obtained by solving a transport equation for the species, computing the source term at each iteration by combining the source terms from the flamelet library and accumulating the results, e.g.: P-^ = pu • VY NO Z^> = J -al - V • [pDVY ) NO (Y \c)p(QdC. NO = (^NO~ (2.8) (2.9) CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 40 The solid line, in contrast, is obtained by evaluation of: Y = NO I'a Jo • (Y \6p(0d<; NO (2.10) While it should be necessary to solve Eq. 2.8 for a scalar such as temperature which is used for the inversion process, as long as the feedback to the flow solution is minimal, as for nitric oxide, it should be valid to instead evaluate Eq. 2.10. That the two solutions differ so dramatically is cause for some concern. Specifically, it suggests that there is a mis-match between the mean diffusion of the flow and in the selected flamelets. If we were to compare the scalar dissipation in the flow, that is: Xa = 2D (VZ • VZ) (2.11) with the scalar dissipation resulting from combining the flamelets chosen: (2.12) we might expect to find that Xb would exceed Xa for this circumstance. If the results of Figure 2.5 were reversed, then we might expect to find that Xb < Xa- This is clearly a cause for concern; it is possible that the flamelet library developed did not in fact contain flamelets corresponding to sufficiently low enough dissipation rates to match the flow field. It is also possible that the library did contain the appropriate flamelets, but that they were not selected. In this case, it might seem appropriate to enforce the condition given by equating Eq. 2.11 and Eq. 2.12, and perhaps this should be done in the future. Another interpretation is also possible. Due to the different combinations of flamelets present in the field, the residence time of NO in the flamelet library at mixture fractions with negative source terms is longer for the solutions chosen than its actual residence time at these conditions in the flow. Thus, less of the NO formed is destroyed after formation in the actual flow than in the laminar flamelet solutions. By this interpretation, although the result of solving a transport equation deviated from the experiment in this case, there is not an indication that something is broken in the model that needs to be fixed by enforcing a scalar dissipation based constraint. Instead, the most likely source of the problem is errors in the solution for mixing and transport in the flow solution. CHAPTER 2.6 2. SIMULATION OF THE SANDIA FLAME 'D' 41 Conclusion A n adaptation of conditional source term estimation has been used to predict temperature and scalar fields for the Sandia Flame ' D ' . Despite an inability to incorporate extinction/reignition and forceful methods needed to address the near field, which results in an early reaction zone that experiences more complete combustion than the experimental data, the predictions of temperature and CO are promising. Further work is necessary i n order to obtain a stable solution when the effects of extinction, which can be expected to extend the flame length and produce higher levels of CO, are included. A s well, basis functions which are representative of the edge flame expected to provide (re)ignition in this flame will be necessary i n the future. W h i l e the prediction of nitric oxide concentration is poor, that trends are predicted and the result is consistent w i t h what would be expected in light of the temperature over-prediction is encouraging. CHAPTER 2. SIMULATION OF THE SANDIA CSE/LFD RANS Sandia Flame D Validation - Centeriine 2500 FLAME 'D' 42 CSE/LFD RANS Sandia Flame D Validation - Centeriine 0.08 r 2000 a % 1500 8o.04 1000 500 0 20 40 x/d 60 80 0 CSE/LFD RANS Sandia Flame D Validation - x/D=30 20 40 X/d 60 80 CSE/LFD RANS Sandia Flame D Validation - x/D=30 2500 2000 1500 1000 500 0 2500 2 4 r/d 6 8 CSE/LFD RANS Sandia Flame D Validation - x/D=45 CSE/LFD RANS Sandia Flame D Validation - x/D=45 0.04 2000, 1500 1000 500 Figure 2.4: C S E / L F D Validation - Comparison of A x i a l and R a d i a l Profiles to Experimental Results (Solid lines predictions, symbols experimental measurements) CHAPTER 2. SIMULATION OF THE SANDIA FLAME D' l 43 REFERENCES 44 REFERENCES [47] Proceedings of the Third International Workshop on Measurement and Computation of Turbulent Nonpremixed Flames, Boulder, Colorado, 1998. URL http://www.ca.sandia.gov/TNF/3rdWorkshop/TNF3.html. Section 2 - Sandia / Darmstadt Piloted CH^/Avc F l a m e D . [48] R. Barlow and J . Frank. Effects of turbulence on species mass fractions i n methane/air jet flames. In Proceedings of the 1998 27th International Symposium on Combustion, Aug 2-Aug 7 1998, volume 1, pages 10871095, Sandia N a t l L a b , Livermore, C A , U S A , 1998. Combustion Inst, Pittsburg, P A , U S A . I S B N 0082-0784. [49] C . B l a i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of B r i t i s h C o l u m b i a , 2003. [50] W . K . Bushe and H . Steiner. Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, J u l y 1999. [51] W . K . Bushe and H . Steiner. Laminar flamelet decomposition for conditional source-term estimation. Physics of Fluids, 15(6):1564-1575, June 2003. [52] Fluent 6.1. Fluent Inc., 2003. [53] D . M . G . Gregory P. Smith et al. Gri-mech 3.0. U R L h t t p : //www. me. b e r k e l e y . e d u / g r i - m e c h / . [54] A . K l i m e n k o and R. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [55] N . Peters. L a m i n a r diffusion flamelet models i n non-premixed turbulent combustion. Progress in Energy and Combustion Science, 10:319-339, 1984. [56] W . Press, S. Teukolsky, W . Vetterling, a n d B . Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, U K , 1992. REFERENCES 45 [57] G . Ruetsch, L. Vervisch, and A . L i n a n . Effects of heat release on triple flames. Physics of Fluids, 7(6): 1447-1454, June 1995. [58] H . Steiner and W . Bushe. Large eddy simulation of a turbulent reacting jet w i t h conditional source-term estimation. Physics of Fluids, 13(3): 754-769, M a r c h 2001. [59] A . Tikhonov and V . I. Arsenin. Solutions of ill-posed problems. Halsted Press, Washington, W i n s t o n , New York, 1977. Translation of M e t o d y resheniia nekorrektnykh zadach. 46 3. N O N - P R E M I X E D M E T H A N E AUTO-IGNITION CHAPTER 3.1 3. NON-PREMIXED METHANE AUTO-IGNITION 47 Introduction A predictive simulation of the autoignition process of non-premixed methane introduced into a turbulent environment was performed. Closure for the chemical source term was obtained by an extension to the conditional sourceterm estimation w i t h laminar flamelet decomposition approach proposed i n [62]. T h e ambient oxidizer conditions — the high pressure and moderate temperatures characteristic of compression ignition engines — were chosen w i t h the intent to validate the combustion model used under engine relevant conditions. T h i s was obtained by comparison of the predicted ignition delay to experimental results obtained from the shock-tube facility at U B C at several initial temperatures. E a r l y experimental measurements of the ignition delay time were obtained by Fraser et a l . [63], who injected methane into a constant volume combustion bomb following the premixed combustion of a hydrogen-ethylene-oxygennitrogen mixture to produce the necessary pressure and temperature of the oxidizer and retain an 'air-like' mixture. Subsequently, other studies have used the same facility to explore several parameters potentially affecting the ignition delay, including the ambient pressure, temperature, and composition of the fuel mixture [73]. M o r e recently, Ishiyama et a l . [67] used a similar procedure to repeat Fraser et al.'s work i n their own facility. B o t h studies used pressure measurements to identify ignition; either an increase past a threshold based on pre-ignition measurements [63, 73] or a sudden change i n the rate of pressure rise [67]. Very recent experiments still underway at the University of B r i t i s h C o l u m b i a have sought to measure the ignition delay by injection of methane into the test section of a shock tube. A description of the apparatus can be found i n [66]. T h e criterion used to define an ignition event i n the U B C study differed from the earlier work using pressure based criteria i n that the U B C study identified the onset of ignition w i t h optical observations. T h r o u g h the use of a high-speed camera, the U B C study defined the end of the ignition delay as the time at which optical emission leading to a sustained optical kernel first appeared [77]. T h e results of these independent efforts show some degree of inconsistency, especially at low oxidizer temperatures and correspondingly longer ignition delays. However, the most recent results [67, 78] agree well w i t h one another. G i v e n the dramatic differences i n the experimental apparatus and techniques used to detect ignition, it is remarkable that such agreement exists, and inspires substantial confidence i n the validity of the results obtained. CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 48 O n the simulation front as well, previous efforts have been made to develop a computational method to model ignition. B o t h Conditional Moment Closure [68] and Flamelet models [69] have been applied to the simulation of the experiment described by Fraser et al. [63]. B o t h methods produced results which agreed remarkably w i t h the experimental data presented by Fraser et al. [63]. However, the subsequent introduction into the literature of experimental data disagreeing w i t h Fraser et al.'s results by Ishiyama et al. and the anticipated publishing of supporting data from the U B C shock tube give reason to be somewhat cautious regarding these simulations. In the later discussion, some speculation w i l l be made regarding the possible explanation for the discrepancy, but further careful study b o t h experimentally and through modeling exercises is necessary. In addition to providing an assessment of the performance of the C S E - L F D combustion model, this work seeks to investigate i n some detail the process leading up to the detection of an ignition event, w i t h the goal to produce an ignition delay prediction that is directly comparable to that measured in the U B C shock tube facility by optical means. Further, as the method utilized to measure the ignition experimentally provided a spatial location for ignition, the current effort was also concerned w i t h predicting the location of ignition. 3.2 3.2.1 Formulation CSE Combustion Model Closure of the chemical source term was obtained by an extension to conditional source-term estimation [61] whereby the conditional averages of key scalars (T, CO &c CH OH mass fraction) are decomposed into a collection of solutions to the unsteady laminar flamelet equations. T h e details of the method have been described i n Chapter 1; however, it may be useful to review the general approach taken before considering the specifics of this implementation. 3 In general, conditional source-term estimation seeks to utilize the (first moment) C M C hypothesis [70] to treat the chemical source terms, but without solving directly for the conditional means. Instead, it is noted that the conditional and unconditional mean quantities can be related by virtue of CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 49 the following integral equation: Yi= /VilCMCR, (3-1) Jo where the same also holds for temperature. If the conditional average i n E q . 3.1 can be assumed constant for some ensemble of cells i n a simulation where the other quantities appearing (Yi,p(C.)) differ, a system of integral equations can be formed which can i n theory be solved for a discreet approximation to the conditional mean i n question. T h i s has been done i n the L E S paradigm and produced promising results using strongly reduced chemistry and first moment closure based on the conditional means to close the chemical source terms [76]. T h e unconditional mean source terms are then obtained by evaluating the forward problem offered by E q . 3.1. Incorporating the effects of complex chemistry necessary to treat the autoignition process [79] using the method described above would require solving a transport equation for every species involved, as well as evaluating the reaction rate expressions for each step i n the mechanism. T h i s imposes prohibitive computational costs if a realistic mechanism is used. Instead, the method proposed i n [62], whereby the conditional mean is decomposed into an ensemble of solutions to the laminar flamelet equations is utilized. In this formulation, solutions to the laminar flamelet equations are generated a priori, and E q . 3.1 becomes: _ T « f 1 / J o N Y, n(T\0 p{Z)dZ, (3.2) a n n=l where the approximation: <T|C> « f> n (T|C) n n=l has been used. T o obtain the library of basis functions (T\Qn, we solve the unsteady laminar flamelet equations i n mixture fraction space and store, for each basis function, the mean and rate of change due to chemical reaction for every scalar (mass fractions for a l l species and temperature) present at each point i n mixture fraction space. If we are w i t h i n the flamelet regime, it would be reasonable to suppose that the coefficients a which best satisfy E q . 3.2 would be consistent regardless of which scalar was considered. T h i s implicitly n CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 50 assumes that there is a unique solution to Eq. 3.2, which is not necessarily the case for all scalars. Noting the opportunity for tremendous computational savings, Bushe and Steiner [62] proposed to solve for the coefficients using only one representative scalar, such as temperature, and to then use those coefficients to reconstruct the conditional averages of all other scalars of interest and, further, the chemical source terms. If the first approximation were strictly valid, i.e., determination of a led to a unique solution for the mean mass fraction of every species at every point in mixture fraction space as well as temperature, then the closure approximation reduces to a first-moment closure in conditioning variable space. Considering a set of basis functions generated by solving the unsteady flamelet equations with varied scalar dissipation rate, it becomes apparent that there are multiple basis functions possible with very similar temperature conditional means, and that no unique solution exists which is distinct enough from the others to be distinguishable in relation to the mean field through Eq. 3.2 given the expected level of noise in the mean field. The intermediate concentrations and source terms for these solutions with similar temperatures are widely varying, as will be shown in Figures 3.1 and 3.2. W i t h this observation, it is obvious that determination of a based on temperature alone has some serious limitations. Instead, it is proposed here that transport equations be solved for, and the determination of a be based on, sufficient scalars to discriminate between all relevant flamelet solutions. In the present case temperature, and the species mass fractions of carbon monoxide (Yco) d methanol (YCH OH), were all used to form the set of integral equations for inversion. The system of integral equations written with the intention of solving for a is then written: n n n a n 3 n / \ fj Y C O ( fo«n(Y \OM)d( ^ J co Jo « n YCH 0 3 V : 2 / \ Ji^iTlO^QdC V (Y H OH\0iP(Q<K C 3 / Writing the above system of equations for each cell in an ensemble over which the conditional average is considered constant, a realistic determination of the combination of basis functions that best describes the character CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 51 of the flame can be made. In this case, it was expected that planes normal to the jet centeriine would have a similar structure, and it is o n these planes that the conditional average is assumed to be constant. A practical difficulty is encountered i n that solving this system is not trivial. Once the integrals are replaced by numerical approximations, the resulting system is extremely poorly posed. T o remedy the ill-posedness, the problem is reformulated as a minimization problem and solved i n a least squares sense, balancing the solution between 'goodness of fit' to the data appearing on the left hand side of E q . 3.3 and some a priori knowledge regarding the solution. W h i l e constraints on the first derivative of the solution for a w i t h respect to basis function order i n the library have been used i n earlier work [62], this approach has several drawbacks. T h e relationship between flamelets which reside next to each other i n the library is not necessarily unique, as it depends on the manner i n which the library was constructed. Simple sorting based on scalar dissipation rates becomes unclear when the library is built from solutions to the unsteady flamelet equations where the mean scalar dissipation rate is not varied monotonically i n time, and the time elapsed i n the calculation between the extraction of basis functions may vary. A s a matter of convenience, it is desirable that the order i n which the basis function are stored i n the library have negligible effect on the solution, which can never be the case when some functional form for a [n) is required. In the present formulation, where the system is evolving i n time i n a physical manner, the a priori information is provided through simultaneous minimization of the difference between the left and right hand sides of E q . 3.3 and the difference between the coefficient for each basis function at the current time step and the coefficient selected for the same basis function at the previous time step. T h e actual equation solved is that given i n E q . 1.37 (Chapter 1, page 22). T h e relative penalty for errors i n matching the unconditional means and changing the coefficients between time steps (A i n E q . 1.37 ) is chosen such that the solution w i l l only be pulled toward the a priori guess the m i n i m u m amount necessary to produce a well behaved solution. In addition to a more plausible physical i n terpretation, and allowing for the libraries of basis functions to be built much more arbitrarily, this method seems to have nearly eliminated the difficulties encountered by previous work [60] regarding negative coefficients. Once the appropriate coefficient vector a has been determined, the u n conditional mean of any scalar of interest present i n the library can be determined by evaluation of the forward problem given by E q . 3.2, replacing temperature w i t h the scalar of interest. Closure for the chemical source t e r m n n CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 52 is also recoverable by convolution of the appropriate conditional mean reaction rate constructed by linearly combining the source terms from each basis function weighted by the corresponding a and the pdf of the conditioning variable i n that instance. T h e reader may have noticed a switch i n E q . 3.3 to Favre-averaged quantities; this implies use of the Favre density function, p(£; Z, Z" ), i n place of a probability density function, Z, Z" ). T o recover the non-Favre averaged source terms necessary for the R A N S transport equations, the forward problem to evaluate becomes: n 2 2 (3.4) where the relation between the Favre and non-Favre pdf suggested by K l i menko and Bilger [70] has been utilized. 3.2.2 Flow Field A n axisymmetric computational domain, consisting of 55 cells i n the radial direction and 100 cells in the axial dimension, and covering physical dimensions of 0.05m x 0.15m, was used for the computation. T h e mesh was stretched w i t h a successive ratio of 1.05 and 1.01 i n the radial and axial directions, respectively, to concentrate the bulk of the cells in the region occupied by the jet prior to ignition. A structured mesh was used to facilitate the grouping of cells into ensembles on planes perpendicular to the centerline for the C S E inversion process. T h e stretching in the axial direction allowed for an especially dense grid i n the region containing the plume; this is especially important i n this type of simulation as the conditional average is inverted for independently on each plane of cells normal to the jet centerline. T h e implicit assumption is that the amount of local extinction or ignition is reasonably consistent on planes normal to the jet axis. A n overly coarse grid not only results in a less accurate solution to the scalar transport equations, but restrains the variation of the conditional averages i n physical space and makes the model more susceptible to errors resulting from spatial changes i n the flame structure. It should be noted that this grid is overwhelmingly more fine than those used i n early studies employing C S E for simulation of autoignition process [60]. In addition to properly resolving the flow, it was desired to make more data points w i t h mixture fraction between pure fuel and pure oxidizer available to the C S E inversion process. CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 53 T h e domain was initialized w i t h quiescent air at 30 bar and zero mixture fraction and variance; a pure methane jet entered from the b o t t o m left. T h e boundary condition for the methane jet was a pressure boundary condition resulting i n a moderately under-expanded jet. A l t h o u g h the actual pressure at the nozzle outlet i n the experiments was not known, the calculations given by Iaconis [66] suggested that the nozzle exit pressure be on the order of approximately 75% of the upstream pressure, and this was used as a rough guide. Some experimentation w i t h the pressure boundary i n the simulation was performed; i n general, the results are insensitive to the pressure a n d turbulence parameters at the boundary as the jet quickly expands itself and generates its own turbulence much nearer the nozzle than where significant reaction begins to occur. After 1.5ms, as i n the experiment, the inflow was terminated. T h e governing equations were discretized using a first order upwind scheme, and an implicit time step of 10/is was chosen to advance the simulation i n time. T h e chemical source term, however, was computed only at the start of iteration for every time step, and so was effectively advanced explicitly i n time. Inversion for the discrete conditional averages was performed on a 51 point grid i n mixture fraction space, w i t h the majority of the points clustered around the stoichiometric mixture fraction. 3.2.3 Basis Functions P r i o r to undertaking the flow simulation for each library of basis functions was constructed by evolving the unsteady laminar flamelet equations (described i n Section 1.2) for temperature and species mass fraction i n time using realistic chemistry developed by Huang et a l . [65]. T h i s mechanism was used without adjustment of the various constants tuned to premixed autoignition experiments. . A n internally adiabatic mixing solution between boundary conditions matching the fuel and oxidizer initial states was used as an initial condition. T h e coupled partial differential equations were then reduced to a system of ordinary differential equations using a method of lines approach ; the system of O D E ' s was solved using L S O D E [64]. E v e r y 50 [is i n flamelet time, the scalar values and their time rate of change due to chemical reaction only were output. T h e calculation was terminated at 10ms, well after ignition had occurred. T h e scalar dissipation (x) was presumed to have 1 The author acknowledges that theflameletcode modified for this purpose was provided by W.K. Bushe 1 CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 54 the functional form: (3.5) (Z,t) = xo(t)x(Z) X where the functional dependence on mixture fraction, x(Z) was taken to be the solution for a one-dimensional mixing layer offered by Peters [74]: x(Z) = e - ( 2 e r f c " 1 ( 2 Z ) ) . (3.6) 2 T h e temporal portion, Xo(^), was chosen to have the functional form: Xo(t) =c e^ +c . x (3.7) 3 M a t c h i n g to the approximated behaviour of mean scalar dissipation along the stoichiometric interface for the flow field in question yielded constants ci w 450, c = 7 x 10~ . A collection of solutions for c = 0 w i t h an oxidizer temperature of 1150k are shown in Figures 3.1 - 3.3. 6 2 3 x 10 1200 | 1100 1000 S E a 1 — — Flamelet #5 Flamelet #15 Flamelet #20 Flamelet #25 900 BOO i S. 700 E a. ro I — I 600 SOO I 400 300 0.2 0.4 0.6 Mixture Fraction 0.8 0.4 0.6 Mixture Fraction Figure 3.1: Flamelet Solutions for Temperature and Enthalpy Source Term T h e different scalar conditional means for four solutions shown, labeled as Flamelets 5, 15, 20, and 25, were all saved during the solution of the flamelet equations at the same time. These four scalars are only a part of each of the basis functions; the full basis function is comprised of conditional averages for all 41 scalars in mechanism (40 species and temperature) as well as a rate of change for each due to chemical reaction. A s this selection of solutions are all taken prior to ignition in the laminar library, there is no noticeable difference in the conditional average of temperature between the solutions. However, from the evolution of the other parts of the solution it is apparent that changes are in progress; the enthalpy source term increases as we approach CHAPTER 3. 4 c3.5 2 3 NON-PREMIXED — | — u_ 4 Flamelet #5 Flamelet #15 Flamelet #20 Flamelet #25 2 S 2 o 1 2 5 o1.5 a I AUTO-IGNITION 55 x10~ — 3.5 las §2.5 I'« METHANE 1 1 ! ! 1 Flamelet #5 Flamelet #15 Flamelet #20 Flamelet #25 1 \ 1 0.5 °0.5 0.2 0 Figure 3.2: 0.4 0.6 Mixture Fraction 0.E 0 0.4 0.6 Mixture Fraction Flamelet Solutions for Species Mass Fraction - Yco and YCH OH 3 ignition, as do the minor species. T h e two species mass fractions considered in Figure 3.2 are used to help determine the appropriate combination of basis functions. T h e conditional mass fraction of CH3OH was chosen as it evolves gradually in time prior to ignition, while CO is included as there is a sudden change in its behaviour near ignition, but before the sudden change in temperature. T h e combination provides a means of distinguishing between flamelets i n all stages of the ignition process; as well, the inclusion of CO was seen in Chapter 2 to be necessary to choose between flamelets i n the presence of extinction, so its use here foreshadows the extension of this work to incorporate the flame propagation and burning subsequent to ignition. In Figure 3.3, we see the same conditional averaged quantities shown i n the early figures, but here we are looking at a contour plot of all the solutions forming the first scalar dissipation history in the library. Ignition occurs near flamelet basis function number 50; the gradual evolution of YCH OH prior to that point is clear, whilst the scale chosen to show the changes i n temperature and Yco shortly following ignition obscures any changes before ignition. 3 In the exact flow, although the mean scalar dissipation decays exponentially w i t h distance from the nozzle exit, one can expect unresolved excursions to much higher local scalar dissipation rates. T h e effect of these can be substantial; M a s o n et al. [71], observed that during the ignition process temporary excursions in scalar dissipation, far in excess of the critical value that would normally prevent ignition, do not necessarily prevent ignition but may alter the delay, thus the history of \ is relevant. In order for the method to be able to capture these effects, it is necessary to incorporate flamelets w i t h much higher scalar dissipation rates in the library. T h i s was accom- CHAPTER 0 3. NON-PREMIXED 0.2 0.4 0.6 0.8 1 Temperature [K] METHANE 0.2 0 AUTO-IGNITION 0.4 0.6 56 0.8 Natural Log of Enthalpy Source Term [ln(erg/g/s)] x10 0 0.2 0.4 0.6 0.8 1 CO Mass Fraction 0.2 0 CH OH 3 0.4 0.6 0.8 Mass Fraction Figure 3.3: Flamelet Solutions for Species Mass Fraction - Contour Plots CHAPTER 3. NON-PREMIXED METHANE 57 AUTO-IGNITION plished by rebuilding the library of basis functions for additional values of the parameter c i n E q . 3.7. T h e final library used incorporated sub-libraries generated w i t h c = 0,0.5,10,100 concatenated as illustrated for the temperature solutions in Figure 3.4. T h e portion of the library w i t h c = 100 showed approximately a 10% increase in the magnitude of the temperature source term over the portion w i t h c = 0 for the pre-ignition solutions; the change was much more pronounced (exceeding of 100%) after ignition occurs in the flamelet calculation. 3 3 3 3 Figure 3.4: Flamelet Temperature [K] Contours for Complete L i b r a r y Here, we see a practical benefit of not enforcing dependence of a on n. A t a given time, it is conceivable that if the ideal combination of basis functions included a large amount of solution #50, we can see from Figure 3.4 that it might also include some amount of solution #125, which is very similar in terms of its conditional mean of temperature. In fact, since 75 flamelets were included using each scalar dissipation history and the solutions were extracted at the same time interval, solutions 50, 125, 200, and 275 all represent the same flamelet time, but the solutions had been subjected to differing strain rates. Insisting on smoothness i n n space would discourage this type of combination, i n contrast to the evolutionary stabilization used here which imposes no such constraint. n CHAPTER 3.3 3.3.1 3. NON-PREMIXED METHANE AUTO-IGNITION 58 Results and Discussion Predictions of Ignition Delay T h e predicted variation of the non-premixed ignition delay time w i t h temperature is shown i n Figure 3.5 along w i t h the experimental data from the U B C shock tube facility [78]. T h e results shown were produced using a n identical procedure; only the initial temperature field i n the flow solver and the initial oxidizer temperature used to build the library of flamelet solutions used was changed to match the ambient conditions. Despite the substantial scatter 1 • 0.8 • U B C Shock Tube C S E - L F D Temperature Criteria - 2.23 • -o- Delay [ms] 0.6 •o 0.4 • OT E -0.2 • ~ b c ~c • 0 • >* • a a a -0.2 -0.4 -0.6 Ipjnitk „ 0.67 , 0.75 g 1 , 0.8 1000/T[1/k] 0.85 0.9 Figure 3.5: Effect of Ambient Temperature on Ignition Delay i n the experimental data, there appears to be reasonable agreement between the simulation and the experimental results. A s a matter of interest, the fit of the simulation predictions to the experimental data was compared to the fit of a linear regression based on the available data. A s the object of this simulation exercise was to capture the effect of temperature on ignition delay, the least squares regression line for the experimental data provides an indication of the best agreement we could hope for given the scatter i n the experimental data. T h e results of the allocation of variance computations, and calculated R values are shown i n Table 3.1. W h i l e the coefficient of determination is low for b o t h the simulation and the regression line, the comparison indicates that the simulation predictions fit the experimental data nearly as well as a 2 CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 59 regression line, which is heartening for a predictive simulation not tuned to match the data. Simulation Prediction 0.288 0.52 SSE K 2 (SS yy = 0.596, fff^ Linear Regression 0.223 0.63 = 1.29 ) Table 3.1: A l l o c a t i o n of Variance for Ignition Delay In addition to predicting the time at which ignition occurs, the model used allows spatial variation along the jet axis so that it is also able to predict the location at which the identified ignition site appeared. Table 3.2 summarizes the time and location at which ignition sites are predicted. T h e laminar ignition delay, which is extracted from the library of basis functions used for each case, is included here to facilitate later discussion. T h i s laminar delay is obtained by averaging of the ignition delays for the various histories i n the appropriate library of basis functions, using the conditionally averaged temperature increase criterion. Returning to the location data, a scatter plot of the experimentally observed ignition location against the predicted ignition location is inconclusive. W h i l s t some of the points coincide either exactly or nearly so, the majority of the experimental points fall between the simulation prediction and twice the distance downstream. T h i s is not entirely unpromising; given the substantial scatter i n the experimental d a t a it is impossible to claim that the test casts significant doubt on the validity of the simulation. Additionally, interpretation of the experimental data is difficult due to many simultaneous or nearly simultaneous ignition sites. W h i l e it does appear that many of the experimentally observed ignition sites are farther downstream than the prediction, this could be explained by differences i n the jet penetration at the time of injection, arising from errors i n the momentum flux of the jet between the two due to the crude approximations made i n modeling the injector. Disappointingly, the visual methods used to observe the ignition in the shock tube do not allow for measurement of the jet penetration, so confirming or disproving this significant possibility is arduous at best, and beyond the scope of this work. In any case, it is encouraging that i n no case did the simulation suggest that ignition occurred closer t h a n CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 60 29 diameters downstream of the nozzle, and nearly all of the experimentally observed ignition sites were at least 20 diameters downstream. T[k] 1150 1200 1250 1300 1400 t ( ign) x d 1.49 1.35 1.22 1.02 0.68 0.0411 0.0373 0.0336 0.0360 0.0289 t (lam) d 2.39 1.42 0.92 0.62 0.27 Table 3.2: Location and M i x t u r e Fraction at td 3.3.2 Relation to Prior Studies To put the results of the current work into context i n the literature, Figure 3.6 compares the current results to experimental results published by Fraser et a l . [63] a n d Ishiyama et a l . [67] as well as previous simulation results. T h e earlier simulation results are those presented by K i m et a l . [68] using conditional moment closure, a n d an earlier implementation of conditional source-term estimation reported by B l a i r [60]. A n t h i r d previous simulation undertaken by K i m et al. [69] also attempted to predict Fraser et al.'s results using a flamelet model. Those results very nearly coincide w i t h those i n [68], and so are omitted from the figure for clarity. It would appear prima facie that b o t h the experimental results used for validation i n this work as well as the current simulation predictions deviate from the results presented by Fraser i n [63] (as well as the simulation results presented i n [68, 69]), although the comparison to the d a t a from Ishiyama et a l . is much more promising. There are, of course, substantial differences i n the physical configuration used to measure the delay; b o t h Fraser et al. and Ishiyama et al. were injecting into a combustion bomb as opposed to a shock tube. A s well, the injector nozzle diameter, pressure ratio, pressure of the oxidizer prior to injection, control characteristics, a n d pre-combustion composition of the oxidizer varied substantially. However, the disagreement between the two earlier experimental d a t a sets, combined w i t h the strong resemblance between Ishimyama et al.'s d a t a a n d U B C ' s shock tube d a t a suggest that these factors may not be solely responsible. T h e pressure of the oxidizer prior to injection was investigated by Fraser et al., a n d the dependence of the ignition delay on that pressure was found CHAPTER 2 3. NON-PREMIXED METHANE + Fraser et al. CSE-LFD 1.5 -o• Kim et al. • U B C Shock Tube 1 * Ishiyama et al. Blair 4.48 *% 00 1.65^ >. E I 61 AUTO-IGNITION Q 0 0.61 S-0.5 -1 + + ++ r+ -1.5 -2 0.6 + ++ c c n + 0.22 + 0.65 o *^ f 0.7 0.75 0.8 1000/T[1/k] 0.85 0.9 0.95 Figure 3.6: Effect of Ambient Temperature on Ignition Delay - Comparison P r i o r Works CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 62 to be minor, certainly insufficient to account for differences of this magnitude. In order to elucidate the reasons for the startling discrepancy, we must consider once again the question of what defines ignition. In obtaining the experimental data, Eraser et al. [63] utilized a measure of the pressure rise i n the vessel to identify ignition. T h e more recent investigation discussed by Ishiyama et al. [67] reconsidered the question of determining ignition by pressure measurement i n a constant volume chamber and reached the conclusion that, although the pressure measure used by Fraser et al. is well suited to liquid fuels, for gaseous fuels detection of ignition based on a rate of change of pressure i n the vessel is more appropriate. Before the pressure rises noticeably, a significant heat release and corresponding temperature increase must have occurred. T h i s means that the actual start of significant reaction must have occurred some time earlier. Further, the rate of heat release during this period following ignition is not necessarily constant between different oxidizer temperatures; specifically the delay between the true onset of ignition and that detected by the pressure rise will vary w i t h the amount of fuel mixed before ignition. In particular, by inspection of Figure 4 i n [63], for longer ignition delays, the method of analysis used to detect ignition from the pressure rise i n the vessel employed by Fraser et al. appears to overestimate the ignition delay by nearly 2 ms for long ignition delay times, and by perhaps 0.3 ms for short ignition delays times. T h e measure used by Ishiyama et al. is much more likely to coincide w i t h the onset of ignition determined by optical measures at U B C , as it defines ignition as occurring much earlier — essentially as soon as the combustion causes a change in pressure in the bomb, rather than after combustion has raised the pressure i n the bomb significantly. However, it is still expected to be longer than the delay detected at U B C as light emissions can be expected to be observed before the pressure responds. T h e earlier C M C and flamelet simulations are i n reasonable agreement w i t h Fraser et al.'s results, but not w i t h the more recent experiments. A s discussed above, the criteria used there for ignition is one which w i l l identify ignition after the system has been burning for some significant amount of time. A s such, these results are no more directly comparable to the current effort t h a n Fraser et al.'s data. T h a t these simulations agree reasonably well w i t h the matching data does not necessarily guarantee that they would be able to identify the onset of ignition using the more stringent criteria used i n the latter experiments. After ignition (by the new definition) has occurred, the rate of pressure rise will be dominated by the rate of heat release from CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 63 the reacting system, and may either exaggerate or mask errors i n the time at which ignition actually occurred. W h i l e it would be interesting to compare the current simulation to these previous results, it is first necessary to develop a model for the basis functions which incorporates edge flame effects [75], which are responsible for the initial flame propagation following ignition. W i t h o u t such effects, it is unrealistic to continue the current simulation beyond ignition to the time at which a criterion corresponding to the earlier work could be applied. Despite this limitation, which applied to the basis functions used by B l a i r [60] as well, B l a i r used a criteria to define ignition which would also occur much later than ignition as defined here — that the mean control volume average temperature exceeded 2000K. T h i s can only occur well after the ignition prediction based on the current criteria. T h e mechanism by which B l a i r obtained ignition delays shorter by nearly an order of magnitude based on a criteria that should have delayed the identification of ignition by a substantial margin relative to the current study, while using a method similar i n principle to that employed in this paper is unclear. There are, however, several differences between the current method and B l a i r ' s implementation which may be responsible for the discrepancy. F i r s t , B l a i r used only temperature to perform the inversion for ~a; discussion elsewhere presents the difficulty distinguishing between flamelet basis functions on temperature alone. However, B l a i r d i d not appear to suffer from a uniqueness problem problem related to using only temperature; Figure 4.3 i n [60] indicates that i n the library of basis functions used, the conditional temperature appears to be changing substantially between flamelets. However, the lack of flamelets w i t h similar temperatures highlights a larger issue: i n the library used by B l a i r , there are at best two flamelets (one of which is the initial condition) present which have a peak temperature less than 1800K. T h i s suggests that the ignition process is almost entirely unrepresented in the library of basis functions. A s well, the lack of 'ignited' solutions w i t h similar temperatures is the result of using only one scalar dissipation history, which was increasing i n time as opposed to the expected decrease i n time i n the problem under consideration, to b u i l d the library. T h i s suggests that it is quite possible that if an additional scalar were included, the minimization of the least squares problem necessary to determine ~a may have had a very large residual indeed w i t h this library, as it is unlikely that any combination of the solutions present accurately matches the physics in the flow. A s a result, the physically unjustifiable smoothing imposed on ct (n) to stabilize the inversion process would likely have h a d a n CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 64 significant impact on the chosen solution. A s any of these major issues could potentially explain the differences between the current approach and B l a i r ' s implementation, it is difficult to ascertain to what extent each is responsible. T h e effect of other potentially less substantial differences, such as a much finer grid i n the current study to allow for more rapid variation of the flame structure i n space a n d more data points for the inversion process, may be entirely masked, making further comparison unproductive. 3.3.3 Behaviour Prior to Ignition and Identification of Ignition One of the advantages of the closure chosen is that despite solving for a n (unconditionally) ensemble averaged solution, we can obtain the conditional averages of any scalar up to the validity of the closure. A key feature of the ignition process is that some finite amount of time elapses without a noticeable increase i n temperature and a slow build up i n reaction rate, followed by a rapid increase i n temperature and reaction rate [72]. W h i l e early implementations of C S E + L F D i n K I V A - 3 V produced a steady increase i n temperature from the start of the simulation [60], the flavour presented i n this paper follows the physical behaviour much more closely. T h e conditional averages of the three scalars used for the inversion process are shown i n F i g ure 3.7 along w i t h the conditional average of the creation rate of the CH radical at several times prior to ignition. For clarity, only the conditional average from the plane showing the largest increase over the initial condition is shown. Here we see that there is relatively minor change i n the conditional average of temperature prior to ignition, followed by the development of a local increase i n temperature near the stoichiometric mixture fraction. Following ignition, the increase i n temperature spreads outward from the stoichiometric mixture fraction. T h e lack of change i n the conditional average for temperature prior to ignition suggests, once again, that this is a poor scalar to base a determination of the state of the system on prior to ignition. T h e concentration of carbon monoxide also increases dramatically, but there is a more apparent increase i n the concentration prior to ignition. T h e evolution of the conditional average of the intermediate CH OH is more interesting. T h e peak conditional magnitude of this species is substantial well before ignition takes place, and changes location i n mixture fraction space as ignition 3 CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 65 Figure 3.7: Evolution of Conditional Averages in Mixture Fraction Space, T ox = 1150 A" CHAPTER 3. NON-PREMIXED METHANE 66 AUTO-IGNITION occurs. C o m p a r i n g to Figure 3.3, we see that there is a dramatic shift to richer mixture fractions for the residence of this species near ignition. T h i s shift, as well as the non-negligible concentration of the scalar well prior to ignition, justifies the choice of this scalar to aid i n selecting the appropriate combination of flamelets prior to ignition. Finally, the behaviour of the CH creation rate is shown. T h i s quantity, which is obtained from reactions forming CH only (as opposed to the net creation rate), was investigated i n order to explore a quantity directly related to the optical emissions recorded i n the U B C shock tube experiments. In the absence of soot, the hypothesis is that the light emissions observed are primarily the result of CH radicals decaying from an excited state, and that the quantity of excited molecules is proportional to the creation rate. T h e behaviour is promising for the detection of ignition; as w i t h temperature, there is a very dramatic increase just as ignition occurs. Setting the absolute level at which to declare ignition has occurred is not obvious; i n the U B C shock tube, the threshold was determined by the sensitivity of the camera used, but the relationship between the CH creation rate and the luminosity seen by the camera is unclear. In Figure 3.8, the m a x i m u m increase i n the conditional average of the various scalars at any mixture fraction and on any plane is plotted against time for a selection of oxidizer temperatures. T h e various quantities are normalized by a threshold chosen for the quantity, such that ignition was declared based o n the various quantities when the quantity plotted exceeds a value of one. A l t h o u g h the choice of where to set the thresholds was somewhat arbitrary, they were selected to be such that they were the lowest value that could be selected to not be confused by fluctuations prior to the dramatic increase, and to occur very shortly after the visually apparent 'knee' i n the various plots. W i t h the thresholds chosen, there is very little difference i n the ignition delay identified by consideration of any of the three scalars, as is apparent from Figure 3.9. G i v e n this comparison, the ignition delay identified by a n increase of 75 degrees i n the conditional average of temperature at any mixture fraction, anywhere i n the flow, was deemed to be well correlated w i t h the optical criteria used i n the experiments. T h e data presented earlier (Figures 3.6, 3.5) was generated using this criterion. In addition to consistency i n the predicted ignition delay, the location along the centreline at which ignition occurred was within one plane regardless if the temperature or CH based criteria was used. T h e CH detection, as expected, displayed its maximum increase at slightly richer mixture fractions (Z\{T\c,) 0.083) than the temperature criterion (Z\(T\Q — w max i n c max i n c CHAPTER 3. NON-PREMIXED g>10 AUTO-IGNITION — Temperature CH Production Rate CO Source Term j S 6 J E 1• a 2 1 METHANE 2 1 l[ms] T O I T m 67 — Temperature CH Production Rate — CO Source Term /II /' e E ill 'S 2 1 jjy D = 1150k • ! 2 0.5 T ox = 1300k 1 t[ms] 1.5 2 = 1200k 1400k Figure 3.8: Evolution of Peak Increase i n Conditional Averages 0.5 r E %. 0 o - * Temperature Increase C H Production Rate C O Source Term -°0.7 0.75 0.8 1000/T[1/k] 0.85 0.9 Figure 3.9: Comparison of Results for Different Ignition C r i t e r i a CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 68 0.064). 3.3.4 Effect of Turbulence on Ignition Delay- In considering the autoignition of a non-premixed system, the effect of turbulence o n the delay before ignition is not immediately obvious. There are competing effects; increased turbulence levels will increase the rate of mixing. W h i l e this means that the fuel and oxidizer will m i x sooner, it also means that any radicals produced or temperature increase will be dissipated quickly. Further, i n a turbulent flow, there will be a distribution of strain rates; w i t h i n any given region, there will be a variety of strain histories present [72]. T h e laminar and turbulent ignition delays from Table 3.2 are plotted against each other i n Figure 3.10. W e see that at high temperatures, where the delay is short, the effect of turbulence is to increase the ignition delay. A t low temperatures, where the delay is longer, the effect of turbulence is to decrease the delay. Interestingly, the ratio of turbulent to laminar ignition delays is almost linear w i t h 1000/T; a linear regression gives a coefficient of determination of R = 0.98. T h i s effect can be related to the earlier framework by the following speculation: when the ignition delay is short, the effect of turbulence dissipating temperature and radical build up is dominant; although the presence of turbulence also suggests that there will be regions of low strain, the likelihood of one of these regions having been already sufficiently mixed for ignition to occur is low. A t longer ignition delays, the majority of the field is well mixed prior to ignition. T h e n , the distribution of strain rates suggests that there is a strong likelihood that one of the well mixed regions w i l l encounter low strain, resulting i n ignition promotion. 2 3.4 Conclusion A simulation has been undertaken to predict the ignition delay of nonpremixed methane injected into hot air under high pressure at high velocity. T h e technique used to close the source term i n the R A N S paradigm takes advantage of concepts from b o t h the C M C and flamelet models. A separate conditional average is permitted (and encouraged) on each plane of cells axially equidistant from the nozzle. T h i s allows for a separate combination of the basis functions on each plane; as the distance from the nozzle increases, the mean strain rate can be expected to decay, so a different combination CHAPTER 3. NON-PREMIXED 0.75 METHANE 0.8 1000/T[1/k] AUTO-IGNITION 0.85 0.9 Figure 3.10: Ratio Between Turbulent and Laminar Ignition Delays 69 CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 70 of basis functions w i l l be necessary. T h e predicted ignition delay was i n excellent agreement w i t h the experimental data at a variety of oxidizer temperatures. Despite the significant scatter in the experimental d a t a used for direct validation, the agreement i n trends w i t h an additional experimental d a t a set obtained from an independent source lends to significant confidence i n the results. Scatter i n experimental measurements of the spatial location of the ignition event, as well as the occurrence of multiple, nearly simultaneous ignition events in the experiment, prevented a meaningful validation of the predicted ignition location. However, the limited comparison that was possible did not present any evidence that we should be skeptical of the simulation results. In order to produce simulations valid beyond ignition, it is necessary to develop libraries of basis functions which encompass the appropriate physics to address the flame propagation following ignition, which are not captured in the solutions to the laminar flamelet equations. REFERENCES 71 REFERENCES [60] C . B l a i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of B r i t i s h C o l u m b i a , 2003. [61] W . K . Bushe and H . Steiner. Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, J u l y 1999. [62] W . K . Bushe and H . Steiner. Laminar flamelet decomposition for conditional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [63] R. A . Fraser, D . L. Siebers, and C . F . Edwards. A u t o i g n i t i o n of methane and natural gas i n a simulated diesel environment. SAE Technical Paper Series, 1991. [64] A . C . Hindmarsh. O D E P A C K , A systematized collection of ode solvers. In R. Stepleman et al., editors, Scientific Computing, pages 55-64, A m sterdam, 1983. North-Holland. [65] J . Huang, P. H i l l , W . Bushe, a n d S. M u n s h i . Shock-tube study of methane ignition under engine-relevant conditions: experiments a n d modeling. Combustion and Flame, 136(1-2):25-42, 2004. [66] J . Iaconis. A n investigation of methane autoignition behaviour under diesel engine-relevant conditions. Master's thesis, University of B r i t i s h C o l u m b i a , 2003. [67] T . Ishiyama, M . Shioji, et a l . Characteristics of spontaneous ignition and combustion i n unsteady high-speed gaseous fuel jets. In 2003 JSAE/SAE International Spring Fuels & Libricants Meeting, number J S A E 20030225 / S A E 2003-01-1922, Yokohama, Japan, M a y 19-22 2003. [68] S. H . K i m , K . Y . H u h , and R. A . Fraser. Modeling autoignition of a turbulent methane jet by the conditional moment closure model. Proceedings of the Combustion Institute, 28:185-191, 2000. REFERENCES 72 [69] S . - K . K i m , Y . Y u , et al. Numerical investigation of the autoignition of turbulent gaseous jets i n a high-pressure environment using the multiplerif model. Fuel, 83:375-386, 2004. [70] A . K l i m e n k o and R. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [71] S. D . M a s o n , J . H . C h e n , and H . G . Im. Effects of unsteady scalar dissipation rate on ignition of non-premixed hydrogen/air mixtures i n counterflow. Procedings of the Combustion Institute, 29:1629-1636, 2003. [72] E . Mastorakos, T . B a r i t a u d , and T . Poinsot. Numerical simulations of autoignition i n turbulent mixing flows. Combustion and Flame, 109: 198-223, 1997. [73] J . D. Naber, D. L. Siebers, S. S. Dijulio, and C . K . Westbrook. Effects of natural-gas composition on ignition delay under diesel conditions. Combustion and Flame, 99:192-200, 1994. [74] N . Peters. L a m i n a r diffusion flamelet models i n non-premixed turbulent combustion. Progress in Energy and Combustion Science, 10:319-339, 1984. [75] G . Ruetsch, L. Vervisch, and A . L i n a n . Effects of heat release on triple flames. Physics of Fluids, 7(6): 1447-1454, June 1995. [76] H . Steiner and W . Bushe. Large eddy simulation of a turbulent reacting jet w i t h conditional source-term estimation. Physics of Fluids, 13(3): 754-769, M a r c h 2001. [77] G . Sullivan. Personal communication. 2004. [78] G . Sullivan. Personal communication. M a y 2004. [79] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress in Energy and Combustion Science, 28:193-266, 2002. 4. DISCUSSION AND CONCLUSION CHAPTER 4. DISCUSSION AND CONCLUSION 74 In the preceding chapters, two implementations of a novel method of capturing the physics of non-premixed turbulent combustion proposed by Bushe and Steiner [82] have been discussed. Chapter 2 presented the first application of the method to a reacting flow w i t h a steady ensemble averaged solution. W h i l e not entirely successful, that a converged solution was obtained sets the framework for future development. T h e second paper, presented i n Chapter 3, discusses a successful simulation of non-premixed methane autoignition, and demonstrates that where the basic requirements of the model are satisfied its performance can be excellent. W h e n accounting for the effects of chemical reaction w i t h conditional source-term estimation using laminar flamelet decomposition, the effects of arbitrarily complex chemistry can be accounted for w i t h m i n i m a l computational cost. T h i s is because the majority of the computational effort associated w i t h the chemistry is encountered during the library generation, which can be performed independently from the flow simulation. Further, as the flamelet solutions for the varied strain histories necessary for a sufficiently general library are entirely independent, the calculations for library generation can be easily performed in parallel. W h e n the actual simulation of the turbulent flow is performed, only the smallest timescale on which the scalars used in the inversion process evolve must be resolved, as transport equations are solved for only these species. Events which prove difficult to capture w i t h conditional moment closure using a single conditioning variable, such as extinction and ignition, can be captured i n an unsteady simulation, a l though there are yet to be resolved complications associated w i t h capturing unsteady effects in a steady solver. In general, the model has broad applicability and should be useful i n any circumstance where the necessary physics can be captured i n the library of basis functions. In early formulations of conditional source-term estimation where the conditional averages were inverted for directly and conditional reaction rates computed from the first conditional moment, it was necessary to assume some degree of spatial homogeneity of the conditional averages, and that the fluctuations about the conditional means were negligible [81]. In this formulation using decomposition, we must tighten the assumption of spatial homogeneity but are able to relax the assumption regarding the fluctuations. W h e n closing the source term directly from conditional means determined by inversion, it is necessary to assume spatial homogeneity only across cells where there is an overlap of the pdf of the conditioning variable. Here, it is necessary to assume spatial homogeneity of b o t h the conditional averages CHAPTER 4. DISCUSSION AND CONCLUSION 75 (even i n the absence of overlap of the pdf of mixture fraction) and of the flame structure (as captured by the vector of coefficients 7?) over the ensemble of cells used for inversion. A s well, it is necessary that any fluctuations about the conditional mean modify the conditional mean of the scalars used for inversion sufficiently to make the change detectable. It is not, however, essential that the fluctuations be negligible — i n fact, it is expected that they will not be. To see why this is so, and how the method is able to escape without the use of a second condition, it is necessary for a moment to consider how the model functions in a hypothetical circumstance. In treating variation around the mean conditioned on a single variable, it is noteworthy that there are two types of variation about the conditional averages to be concerned with. T h e first is the approximately Gaussian variation about the conditional mean w i t h small standard deviation due to small local fluctuations i n the strain and other effects not related to mixture fraction. T h e second is the variation observed when the local fluctuations i n the strain become large enough to cause major phenomenological changes in the flame structure, or external effects such as a passing edge flame modify the structure of the reaction zone. A l l first moment closures must assume that the former variation is small; second moment closures have been developed to address larger variations [83, 84]. T h e latter form of variation is more problematic, and causes difficulty for second moment C M C methods as well. To see why the latter form is so detrimental to first moment closures, suppose that w i t h i n a given region, conditionally averaging individual points obtained from a fully resolved sample produces a conditional mean of temperature coinciding w i t h the dashed line on the left of Figure 4.1. Further suppose that approximately half the points are clustered about the b o t t o m curve (solution A ) , and the remainder are clustered near the top curve (solution C ) . T h i s would represent a situation where a significant amount of local extinction was present. In this case, the conditional mean for the mass fraction of carbon monoxide could be expected to be the dashed curve on the right of Figure 4.1, formed again from approximately half the sample points clustered about each of solutions A and C . T h e solid curves i n the figures represent solutions to the flamelet equations for various levels of scalar dissipation — 'Solution C is a solution for low x, and 'Solution A ' is a solution for high x, near extinction. 'Solution B ' is a flamelet solution for intermediate x- T h e dashed lines are a conditional average of Solutions A and C . A l t h o u g h it necessitates making a flamelet assumption to relate the two scalars, it is certainly possible to conceive of a turbulent flame where this CHAPTER 4. DISCUSSION AND CONCLUSION 76 — — Solution A Solution 6 Solution C <A,C|(> §0.02 0.2 0.4 0.6 Mixture Fraction 0.2 0.4 0.6 Mixture Fraction Figure 4.1: Flamelet Solutions for T and Yco', {A,C\Q average of Solutions A and C o.e indicates conditional situation might exist — e.g., one which falls i n the flamelet regime, and i n which approximately half of the realizations are close to extinction and half are not. A first order C M C method would compute the conditional reaction rates based on the conditional means indicated by the dashed lines, whilst a strictly valid closure would compute the reaction rate for each T, Yco pair of points and then conditionally average the result. For this situation, the first order C M C closure, as well as second order closures which presume a Gaussian distribution about the conditional mean would be highly inaccurate. In C S E + L F D , faced w i t h this situation, we have to choose the coefficients aA,<*B,ctc- Choosing = 0 , Q B = 1.0, ac = 0 would agree well w i t h the temperature data. If the inversion was performed based on temperature alone, as proposed i n [82] and implemented i n [80], this would be the expected result. T h e conditional source term used would then be that calculated based on {T\Q , (Yi\Q . T h i s would obviously be error prone — worse i n fact than the first order C M C closure. However, using carbon monoxide also for inversion, the method would find a poor fit to the data for Yco- If the scalars chosen are sufficient to discriminate effectively between the basis solutions offered, the only way to obtain simultaneous agreement w i t h both the temperature and carbon monoxide fields is to choose = 0.5, O B = 0.0, ac = 0.5. T h e n , the source term used will be an average of the source term calculated based o n (T\Q , ( ^ I O A and (T\Q , ( F j | 0 , which is much closer to the ideal behaviour. W i t h this type of closure, the use of a second conditioning variable is unnecessary; the physics which would normally be captured by using a second conditioning variable are implicitly captured by the combination of flamelets chosen to match the behaviour of the inversion B B A C C CHAPTER 4. DISCUSSION 77 AND CONCLUSION scalars. However, the conditional average we are seeking must be constant across the cells used for inversion when the kernel varies, so that each cell w i t h a different kernel provides information about a different aspect of the behaviour of the same underlying function. The above discussion and claims are predicated on the presence of solutions representing all possible physical conditions i n the library used for decomposition, and on choosing the appropriate combination of basis functions based on the mean field supplied being possible. In order to discriminate between the basis functions, the ensemble of cells used for inversion must present data (unconditional means) which show the effect of filtering the u n derlying function in sufficiently varied ways that we choose between basis functions appropriately to reflect the behaviour of the conditional mean. Distinguishing between basis functions appropriately also entails the use of scalars for the inversion process which effectively indicate different aspects of the solution's behaviour. Rather than think of the function that we are inverting for as a one-dimensional function of mixture fraction, we are actually searching to minimize the sum of the components of a vector function, / , where each component of / represents the error in the prediction for a given scalar. T h i s statement is similar to E q . 1.38, and can be written as: m i n { | | / i | | + ||/ || + ||/s|| + --- + ||// i||} 2 where (4.1) + are the individual components of / . In the continuous exten- sion of a discreet basis function space (n*), each component of / can be written as a product of the appropriate coefficient, and the corresponding scalar value for that position i n basis function space; we are then t r y i n g to minimize a functional of the coefficient function and the scalar functions: m i n {||-?||} ; ^ = ( a ( n ) ( T | 0 (n), a ( n ) (Yi|C>, • • •, «(") <Xi\0 ) (4-2) Considering E q . 4.2, it is apparent that the solution of E q . 4.1 for a(n*) is actually a functional dependent on (T\() (ri*), (Yi|C) (n*), • • •, (Yi\0 (ri*). T h i s is clearly a difficult situation; it might at first appear that inversion i n this pan-dimensional space is hopeless. Fortunately, we can make assumptions that simplify the problem. F i r s t , for a practical approximation the basis function space ri* is not infinite. O n l y the physics relevant to the problem at hand must be included; for example, i n unbounded flows wall effects need not be included. O u r a priori knowledge of the important physics CHAPTER 4. DISCUSSION AND CONCLUSION 78 and modeling assumptions can reduce the domain of n* to a reasonable size. Further, we seek only a discreet approximation i n n*. Lastly, by inspection of the solutions generated for the discreet subset of n*, the functional is frequently homogeneous i n many of the dimensions, and there are only a limited number of approximately orthogonal coordinates i n the space. These important coordinates which must be retained as scalars used for inversion must be chosen to reflect the physics of interest in each circumstance, but doing so is relatively straightforward once the library of basis functions has been assembled. T h e presence of all relevant physics in the library of basis functions is of paramount importance, as is the use of sufficient scalars for inversion to discriminate between those physics. In the first paper presented, the poor performance of the model is due in large part to the inability of the model to account for the premixed pilot flame. T h i s is a classic example, where the physics necessary to describe this process were not present i n the library. Other physics of relevance to the simulation of the Sandia F l a m e ' D ' , as well as continuing the autoignition simulations beyond ignition include edge flames [85], which can never be captured i n a library built from simple solutions to the laminar flamelet equations. T h i s highlights that although a reasonable place to start, the use of laminar flamelets as basis functions has little future. Heat transfer between the flow and external objects is also of relevance in many practical flows, as are buoyancy effects — these also w i l l necessitate more comprehensive solutions for use as basis functions. Fortunately, the formulation is independent of the method used to construct the basis functions — they could be as well built from solutions constructed from C M C , or D N S for that matter. In general, C S E w i t h decomposition can be viewed as a way to combine easily obtainable elementary solutions to b u i l d a more realistic simulation. A s more complex physics are added, it is likely that the number of scalars necessary for inversion may increase substantially to discriminate between the additional physics, as it is necessary to use a scalar for inversion that can discriminate between basis functions across all phenomenological transitions. Once the appropriate combination of basis functions has been chosen, the mean reaction rate is obtained by virtue of an expression similar to E q . 1.29; that is: (4.3) 4. DISCUSSION CHAPTER AND CONCLUSION 79 Since the p d f of mixture fraction is approximated throughout this work by a presumed functional form, it is of course necessary to assume such a form is physically correct. A s mentioned i n [82], however, this is a common assumption not specific to the C S E methods, and it would be trivial to adjust the formulation to utilize more complex representations of the pdf. C o m p a r i n g E q . 4.3 (written for temperature) w i t h the corresponding closure using double-condition C M C (using scalar dissipation, w i t h sample space X as the second conditioning variable): 1 w « / I < w | C , X > ( C , x W x ' « f f (viCx'jpiOpMdCdx' (4.4) Jc J ' Jc J > x X Bushe and Steiner [82] noted that coefficients a take the place of p(x')dx' when a numerical quadrature is used to perform the integration. In the present formulation, it is clear that defining a specific value of x' f ° the basis functions used, where the history of x as well as the instantaneous value of x when the solution is extracted from the time dependent solution is relevant, is difficult. In general, as the basis functions could vary independently of x or Z, the quantity a is more appropriately likened to p(a')da', where a is a n arbitrary scalar which must satisfy only the condition of being independent of Z, so that it captures variations i n the flame structure not correlated w i t h mixture fraction. T h i s poses a minor conceptual difficulty, as it implies that we are i n fact performing a numerical integration over a', the sample space for a. A requirement, then, is that the library must be sufficiently fine grained such that a numerical integration over a' is not error-prone. T h i s could be a problem where the flame structure varies rapidly w i t h n , e.g., for basis functions near extinction. For this phenomenon, the scalar a would be strongly related to x\ at high x the log-normal p d f of x approximately linear; however, the reaction rate varies i n a non-linear way and the product (u\x')p(x') the integral of E q . 4.3 must be highly non-linear, so the constraint may be arduous. In Chapter 2, it was mentioned that when flamelets w i t h high scalar dissipation were included i n the library the solution became unstable, and it was suggested that this was due to the transient nature of these flamelets being included i n a steady solution. Another possibility is now apparent: the errors arising from the effective discretization i n scalar dissipation sample space may be large. A s the flamelets used for basis functions were extracted at constant time intervals, and the structure was changing more rapidly i n flamelet time near extinction, the effective discretization became much coarser when the n r n 1S s e e n m CHAPTER 4. DISCUSSION AND CONCLUSION 80 extinguishing flamelets were included. A l t h o u g h not definitive, this line of reasoning suggests that it may be worthwhile to explore the use of a much finer library, including the extinction process. T h e drawback, of course, is that the library becomes much larger, the inversion process more timeconsuming, and the uniqueness of the solution suffers. G i v e n that i n the flamelets used for the decomposition i n this work there is a strong correlation between the scalar dissipation and the structure of the flamelet solutions, and it is possible to solve a modeled transport equation for the mean scalar dissipation, it might seem that we should add x f ° the list of scalars used for inversion, and solve: X« f p ( C ) [ « • <xK>] ~ f p(C) • XofiO] d{ (4.5) However, when extinction is included, it is also necessary to include an extinguished flamelet, i.e., a pure mixing solution. For this solution, the flamelet's behaviour i n terms of scalar concentrations and source terms is independent of the scalar dissipation, so the value of (xo) for that basis function is undefined — it could take any value. If we let the extinguished flamelet be the last flamelet in the library, it then makes more sense to impose the constraint: TAf-i X> /%(C) Jo E («n(Xo)n) /(CR (4-6) _n=l It is not yet clear if this additional constraint would produce an improvement — hopefully, given the assumptions of the model discussed so far, it should be implicitly satisfied based on the scalars already employed for the inversion. It may be useful as a check to evaluate E q . 4.6 and test it's validity; if it is violated, then likely there is a violation of the earlier assumptions. Finally, some comments are warranted regarding the choice of a priori information used to stabilize the inversion process. In choosing to guide the solution chosen by the inversion process toward the solution from the previous inversion, we are making a physical assumption about the expected behaviour. A l t h o u g h i n general the penalty for deviating from the previous solution is set only just high enough to ensure a unique solution, the penalty exists nevertheless, so the solution will always incorporate the a priori assumption to some extent. In searching for a solution which is steady in the Eulerian sense, our choice seems logical. It will provide stability for the inversion as the solution evolves, limit the rate at which it can change CHAPTER 4. DISCUSSION AND CONCLUSION 81 (effectively under-relaxing the solution), and have no effect once a steady solution is obtained. In an unsteady calculation, it might be suggested that a more physical choice would the the solution from an earlier time i n the Lagrangian sense. T h i s claim would treat the flame structure as a convected entity, which is not entirely unreasonable. T h e difficulty is found i n that it would also insist that the flame structure across all cells over which that structure is considered constant must somehow be transported as a group, which is implausible. Interestingly, introducing a bias toward a past solution i n the Lagrangian sense would have to produce a solution equal to the past solution i n the Eulerian sense for a problem w i t h an E u l e r i a n steady solution. A n elegant method of incorporating influence from the previous structure i n the Lagrangian sense was not explored i n this work, but may be an interesting avenue for future development. A revised formulation of conditional source-term estimation w i t h l a m i nar flamelet decomposition has been presented and used successfully i n the simulation of turbulent non-premixed methane ignition. T h e framework has also been used to develop a simulation of the Sandia F l a m e ' D ' . T h i s formulation differs from previous implementations i n several ways. Transport equations are solved for multiple species which are used i n the inversion process to determine the flamelet structure; using multiple scalars allows the method to discriminate between the various physical conditions represented i n the library of basis functions. T h e library of basis functions used has been extended to incorporate more varied solutions, i n order to more fully represent the relevant physics. Extending the library of basis functions was made possible by the use of a quasi-physical method of introducing a priori information into the inversion process, which forces the solution to evolve between successive iterations only so far as supported by changes i n the flow field. A slightly revised interpretation for the method was also suggested. T h e successful prediction of the temperature dependence of the ignition delay for an igniting methane jet, as well as obtaining a steady solution for the Sandia Flame ' D ' i n the R A N S paradigm effectively demonstrates the potential of the model, and provides a framework for future development. REFERENCES 82 REFERENCES [80] C . B l a i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of B r i t i s h C o l u m b i a , 2003. [81] W . K . Bushe and H . Steiner. Conditional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, J u l y 1999. [82] W . K . Bushe and H . Steiner. L a m i n a r flamelet decomposition for conditional source-term estimation. Physics of Fluids, 15(6):1564-1575, June 2003. [83] S. H . K i m , K . Y . H u h , and R . W . Bilger. Second-order conditional moment closure modeling of local extinction and reignition i n turbulent non-premixed hydrocarbon flames. Proceedings of the Combustion Institute, 29:2131-2137, 2002. [84] E . Mastorakos and R. Bilger. Second order conditional moment closure for the autoignition of turbulent flows. Physics of Fluids, 10(6), June 1998. Letter. [85] G . Ruetsch, L. Vervisch, and A . L i n a n . Effects of heat release on triple flames. Physics of Fluids, 7(6): 1447-1454, June 1995. RANS FIELDS APPENDIX A. RANS FIELDS 84 In Figure A . l , the R A N S fields showing contours of the unconditional mean temperature and CH creation rate are given at the time of ignition. It is important to note that although there is a postulated correlation between the reaction rate of the CH radical and light emissions, the specifics of the relation are not known in sufficient detail to relate the figure to images captured experimentally. Further, these contours are ensemble averages, which need not necessarily have any direct relation to any specific realization. R A N S Temperature Field, 1150K R A N S C H Creation Field, 1150K Oxidizer Oxidizer 0.015 _ 0.01 *0.005 0 0.04 400 500 600 700 800 900 0.05 1000 1100 R A N S Temperature Field. 12O0K Oxidizer 500 600 700 8O0 900 1000 05 1 1.5 2 25 3 3JS 4 R A N S C H Creation Field, 1200K Oxidizer 0.04 400 0 0.05 1100 1200 R A N S Temperature FieW, 1300K Oxidizer 400 500 600 700 800 900 1000 1100 1200 1300 400 600 1000 12O0 2 3 4 5 R A N S C H Creation Field, 1400K Oxidizer R A N S Temperature Field, 1400K Oxidizer 1400 Figure A . l : R A N S Temperature and C H Creation Contours at td 4.5
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Combustion modeling with conditional source-term estimation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Combustion modeling with conditional source-term estimation and laminar flamelet decomposition Grout, Ray W. S. 2004
pdf
Page Metadata
Item Metadata
Title | Combustion modeling with conditional source-term estimation and laminar flamelet decomposition |
Creator |
Grout, Ray W. S. |
Date Issued | 2004 |
Description | Conditional source-term estimation with laminar flamelet decomposition has been utilized to model the mean chemical source term in a predictive RANS simulation of two different problems. With this model, integral equations relating the unconditional mean temperature and species fields to the conditional means are used to determine the flame structure as a combination of basis functions formed from solutions to the unsteady laminar flamelet equations. In the simulation of the Sandia Flame 'D', a well studied co-flowing piloted jet flame with a steady average solution, a converged solution was obtained which captured the trends in the temperature and major species, although the nitric oxide prediction overestimated the peak concentration by a substantial margin. Simulation of the autoignition process of non-premixed methane produced simulation predictions in excellent agreement with the experimental data. Using a new, more stringent, criteria to define ignition than earlier studies, the effect of ambient temperature on the ignition delay was captured, as was the expected physical behaviour prior to ignition. The errors in the simulation of the co-flowing piloted jet can be attributed to a large degree to the lack of ability of the basis functions used to account for the effect of the pilot flame. That the autoignition simulation was much more successful highlights the importance of including all relevant physics in the library of basis functions. Earlier formulations of the model were extended by expanding the sample space for the basis functions to include a wide variety of solutions to the laminar flamelet equations; in order to distinguish between the larger set, a new method of stabilizing the numerics was found to be necessary and the number of scalars used to determine the flame structure was increased. The formulation is easily extensible to libraries of basis functions which are capable of including effects such as pilot and edge flames not captured by the laminar flamelet solutions. |
Extent | 8842098 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-11-24 |
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.0080709 |
URI | http://hdl.handle.net/2429/15580 |
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-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2004-0470.pdf [ 8.43MB ]
- Metadata
- JSON: 831-1.0080709.json
- JSON-LD: 831-1.0080709-ld.json
- RDF/XML (Pretty): 831-1.0080709-rdf.xml
- RDF/JSON: 831-1.0080709-rdf.json
- Turtle: 831-1.0080709-turtle.txt
- N-Triples: 831-1.0080709-rdf-ntriples.txt
- Original Record: 831-1.0080709-source.json
- Full Text
- 831-1.0080709-fulltext.txt
- Citation
- 831-1.0080709.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-0080709/manifest