Combustion Modeling with Conditional Source-Term Estimation and Laminar Flamelet Decomposition by Ray W.S. Grout B.A.Sc, The University of British Columbia, 2002 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Department of Mechanical Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 8, 2004 © Ray W.S. Grout, 2004 THE UNIVERSITY OF BRITISH COLUMBIA FACULTY OF GRADUATE 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 print) <^3/o Sr/2JZX=><J-. / _( ' Date (dd/mm/yyyy Title of Thesis: Degree: Year: Department of The University of British Columbia Vancouver, BC Canada grad.ubc.ca/forms/?formlD=THS page 1 of 1 last updated: 20-Jul-04 11 ABSTRACT 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 condi tional means are used to determine the flame structure as a combination of basis functions formed from solutions to the unsteady laminar flamelet equa tions. 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 ex perimental 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 er rors 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 va riety 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 struc ture 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. iii CONTENTS Abstract ii Contents iiList of Tables v List of Figures vi Nomenclature viiAcknowledgements x Introduction xi 1 Background 1 1.1 Turbulent Combustion 2 1.2 Flamelet Models 3 1.3 Conditional Moment Closure 12 1.4 Conditional Source-Term Estimation 17 References 25 2 Simulation of the Sandia Flame 'D' 30 2.1 Introduction 31 2.2 Problem Description 2 2.3 Background2.4 Formulation 5 2.4.1 Chemical Source Term 32.4.2 Flow Solution 36 2.5 Results and Discussion 7 2.6 Conclusion 41 References 4 CONTENTS iv 3 Non-Premixed Methane Auto-Ignition 46 3.1 Introduction 47 3.2 Formulation 8 3.2.1 CSE Combustion Model 43.2.2 Flow Field 52 3.2.3 Basis Functions 3 3.3 Results and Discussion 8 3.3.1 Predictions of Ignition Delay 53.3.2 Relation to Prior Studies 60 3.3.3 Behaviour Prior to Ignition and Identification of Ignition 64 3.3.4 Effect of Turbulence on Ignition Delay 68 3.4 Conclusion 6References 71 4 Discussion and Conclusion 73 References 82 A RANS Fields\ V LIST OF TABLES 3.1 Allocation of Variance for Ignition Delay 59 3.2 Location and Mixture Fraction at tj, 60 vi LIST OF FIGURES 1.1 Borghi Diagram 8 1.2 Scatter Plot of Temperature vs. Mixture Fraction for Sandia Flame 'D\ x/D = 30 [1] 13 1.3 Relation Between Temperature and Mixture Fraction, Log of Scalar Dissipation. Portion of Figures 1 and 2 in [9] 17 2.1 Selected Basis Function Solutions 34 2.2 Predicted Temperature Field 8 2.3 CSE/LFD Validation - Comparison of Conditional Averages to Experimental Results (Solid lines predictions, symbols ex perimental measurements) 39 2.4 CSE/LFD Validation - Comparison of Axial and Radial Pro files to Experimental Results (Solid lines predictions, symbols experimental measurements) 42 2.5 Nitric Oxide Prediction 3 3.1 Flamelet Solutions for Temperature and Enthalpy Source Term 54 3.2 Flamelet Solutions for Species Mass Fraction - Yco and YCH3OH 55 3.3 Flamelet Solutions for Species Mass Fraction - Contour Plots 56 3.4 Flamelet Temperature [K] Contours for Complete Library ... 57 3.5 Effect of Ambient Temperature on Ignition Delay 58 3.6 Effect of Ambient Temperature on Ignition Delay - Compari son Prior Works 61 3.7 Evolution of Conditional Averages in Mixture Fraction Space, Tox = U50K 5 3.8 Evolution of Peak Increase in Conditional Averages 67 3.9 Comparison of Results for Different Ignition Criteria 67 3.10 Ratio Between Turbulent and Laminar Ignition Delays .... 69 4.1 Flamelet Solutions for T and Yco 76 LIST OF FIGURES A.l RANS Temperature and CH Creation Contours at td 84 Vlll NOMENCLATURE Symbols ~u Velocity ~x Spatial Coordinate T Temperature Y Species Mass Fraction Z Mixture Fraction £ Sample space for mixture fraction Z"2 Mixture Fraction Fluctuations (Variance) Le Lewis Number k Thermal Conductivity D Diffusivity Cp Heat capacity Ka Karlovitz Number Da Damkohler Number tk Kolmogorov time scale tc Chemical time scale tj Integral time scale tref Non-premixed reference time scale lref Non-premixed reference length scale uTef Non-premixed reference velocity scale v Kinematic Viscosity e Dissipation rate of turbulent kinetic energy p Density cjj Source term for species i X Scalar Dissipation Rate a Arbitrary Scalar an Coefficient for flamelet n ~a Vector of coefficients (ai,...,aN) ipjtk Integral of p(C') across discretization interval Notation ~v Vector M Matrix (T\Q Conditional average of quantity to left of |, conditioned on 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 in minimization NOMENCLATURE ix Subscripts i Species index n Basis function index j Computational cell (physical space) k Finite difference point (mixture fraction space) Abbreviations Pdf Probability density function CMC Conditional Moment Closure CSE Conditional Source-term Estimation RANS Reynolds Averaged Navier-Stokes LES Large Eddy Simulation DNS Direct Numerical Simulation X ACKNOWLEDGEMENTS I would like to thank the many people who supported and encouraged me throughout this project. First, I am grateful to my supervisor, Dr. W. Kendal 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 "Six pounds of cool in a five pound sack". Also, thanks are due to my friends and family (especially my parents) who kept my morale high, and life in the office interesting. Among the many people who have helped me out along the way are Colin, Chad, and Jorge, who were kind enough to proofread this thesis. I am also indebted to Colin for fighting the good fight with KIVA and sparing me the trials of adding transport equations to that particular code. As well, I would perhaps not have undertaken this project had it not been for the encouragement of Dave Mumford and Dr. Hill at Westport Innovations while I was an engineering undergrad. The financial assistance of Westport Innovations, the Natural Science and Engineering Research Council of Canada, and the BC Advanced Systems Institute was gratefully appreciated. Finally, I would like to express my gratitude to Amanda, 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 July 2004 INTRODUCTION INTRODUCTION xii There is significant motivation to improve the performance of combustion-based energy conversion devices. As 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 in eastern North America, 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 in tightening EPA and European acceptable limits for automotive emissions, or broad base policy initiatives such as the Kyoto Accord on greenhouse gas reduction. At the same time, the demand for affordable energy in developing economies has stretched the world's energy resources. On all fronts, it is obvious that we must improve our primary means of energy conversion — controlled combustion — in terms of both efficiency and minimizing the release of harmful emissions. Realizing this objective necessitates accurate tools for modeling the com bustion process. The design of devices which control the combustion process for greatest efficiency and least environmental impact can only be done with 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. While it is possible to build a prototype of, for example, a proposed combustion environment's geometry, obtaining the precise tem perature 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 main areas of interest which are tightly coupled. The macro-mixing behaviour of the flow can be treated using conventional computational fluid dynam ics (CFD) 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. Micro-mixing phenomena, and the resulting turbulence-chemistry interactions, are not directly captured by a CFD solution for en semble or spatial averages. Since direct numerical simulation (DNS) (which INTRODUCTION xiii does capture these effects) is beyond current computing limits for either real istic geometries or chemical kinetic mechanisms, models to account for these effects must be provided. The subject of this thesis is one such model, conditional source-term esti mation (CSE). After a general formulation of the model in the latter part of the next chapter, two papers are presented which describe the application of the model with little adjustment to two diverse problems using commercial CFD tools commonly employed by industry. The first, an extended version of that presented at the 2004 Canadian Section of the Combustion Institute Spring Technical Meeting1, addresses the simulation of an atmospheric pres sure co-flowing piloted methane jet. The second considers the auto-ignition of non-premixed methane under high pressure and moderate temperature such as found in a diesel engine, and is scheduled for submission in the near future. These two applications illustrate the flexibility of the model, and 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 to combustion modeling, will follow. 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. BACKGROUND 2 1.1 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 in 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. Typical reaction mechanisms for the combustion of methane in air may involve 40-50 species, and 200-300 elementary reactions [14, 15]. Solution of such a reac tion mechanism for a fixed composition requires solving a system of ordinary differential equations (ODEs) involving one equation for each species, plus temperature, which may evolve on timescales varying by several orders of magnitude. When transport is included, a numerical solution to the Navier-Stokes and species evolution equations is physically accurate and theoreti cally 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. The 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 in a non-reacting flow are inadequate for turbulent combustion. Numerical simulation of practical turbulent flows solve the governing equations for average values, in either an ensemble averaged (RANS) or spatially filtered (LES) sense, and provide a model for the mean effect of the unresolved (fluctuating or sub-grid scale) components based on the available mean values. While this approach is reasonably effective at approximating these effects in the eddy viscosity, attempting to estimate chemical source terms directly from mean quantities is hopeless — in a highly non-linear Arrhenius reaction rate expression, the fluctuations about the mean for temperature and species concentration will dominate the result. Assuming that the unresolved variation for all dependent variables of in terest is correlated, it is possible to relate the unresolved variation in all scalars to a representative scalar. It is then necessary to know the proba bility density function of only that representative scalar. In the literature, this is done in different ways to obtain two of the most popular modeling ap proaches for non-premixed turbulent combustion — Laminar Flamelet mod eling (LFM) and Conditional Moment Closure (CMC). The representative scalar most appropriate for non-premixed combustion is termed the mixture fraction, Z. It is a conserved scalar in the flow, often assumed to be trans-CHAPTER 1. BACKGROUND 3 ported with unity Lewis number, and represents the degree of mixedness. Neither created nor destroyed, conserved scalars are free from the encum brance of a source term. This feature allows accurate transport equations to be written for the mean and higher moments of conserved scalars. Express ing the mean values of unconserved dependent scalars (temperature, species mass fractions) in 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. The probability density function for the mixture fraction is commonly assumed to be a Beta pdf. With such a presumption, only the first two moments need to be estimated from the flow. Laminar Flamelet models and Conditional Moment Closure, while sharing the use of a conserved scalar, differ dramatically in their conceptual origins. As always, there is a trade-off between the performance of the model in terms of accuracy and economy. Laminar Flamelet models are conceptually straightforward, simple to implement, and easily incorporate the effects of complex chemistry with 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 in 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 in the resulting equations introduces modeling assumptions of questionable accuracy. As an alternative to these two presumed-pdf / conserved scalar models, a class of models known as 'pdf models exist, for which the evolution of the joint probability density function (pdf) of all 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. As there is little from this model which is essential to the body of this thesis, it will not be further discussed. The balance of this chapter will first review flamelet modeling and CMC, then finish with the development of the conditional source-term estimation model which is the central topic of this thesis. 1.2 Flamelet Models Pioneered for use in 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. By assuming a laminar structure, it is possible to close the transformed equations in a form decoupled from the original (~x,t) coordinates. The chemistry problem can then be solved in conserved scalar space only. This results in a highly cost-effective method, but unfortunately it is only strictly valid where it is justifiable to make the assumption of the laminar structure. Still, flamelet modeling is conceptually attractive and serves to provide an introduction to the use of a conserved scalar and presumed pdf in combustion modeling. As well, the flavour of conditional source-term estimation explored in this work makes use of solutions to the flamelet equations. In a general turbulent reacting flow, the key scalars (velocity, temper ature, 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) The mixture fraction (Z) is also a dependent variable naturally a function of (~x,t). Although an expression for the mixture fraction can be fabricated from the individual species concentrations by combining them so as to pro duce a conserved scalar and normalizing such that the resulting scalar takes a value of 1 in the fuel stream and 0 in pure oxidizer, convention in the simulation literature now follows a more convenient definition. This defini tion is simply that the mixture fraction is transported as a conserved scalar with unity Lewis number [4, 44]; that is, it satisfies Eq. 1.5 with Co = 0 and D = -J^r. Boundary conditions are Z = 0 in the oxidizer stream and Z = 1 in the fuel stream. Although some would argue that this definition lacks a clear physical interpretation, there is little point in debating the issue. First, for the purposes of the flamelet models (and, as we shall see later, CMC), the physical interpretation is irrelevant beyond the existence of a strong cor relation 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 depen dent 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 in terms of this now widely accepted definition. The 'flamelet equations' are developed in detail by Peters [29] for a glob ally unity Lewis number. An extended formulation incorporating the effects of differential diffusion was evaluated by Pitsch [32] for a test case involving a CHAPTER 1. BACKGROUND 5 turbulent jet flame where molecular diffusivity was of the same order of mag nitude as turbulent diffusivity in the near field. As 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. Noting the correlation between the mixture fraction and the other scalars of interest, Peters observes that it should be possible to perform a Crocco-transformation for the governing equations, which expresses dependent vari ables of interest in 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 depen dent 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. To 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+ dtdZ { ' ] d 92 3 (1.3) dxi dxi dZ d _ d dZ d dxk ~ dZk + dxkdZ' { " } He then applies them to the species transport equation (the diffusion term has been modeled using Fick's Law i.e.: V • ji « — pDVYi) : dY _^ + pH • VYi - V • (pDVYi) = pui (1.5) resulting in the transformed species transport equation: The operator R is defined in [29], and involves derivatives with respect to the higher dimensions (zT2,zT3). The first observation one may make is that the result, Eq. 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). This is obviously a mathematical necessity — we be gan with a 4-dimensional coordinate system, so certainly 4 new orthogonal 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 frac tion is negligible, Peters discards the terms included in R(Y{), and hence the dependence on Z2, Z3. In doing so, it is necessary to assume that the struc ture 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. This 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. The 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. Early 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, with boundary conditions Z(y —»• —oo) = 0, Z(y —> oo) = 1. The resulting functional form is: -2[erfc"1(2Z)]2 X*(Z) = —r -j. (1.7) ^ ' e-2[erfc-l(2Z0)]2 V ' Alternatively, from an analytical solution to a symmetric one-dimensional mixing layer where the fuel originates from the plane of symmetry, Pitsch et al. [34] obtain: *•<*>-f£i ™ Despite the observation by Kim et al. [19] that the latter form is incon sistent with the transport equation for the pdf of Z, it now commonly used, e.g., in [32]. Other forms for this relationship are typically used in CMC formulations, and will 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. At this point, flamelet models spread out into several different direc tions with varying success. Before considering some of the approaches, it is worthwhile to note the implication of the assumptions made thus far. In dis regarding 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. The 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. An alternative state ment is that the chemical reactions must occur on a time scale shorter than the Kolmogorov time scale. This can be visualized on the Borghi 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]. Although laminar burning velocities are reasonable for premixed combustion, for non-premixed combustion Lentini proposes using revised definitions for the reference scales; specifically: Iref = \/vtc Uref = yj^ (1.9) Constant values of two dimensionless numbers, the Karlovitz number (Ka) and the Damkohler number (Da) appear as straight lines on such a plot. The Karlovitz number is defined as the ratio of chemical time scale of the laminar flame to the turbulent (Kolmogorov) time scale [46]: Ka=^ tk The Damkohler number is defined as the ratio between the (integral) flow time scale and the chemical time scale [46]: Da^. tc Using Lentini's revised definitions of the non-premixed length and time scales, tc in the above definitions is replaced by trer = defined as in Eq. 1.9. In keeping with the earlier criteria, strictly speaking the validity of the flamelet models require tc < tk, or Ka < 1. This region is the portion of Figure 1.1 below the line shown for Ka = 1. There is also a regime where the turbulence interacts with the flame, but the flow does not completely overwhelm the flame. This region, commonly known as the 'Thickened flame' or 'perturbed flamelet' regime falls in the region of Figure 1.1 between the CHAPTER 1. BACKGROUND 8 togfy Figure 1.1: Borghi Diagram 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. Within the flamelet regime, he obtained reasonable results. However, just outside the flamelet regime in the perturbed flamelet regime Lentini noticed under-prediction in minor species such as OH of more than 50%. The implication is that while conceptually useful and convenient for providing engineering level solutions, basic flamelet models should be treated with caution due to a limited regime of applicability. With this in mind, let us now turn our attention to the various types of flamelet models. The most straightforward implementation of flamelet models is to neglect transient effects, such as extinction and re-ignition, and solve for steady so lutions to Eq. 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, in which the scalar dissipation (xo) appears as a parameter. Solutions to the system rep resent a balance between mixing (associated with 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 so lutions in a 2-D lookup table as functions of Z, x- To then recover the mean value of an arbitrary scalar, a, the problem is simply to evaluate: A common assumption (e.g., [23, 29]) is to assume statistical indepen dence of Z and xo- Then, 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 Beta pdf (e.g., [24, 44]); although early (limited) use was made of a 'clipped-gaussian' pdf [25] it has lost favour. The Beta pdf choice has been tested by Wall et al. [45]; it has also been noted by Liew 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 Beta pdf, transport equa-(1.10) CHAPTER 1. BACKGROUND 10 tions are typically solved in physical space for the mean and variance of mixture fraction. The 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 discus sion of models available for the mean scalar dissipation rate are presented by Sanders and Gokalp [39]. The 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 cx has a value of approximately 2. The second parameter for the log-normal pdf is often quite poorly estimated; based on observations made by Sreenivasan et al. (1977), Lentini [23] assigns the value a = 2, while Liew 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 CFD code Fluent [13], the Laminar 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 in this way has some serious drawbacks. As 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 in 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 in Eq. 1.6 and attempting to relate the Lagrangian r to an Eulerian time/space coordinate system. This 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. At 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 CHAPTER 1. BACKGROUND 11 jointly by Barths, Pitsch, and Peters [3] as the 'Representative Interactive Flamelet' (RIF), a good description of what is now known as the 'Eulerian Particle Flamelet Model' can be found in [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. The 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. This enthalpy is then used to de termine a new temperature to be fed back to the flow calculation. Transport equations for energy or individual concentrations are not solved in physical space; aside from the temperature interaction, there is no feedback from the flamelet code to the solution in physical space. These methods have been applied to problems such as the treatment of internal combustion engines [2] and mild combustion burners in furnaces [10] with reasonable success at predicting temperature and major species fields. While the majority of the use of flamelet models is in RANS simula tions, the model has also been applied to LES computations. Cook et al. [11] describe an application of steady flamelets in LES, using much the same approach as the earlier description; the various species concentrations are de termined through a lookup table based on mixture fraction and mean scalar dissipation. An implementation making use of an unsteady formulation is given by Pitsch [33], Pitsch and Steiner [35]. In the latter work, the formu lation and solution procedure closely resembles that commonly employed for the modeled CMC equations which will be described in the next section. The flamelet formulation, however, leaves the specific omissions from the equa tions unclear. Broad 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 CMC counterparts do produce improve ment, but it is unclear whether any additional unclosed terms may have been overlooked. CHAPTER 1. BACKGROUND 12 1.3 Conditional Moment Closure Conditional Moment Closure (CMC) methods seek to reduce the magnitude of the fluctuations about the mean quantities by conditioning the means on one or more closely correlated variables. Proposed independently by Bilger and Klimenko and unified in [22], the reward for reduced fluctuations about the conditional mean quantities is that simple first order closure for the chem ical source term is often sufficient. The dramatically different derivations suggested by Klimenko [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 in Figure 1.2, obtained from experimental measurements of the Sandia Flame 'D' [1] at 30 diameters downstream from the nozzle; this flame will be described in detail when it is considered in Chapter 2. At 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 in black) when mixture fraction is used as a conditioning variable will be much smaller than that about the unconditional mean. It is worthwhile to briefly consider the origin of the CMC equations, so that the modeling assumptions are clearly visible. To derive the CMC equations using the decomposition approach, Bilger makes use of the decom position: 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 in the use of Q in this section only for brevity and consistency with Bilger's notation. The 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 in a gen eral species transport equation, as in Eq. 1.5, Eq. 1.11 is differentiated as CHAPTER 1. BACKGROUND 13 2200 Mixture Fraction Figure 1.2: Scatter Plot of Temperature vs. Mixture Fraction for Sandia Flame 'D\ x/D = 30 [1] CHAPTER 1. BACKGROUND 14 necessary: dY = dQ dQdZ dY^ dt dt dC dt dt K ' Vr = VQ + ^ VZ + VY" (1.13) V • (pDVY) = V • (pDVQ) + ||v • (pDVZ) + pD(VZfd^ + PDVZ.V^ (I-") + V • (pDVY") Substituting the above into Eq. 1.5, using that Z also satisfies Eq. 1.5 with LJ = 0 to simplify, and taking the conditional average, Bilger arrives at the 'unclosed' CMC evolution equation for species concentration: <p|C> ?r + (p\0 0*10 • VQ - <P|C> (xlO ^Q where: ' \ i-s/ * ~e xrib/ \/vi\>/ _Qr2 = (p\Q(u\0 + eQ + eY (1.15) eQ = ( V • (pDVQ) + pDVZ • V^|C ) (1.16) / dY" _^ \ {p— + pH • VY" - V • (DpVY")\(^ (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• «PIC>(u"Y"\QP(0).. v• «PIC>P(0) ,„„ Where a gradient transport model has been used to close the term ( u "Y"\Q, similar to that used in the development of the transport equation for tur bulent kinetic energy, where the term u'u' is closed by the approximation u'u' « —DtVu. This approximation is commonly used to close terms involv ing a correlation of fluctuations in velocity and a scalar, e.g., in [31]. The CHAPTER 1. BACKGROUND 15 assumptions used for closure up to this point, along with a linear approxi mation for the conditional average of velocity, have been broadly adopted, e.g.,in [12, 16, 28]. The final closure approximation is the functional form for the conditional scalar dissipation. Kim et al. [16, 18] and Devaud and Bray [12] all make reference to an expression obtained by Girimaji (1992) by consideration of the pdf transport equation. Kim and Huh [16] also make reference to the counterflow mixing model used in flamelet calculations. In their work on spray autoignition, Kim and Huh [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. The computational cost associated with solving Eq. 1.15 using the de scribed approaches to obtain closure is significant. However, the conditional averages tend to vary slowly in 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! This has been done, e.g., by Swaminathan and Bilger [42], and is known as 0-dimensional CMC. In a homogeneous field, the conditional transport terms vanish, and the CMC species evolution equation reduces to a form reminiscent of the flamelet equa tions: d~i + (xlO = WO (i-w) Although this form of the CMC equations bears resemblance to the flamelet equations, there is an important difference in proceeding from the current state. While the assumptions made early in the derivation of the flamelet models disregard much of the physical behaviour, the rigorously justifiable derivation of the CMC equations, up to the point where closure assumptions are made, requires only the fundamental assumption that the variations in the various scalars are well correlated with the conditioning variable. This means that the framework to reconsider or improve the clo sure for a given term is available. All of the necessary physics are present, and we are consciously aware of the assumptions made to neglect or model various terms. This makes the model ideally suited to improvements not only in modeling, but in basic understanding as well. While quite successful in some circumstances, variations about the con ditional averages in violation of the fundamental assumption can lead to errors in first moment closure for the conditional reaction rates. The reme-CHAPTER 1. BACKGROUND 16 dies for this vary; allowing for spatial variation of the conditional averages, in 1 or more dimensions, can alleviate the difficulty. This is successful as large variations about the conditional means are often the result of local ized 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 Kim et al. [17], Mastorakos and Bilger [27]. However, this requires modeling of several unclosed terms in the transport equation for the conditional fluctuations in 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. mix ture fraction from a direct numerical simulation of a flame with substantial levels of local extinction performed by Sripakagorn given as Figure 1 in [9]. In a case where the variation about the conditional average is significant, as in 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 with the scalar variations not correlated with mixture fraction. Work 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 with promising results [26]. In a realistic combustion simulation, the choice of second conditioning variable is unclear. Cha 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 will ultimately prove unsuitable. As the strain on a flame increases, the local gradients become steeper and the rate of mixing increases. For a mild increase in mixing rates a corresponding increase in reaction rates, and a moderate change in the species conditional averages is observed. When the rate of mixing becomes too great, the tem perature and intermediates are mixed away too quickly to be replaced by chemical reaction and the flame is extinguished. Once a flame has been ex tinguished, there is no longer a relation between the temperature and scalar dissipation — regardless of how the flame is strained, the temperature will 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. An increase in 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 concen trations have dropped too far. The various scalar dissipation perturbation possibilities give rise to a plethora of possible scalar fields for the same scalar dissipation magnitude. The right side of Figure 1.3, which forms a portion of Figure 2 in [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 with temperature. Obviously, this is not so for the scalar dissipation. Z 'Og Xst Figure 1.3: Relation Between Temperature and Mixture Fraction, Log of Scalar Dissipation. Portion of Figures 1 and 2 in [9] The preceding discussion suggests that there are two factors to consider when evaluating reaction rates: the mixture fraction, and the elusive addi tional quantity that somehow represents the cumulative effects of the mixing field on the flame throughout its recent history up to and including the present. In the next section the detailed derivation of an approach which has the potential to incorporate these two considerations will be considered. 1.4 Conditional Source-Term Estimation Conditional source-term estimation (CSE) ultimately uses the same first or der closure based on conditional means to evaluate the chemical source terms CHAPTER 1. BACKGROUND 18 employed in CMC. Rather than solving directly for the conditional averages, the flow is formally solved for the unconditional means in physical space. The conditional means, which are assumed constant in 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. This method was tested in an a priori sense against a DNS database [7] and was found to be effective when little extinction was present using only the mixture fraction as a conditioning variable. The 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 LES calculation of the Sandia Flame 'D' using 2-step chemistry with substantial success [40]. More recently, Bushe and Steiner [8] proposed that the conditional means and source terms be decomposed into a linear combination of solutions to the Laminar 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 with the flamelets are then combined and integrated over the pdf of mixture fraction to close the enthalpy source term. This method produced promising results in an a priori test in [8]. Later, Blair [6] used this formulation in a KIVA-II simulation of a turbu lent methane jet injected into hot air. However, Blair's implementation was unable to successfully predict the ignition delay when compared to exper imental data. In the remainder of this section, the general formulation of CSE with Laminar Flamelet Decomposition used throughout this thesis will be presented. Following Bushe and Steiner [7], first note that the unconditional aver age 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"2^j [22], the relationship then holds for Favre averages: (1.20) As per [8], decomposition of the conditional average in question into a sum-CHAPTER 1. BACKGROUND 19 mation of basis functions, e.g., <T|C> « 5>n <T|0n N (1.21) n=l can be expected to be valid for finite N provided the basis functions (T\Qn capture characteristics of the true conditional average. Notably, it must be possible to reconstruct the potentially sharp gradients with respect to mixture fraction. Substitution of Eq. 1.21 into Eq. 1.20 produces: Jo N .n=l (1.22) At this point, we can see that Eq. 1.22 is a Fredholm Equation of the first kind: with g(t)= j K(t,s)f(s)ds, J a s = ( a = 0 b = 1 9 = f K=p(Q N / = 5>n(T|C>n. (1.23) (1.24) (1.25) (1.26) (1.27) Choosing an ensemble of J cells over which an=i..JVj=i..j is constant, we can write Eq. 1.22 numerous times, and form a system of equations resem bling (18.4.5) in [37] (adapting to match notation): 9j Jo N n=l (1.28) Given J distinct 'observations' of Qj (possibly contaminated with noise), each with 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. BACKGROUND 20 sense to obtain an approximation of the underlying function at the quadra ture points. As the resulting system of equations is typically ill-posed (it could only be well posed if each of the N observations came conveniently supplied with a response kernel that was a delta function at each of the quadrature points), the recommendation in [37] is to add a priori informa tion in the form of smoothing to stabilize the solution. In our case, we do not actually wish to solve for the function /(£) = ^2n=i a™ CHOn a^ quadrature points in £-space, but rather for the basis function coefficients an. As will be discussed in the following chapters, that a„ should have any elementary func tional dependence on the position in the library, n, is unphysical. Instead, a more physical a priori belief is that, if started from a physically realiz able condition, the solution should evolve smoothly in time. Turning now to the algebra necessary to accomplish this, we revisit a rearranged version of Eq. 1.22 and switch to vector/matrix notation: Jo (1.29) Forming a if-point quadrature in £-space, / P(C)/(CK « 5>(C*)/(CkKAfc Jo fc=i K E fc=i 1 /"Cfc+1 5 / p(C)dC K /(a) = £>/(&)• (i-3o) fe=i A one-sided approximation (e.g., tbx = /^_1p(C')^C) is substituted for il>i,4>K- Using Eq. 1.30 to discretize Eq. 1.29 and writing it for each of J computational cells in an ensemble over which we assume ~a is spatially homogeneous in the fashion of Eq. 1.28, neglecting rij, we obtain: f(T\Ok=l,n=l ^j=l,fc=l • • • ^1,K \TjJ \ ... jjjcj V (T\C)K,i or, where A is defined by comparison to Eq. 1.31, (T\C)K,N) a (1.31) Aa. (1.32) CHAPTER 1. BACKGROUND 21 We can also write an equation resembling Eq. 1.20 for any scalar of in terest. Repeating the above procedure for an arbitrary species mass fraction Yi, and defining: Bi = /^I0fc=1,n=1 ••• (Yi\OhN we can also form the system of equations: (1.33) Yi « Bi~a. (1.34) Following Bushe and Steiner [8], we suppose that it is possible to recon struct the flame of interest from a linear combination of N basis functions, where the nth basis function solution for (T\()n, (li|C)„ are corresponding scalar profiles from the same quasi-physical flame. In this situation, solution of either Eq. 1.32 or Eq. 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 ques tion of existence, there is a significant question regarding uniqueness. Later, in 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\Qn 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 to be solved be expanded by solving, simultaneously, Eq. 1.32 and Eq. 1.34, using as many species mass fractions as necessary to differentiate between the basis functions in the library. We now have: M~a ?s R, (1.35) with: M = ( A \ \BIJ R = \YiJ (1.36) Despite our hope that Eq. 1.35 contains sufficient information to select the most suitable combination of basis functions, it is still potentially ill-posed. We have doubts both about the uniqueness of the solution and the CHAPTER 1. BACKGROUND 22 existence of a suitable solution. The 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 Eq. 1.30, etc.), and in 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 turn here to the more general Tikhonov Reg-ularization [43], rather than the Linear Regularization recommended in [37] and used by Bushe and Steiner [8]. While both of these regularization ap proaches include introducing a priori information into the solution, the gen eral case allows us more freedom in choosing that information. This 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 con straints. Tikhonov suggests reformulation of the problem as a minimization problem, where we simultaneously minimize the residual in the problem we wish to solve as well as the deviation of the solution from an expected value: This reformulation as a minimization problem addresses the question of existence of a solution. However, it means that we must take care in assem bling Eq. 1.35 to normalize both sides of each equation, lest the minimiza tion procedure place greater emphasis on satisfying those equations involving scalars with larger absolute magnitudes. The alternative formulation of the minimization which perhaps makes this more clear is: Considered in this form, it is obvious that we can adjust the influence of each scalar by adjusting the coefficients A0.../+i. However, since we wish to satisfy the equations for all scalars with the same weight, we can accomplish our objective simply by normalizing Eq. 1.32 and Eq. 1.34 by the maximum present so that all of the scalars fall in the range [0,1]. Returning to the original minimization problem Eq. 1.37, the solution in a least squares sense (1.37) (1.38) CHAPTER 1. BACKGROUND 23 is given by the solution to the linear system of equations: (M*M + AJ) ~a = M*~R + Ac*0. (1.39) We now have an opportunity to provide an arbitrary initial solution of0, 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. While the smoothing used for stabilization in earlier work by Bushe and Steiner [8] implies a functional dependence of an on n, the claim in this work is that such functional dependence is unphysical. Insistence on such dependence implies that if, for example, the order in which the basis functions are stored or the spacing between the basis functions in terms of flamelet time is altered, the physics captured will be altered. As will be discussed later, this is clearly unacceptable. In the two applications of this method discussed in 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. The above formulation has a great deal in common with that proposed by Bushe and Steiner [7, 8] and used by Blair [6], Steiner and Bushe [40]. However, there are significant differences, which are believed to be responsible in part for the relative success of the implementations presented in this thesis. The two most notable which have appeared so far include: • Regularization Method. All of the previous formulations imposed smooth ing on the quantity sought by the integral inversion. In [40], progress variables were solved for, which do tend to be smooth in mixture frac tion space. However, this necessitates inversion for a separate progress variable for each reaction in the mechanism used, which is a pro hibitively large set for realistic chemistry. In [6, 8], the smoothing was for the coefficient in basis function index space; the drawbacks to this have already been touched on, and will be discussed further in the following chapters. • Inversion Based on Multiple Species. While 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 with the mean fields for minor species, which discriminate much more clearly between the flamelets. CHAPTER 1. BACKGROUND 24 Additional 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. REFERENCES 25 REFERENCES [1] 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 CH4/Air Flame D. [2] H. Barths, C. Hasse, et al. Simulation of combustion in direct injection diesel engines using an eulerian particle flamelet model. In Twenty-Eighty Symposium (International) on Combustion, pages 1161-1168. The Combustion Institute, 2000. [3] H. Barths, H. Pitsch, and N. Peters. 3d simulation of di diesel combus tion 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 Com bustion Science, l(2-3):87-109, 1976. [5] R. Bilger. Conditional moment clsoure for turbulent reacting flow. Physics of Fluids A - Fluid Dynamics, 5(2):436-444, 1993. [6] C. Blair. Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of British Columbia, 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, July 1999. [8] W. K. Bushe and H. Steiner. Laminar flamelet decomposition for condi tional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [9] CM. Cha, G. Kosaly, and H. Pitsch. Modeling extinction and reigni-tion in 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 mild 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 in 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. Combus tion and Flame, 132:102-114, 2003. [13] Fluent 6.1. Fluent Inc., 2003. [14] D. M. G. Gregory P. Smith et al. Gri-mech 3.0. URL http: //www. me. berkeley. edu/gri-mech/. [15] J. Huang, P. Hill, W. Bushe, and S. Munshi. Shock-tube study of methane ignition under engine-relevant conditions: experiments and modeling. Combustion and Flame, 136(l-2):25-42, 2004. [16] S. H. Kim and K. Y. Huh. Use of the conditional moment closure model to predict NO formation in a turbulent CH4/H2 flame over a bluff-body. Combustion and Flame, 130:94-111, 2002. [17] S. H. Kim, K. Y. Huh, and R. W. Bilger. Second-order conditional moment closure modeling of local extinction and reignition in turbu lent non-premixed hydrocarbon flames. Proceedings of the Combustion Institute, 29:2131-2137, 2002. [18] S. H. Kim, K. Y. Huh, and R. A. Fraser. Modeling autoignition of a turbulent methane jet by the conditional moment closure model. Pro ceedings of the Combustion Institute, 28:185-191, 2000. [19] S.-K. Kim, S.-M. Kang, and Y.-M. Kim. Flamelet modeling for com bustion processes and NOx formation in the turbulent nonpremixed CO/H2/N2 jet flames. Combustion Science and Technology, 168:47-83, 2001. [20] W. T. Kim and K. Y. Huh. 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 in turbu lent flow. Fluid Dynamics, 25(3):327-334, 1990. [22] A. Klimenko 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 Technol ogy, 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:199-213, 1984. [25] F. Lockwood and A. Naguib. The prediction of the fluctuations in 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 British Columbia, 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 and Y. Wright. Simulations of turbulent spray auto ignition with elliptic conditional moment closure. In Proceeding of the European Combustion Meeting, 2003. [29] N. Peters. Laminar diffusion flamelet models in 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". Com bustion and Flame, 116:675-676, 1999. [31] N. Peters. Turbulent Combustion. Cambridge University Press, New York, 2000. REFERENCES 28 [32] H. Pitsch. Unsteady flamelet modeling of differential diffusion in turbu lent jet diffusion flames. Combustion and Flame, 123(3):358-374, 2000. [33] H. Pitsch. Improved pollutant predictions in 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.Chen, and N. Peters. Unsteady flamelet modeling of tur bulent hydrogen-air diffusion flames. In Twenty-Seventh Symposium (International) on Combustion, pages 1057-1064. The Combustion In stitute, 1998. [35] H. Pitsch and H. Steiner. Scalar mixing and dissipation rate in large-eddy simulations of non-premixed turbulent combustion. Proceedings of the Combustion Institute, 28:41-49, 2000. [36] S. Pope. Pdf methods for turbulent reactive flows. Progress in Energy and Combustion Science, 11:119-192, 1985. [37] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, UK, 1992. [38] G. Ruetsch, L. Vervisch, and A. Linan. 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 in variable density turbulent axisymmetric jets and diffusion flames. Physics of Fluids, 10(4):938-948, April 1998. [40] H. Steiner and W. Bushe. Large eddy simulation of a turbulent reacting jet with conditional source-term estimation. Physics of Fluids, 13(3): 754-769, March 2001. [41] N. Swaminathan. Flamelet regime in non-premixed combustion. Com bustion and Flame, 129:217-219, 2002. [42] N. Swaminathan and R. Bilger. Assessment of combustion and submod els 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, Winston, New York, 1977. Translation of Metody resheniia nekorrektnykh zadach. [44] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress in Energy and Combustion Science, 28:193-266, 2002. [45] C. Wall, B. J. Boersma, and P. Moin. An evaluation of the assumed beta probability density function subgrid-scale model for large eddy simula tion of nonpremixed, turbulent combustion with heat release. Physics of Fluids, 12(10):2522-2529, October 2000. [46] J. Warnatz, U. Mass, and R. Dibble. Combustion: physical and chem ical fundmentals, modeling and simulation, experiments, pollutant for mation. Springer-Verlag, New York, 2nd edition, 1999. 2. SIMULATION OF THE SANDIA FLAME 'D' CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 31 2.1 Introduction Computer simulations of turbulent reacting flows seeking to solve for mean quantities in either a filtered (Large Eddy Simulation (LES)) or Reynolds Average (RANS) sense are hampered by the closure problem in turbulent combustion, that is, determining the mean source-terms due to chemical reaction. The difficulty arises because unresolved fluctuations about mean quantities dominate the reaction rate expressions. Conditional moment closure (CMC) [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 fluctua tions about the mean quantities. This allows for much improved closure of the source-term compared to evaluating the source-term with uncondi tional means. However, this method requires adding another dimension to the solution for each conditioning variable leading to a substantial increase in computational cost. Further, although applicable to arbitrarily complex chemistry in principle, an additional transport equation must be solved for each species added. Conditional source-term estimation (CSE) [50, 58] al lows use of the CMC hypothesis in the LES paradigm by solving for only unconditional means and then inverting integral equations to determine the conditional means necessary to close the source-terms. This method can also be applied in the RANS paradigm resulting in substantial computational sav ings compared to traditional CMC. However, using realistic chemistry still results in a prohibitively large number of transport equations. Laminar flamelet models [55] have a substantial advantage in terms of computational efficiency in that the chemistry is solved a priori and tabu lated in terms of representative scalars such as the mixture fraction used as a conditioning variable in CMC. However, these models are, strictly speak ing, valid only within the flamelet regime and there is considerable concern regarding their applicability outside this regime. As well, except in limited circumstances, it is impossible to represent a flame of practical interest with a single flamelet solution. A recently proposed extension to CSE, where the conditional averages are solved for in 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. This proposal has been tested by comparison to results of a DNS simulation in the a priori sense and the RANS paradigm. CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 32 In this chapter, the method has been adapted for use in a commercial RANS code and used in a predictive simulation. 2.2 Problem Description The arrangement considered in this paper is a piloted methane-air jet flame. The 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. This geometry has been utilized to generate extensive experimental data for a variety of flames by Sandia National Laboratories [47, 48]. The specific conditions used here are those of Sandia Flame 'D'. In this flame, the fuel is 25% CH4 and 75% air, and the jet Reynolds number is 22,400. Under these conditions, there is only a small probability of local extinction occurring [47]. The pilot flame is a premixed C2H2 / H2 / CO2 / N2 / air flame used to stabilize the diffusion flame on the burner. Conditional source-term estimation with 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. The conditional source-term is then obtained from the corresponding flamelet solutions, and integrated over the Favre density function (probabil ity density function with 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 in terms of the conditional averages and Favre density function of mixture fraction: Given 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. The choice of the ensemble of cells over which the structure of the flame is constant in a CFD computation is left to 2.3 Background (2.1) 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)B , (2.2) n=l where N need not be large if the basis functions Tn(£) 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 an: f>\_ (£zLi«n(TK)npi(<;)<K\ (23) 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 arbi trarily complex chemistry. If a variety of scalar dissipation histories are incorporated, the likelihood of obtaining solutions which resemble the phys ical 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 nec essarily 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 Eq. 2.3, we note that the equation is a Fredholm equation of the first kind and replace the integral with a numerical approximation [56]. With this observation, it is straightforward to cast Eq. 2.3 into a system of approximate linear equations of the form M~a « R, where the number of rows in M need not equal the number of rows in "a , and the problem may be ill-posed (implying possible problems related to existence, uniqueness, and continuity). The approximate nature of the equality adds further difficulties in that the data on the right hand side are contaminated with noise due to the accumulated errors in 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 Eq. 2.3 for an, the mean chemical source-term can then be recovered in the manner of Eq. 2.1. Once a suitable choice for ~a has been made it is possible to obtain, using the same 7* in the forward problem, the conditional averages for any arbitrary species. The conditional source terms for all species can then be CHAPTER 2. SIMULATION OF THE SANDIA FLAME lD' 35 computed using first conditional moment closure. For efficiency, these source terms can also be determined in advance and obtained from a table. 2.4 Formulation 2.4.1 Chemical Source Term In this work, closure of the chemical source term was found by extending CSE with 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 in the context of a realistic simulation covering a wide range of conditions, in order to represent the physics present in the library of flamelet solutions, the conditional average of temperature alone is insufficient to discriminate between flamelets. We therefore augment our system of integral equations with similar relations for the mass fraction of CO: ( f, \ YcOi V / /o«»<r|C>„p(C)dC \ Slan(YCo\Qnp{QdC 1 (2.5) ) For this case, where the configuration under consideration is the Sandia Flame '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 in our system for each cell in these ensembles. Our hypothesis, supported by inspection of the library used and noting that these two scalars are suffi cient 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. We employ Tikhonov Regularization [59], where the functional to be min imized takes the form: min{||M~a - R\\2 + A||c? - 7^°||2} which has the least-squares solution: (M*M + XI) 7^ = M*~R + Aa°. (2.6) (2.7) CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 36 In Eq. 2.6, there is a balance between minimization of the error in 'fitting' the unknown function to the data and incorporating a priori information about the function. The 'regularization parameter' A was chosen to add just enough a priori information to produce a well-behaved solution. We 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. This 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 in the mean sense such as this the choice is even more logical. As well, beginning with a physically realizable state, even in the absence of any useful data for the inversion process a physical solution will result. Other regularization possibilities (e.g., smoothing) lack this justification and may not lead to a physically realizable solution. For example, Blair [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 RANS simulation of the flame was performed using the k — e turbulence model with standard model constants, except for the value of Ce2 = 1.8, in preference to the standard value of C£2 = 1.94 to obtain better agreement with the experimentally measured mixing field near the nozzle exit. In ad 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. The 2-d axisymmetric mesh used encom passed 130 diameters in the axial direction, and 90 diameters in the radial direction. 120 cells in the axial dimension, stretched with a consecutive ratio of 1.04, and 78 cells in the radial direction, stretched a variable consecutive ratio ranging from 1.05 — 1.09, were used. The size and stretching of the mesh was chosen to closely match that used by Hinz 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. The solution was obtained using the iterative solver in Fluent [52]. CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 37 2.5 Results and Discussion Solutions to the laminar flamelet equations with boundary conditions match ing 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 co-flow, at Z — 1 matching the fuel stream, and at Z — Zpu 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 Sec tion 1.2 to extract the mean scalar dissipation at the stoichiometric interface as a function of distance from the nozzle. Assuming that the stoichiomet ric 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 correspond ing to higher scalar dissipation rates. The solution was sampled at regular intervals to build the library for use in the following RANS simulation. How ever, when solutions corresponding to flamelets near extinction events were included in the library, the RANS 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 ex tinguished 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 sin-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 dramati cally. 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 phys ical and mixture fraction space despite the severe restriction of having only burning basis functions available for the decomposition. However, it is ap parent 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 0.4 0.6 Mixture fraction x/d=45 2500 2000 1500 1000 500 2500 2000 1500 1000 0.4 0.6 Mixture fraction x/d=45 0.2 0.4 0.6 0.8 Mixture fraction Figure 2.3: CSE/LFD Validation - Comparison of Conditional Averages to Experimental Results (Solid lines predictions, symbols experi mental 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 • VYNO - V • [pDVYNO) = (^NO~ (2.8) Z^> = Jl-a- (YNO\c)p(QdC. (2.9) CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 40 The solid line, in contrast, is obtained by evaluation of: YNO= I'a • (YNO\6p(0d<; (2.10) Jo 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: with the scalar dissipation resulting from combining the flamelets chosen: 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. Xa = 2D (VZ • VZ) (2.11) (2.12) CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 41 2.6 Conclusion An 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 temper ature and CO are promising. Further work is necessary in 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. As well, basis functions which are representative of the edge flame expected to provide (re)ignition in this flame will be necessary in the future. While the prediction of nitric oxide concentration is poor, that trends are predicted and the result is consistent with what would be expected in light of the temperature over-prediction is encouraging. CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 42 CSE/LFD RANS Sandia Flame D Validation - Centeriine 2500 2000 1500 1000 500 a % CSE/LFD RANS Sandia Flame D Validation - Centeriine 0.08 r 0 20 40 60 80 x/d CSE/LFD RANS Sandia Flame D Validation - x/D=30 2500 2000 1500 1000 500 0 2 4 6 8 r/d CSE/LFD RANS Sandia Flame D Validation - x/D=45 2500 2000, 1500 1000 500 8o.04 0 20 40 60 80 X/d CSE/LFD RANS Sandia Flame D Validation - x/D=30 CSE/LFD RANS Sandia Flame D Validation - x/D=45 0.04 Figure 2.4: CSE/LFD Validation - Comparison of Axial and Radial Profiles to Experimental Results (Solid lines predictions, symbols exper imental measurements) CHAPTER 2. SIMULATION OF THE SANDIA FLAME lD' 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 Flame D. [48] R. Barlow and J. Frank. Effects of turbulence on species mass fractions in methane/air jet flames. In Proceedings of the 1998 27th International Symposium on Combustion, Aug 2-Aug 7 1998, volume 1, pages 1087-1095, Sandia Natl Lab, Livermore, CA, USA, 1998. Combustion Inst, Pittsburg, PA, USA. ISBN 0082-0784. [49] C. Blair. Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of British Columbia, 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, July 1999. [51] W. K. Bushe and H. Steiner. Laminar flamelet decomposition for condi tional 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. URL http: //www. me. berkeley. edu/gri-mech/. [54] A. Klimenko and R. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [55] N. Peters. Laminar diffusion flamelet models in non-premixed turbulent combustion. Progress in Energy and Combustion Science, 10:319-339, 1984. [56] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, UK, 1992. REFERENCES 45 [57] G. Ruetsch, L. Vervisch, and A. Linan. 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 with conditional source-term estimation. Physics of Fluids, 13(3): 754-769, March 2001. [59] A. Tikhonov and V. I. Arsenin. Solutions of ill-posed problems. Halsted Press, Washington, Winston, New York, 1977. Translation of Metody resheniia nekorrektnykh zadach. 46 3. NON-PREMIXED METHANE AUTO-IGNITION CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 47 3.1 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 source-term estimation with laminar flamelet decomposition approach proposed in [62]. The ambient oxidizer conditions — the high pressure and moderate temperatures characteristic of compression ignition engines — were chosen with the intent to validate the combustion model used under engine relevant conditions. This was obtained by comparison of the predicted ignition delay to experimental results obtained from the shock-tube facility at UBC at several initial temperatures. Early experimental measurements of the ignition delay time were obtained by Fraser et al. [63], who injected methane into a constant volume combus tion bomb following the premixed combustion of a hydrogen-ethylene-oxygen-nitrogen 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]. More recently, Ishiyama et al. [67] used a similar procedure to repeat Fraser et al.'s work in their own facility. Both studies used pressure measurements to identify ignition; either an increase past a threshold based on pre-ignition measurements [63, 73] or a sudden change in the rate of pressure rise [67]. Very recent experiments still underway at the University of British Columbia 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 in [66]. The criterion used to define an ignition event in the UBC study differed from the earlier work using pressure based criteria in that the UBC study identified the onset of ignition with optical ob servations. Through the use of a high-speed camera, the UBC study defined the end of the ignition delay as the time at which optical emission leading to a sustained optical kernel first appeared [77]. The results of these independent efforts show some degree of inconsistency, especially at low oxidizer temper atures and correspondingly longer ignition delays. However, the most recent results [67, 78] agree well with one another. Given the dramatic differences in the experimental apparatus and techniques used to detect ignition, it is remarkable that such agreement exists, and inspires substantial confidence in the validity of the results obtained. CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 48 On the simulation front as well, previous efforts have been made to de velop a computational method to model ignition. Both Conditional Moment Closure [68] and Flamelet models [69] have been applied to the simulation of the experiment described by Fraser et al. [63]. Both methods produced results which agreed remarkably with the experimental data presented by Fraser et al. [63]. However, the subsequent introduction into the literature of experimental data disagreeing with Fraser et al.'s results by Ishiyama et al. and the anticipated publishing of supporting data from the UBC shock tube give reason to be somewhat cautious regarding these simulations. In the later discussion, some speculation will be made regarding the possible ex planation for the discrepancy, but further careful study both experimentally and through modeling exercises is necessary. In addition to providing an as sessment of the performance of the CSE-LFD combustion model, this work seeks to investigate in some detail the process leading up to the detection of an ignition event, with the goal to produce an ignition delay prediction that is directly comparable to that measured in the UBC 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 with predicting the location of ignition. 3.2 Formulation 3.2.1 CSE Combustion Model Closure of the chemical source term was obtained by an extension to con ditional source-term estimation [61] whereby the conditional averages of key scalars (T, CO &c CH3OH mass fraction) are decomposed into a collection of solutions to the unsteady laminar flamelet equations. The details of the method have been described in Chapter 1; however, it may be useful to review the general approach taken before considering the specifics of this implementation. In general, conditional source-term estimation seeks to utilize the (first moment) CMC hypothesis [70] to treat the chemical source terms, but with out 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 in Eq. 3.1 can be assumed constant for some ensemble of cells in a simulation where the other quantities appearing (Yi,p(C.)) differ, a system of integral equations can be formed which can in theory be solved for a discreet ap proximation to the conditional mean in question. This has been done in the LES paradigm and produced promising results using strongly reduced chem istry and first moment closure based on the conditional means to close the chemical source terms [76]. The unconditional mean source terms are then obtained by evaluating the forward problem offered by Eq. 3.1. Incorporating the effects of complex chemistry necessary to treat the au toignition process [79] using the method described above would require solv ing a transport equation for every species involved, as well as evaluating the reaction rate expressions for each step in the mechanism. This imposes pro hibitive computational costs if a realistic mechanism is used. Instead, the method proposed in [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 Eq. 3.1 becomes: _ f1 N T« / Y,an(T\0np{Z)dZ, (3.2) Jo n=l where the approximation: <T|C> « f>n (T|C)n n=l has been used. To obtain the library of basis functions (T\Qn, we solve the unsteady laminar flamelet equations in 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 all species and temperature) present at each point in mixture fraction space. If we are within the flamelet regime, it would be reasonable to suppose that the coefficients an which best satisfy Eq. 3.2 would be consistent regardless of which scalar was considered. This implicitly 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 an 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. With this observation, it is obvious that determination of an based on temperature alone has some serious limitations. Instead, it is proposed here that transport equations be solved for, and the determination of an 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) and methanol (YCH3OH), were all used to form the set of integral equations for inversion. The system of integral equations written with the intention of solving for an is then written: / fj \ ( Ji^iTlO^QdC \ YCOJ ^ fo«n(Yco\OM)d( YCH302 Jo «n (YCH3OH\0iP(Q<K V : / V / Writing the above system of equations for each cell in an ensemble over which the conditional average is considered constant, a realistic determina tion 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 on these planes that the conditional average is assumed to be constant. A practical difficulty is encountered in that solving this system is not trivial. Once the integrals are replaced by numerical approximations, the resulting system is extremely poorly posed. To remedy the ill-posedness, the problem is reformulated as a minimization problem and solved in a least squares sense, balancing the solution between 'goodness of fit' to the data appearing on the left hand side of Eq. 3.3 and some a priori knowledge regarding the solution. While constraints on the first derivative of the solution for an with respect to basis function order in the library have been used in earlier work [62], this approach has several drawbacks. The relationship between flamelets which reside next to each other in the library is not necessarily unique, as it depends on the manner in 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 in time, and the time elapsed in the calculation between the extraction of basis functions may vary. As a matter of convenience, it is desirable that the order in which the basis function are stored in the library have negligible effect on the solution, which can never be the case when some functional form for an[n) is required. In the present formulation, where the system is evolving in time in a physical manner, the a priori information is provided through simultaneous minimization of the difference between the left and right hand sides of Eq. 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. The actual equation solved is that given in Eq. 1.37 (Chapter 1, page 22). The relative penalty for errors in matching the unconditional means and changing the coefficients between time steps (A in Eq. 1.37 ) is chosen such that the solution will only be pulled toward the a priori guess the minimum amount necessary to produce a well behaved solution. In addition to a more plausible physical in 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 un conditional mean of any scalar of interest present in the library can be de termined by evaluation of the forward problem given by Eq. 3.2, replacing temperature with the scalar of interest. Closure for the chemical source term CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 52 is also recoverable by convolution of the appropriate conditional mean reac tion rate constructed by linearly combining the source terms from each basis function weighted by the corresponding an and the pdf of the conditioning variable in that instance. The reader may have noticed a switch in Eq. 3.3 to Favre-averaged quantities; this implies use of the Favre density function, p(£; Z, Z"2), in place of a probability density function, Z, Z"2). To re cover the non-Favre averaged source terms necessary for the RANS transport equations, the forward problem to evaluate becomes: where the relation between the Favre and non-Favre pdf suggested by Kli-menko and Bilger [70] has been utilized. An axisymmetric computational domain, consisting of 55 cells in the ra dial direction and 100 cells in the axial dimension, and covering physical dimensions of 0.05m x 0.15m, was used for the computation. The mesh was stretched with a successive ratio of 1.05 and 1.01 in the radial and axial direc tions, 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 CSE inversion process. The stretching in the axial direction allowed for an especially dense grid in the region containing the plume; this is especially important in this type of simulation as the conditional average is inverted for independently on each plane of cells normal to the jet centerline. The implicit assumption is that the amount of local extinction or ignition is rea sonably consistent on planes normal to the jet axis. An overly coarse grid not only results in a less accurate solution to the scalar transport equations, but restrains the variation of the conditional averages in physical space and makes the model more susceptible to errors resulting from spatial changes in the flame structure. It should be noted that this grid is overwhelmingly more fine than those used in early studies employing CSE for simulation of autoignition process [60]. In addition to properly resolving the flow, it was desired to make more data points with mixture fraction between pure fuel and pure oxidizer available to the CSE inversion process. (3.4) 3.2.2 Flow Field CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 53 The domain was initialized with quiescent air at 30 bar and zero mixture fraction and variance; a pure methane jet entered from the bottom left. The boundary condition for the methane jet was a pressure boundary condition resulting in a moderately under-expanded jet. Although the actual pressure at the nozzle outlet in 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 with the pressure boundary in the simulation was performed; in general, the results are insensitive to the pressure and 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 in the experiment, the inflow was terminated. The governing equations were discretized using a first order upwind scheme, and an implicit time step of 10/is was chosen to advance the simulation in time. The chemical source term, however, was computed only at the start of iteration for every time step, and so was effectively advanced explicitly in time. Inversion for the discrete conditional averages was performed on a 51 point grid in mixture fraction space, with the majority of the points clustered around the stoichiometric mixture fraction. 3.2.3 Basis Functions Prior to undertaking the flow simulation for each library of basis func tions was constructed by evolving the unsteady laminar flamelet equations (described in Section 1.2) for temperature and species mass fraction in time using realistic chemistry developed by Huang et al. [65]. This mechanism was used without adjustment of the various constants tuned to premixed autoignition experiments. . An internally adiabatic mixing solution between boundary conditions matching the fuel and oxidizer initial states was used as an initial condition. The coupled partial differential equations were then reduced to a system of ordinary differential equations using a method of lines approach1; the system of ODE's was solved using LSODE [64]. Every 50 [is in flamelet time, the scalar values and their time rate of change due to chem ical reaction only were output. The calculation was terminated at 10ms, well after ignition had occurred. The scalar dissipation (x) was presumed to have 1The author acknowledges that the flamelet code modified for this purpose was provided by W.K. Bushe CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 54 the functional form: X(Z,t) = xo(t)x(Z) (3.5) 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(erfc"1(2Z))2. (3.6) The temporal portion, Xo(^), was chosen to have the functional form: Xo(t) =cxe^ +c3. (3.7) Matching to the approximated behaviour of mean scalar dissipation along the stoichiometric interface for the flow field in question yielded constants ci w 450, c2 = 7 x 10~6. A collection of solutions for c3 = 0 with an oxidizer temperature of 1150k are shown in Figures 3.1 - 3.3. 1200 1100 1000 S 900 E a BOO 1 S. 700 E I— 600 SOO 400 300 a. ro I x 10 | — Flamelet #5 Flamelet #15 — Flamelet #20 Flamelet #25 i I 0.2 0.4 0.6 0.8 Mixture Fraction 0.4 0.6 Mixture Fraction Figure 3.1: Flamelet Solutions for Temperature and Enthalpy Source Term The 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. As 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. NON-PREMIXED METHANE AUTO-IGNITION 55 4 c3.5 2 3 u_ §2.5 I 2 '« 5 o1.5 a I 1 °0.5 0 — Flamelet #5 Flamelet #15 | — Flamelet #20 Flamelet #25 0.2 0.4 0.6 0.E Mixture Fraction 4 3.5 las S 2 2 o 1 0.5 0 x10~ — Flamelet #5 Flamelet #15 1 Flamelet #20 1 1 Flamelet #25 ! ! 1 \ 0.4 0.6 Mixture Fraction Figure 3.2: Flamelet Solutions for Species Mass Fraction - Yco and YCH3OH ignition, as do the minor species. The two species mass fractions considered in Figure 3.2 are used to help determine the appropriate combination of basis functions. The 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. The combination provides a means of distinguishing between flamelets in all stages of the ignition process; as well, the inclusion of CO was seen in Chapter 2 to be necessary to choose between flamelets in 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 in 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 YCH3OH prior to that point is clear, whilst the scale chosen to show the changes in temperature and Yco shortly following ignition obscures any changes before ignition. In the exact flow, although the mean scalar dissipation decays exponen tially with distance from the nozzle exit, one can expect unresolved excur sions to much higher local scalar dissipation rates. The effect of these can be substantial; Mason 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 with much higher scalar dissipation rates in the library. This was accom-CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 56 0 0.2 0.4 0.6 0.8 1 Temperature [K] 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 Natural Log of Enthalpy Source Term [ln(erg/g/s)] x10 0 0.2 0.4 0.6 0.8 CO Mass Fraction CH3OH Mass Fraction Figure 3.3: Flamelet Solutions for Species Mass Fraction - Contour Plots CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 57 plished by rebuilding the library of basis functions for additional values of the parameter c3 in Eq. 3.7. The final library used incorporated sub-libraries generated with c3 = 0,0.5,10,100 concatenated as illustrated for the tem perature solutions in Figure 3.4. The portion of the library with c3 = 100 showed approximately a 10% increase in the magnitude of the temperature source term over the portion with c3 = 0 for the pre-ignition solutions; the change was much more pronounced (exceeding of 100%) after ignition occurs in the flamelet calculation. Figure 3.4: Flamelet Temperature [K] Contours for Complete Library Here, we see a practical benefit of not enforcing dependence of an on n. At 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 in n space would discourage this type of combination, in contrast to the evolutionary stabilization used here which imposes no such constraint. CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 58 3.3 Results and Discussion 3.3.1 Predictions of Ignition Delay The predicted variation of the non-premixed ignition delay time with temper ature is shown in Figure 3.5 along with the experimental data from the UBC shock tube facility [78]. The results shown were produced using an identical procedure; only the initial temperature field in 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 0.6 „ 0.4 OT E -0.2 ~c ~ 0 -0.2 -0.4 -0.6 • UBC Shock Tube • -o- CSE-LFD Temperature Criteria - 2.23 • • •o Delay [ms] • Delay [ms] c • b Ipjnitk • > * a b Ipjnitk • a a , g 1 , 0.67 0.75 0.8 0.85 1000/T[1/k] 0.9 Figure 3.5: Effect of Ambient Temperature on Ignition Delay in the experimental data, there appears to be reasonable agreement between the simulation and the experimental results. As 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. As the object of this simu lation 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 in the experimental data. The results of the allocation of variance computations, and calculated R2 values are shown in Table 3.1. While the coefficient of determination is low for both the simulation and the regression line, the comparison indicates that the simulation predictions fit the experimental data nearly as well as a CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 59 regression line, which is heartening for a predictive simulation not tuned to match the data. 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. The laminar ig nition delay, which is extracted from the library of basis functions used for each case, is included here to facilitate later discussion. This laminar delay is obtained by averaging of the ignition delays for the various histories in 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. Whilst 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. This is not entirely unpromising; given the substantial scatter in the experimental data it is impossible to claim that the test casts significant doubt on the validity of the simulation. Additionally, interpretation of the experimental data is diffi cult due to many simultaneous or nearly simultaneous ignition sites. While it does appear that many of the experimentally observed ignition sites are farther downstream than the prediction, this could be explained by differ ences in the jet penetration at the time of injection, arising from errors in the momentum flux of the jet between the two due to the crude approximations made in 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 in no case did the simulation suggest that ignition occurred closer than Simulation Linear Prediction Regression SSE K2 0.288 0.223 0.52 0.63 (SSyy = 0.596, fff^ = 1.29 ) Table 3.1: Allocation of Variance for Ignition Delay 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] td (xign) td (lam) 1150 1.49 0.0411 2.39 1200 1.35 0.0373 1.42 1250 1.22 0.0336 0.92 1300 1.02 0.0360 0.62 1400 0.68 0.0289 0.27 Table 3.2: Location and Mixture Fraction at td 3.3.2 Relation to Prior Studies To put the results of the current work into context in the literature, Figure 3.6 compares the current results to experimental results published by Fraser et al. [63] and Ishiyama et al. [67] as well as previous simulation results. The earlier simulation results are those presented by Kim et al. [68] using conditional moment closure, and an earlier implementation of conditional source-term estimation reported by Blair [60]. An third previous simulation undertaken by Kim et al. [69] also attempted to predict Fraser et al.'s results using a flamelet model. Those results very nearly coincide with those in [68], and so are omitted from the figure for clarity. It would appear prima facie that both the experimental results used for validation in this work as well as the current simulation predictions devi ate from the results presented by Fraser in [63] (as well as the simulation results presented in [68, 69]), although the comparison to the data from Ishiyama et al. is much more promising. There are, of course, substantial differences in the physical configuration used to measure the delay; both Fraser et al. and Ishiyama et al. were injecting into a combustion bomb as opposed to a shock tube. As well, the injector nozzle diameter, pressure ratio, pressure of the oxidizer prior to injection, control characteristics, and pre-combustion composition of the oxidizer varied substantially. However, the disagreement between the two earlier experimental data sets, combined with the strong resemblance between Ishimyama et al.'s data and UBC's shock tube data suggest that these factors may not be solely responsible. The pressure of the oxidizer prior to injection was investigated by Fraser et al., and the dependence of the ignition delay on that pressure was found CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 61 + Fraser et al. -o- CSE-LFD • Kim et al. • UBC Shock Tube * Ishiyama et al. Blair *% r+ + ++ +++ ++ + 4.48 00 1.65^ >. Q 2 1.5 1 E I 0 S-0.5 -1 -1.5 -2 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1000/T[1/k] Figure 3.6: Effect of Ambient Temperature on Ignition Delay - Comparison Prior Works o *^ c cn 0.61 f 0.22 CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 62 to be minor, certainly insufficient to account for differences of this magni tude. 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 in the vessel to identify ignition. The more recent investigation discussed by Ishiyama et al. [67] reconsidered the question of determining ignition by pres sure measurement in 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 in the vessel is more appropriate. Before the pressure rises noticeably, a significant heat release and cor responding temperature increase must have occurred. This 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 with the amount of fuel mixed before ignition. In particular, by inspection of Figure 4 in [63], for longer ignition delays, the method of anal ysis used to detect ignition from the pressure rise in 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. The measure used by Ishiyama et al. is much more likely to coincide with the onset of ignition determined by optical measures at UBC, 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 in the bomb significantly. However, it is still expected to be longer than the delay detected at UBC as light emissions can be expected to be observed before the pressure responds. The earlier CMC and flamelet simulations are in reasonable agreement with Fraser et al.'s results, but not with the more recent experiments. As discussed above, the criteria used there for ignition is one which will identify ignition after the system has been burning for some significant amount of time. As such, these results are no more directly comparable to the current effort than Fraser et al.'s data. That these simulations agree reasonably well with the matching data does not necessarily guarantee that they would be able to identify the onset of ignition using the more stringent criteria used in 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 in the time at which ignition actually occurred. While it would be interesting to com pare 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 igni tion. Without 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 Blair [60] as well, Blair 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. This can only occur well after the igni tion prediction based on the current criteria. The mechanism by which Blair 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 in principle to that employed in this paper is unclear. There are, however, several dif ferences between the current method and Blair's implementation which may be responsible for the discrepancy. First, Blair used only temperature to perform the inversion for ~a; discussion elsewhere presents the difficulty dis tinguishing between flamelet basis functions on temperature alone. However, Blair did not appear to suffer from a uniqueness problem problem related to using only temperature; Figure 4.3 in [60] indicates that in the library of basis functions used, the conditional temperature appears to be changing substantially between flamelets. However, the lack of flamelets with similar temperatures highlights a larger issue: in the library used by Blair, there are at best two flamelets (one of which is the initial condition) present which have a peak temperature less than 1800K. This suggests that the ignition process is almost entirely unrepresented in the library of basis functions. As well, the lack of 'ignited' solutions with similar temperatures is the result of using only one scalar dissipation history, which was increasing in time as opposed to the expected decrease in time in the problem under consideration, to build the library. This suggests that it is quite possible that if an additional scalar were included, the minimization of the least squares problem necessary to de termine ~a may have had a very large residual indeed with this library, as it is unlikely that any combination of the solutions present accurately matches the physics in the flow. As a result, the physically unjustifiable smoothing imposed on ctn(n) to stabilize the inversion process would likely have had a CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 64 significant impact on the chosen solution. As any of these major issues could potentially explain the differences between the current approach and Blair's implementation, it is difficult to ascertain to what extent each is responsible. The effect of other potentially less substantial differences, such as a much finer grid in the current study to allow for more rapid variation of the flame structure in space and 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 an (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 in temperature and a slow build up in reaction rate, fol lowed by a rapid increase in temperature and reaction rate [72]. While early implementations of CSE+LFD in KIVA-3V produced a steady increase in temperature from the start of the simulation [60], the flavour presented in this paper follows the physical behaviour much more closely. The conditional averages of the three scalars used for the inversion process are shown in Fig ure 3.7 along with 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 in the conditional aver age of temperature prior to ignition, followed by the development of a local increase in temperature near the stoichiometric mixture fraction. Following ignition, the increase in temperature spreads outward from the stoichiometric mixture fraction. The lack of change in the conditional average for temper ature 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. The concentration of carbon monoxide also increases dramatically, but there is a more apparent increase in the concentration prior to ignition. The evolution of the conditional average of the intermediate CH3OH is more interesting. The peak conditional magnitude of this species is substantial well before ig nition takes place, and changes location in mixture fraction space as ignition CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 65 Figure 3.7: Evolution of Conditional Averages in Mixture Fraction Space, Tox = 1150 A" CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 66 occurs. Comparing to Figure 3.3, we see that there is a dramatic shift to richer mixture fractions for the residence of this species near ignition. This shift, as well as the non-negligible concentration of the scalar well prior to ignition, justifies the choice of this scalar to aid in selecting the appropri ate combination of flamelets prior to ignition. Finally, the behaviour of the CH creation rate is shown. This quantity, which is obtained from reactions forming CH only (as opposed to the net creation rate), was investigated in order to explore a quantity directly related to the optical emissions recorded in the UBC 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. The behaviour is promising for the detec tion of ignition; as with 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; in the UBC shock tube, the threshold was deter mined 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 maximum increase in the conditional average of the var ious scalars at any mixture fraction and on any plane is plotted against time for a selection of oxidizer temperatures. The various quantities are normal ized by a threshold chosen for the quantity, such that ignition was declared based on the various quantities when the quantity plotted exceeds a value of one. Although 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' in the various plots. With the thresholds chosen, there is very little difference in the ignition delay identified by consideration of any of the three scalars, as is apparent from Figure 3.9. Given this comparison, the ignition delay identi fied by an increase of 75 degrees in the conditional average of temperature at any mixture fraction, anywhere in the flow, was deemed to be well correlated with the optical criteria used in the experiments. The data presented earlier (Figures 3.6, 3.5) was generated using this criterion. In addition to consistency in the predicted ignition delay, the location along the centreline at which ignition occurred was within one plane regard less if the temperature or CH based criteria was used. The CH detection, as expected, displayed its maximum increase at slightly richer mixture frac tions (Z\{T\c,)max inc w 0.083) than the temperature criterion (Z\(T\Qmax inc — CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 67 g>10 S 6 J E 1 • a 2 1 2 — Temperature CH Production Rate CO Source Term j 1 l[ms] e E 'S 2 1 2 TOI = 1150k — Temperature CH Production Rate — CO Source Term /II /' ill jjy! • D 0.5 1 1.5 2 t[ms] Tox = 1200k Tm = 1300k 1400k Figure 3.8: Evolution of Peak Increase in Conditional Averages 0.5 r E %. 0 o Temperature Increase -* CH Production Rate CO Source Term -°0.7 0.75 0.8 0.85 1000/T[1/k] 0.9 Figure 3.9: Comparison of Results for Different Ignition Criteria 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 tur bulence on the delay before ignition is not immediately obvious. There are competing effects; increased turbulence levels will increase the rate of mixing. While this means that the fuel and oxidizer will mix sooner, it also means that any radicals produced or temperature increase will be dissipated quickly. Further, in a turbulent flow, there will be a distribution of strain rates; within any given region, there will be a variety of strain histories present [72]. The laminar and turbulent ignition delays from Table 3.2 are plotted against each other in Figure 3.10. We see that at high temperatures, where the delay is short, the effect of turbulence is to increase the ignition delay. At low tem peratures, 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 with 1000/T; a linear regression gives a coefficient of determi nation of R2 = 0.98. This effect can be related to the earlier framework by the following speculation: when the ignition delay is short, the effect of tur bulence 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. At longer ignition delays, the majority of the field is well mixed prior to ignition. Then, the distribution of strain rates suggests that there is a strong likelihood that one of the well mixed regions will encounter low strain, resulting in ignition promotion. 3.4 Conclusion A simulation has been undertaken to predict the ignition delay of non premixed methane injected into hot air under high pressure at high velocity. The technique used to close the source term in the RANS paradigm takes advantage of concepts from both the CMC and flamelet models. A separate conditional average is permitted (and encouraged) on each plane of cells ax-ially equidistant from the nozzle. This 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 METHANE AUTO-IGNITION 69 0.75 0.8 0.85 1000/T[1/k] 0.9 Figure 3.10: Ratio Between Turbulent and Laminar Ignition Delays CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 70 of basis functions will be necessary. The predicted ignition delay was in excellent agreement with the experimental data at a variety of oxidizer tem peratures. Despite the significant scatter in the experimental data used for direct validation, the agreement in trends with an additional experimental data set obtained from an independent source lends to significant confidence in the results. Scatter in experimental measurements of the spatial location of the ignition event, as well as the occurrence of multiple, nearly simultane ous 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 sim ulation results. In order to produce simulations valid beyond ignition, it is necessary to develop libraries of basis functions which encompass the appro priate 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. Blair. Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of British Columbia, 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, July 1999. [62] W. K. Bushe and H. Steiner. Laminar flamelet decomposition for condi tional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [63] R. A. Fraser, D. L. Siebers, and C. F. Edwards. Autoignition of methane and natural gas in a simulated diesel environment. SAE Technical Paper Series, 1991. [64] A. C. Hindmarsh. ODEPACK, A systematized collection of ode solvers. In R. Stepleman et al., editors, Scientific Computing, pages 55-64, Am sterdam, 1983. North-Holland. [65] J. Huang, P. Hill, W. Bushe, and S. Munshi. Shock-tube study of methane ignition under engine-relevant conditions: experiments and modeling. Combustion and Flame, 136(1-2):25-42, 2004. [66] J. Iaconis. An investigation of methane autoignition behaviour under diesel engine-relevant conditions. Master's thesis, University of British Columbia, 2003. [67] T. Ishiyama, M. Shioji, et al. Characteristics of spontaneous igni tion and combustion in unsteady high-speed gaseous fuel jets. In 2003 JSAE/SAE International Spring Fuels & Libricants Meeting, number JSAE 20030225 / SAE 2003-01-1922, Yokohama, Japan, May 19-22 2003. [68] S. H. Kim, K. Y. Huh, and R. A. Fraser. Modeling autoignition of a turbulent methane jet by the conditional moment closure model. Pro ceedings of the Combustion Institute, 28:185-191, 2000. REFERENCES 72 [69] S.-K. Kim, Y. Yu, et al. Numerical investigation of the autoignition of turbulent gaseous jets in a high-pressure environment using the multiple-rif model. Fuel, 83:375-386, 2004. [70] A. Klimenko and R. Bilger. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [71] S. D. Mason, J. H. Chen, and H. G. Im. Effects of unsteady scalar dissi pation rate on ignition of non-premixed hydrogen/air mixtures in coun-terflow. Procedings of the Combustion Institute, 29:1629-1636, 2003. [72] E. Mastorakos, T. Baritaud, and T. Poinsot. Numerical simulations of autoignition in 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. Laminar diffusion flamelet models in non-premixed turbulent combustion. Progress in Energy and Combustion Science, 10:319-339, 1984. [75] G. Ruetsch, L. Vervisch, and A. Linan. 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 with conditional source-term estimation. Physics of Fluids, 13(3): 754-769, March 2001. [77] G. Sullivan. Personal communication. 2004. [78] G. Sullivan. Personal communication. May 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 cap turing the physics of non-premixed turbulent combustion proposed by Bushe and Steiner [82] have been discussed. Chapter 2 presented the first applica tion of the method to a reacting flow with a steady ensemble averaged solu tion. While not entirely successful, that a converged solution was obtained sets the framework for future development. The second paper, presented in Chapter 3, discusses a successful simulation of non-premixed methane au toignition, and demonstrates that where the basic requirements of the model are satisfied its performance can be excellent. When accounting for the effects of chemical reaction with conditional source-term estimation using laminar flamelet decomposition, the effects of arbitrarily complex chemistry can be accounted for with minimal computa tional cost. This is because the majority of the computational effort associ ated with 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 genera tion can be easily performed in parallel. When 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 with conditional moment closure using a single conditioning variable, such as extinction and ignition, can be captured in an unsteady simulation, al though there are yet to be resolved complications associated with capturing unsteady effects in a steady solver. In general, the model has broad applica bility and should be useful in any circumstance where the necessary physics can be captured in 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. When 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 both the conditional averages CHAPTER 4. DISCUSSION AND CONCLUSION 75 (even in 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. As 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 — in 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. The first is the approximately Gaussian variation about the conditional mean with small standard deviation due to small local fluctuations in the strain and other effects not related to mixture fraction. The second is the variation observed when the local fluctuations in 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. All first moment closures must assume that the former variation is small; second moment closures have been developed to address larger variations [83, 84]. The latter form of variation is more problematic, and causes difficulty for second moment CMC methods as well. To see why the latter form is so detrimental to first moment closures, suppose that within a given region, conditionally averaging individual points obtained from a fully resolved sample produces a conditional mean of tem perature coinciding with the dashed line on the left of Figure 4.1. Further suppose that approximately half the points are clustered about the bottom curve (solution A), and the remainder are clustered near the top curve (so lution C). This 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. The solid curves in the fig ures 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- The dashed lines are a conditional average of Solutions A and C. Although 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 §0.02 — Solution A Solution 6 — Solution C <A,C|(> 0.2 0.4 0.6 Mixture Fraction 0.2 0.4 0.6 o.e Mixture Fraction Figure 4.1: Flamelet Solutions for T and Yco', {A,C\Q indicates conditional average of Solutions A and C situation might exist — e.g., one which falls in the flamelet regime, and in which approximately half of the realizations are close to extinction and half are not. A first order CMC 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 CMC closure, as well as second order closures which presume a Gaussian distribution about the conditional mean would be highly inaccurate. In CSE+LFD, faced with this situation, we have to choose the coefficients aA,<*B,ctc- Choosing = 0,QB = 1.0, ac = 0 would agree well with the temperature data. If the inversion was performed based on temperature alone, as proposed in [82] and implemented in [80], this would be the expected result. The conditional source term used would then be that calculated based on {T\QB, (Yi\QB. This would obviously be error prone — worse in fact than the first order CMC 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 with both the temperature and carbon monoxide fields is to choose = 0.5, OB = 0.0, ac = 0.5. Then, the source term used will be an average of the source term calculated based on (T\QA , (^IOA and (T\QC , (Fj|0C, which is much closer to the ideal behaviour. With 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 CHAPTER 4. DISCUSSION AND CONCLUSION 77 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 with 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 solu tions representing all possible physical conditions in the library used for de composition, 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 un 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 actu ally 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. This statement is similar to Eq. 1.38, and can be written as: min{||/i|| + ||/2|| + ||/s|| + --- + ||//+i||} (4.1) where 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 in basis function space; we are then trying to minimize a functional of the coefficient function and the scalar functions: min {||-?||} ; ^ = (a(n) (T|0 (n), a(n) (Yi|C>, • • •, «(") <Xi\0 ) (4-2) Considering Eq. 4.2, it is apparent that the solution of Eq. 4.1 for a(n*) is actually a functional dependent on (T\() (ri*), (Yi|C) (n*), • • •, (Yi\0 (ri*). This is clearly a difficult situation; it might at first appear that inversion in this pan-dimensional space is hopeless. Fortunately, we can make as sumptions that simplify the problem. First, for a practical approximation the basis function space ri* is not infinite. Only the physics relevant to the problem at hand must be included; for example, in unbounded flows wall ef fects need not be included. Our 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 in n*. Lastly, by inspec tion of the solutions generated for the discreet subset of n*, the functional is frequently homogeneous in many of the dimensions, and there are only a limited number of approximately orthogonal coordinates in 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. The 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. This is a classic example, where the physics necessary to describe this process were not present in the library. Other physics of relevance to the simulation of the Sandia Flame 'D', as well as continuing the autoignition simulations beyond ignition include edge flames [85], which can never be captured in a library built from simple solutions to the laminar flamelet equations. This 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 will 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 CMC, or DNS for that matter. In general, CSE with decomposition can be viewed as a way to combine easily obtainable elementary solutions to build a more realistic simulation. As 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 Eq. 1.29; that is: (4.3) CHAPTER 4. DISCUSSION AND CONCLUSION 79 Since the pdf 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. As mentioned in [82], however, this is a common assumption not specific to the CSE methods, and it would be trivial to adjust the formulation to utilize more complex representations of the pdf. Comparing Eq. 4.3 (written for temperature) with the corresponding clo sure using double-condition CMC (using scalar dissipation, with sample space X1 as the second conditioning variable): w« / I <w|C,X>(C,xWx'« f f (viCx'jpiOpMdCdx' (4.4) Jc Jx> Jc JX' Bushe and Steiner [82] noted that coefficients an 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°r 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 an is more appropriately likened to p(a')da', where a is an arbitrary scalar which must satisfy only the condition of being independent of Z, so that it captures variations in the flame structure not correlated with mixture fraction. This poses a minor conceptual difficulty, as it implies that we are in 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. This could be a problem where the flame structure varies rapidly with 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 pdf of x 1S approximately linear; however, the reaction rate varies in a non-linear way and the product (u\x')p(x') seen m the integral of Eq. 4.3 must be highly non-linear, so the constraint may be arduous. In Chapter 2, it was mentioned that when flamelets with high scalar dissi pation were included in the library the solution became unstable, and it was suggested that this was due to the transient nature of these flamelets being included in a steady solution. Another possibility is now apparent: the errors arising from the effective discretization in scalar dissipation sample space may be large. As the flamelets used for basis functions were extracted at constant time intervals, and the structure was changing more rapidly in flamelet time near extinction, the effective discretization became much coarser when the CHAPTER 4. DISCUSSION AND CONCLUSION 80 extinguishing flamelets were included. Although 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. The drawback, of course, is that the library becomes much larger, the inversion process more time-consuming, and the uniqueness of the solution suffers. Given that in the flamelets used for the decomposition in 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 in 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 E («n(Xo)n) _n=l It is not yet clear if this additional constraint would produce an improve ment — 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 Eq. 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 pre vious inversion, we are making a physical assumption about the expected behaviour. Although in general the penalty for deviating from the previ ous solution is set only just high enough to ensure a unique solution, the penalty exists nevertheless, so the solution will always incorporate the a pri ori 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 X> /%(C) Jo /(CR (4-6) 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 in the Lagrangian sense. This claim would treat the flame structure as a convected entity, which is not entirely unreasonable. The difficulty is found in 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 solu tion in the Lagrangian sense would have to produce a solution equal to the past solution in the Eulerian sense for a problem with an Eulerian steady solution. An elegant method of incorporating influence from the previous structure in the Lagrangian sense was not explored in this work, but may be an interesting avenue for future development. A revised formulation of conditional source-term estimation with lami nar flamelet decomposition has been presented and used successfully in the simulation of turbulent non-premixed methane ignition. The framework has also been used to develop a simulation of the Sandia Flame 'D'. This for mulation differs from previous implementations in several ways. Transport equations are solved for multiple species which are used in the inversion pro cess to determine the flamelet structure; using multiple scalars allows the method to discriminate between the various physical conditions represented in the library of basis functions. The library of basis functions used has been extended to incorporate more varied solutions, in 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 in the flow field. A slightly revised interpretation for the method was also suggested. The successful prediction of the temperature dependence of the ignition de lay for an igniting methane jet, as well as obtaining a steady solution for the Sandia Flame 'D' in the RANS paradigm effectively demonstrates the potential of the model, and provides a framework for future development. REFERENCES 82 REFERENCES [80] C. Blair. Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of British Columbia, 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, July 1999. [82] W. K. Bushe and H. Steiner. Laminar flamelet decomposition for condi tional source-term estimation. Physics of Fluids, 15(6):1564-1575, June 2003. [83] S. H. Kim, K. Y. Huh, and R. W. Bilger. Second-order conditional moment closure modeling of local extinction and reignition in turbulent non-premixed hydrocarbon flames. Proceedings of the Combustion Insti tute, 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. Linan. 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 RANS 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 cap tured experimentally. Further, these contours are ensemble averages, which need not necessarily have any direct relation to any specific realization. RANS Temperature Field, 1150K Oxidizer RANS CH Creation Field, 1150K Oxidizer 0.04 0.05 400 500 600 700 800 900 1000 1100 RANS Temperature Field. 12O0K Oxidizer 0.04 0.05 400 500 600 700 8O0 900 1000 1100 1200 RANS Temperature FieW, 1300K Oxidizer 400 500 600 700 800 900 1000 1100 1200 1300 RANS Temperature Field, 1400K Oxidizer 0.015 _ 0.01 *0.005 0 0 05 1 1.5 2 25 3 3JS 4 4.5 RANS CH Creation Field, 1200K Oxidizer 2 3 4 5 RANS CH Creation Field, 1400K Oxidizer 400 600 1000 12O0 1400 Figure A.l: RANS Temperature and CH Creation Contours at td
- 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 | 2004 |
Date Issued | 2009-11-24T20:35:25Z |
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 |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
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
- ubc_2004-0470.pdf [ 8.43MB ]
- Metadata
- JSON: 1.0080709.json
- JSON-LD: 1.0080709+ld.json
- RDF/XML (Pretty): 1.0080709.xml
- RDF/JSON: 1.0080709+rdf.json
- Turtle: 1.0080709+rdf-turtle.txt
- N-Triples: 1.0080709+rdf-ntriples.txt
- Original Record: 1.0080709 +original-record.json
- Full Text
- 1.0080709.txt
- Citation
- 1.0080709.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 13 | 18 |
China | 7 | 4 |
France | 5 | 2 |
United Kingdom | 4 | 8 |
India | 4 | 1 |
Republic of Korea | 3 | 0 |
Germany | 2 | 1 |
Sweden | 1 | 0 |
Iran | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 12 | 4 |
Ashburn | 8 | 0 |
Guwahati | 4 | 0 |
Beijing | 3 | 1 |
Oxford | 3 | 3 |
Shenzhen | 2 | 2 |
Shanghai | 2 | 0 |
Los Angeles | 2 | 0 |
Mountain View | 1 | 0 |
Stockholm | 1 | 0 |
Golden | 1 | 0 |
Southampton | 1 | 3 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080709/manifest