- 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 2004
pdf
Page Metadata
Item Metadata
Title | Combustion modeling with conditional source-term estimation and laminar flamelet decomposition |
Creator |
Grout, Ray W. S. |
Date Created | 2009-11-24 |
Date Issued | 2009-11-24 |
Date | 2004 |
Description | Conditional source-term estimation with laminar flamelet decomposition has been utilized to model the mean chemical source term in a predictive RANS simulation of two different problems. With this model, integral equations relating the unconditional mean temperature and species fields to the conditional means are used to determine the flame structure as a combination of basis functions formed from solutions to the unsteady laminar flamelet equations. In the simulation of the Sandia Flame 'D', a well studied co-flowing piloted jet flame with a steady average solution, a converged solution was obtained which captured the trends in the temperature and major species, although the nitric oxide prediction overestimated the peak concentration by a substantial margin. Simulation of the autoignition process of non-premixed methane produced simulation predictions in excellent agreement with the experimental data. Using a new, more stringent, criteria to define ignition than earlier studies, the effect of ambient temperature on the ignition delay was captured, as was the expected physical behaviour prior to ignition. The errors in the simulation of the co-flowing piloted jet can be attributed to a large degree to the lack of ability of the basis functions used to account for the effect of the pilot flame. That the autoignition simulation was much more successful highlights the importance of including all relevant physics in the library of basis functions. Earlier formulations of the model were extended by expanding the sample space for the basis functions to include a wide variety of solutions to the laminar flamelet equations; in order to distinguish between the larger set, a new method of stabilizing the numerics was found to be necessary and the number of scalars used to determine the flame structure was increased. The formulation is easily extensible to libraries of basis functions which are capable of including effects such as pilot and edge flames not captured by the laminar flamelet solutions. |
Extent | 8842098 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-11-24 |
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 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 2004-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/15580 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0080709/source |
Download
- Media
- ubc_2004-0470.pdf [ 8.43MB ]
- ubc_2004-0470.pdf
- 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
- Citation
- 1.0080709.ris
Full Text
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 S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF M A S T E R OF A P P L I E D SCIENCE in The Faculty of Graduate Studies (Department of Mechanical Engineering) We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A 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 Condit ional source-term estimation wi th laminar flamelet decomposition has been uti l ized to model the mean chemical source term in a predictive R A N S simulation of two different problems. W i t h this model, integral equations relating the uncondit ional mean temperature and species fields to the condi- t ional 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 wi th a steady average solution, a converged solution was obtained which captured the trends in the temperature and major species, although the nitr ic 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 wi th the ex- perimental data. Using a new, more stringent, criteria to define ignit ion than earlier studies, the effect of ambient temperature on the ignit ion delay was captured, as was the expected physical behaviour prior to ignit ion. The er- rors in the simulation of the co-flowing piloted jet can be attr ibuted to a large degree to the lack of abil i ty 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 al l relevant physics in the l ibrary of basis functions. Earl ier 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 stabil izing 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. i i i CONTENTS Abstract ii Contents i i i List of Tables v List of Figures v i Nomenclature v i i i Acknowledgements x Introduction xi 1 Background 1 1.1 Turbulent Combustion 2 1.2 Flamelet Models 3 1.3 Condi t ional Moment Closure 12 1.4 Condit ional Source-Term Est imat ion 17 References 25 2 Simulation of the Sandia Flame 'D' 30 2.1 Introduction 31 2.2 Problem Description 32 2.3 Background 32 2.4 Formulat ion 35 2.4.1 Chemical Source Term 35 2.4.2 Flow Solution 36 2.5 Results and Discussion 37 2.6 Conclusion 41 References 44 CONTENTS iv 3 Non-Premixed Methane Auto-Ignition 46 3.1 Introduction 47 3.2 Formulat ion 48 3.2.1 C S E Combustion Model 48 3.2.2 F low F ie ld 52 3.2.3 Basis Functions 53 3.3 Results and Discussion 58 3.3.1 Predict ions of Ignit ion Delay 58 3.3.2 Relat ion to Pr ior Studies 60 3.3.3 Behaviour Pr ior to Ignition and Identification of Ignition 64 3.3.4 Effect of Turbulence on Ignition Delay 68 3.4 Conclusion 68 References 71 4 Discussion and Conclusion 73 References 82 A RANS Fields 83 \ V LIST OF TABLES 3.1 Al locat ion of Variance for Ignition Delay 59 3.2 Locat ion and Mixture Fraction at tj, 60 v i LIST OF FIGURES 1.1 Borghi Diagram 8 1.2 Scatter Plot of Temperature vs. Mix ture Fract ion for Sandia Flame ' D \ x/D = 30 [1] 13 1.3 Relat ion Between Temperature and Mixture Fract ion, Log of Scalar Dissipation. Port ion of Figures 1 and 2 in [9] 17 2.1 Selected Basis Funct ion Solutions 34 2.2 Predicted Temperature F ie ld 38 2.3 C S E / L F D Val idat ion - Comparison of Condi t ional Averages to Experimental Results (Solid lines predictions, symbols ex- perimental measurements) 39 2.4 C S E / L F D Val idat ion - Comparison of Ax ia l and Rad ia l P ro - files to Experimental Results (Solid lines predictions, symbols experimental measurements) 42 2.5 Ni t r ic Oxide Predict ion 43 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 L ibrary . . . 57 3.5 Effect of Ambient Temperature on Ignition Delay 58 3.6 Effect of Ambient Temperature on Ignition Delay - Compar i - son Pr ior Works 61 3.7 Evolut ion of Condit ional Averages in Mixture Fract ion Space, Tox = U50K 65 3.8 Evolut ion of Peak Increase in Condit ional Averages 67 3.9 Comparison of Results for Different Ignition Cr i ter ia 67 3.10 Rat io Between Turbulent and Laminar Ignition Delays . . . . 69 4.1 Flamelet Solutions for T and Yco 76 LIST OF FIGURES A . l R A N S Temperature and C H Creation Contours at td 84 V l l l NOMENCLATURE Symbols ~u Velocity ~x Spat ial Coordinate T Temperature Y Species Mass Fraction Z Mix ture Fract ion £ Sample space for mixture fraction Z"2 Mixture Fraction Fluctuations (Variance) Le Lewis Number k Thermal Conduct iv i ty D Diffusivity Cp Heat capacity Ka Kar lov i tz Number Da Damkohler Number tk Kolmogorov time scale tc Chemical t ime scale tj Integral t ime scale tref Non-premixed reference time scale lref Non-premixed reference length scale uTef Non-premixed reference velocity scale v K inemat ic Viscosity e Dissipation rate of turbulent kinetic energy p Density cjj Source term for species i X Scalar Dissipation Rate a Arb i t rary Scalar an Coefficient for flamelet n ~a Vector of coefficients (ai,...,aN) ipjtk Integral of p(C') across discretization interval Notation ~v Vector M Ma t r i x (T\Q Condit ional average of quantity to left of |, conditioned on condit ioning 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 minimizat ion NOMENCLATURE ix Subscripts i Species index n Basis function index j Computat ional cell (physical space) k Fin i te difference point (mixture fraction space) Abbreviations Pdf Probabi l i ty density function C M C Condit ional Moment Closure C S E Condit ional Source-term Est imat ion R A N S Reynolds Averaged Navier-Stokes L E S Large Eddy Simulat ion D N S Direct Numerical Simulat ion X ACKNOWLEDGEMENTS I would like to thank the many people who supported and encouraged me throughout this project. Fi rst , I am grateful to my supervisor, Dr . W . Kenda l Bushe, who not only guided my research efforts, but shortly after we met convinced me that the best way to describe turbulent combustion simulation is "S ix pounds of cool in a five pound sack". A lso, 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 Co l in , Chad , and Jorge, who were k ind enough to proofread this thesis. I am also indebted to Col in for fighting the good fight wi th K I V A and sparing me the trials of adding transport equations to that particular code. A s well, I would perhaps not have undertaken this project had it not been for the encouragement of Dave Mumford and Dr. H i l l at Westport Innovations while I was an engineering undergrad. The financial assistance of Westport Innovations, the Natura l Science and Engineering Research Counci l of Canada, and the B C Advanced Systems Institute was gratefully appreciated. Final ly, 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 Ju ly 2004 INTRODUCTION INTRODUCTION x i i 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 Nor th Amer ica, or remember the plight of Cal i fornia 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 E P A and European acceptable l imits for automotive emissions, or broad base policy initiatives such as the Kyo to Accord on greenhouse gas reduction. A t the same time, the demand for affordable energy in developing economies has stretched the world's energy resources. O n al l fronts, it is obvious that we must improve our pr imary means of energy conversion — controlled combustion — in terms of both efficiency and minimizing the release of harmful emissions. Real iz ing 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 wi th a thorough understanding of the physics at work and an effective means of evaluating design alternatives at one's disposal. Experimentat ion, valuable for confirming that the final design performs as expected, is l imited as a useful design tool. Whi le it is possible to bui ld a prototype of, for example, a proposed combustion environment's geometry, obtaining the precise tem- perature distr ibution throughout that environment during the combustion event is unfeasible experimentally. Numerical simulations, based on models developed and validated rigorously using fundamental experiments, hold the possibil ity to economically meet the need for detailed information. Development of adequate models involves incorporating the effects of a mult i tude 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 t ightly coupled. The macro-mixing behaviour of the flow can be treated using conventional computational fluid dynam- ics ( C F D ) techniques. Chemical kinetic mechanisms form an entire field of ongoing research; high quality detailed mechanisms are available only for some fuel/oxidizer combinations, tuned by comparison to a wide variety of experimental data. Micro-mixing phenomena, and the resulting turbulence- chemistry interactions, are not directly captured by a C F D solution for en- semble or spatial averages. Since direct numerical simulation (DNS) (which INTRODUCTION xi i i does capture these effects) is beyond current computing l imits 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 wi th l itt le adjustment to two diverse problems using commercial C F D tools commonly employed by industry. The first, an extended version of that presented at the 2004 Canadian Section of the Combust ion Institute Spring Technical Meet ing 1 , addresses the simulation of an atmospheric pres- sure co-flowing piloted methane jet. The second considers the auto-ignit ion 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 il lustrate 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, wi l l 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 product ion/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. Typ ica l 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- t ion 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 l imitations. The typical approach is to not ful ly resolve the flow, and provide a model for the unresolved effects. However, the methods typical ly 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 ( R A N S ) 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. Whi le 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 wi l l dominate the result. Assuming that the unresolved variation for al l dependent variables of in- terest is correlated, it is possible to relate the unresolved variat ion in al l scalars to a representative scalar. It is then necessary to know the proba- bi l i ty 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 ( L F M ) and Condit ional Moment Closure ( C M C ) . 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 wi th 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 wri t ten 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 probabil i ty density function for the mixture fraction is commonly assumed to be a Beta pdf. W i t h such a presumption, only the first two moments need to be estimated from the flow. Laminar Flamelet models and Condit ional Moment Closure, while sharing the use of a conserved scalar, differ dramatical ly 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 wi th minor computational cost. However, the range of applicabil i ty is l imited, and many of the subtleties of the combustion process are not accounted for. Condit ional Moment Closure methods are much more robust in their range of application, have a rigorously justifiable derivation, and are much more computationally expensive, part icularly when realistic chemistry is employed. Also, the necessary modeling of many unclosed terms in the resulting equations introduces modeling assumptions of questionable accuracy. A s an alternative to these two presumed-pdf / conserved scalar models, a class of models known as ' p d f models exist, for which the evolution of the joint probabil i ty density function (pdf) of al l scalars and velocities is solved [36]. In such a formulation, the chemical terms are closed exactly — it is the turbulent transport and molecular mixing which must be modeled. A s there is l i t t le from this model which is essential to the body of this thesis, it wi l l not be further discussed. The balance of this chapter wi l l first review flamelet modeling and C M C , then finish wi th the development of the condit ional 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. B y 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 val id where it is justifiable to make the assumption of the laminar structure. St i l l , flamelet modeling is conceptually attractive and serves to provide an introduction to the use of a conserved scalar and presumed pdf in combustion modeling. A s well, the flavour of conditional source-term estimation explored 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 natural ly a function of (~x,t). A l though 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- t ion is simply that the mixture fraction is transported as a conserved scalar wi th unity Lewis number [4, 44]; that is, it satisfies E q . 1.5 wi th Co = 0 and D = -J^r. Boundary conditions are Z = 0 in the oxidizer stream and Z = 1 in the fuel stream. Al though some would argue that this definition lacks a clear physical interpretation, there is l itt le point in debating the issue. F i rs t , for the purposes of the flamelet models (and, as we shall see later, C M C ) , the physical interpretation is irrelevant beyond the existence of a strong cor- relation between the mixture fraction and the other dependent variables [34]. Second, although it wi l l 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- al ly unity Lewis number. A n extended formulation incorporating the effects of differential diffusion was evaluated by Pi tsch [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 wi l l consider only equal diffusivities. Not ing 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 strict ly val id 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 ( implicit ly normal to the stoichiometric interface), Z2 = X2, Z3 = X3, T = t, and develops the following transformation rules: d_=d_ dZd_ dt~ d r + dtdZ { ' ] d 9 2 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 wi th respect to the higher dimensions (zT2,zT3). The first observation one may make is that the result, E q . 1.6, has not been to transform the species transport equation from a function of (~x,t) to one which is a function of (Z, t), but rather (Z, Z%, Z3, r ) . Th is is obviously a mathematical necessity — we be- gan wi th a 4-dimensional coordinate system, so certainly 4 new orthogonal CHAPTER 1. BACKGROUND 6 coordinates are necessary for a val id transformation. However, claiming that the variation along the coordinates normal to the gradient of mixture frac- t ion 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. Th is 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. Ear ly formulations uti l ize the form obtained by Peters ([29]) from an analytical solution for the spatial variat ion of mixture fraction in an infinite 1-D mixing layer, wi th boundary conditions Z(y —»• —oo) = 0, Z(y —> oo) = 1. The resulting functional form is: - 2 [ e r f c " 1 ( 2 Z ) ] 2 X*(Z) = —r - j . (1.7) ^ ' e - 2 [ e r f c - l ( 2 Z 0 ) ] 2 V ' Alternatively, from an analytical solution to a symmetric one-dimensional mix ing layer where the fuel originates from the plane of symmetry, P i tsch et al . [34] obtain: * • < * > - f £ i ™ Despite the observation by K i m et al . [19] that the latter form is incon- sistent wi th the transport equation for the pdf of Z , it now commonly used, e.g., in [32]. Other forms for this relationship are typical ly used in C M C formulations, and wi l l be discussed later. Once such a functional form has been chosen, xo is the only remaining coupling to the flow field. Final ly, it is important to note that the Euler ian time has been replaced by a Lagrangian 'flamelet t ime', r. A t this point, flamelet models spread out into several different direc- tions wi th 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 val id, the smallest scale of the turbulent structure, the Kolmogorov length, must be larger than the thickness of the reaction zone. A n alternative state- ment is that the chemical reactions must occur on a t ime 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 Lent ini [23]. Al though laminar burning velocities are reasonable for premixed combustion, for non-premixed combustion Lent ini proposes using revised definitions for the reference scales; specifically: Iref = \/vtc Uref = yj^ (1.9) Constant values of two dimensionless numbers, the Kar lov i tz number (Ka) and the Damkohler number (Da) appear as straight lines on such a plot. The Kar lov i tz number is defined as the ratio of chemical t ime 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 t ime scale [46]: Da^. tc Using Lentini 's revised definitions of the non-premixed length and t ime scales, tc in the above definitions is replaced by trer = defined as in E q . 1.9. In keeping wi th the earlier criteria, strictly speaking the validity of the flamelet models require tc < tk, or Ka < 1. This region is the port ion 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. Lent ini evaluated the chemical t ime scale and compared it to the flow time scale during flamelet computations for a turbulent burner. W i th in the flamelet regime, he obtained reasonable results. However, just outside the flamelet regime in the perturbed flamelet regime Lent ini 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 wi th caution due to a l imited regime of applicability. W i t h 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 E q . 1.6. Setting ^ = 0 and neglecting R(Yi) leads to a set of ordinary differential equations in mixture fraction space for the component species mass fractions and (including the corresponding equation obtained from a transformation of the energy equation) temperature, in which the scalar dissipation (xo) appears as a parameter. Solutions to the system rep- resent a balance between mixing (associated wi th 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 Wa l l 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 typical ly 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), Lent ini [23] assigns the value a = 2, while L iew et al . [24] develops an algebraic relationship containing the Kolmogorov and integral scales, as well as various constants (some of which depend on the flow geometry). In the commercial C F D code Fluent [13], the 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. A s pointed out by Barths et al . [3], the history of flamelet solutions is relevant — it is as though the solutions have a 'memory' of their prior experiences, which, although not infinitely long, invalidates the steady state assumption for some cases. There has been some discussion regarding the appropriate criteria on which to assess the val idity 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 ignit ion and extinction — for which an unsteady model is essential. Unsteady flamelet models have the feature of retaining the t ime derivative in E q . 1.6 and attempting to relate the Lagrangian r to an Euler ian t ime/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. A t each point in space, the port ion 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 joint ly by Barths, Pi tsch, and Peters [3] as the 'Representative Interactive Flamelet ' (R IF) , a good description of what is now known as the 'Euler ian Part ic le Flamelet Mode l ' 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 t ime 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 mi ld combustion burners in furnaces [10] wi th reasonable success at predicting temperature and major species fields. Whi le the majority of the use of flamelet models is in R A N S simula- tions, the model has also been applied to L E S computations. Cook et al . [11] describe an application of steady flamelets in L E S , 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. A n implementation making use of an unsteady formulation is given by Pi tsch [33], Pi tsch and Steiner [35]. In the latter work, the formu- lat ion and solution procedure closely resembles that commonly employed for the modeled C M C equations which wi l l 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, l imit the applicabil i ty of the method. Corrections introduced to arrive at equations similar to the C M C counterparts do produce improve- ment, but it is unclear whether any addit ional unclosed terms may have been overlooked. CHAPTER 1. BACKGROUND 12 1.3 Conditional Moment Closure Condit ional Moment Closure ( C M C ) 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 Kl imenko and unified in [22], the reward for reduced fluctuations about the condit ional mean quantities is that simple first order closure for the chem- ical source term is often sufficient. The dramatically different derivations suggested by Kl imenko [21] (based on pdf methods) and Bilger [5] (applying a decomposition to conventional transport equations) producing the same basic equations lends substantial credibil ity to the method. Typical ly, 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 F lame ' D ' [1] at 30 diameters downstream from the nozzle; this flame wi l l be described in detail when it is considered in Chapter 2. A t the axial location shown, only moderate to low levels of extinction are expected. It is apparent that there is a significant correlation between temperature and the mixture fraction, and the variat ion about the conditional mean (shown in black) when mixture fraction is used as a conditioning variable wi l l be much smaller than that about the uncondit ional mean. It is worthwhile to briefly consider the origin of the C M C equations, so that the modeling assumptions are clearly visible. To derive the C M C equations using the decomposition approach, Bilger makes use of the decom- posit ion: Y t) = Q(Z(1t, t), t, t) + Y"Cg, y). (1.11) Y" is the condit ional fluctuation, and Q = (Y\Z = £); we wi l l persist in the use of Q in this section only for brevity and consistency wi th Bilger's notation. The quantity {Y\Z = C) should be interpreted as the condit ional 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 E q . 1.5, E q . 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 ' V r = V Q + ^ V Z + VY" (1.13) V • (pDVY) = V • (pDVQ) + | | v • (pDVZ) + pD(VZfd^ + PDVZ.V^ ( I - " ) + V • (pDVY") Substitut ing the above into E q . 1.5, using that Z also satisfies E q . 1.5 wi th LJ = 0 to simplify, and taking the conditional average, Bilger arrives at the 'unclosed' C M C evolution equation for species concentration: <p|C> ?r + (p\0 0*10 • V Q - <P|C> (xlO ^ Q where: ' \ i -s / * ~e xrib/ \/vi\>/ _Qr2 = (p\Q(u\0 + eQ + eY (1.15) e Q = ( 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• « P I C > ( u " Y " \ Q P ( 0 ) . . v• « P I C > P ( 0 ) , „ „ Where a gradient transport model has been used to close the term ( u "Y"\Q, similar to that used 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 wi th a linear approxi- mation for the condit ional average of velocity, have been broadly adopted, e.g.,in [12, 16, 28]. The final closure approximation is the functional form for the condit ional scalar dissipation. K i m et al . [16, 18] and Devaud and Bray [12] al l make reference to an expression obtained by G i r ima j i (1992) by consideration of the pdf transport equation. K i m and Huh [16] also make reference to the counterflow mixing model used in flamelet calculations. In their work on spray autoignit ion, K i m 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 wi th solving E q . 1.15 using the de- scribed approaches to obtain closure is significant. However, the condit ional 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! Th is has been done, e.g., by Swaminathan and Bilger [42], and is known as 0-dimensional C M C . In a homogeneous field, the conditional transport terms vanish, and the C M C species evolution equation reduces to a form reminiscent of the flamelet equa- tions: d~i + (xlO = WO (i-w) Although this form of the C M C equations bears resemblance to the flamelet equations, there is an important difference in proceeding from the current state. Whi le the assumptions made early in the derivation of the flamelet models disregard much of the physical behaviour, the rigorously justifiable derivation of the C M C equations, up to the point where closure assumptions are made, requires only the fundamental assumption that the variations in the various scalars are well correlated wi th the condit ioning variable. This means that the framework to reconsider or improve the clo- sure for a given term is available. A l l of the necessary physics are present, and we are consciously aware of the assumptions made to neglect or model various terms. Th is makes the model ideally suited to improvements not only in modeling, but in basic understanding as well. Whi le quite successful in some circumstances, variations about the con- dit ional 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 condit ional 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 K i m et al. [17], Mastorakos and Bilger [27]. However, this requires modeling of several unclosed terms in the transport equation for the condit ional 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 condit ional average based on a single conditioning variable may be extreme; the left hand port ion of Figure 1.3 is a scatter plot of temperature vs. mix- ture fraction from a direct numerical simulation of a flame wi th 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 wi th the scalar variations not correlated wi th 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 wi th promising results [26]. In a realistic combustion simulation, the choice of second conditioning variable is unclear. C h a et al. [9] attempted to use the scalar dissipation, after dividing by a presumed functional dependence on mixture fraction to produce a variable somewhat independent of mixture fraction (as shown by Bushe and Steiner [7]). However, as explained earlier by Bushe and Steiner [7], although this quantity has many of the attributes desirable in a second conditioning variable it wi l l ult imately prove unsuitable. A s the strain on a flame increases, the local gradients become steeper and the rate of mixing increases. For a mi ld increase in mixing rates a corresponding increase in reaction rates, and a moderate change in the species condit ional 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 wi l l not increase unti l re-ignition occurs. Re-ignit ion involves either the rate of strain dropping below the level at which autoignition is likely, or an external event such as a passing edge flame [38]. Further complicating matters, there CHAPTER 1. BACKGROUND 17 is a temporal history effect. A n increase 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 port ion of Figure 2 in [9], illustrates the mult ipl ici ty 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 port ion of Figure 1.3, al l of the points plotted have a constant mixture fraction, so an ideal condit ioning variable would have a more or less unique relationship wi th temperature. Obviously, this is not so for the scalar dissipation. Z 'Og Xst Figure 1.3: Relat ion Between Temperature and Mixture Fract ion, Log of Scalar Dissipation. Port ion 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- t ional 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 wi l l be considered. 1.4 Conditional Source-Term Estimation Condi t ional source-term estimation (CSE) ult imately uses the same first or- der closure based on conditional means to evaluate the chemical source terms CHAPTER 1. BACKGROUND 18 employed in C M C . Rather than solving directly for the condit ional averages, the flow is formally solved for the unconditional means in physical space. The condit ional means, which are assumed constant in space over ensembles of cells, are related to the unconditional means and the pdf of the condit ioning variables through integral equations which are then inverted for a discreet approximation to the conditional means. It is these condit ional 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. Th is method was tested in an a priori sense against a D N S database [7] and was found to be effective when litt le extinction was present using only the mixture fraction as a conditioning variable. The addit ion of a modified form of the scalar dissipation as a second conditioning variable produced success when extinction was present. Later, the method was used to obtain closure for the chemistry in a L E S calculation of the Sandia Flame ' D ' using 2-step chemistry wi th 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 condit ional source terms associated wi th 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, B la i r [6] used this formulation in a K IVA- I I simulation of a turbu- lent methane jet injected into hot air. However, Bla i r 's implementation was unable to successfully predict the ignit ion delay when compared to exper- imental data. In the remainder of this section, the general formulation of C S E wi th Laminar Flamelet Decomposit ion used throughout this thesis wi l l be presented. Following Bushe and Steiner [7], first note that the uncondit ional aver- age of any scalar quantity can be related to the condit ional average of that quantity by employing the probabil i ty density function of the condit ioning variable. Using a 'Favre density function' p ((; Z, Z"2^j [22], the relationship then holds for Favre averages: (1.20) A s per [8], decomposition of the conditional average in question into a sum- CHAPTER 1. BACKGROUND 19 mat ion of basis functions, e.g., <T|C> « 5 > n < T | 0 n 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 wi th respect to mixture fraction. Substitut ion of E q . 1.21 into E q . 1.20 produces: Jo N .n=l (1.22) A t this point, we can see that E q . 1.22 is a Fredholm Equat ion of the first k ind: wi th 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 a n=i..JVj=i..j is constant, we can write E q . 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 wi th noise), each wi th 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 typical ly i l l-posed (it could only be well posed if each of the N observations came conveniently supplied wi th a response kernel that was a delta function at each of the quadrature points), the recommendation in [37] is to add a priori informa- t ion 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 wi l l be discussed in the following chapters, that a „ should have any elementary func- t ional 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 condit ion, the solution should evolve smoothly in time. Turning now to the algebra necessary to accomplish this, we revisit a rearranged version of E q . 1.22 and switch to vector/matr ix notation: Jo (1.29) Forming a i f -po in t quadrature in £-space, / P ( C ) / ( C K « 5 > ( C * ) / ( C k K A f c J o 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 E q . 1.30 to discretize E q . 1.29 and writ ing it for each of J computational cells in an ensemble over which we assume ~a is spatial ly homogeneous in the fashion of E q . 1.28, neglecting rij, we obtain: f(T\Ok=l,n=l ^j=l,fc=l • • • ^1,K \ T j J \ . . . jjjcj V (T\C)K,i or, where A is defined by comparison to E q . 1.31, (T\C)K,N) a (1.31) Aa. (1.32) CHAPTER 1. BACKGROUND 21 We can also write an equation resembling E q . 1.20 for any scalar of in- terest. Repeating the above procedure for an arbitrary species mass fraction Yi, and defining: Bi = / ^ I0 f c = 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, ( l i | C ) „ are corresponding scalar profiles from the same quasi-physical flame. In this situation, solution of either E q . 1.32 or E q . 1.34 for any i should yield the same ~a. However, we must recall the ill-posedness of these equations; not only is there a ques- t ion of existence, there is a significant question regarding uniqueness. Later, in Sections 2.3 and 3.2.3, it wi l l be i l lustrated that using a typical l ibrary 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, E q . 1.32 and E q . 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) wi th: M = ( A \ \ B I J R = \YiJ (1.36) Despite our hope that E q . 1.35 contains sufficient information to select the most suitable combination of basis functions, it is st i l l potentially i l l - 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 E q . 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 possibil i ty that it is impossible to reconstruct the conditions in the flow exactly due to the finite size of the l ibrary 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]. Whi le 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. Th is 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 minimizat ion 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: Th is reformulation as a minimizat ion problem addresses the question of existence of a solution. However, it means that we must take care in assem- bl ing E q . 1.35 to normalize both sides of each equation, lest the minimiza- t ion procedure place greater emphasis on satisfying those equations involving scalars wi th larger absolute magnitudes. The alternative formulation of the minimizat ion 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 al l scalars wi th the same weight, we can accomplish our objective simply by normalizing E q . 1.32 and E q . 1.34 by the maximum present so that al l of the scalars fall in the range [0,1]. Returning to the original minimizat ion problem E q . 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 + A c * 0 . (1.39) We now have an opportunity to provide an arbitrary ini t ia l solution of 0, and allow the minimizat ion to select a new solution ~a which deviates from the original only to the extent that the current flow field, incorporated into M, justifies. Whi le the smoothing used for stabil ization 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 t ime is altered, the physics captured wi l l be altered. A s wi l l be discussed later, this is clearly unacceptable. In the two applications of this method discussed in Chapters 2 and 3, the stabil izing solution ~a° wi l l 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 wi th that proposed by Bushe and Steiner [7, 8] and used by Bla i r [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. A l l 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- t ion space. However, this necessitates inversion for a separate progress variable for each reaction in the mechanism used, which is a pro- hibit ively 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 wi l l be discussed further in the following chapters. • Inversion Based on Multiple Species. Whi le 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 wi th the mean fields for minor species, which discriminate much more clearly between the flamelets. CHAPTER 1. BACKGROUND 24 Addi t iona l differences between this and previous implementations are present; however, they are of a problem-specific nature; further discussion wi l l be left unti l the final chapter. REFERENCES 25 REFERENCES [1] Proceedings of the Third International Workshop on Measurement and Computation of Turbulent Nonpremixed Flames, Boulder, Colorado, 1998. U R L h t t p : / / w w w . c a . s a n d i a . g o v / T N F / 3 r d W o r k s h o p / T N F 3 . h t m l . Section 2 - Sandia / Darmstadt Pi loted CH4/Air F lame 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 Combust ion Institute, 2000. [3] H . Barths, H. Pi tsch, and N. Peters. 3d simulation of di diesel combus- t ion 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. Condit ional moment clsoure for turbulent reacting flow. Physics of Fluids A - Fluid Dynamics, 5(2):436-444, 1993. [6] C . Bla i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of Br i t i sh Columbia, 2003. [7] W . K. Bushe and H. Steiner. Condit ional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11 (7): 1896-1906, Ju ly 1999. [8] W . K. Bushe and H. Steiner. Laminar flamelet decomposition for condi- t ional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [9] C M . C h a , G . Kosaly, and H. Pi tsch. Model ing extinction and reigni- t ion in turbulent nonpremixed combustion using a doubly-condit ional moment closure approach. Physics of Fluids, 13(12):3824-3834, 2001. REFERENCES 26 [10] P. Coelho and N . Peters. Numerical simulation of a mi ld 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 applicabil i ty of condit ional 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. U R L h t t p : //www. me. b e r k e l e y . e d u / g r i - m e c h / . [15] J . Huang, P. H i l l , W . Bushe, and S. Munshi . Shock-tube study of methane ignit ion under engine-relevant conditions: experiments and modeling. Combustion and Flame, 136(l-2):25-42, 2004. [16] S. H. K i m and K. Y . Huh. Use of the conditional moment closure model to predict N O formation in a turbulent C H 4 / H 2 flame over a bluff-body. Combustion and Flame, 130:94-111, 2002. [17] S. H. K i m , K. Y . Huh, and R. W . Bilger. Second-order condit ional 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. K i m , K. Y . Huh, and R. A . Fraser. Model ing autoignit ion of a turbulent methane jet by the conditional moment closure model. Pro- ceedings of the Combustion Institute, 28:185-191, 2000. [19] S . -K . K i m , S . - M . Kang , and Y . - M . K i m . Flamelet modeling for com- bustion processes and N O x formation in the turbulent nonpremixed C O / H 2 / N 2 jet flames. Combustion Science and Technology, 168:47-83, 2001. [20] W . T . K i m and K. Y . Huh. Numerical simulation of spray autoignit ion by the first-order conditional moment closure model. Proceedings of the Combustion Institute, 29:569-576, 2002. REFERENCES 27 [21] A . Kl imenko. Mult icomponent diffusion of various admixtures in turbu- lent flow. Fluid Dynamics, 25(3):327-334, 1990. [22] A . Kl imenko and R. Bilger. Condit ional moment closure for turbulent combustion. Progress in Energy and Combustion Science, 25(6):595-687, 1999. [23] D. Lent ini . 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. Condit ional moment closure for methane oxidation using two condit ional variables and stochastic processes. Master 's thesis, Univesity of Br i t ish Columbia, 2004. [27] E . Mastorakos and R. Bilger. Second order condit ional 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- ignit ion wi th ell iptic 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 . Pi tsch. Unsteady flamelet modeling of differential diffusion in turbu- lent jet diffusion flames. Combustion and Flame, 123(3):358-374, 2000. [33] H . Pi tsch. 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. Pi tsch, 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 Combust ion In- stitute, 1998. [35] H. Pi tsch 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. Pd f methods for turbulent reactive flows. Progress in Energy and Combustion Science, 11:119-192, 1985. [37] W . Press, S. Teukolsky, W . Vetterl ing, and B. Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, U K , 1992. [38] G . Ruetsch, L. Vervisch, and A . L inan. 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 modell ing in variable density turbulent axisymmetric jets and diffusion flames. Physics of Fluids, 10(4):938-948, Ap r i l 1998. [40] H. Steiner and W . Bushe. Large eddy simulation of a turbulent reacting jet wi th 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 . Wal l , B. J . Boersma, and P. Mo in . A n evaluation of the assumed beta probabil i ty density function subgrid-scale model for large eddy simula- t ion of nonpremixed, turbulent combustion wi th 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 edit ion, 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 ( R A N S ) 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. Condi t ional moment closure ( C M C ) [54] seeks to address this issue by solving for conditional means where the conditioning variable accounts for much of the unresolved variation, leading to significantly smaller f luctua- tions about the mean quantities. Th is allows for much improved closure of the source-term compared to evaluating the source-term wi th uncondi- t ional 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 arbitrari ly complex chemistry in principle, an addit ional transport equation must be solved for each species added. Condit ional source-term estimation (CSE) [50, 58] al- lows use of the C M C hypothesis in the L E S paradigm by solving for only uncondit ional means and then inverting integral equations to determine the condit ional means necessary to close the source-terms. This method can also be applied in the R A N S paradigm resulting in substantial computational sav- ings compared to tradit ional C M C . However, using realistic chemistry st i l l results in a prohibit ively 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 condit ioning variable in C M C . However, these models are, str ict ly speak- ing, val id only within the flamelet regime and there is considerable concern regarding their applicabil i ty outside this regime. As well, except in l imited circumstances, it is impossible to represent a flame of practical interest wi th a single flamelet solution. A recently proposed extension to C S E , where the condit ional 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 arbitrari ly 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 D N S simulation in the a priori sense and the R A N S paradigm. CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 32 In this chapter, the method has been adapted for use in a commercial R A N S 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. Th is geometry has been uti l ized to generate extensive experimental data for a variety of flames by Sandia Nat ional 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 probabil i ty 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. Condi t ional source-term estimation wi th laminar flamelet decomposition [51] seeks to close the chemical source-term by inversion of integral equations for the condit ional averages of key scalars to determine the combination of flamelet solutions necessary to construct an accurate representation of the flame. The condit ional source-term is then obtained from the corresponding flamelet solutions, and integrated over the Favre density function (probabil- i ty density function wi th 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 condit ional 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 condit ional average (T\Q can be constructed. The choice of the ensemble of cells over which the structure of the flame is constant in a C F D 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\ ( 2 3 ) In searching for appropriate basis functions, an obvious choice is solutions to the unsteady laminar flamelet equations [55]: P——=u(Z,t) + —-—x{Z,t)-Qz2- (2-4) These solutions produce, within the assumptions of the model, conditional averages for each scalar of interest and can be solved a priori using 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 crit ical to choose the correct solution — this is why it is cri t ical that the a priori information not be arbitrary, but rather provide realistic insight into the physical process. For the three solutions shown, the condit ional averages of carbon monoxide are quite different; matching this scalar also to the flow field wi l l dramatically reduce the reliance on the a priori information. Figure 2.1: Selected Basis Function Solutions To solve E q . 2.3, we note that the equation is a Fredholm equation of the first k ind and replace the integral with a numerical approximation [56]. W i t h this observation, it is straightforward to cast E q . 2.3 into a system of approximate linear equations of the form M~a « R, where the number of rows 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 minimizat ion problem, while the addit ion of a priori information encourages a unique solution. If we desired continuity of the solution, it could be encouraged through the a priori information; however, a discontinuous solution does not preclude obtaining a numerical solution as existence and uniqueness issues do. Solving E q . 2.3 for an, the mean chemical source-term can then be recovered in the manner of E q . 2.1. Once a suitable choice for ~a has been made it is possible to obtain, using the same 7* in the forward problem, the conditional averages for any arbitrary species. The conditional source terms for al l 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 C S E wi th laminar flamelet decomposition. Following the earlier method, we formed our system of integral equations relating the uncondit ional 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 l ibrary of flamelet solutions, the condit ional average of temperature alone is insufficient to discriminate between flamelets. We therefore augment our system of integral equations wi th 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 F lame ' 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 l ibrary used and noting that these two scalars are suffi- cient to distinguish between al l 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 + A a ° . (2.6) (2.7) CHAPTER 2. SIMULATION OF THE SANDIA FLAME 'D' 36 In E q . 2.6, there is a balance between minimizat ion of the error in ' f i t t ing' 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 explicit ly 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 wi th a physically realizable state, even in the absence of any useful data for the inversion process a physical solution wi l l result. Other regularization possibilities (e.g., smoothing) lack this justif ication and may not lead to a physically realizable solution. For example, B la i r [49] found a tendency to solutions incorporating a substantial number of negative coefficients using a smoothing operator for the a priori information, which the present choice does not encourage. 2.4.2 Flow Solution A R A N S simulation of the flame was performed using the k — e turbulence model wi th standard model constants, except for the value of C e 2 = 1.8, in preference to the standard value of C £ 2 = 1.94 to obtain better agreement wi th the experimentally measured mixing field near the nozzle exit. In ad- dit ion 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 wi th 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: C S E / L F D 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 A n adaptation of conditional source term estimation has been used to predict temperature and scalar fields for the Sandia Flame ' D ' . Despite an inabi l i ty to incorporate extinction/reignit ion 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. A s well, basis functions which are representative of the edge flame expected to provide (re)ignition in this flame wi l l be necessary in the future. Whi le the prediction of nitr ic 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: C S E / L F D Val idat ion - Comparison of Ax ia l and Radia l 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. U R L http:/ /www.ca.sandia.gov/TNF/3rdWorkshop/TNF3 .html. Section 2 - Sandia / Darmstadt Pi loted CH^/Avc F lame 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 Na t l Lab, Livermore, C A , U S A , 1998. Combust ion Inst, Pi t tsburg, P A , U S A . I S B N 0082-0784. [49] C . Bla i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of Br i t i sh Columbia, 2003. [50] W . K. Bushe and H. Steiner. Condit ional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, Ju ly 1999. [51] W . K. Bushe and H. Steiner. Laminar flamelet decomposition for condi- t ional source-term estimation. Physics of Fluids, 15(6):1564-1575, June 2003. [52] Fluent 6.1. Fluent Inc., 2003. [53] D. M . G . Gregory P. Smith et al . Gri-mech 3.0. U R L h t t p : //www. me. b e r k e l e y . e d u / g r i - m e c h / . [54] A . Kl imenko and R. Bilger. Condit ional 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 . Vetterl ing, and B. Flannery. Numerical Recipes in Fortran. Cambridge University Press, Cambridge, U K , 1992. REFERENCES 45 [57] G . Ruetsch, L. Vervisch, and A . L inan. 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 wi th 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. N O N - P R E M I X E D M E T H A N E A U T O - I G N I T I O N 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 condit ional source- term estimation wi th laminar flamelet decomposition approach proposed in [62]. The ambient oxidizer conditions — the high pressure and moderate temperatures characteristic of compression ignit ion engines — were chosen wi th the intent to validate the combustion model used under engine relevant conditions. Th is was obtained by comparison of the predicted ignit ion delay to experimental results obtained from the shock-tube facility at U B C at several ini t ial temperatures. Ear ly experimental measurements of the ignit ion delay time were obtained by Fraser et al . [63], who injected methane into a constant volume combus- t ion 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-l ike' mixture. Subsequently, other studies have used the same facility to explore several parameters potentially affecting the ignit ion delay, including the ambient pressure, temperature, and composit ion 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. Bo th 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 st i l l underway at the University of Br i t ish Columbia have sought to measure the ignit ion 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 ignit ion event in the U B C study differed from the earlier work using pressure based criteria in that the U B C study identified the onset of ignit ion wi th optical ob- servations. Through the use of a high-speed camera, the U B C study defined the end of the ignit ion 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 ignit ion delays. However, the most recent results [67, 78] agree well wi th one another. Given the dramatic differences in the experimental apparatus and techniques used to detect ignit ion, 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 O n the simulation front as well, previous efforts have been made to de- velop a computational method to model ignit ion. Bo th Condit ional Moment Closure [68] and Flamelet models [69] have been applied to the simulation of the experiment described by Fraser et al . [63]. Bo th methods produced results which agreed remarkably wi th the experimental data presented by Fraser et al . [63]. However, the subsequent introduction into the literature of experimental data disagreeing wi th Fraser et al.'s results by Ishiyama et al . and the anticipated publishing of supporting data from the U B C shock tube give reason to be somewhat cautious regarding these simulations. In the later discussion, some speculation wi l l be made regarding the possible ex- planation for the discrepancy, but further careful study both experimentally and through modeling exercises is necessary. In addit ion to providing an as- sessment of the performance of the C S E - L F D combustion model, this work seeks to investigate in some detail the process leading up to the detection of an ignit ion event, wi th the goal to produce an ignit ion delay prediction that is directly comparable to that measured in the U B C shock tube facil ity by optical means. Further, as the method uti l ized to measure the ignit ion experimentally provided a spatial location for ignit ion, the current effort was also concerned wi th predicting the location of ignit ion. 3.2 Formulation 3.2.1 CSE Combustion Model Closure of the chemical source term was obtained by an extension to con- dit ional 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 uti l ize the (first moment) C M C hypothesis [70] to treat the chemical source terms, but wi th- out solving directly for the conditional means. Instead, it is noted that the condit ional 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 condit ional average in E q . 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 L E S 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 E q . 3.1. Incorporating the effects of complex chemistry necessary to treat the au- toignit ion 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- hibit ive 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 uti l ized. In this formulation, solutions to the laminar flamelet equations are generated a priori, and E q . 3.1 becomes: _ f1 N T « / Y,an(T\0np{Z)dZ, (3.2) J o n=l where the approximation: <T|C> « f>n (T|C)n n=l has been used. To obtain the l ibrary 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 al l species and temperature) present at each point in mixture fraction space. If we are wi th in the flamelet regime, it would be reasonable to suppose that the coefficients an which best satisfy E q . 3.2 would be consistent regardless of which scalar was considered. Th is impl ic i t ly 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) a n d 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 \ Y C O J ^ 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 condit ional average is assumed to be constant. A practical difficulty is encountered in that solving this system is not tr iv ial . 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 minimizat ion 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 E q . 3.3 and some a priori knowledge regarding the solution. Whi le constraints on the first derivative of the solution for an wi th respect to basis function order in the l ibrary have been used in earlier work [62], this approach has several drawbacks. The relationship between flamelets which reside next to each other in the l ibrary is not necessarily unique, as it depends on the manner in which the l ibrary was constructed. Simple sorting based on scalar dissipation rates becomes unclear when the l ibrary 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. A s a matter of convenience, it is desirable that the order in which the basis function are stored in the l ibrary 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 minimizat ion of the difference between the left and right hand sides of E q . 3.3 and the difference between the coefficient for each basis function at the current t ime step and the coefficient selected for the same basis function at the previous time step. The actual equation solved is that given in E q . 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 E q . 1.37 ) is chosen such that the solution wi l l only be pulled toward the a priori guess the minimum amount necessary to produce a well behaved solution. In addit ion to a more plausible physical in- terpretation, and allowing for the libraries of basis functions to be buil t 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- condit ional mean of any scalar of interest present in the l ibrary can be de- termined by evaluation of the forward problem given by E q . 3.2, replacing temperature wi th 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 condit ional mean reac- t ion rate constructed by linearly combining the source terms from each basis function weighted by the corresponding an and the pdf of the condit ioning variable in that instance. The reader may have noticed a switch in E q . 3.3 to Favre-averaged quantities; this implies use of the Favre density function, p(£; Z, Z"2), in place of a probabil i ty density function, Z, Z"2). To re- cover the non-Favre averaged source terms necessary for the R A N S transport equations, the forward problem to evaluate becomes: where the relation between the Favre and non-Favre pdf suggested by K l i - menko and Bilger [70] has been uti l ized. A n 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 wi th 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 ignit ion. A structured mesh was used to facilitate the grouping of cells into ensembles on planes perpendicular to the centerline for the C S E inversion process. 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 ignit ion is rea- sonably consistent on planes normal to the jet axis. A n overly coarse grid not only results in a less accurate solution to the scalar transport equations, but restrains the variation of the conditional averages 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 C S E for simulation of autoignit ion process [60]. In addit ion 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 C S E inversion process. (3.4) 3.2.2 Flow Field CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 53 The domain was init ial ized with quiescent air at 30 bar and zero mixture fraction and variance; a pure methane jet entered from the bot tom left. The boundary condition for the methane jet was a pressure boundary condit ion resulting in a moderately under-expanded jet. A l though 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 wi th 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 t ime 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 explicit ly in time. Inversion for the discrete condit ional averages was performed on a 51 point grid in mixture fraction space, wi th the majori ty of the points clustered around the stoichiometric mixture fraction. 3.2.3 Basis Functions Pr ior to undertaking the flow simulation for each l ibrary 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 t ime using realistic chemistry developed by Huang et al . [65]. Th is mechanism was used without adjustment of the various constants tuned to premixed autoignit ion experiments. . A n internally adiabatic mixing solution between boundary conditions matching the fuel and oxidizer ini t ial states was used as an ini t ial condition. The coupled part ial differential equations were then reduced to a system of ordinary differential equations using a method of lines approach 1 ; the system of O D E ' s was solved using L S O D E [64]. Every 50 [is in flamelet t ime, 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 ignit ion 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 ( e r f c " 1 ( 2 Z ) ) 2 . (3.6) The temporal port ion, Xo(^), was chosen to have the functional form: Xo(t) =cxe^ + c 3 . (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, c 2 = 7 x 10~ 6 . A collection of solutions for c 3 = 0 wi th 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 al l 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 al l taken prior to ignit ion 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 Flamelet #25 ! ! 1 \ 0.4 0.6 Mixture Fraction Figure 3.2: Flamelet Solutions for Species Mass Fract ion - Yco and YCH3OH ignit ion, 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 condit ional mass fraction of CH3OH was chosen as it evolves gradually in time prior to ignit ion, 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 al l stages of the ignit ion 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 ignit ion. 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 al l 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 ignit ion obscures any changes before ignition. In the exact flow, although the mean scalar dissipation decays exponen- t ial ly wi th 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 ignit ion process temporary excursions in scalar dissipation, far in excess of the crit ical value that would normally prevent ignit ion, do not necessarily prevent ignit ion 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 wi th much higher scalar dissipation rates in the library. Th is 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 l ibrary of basis functions for addit ional values of the parameter c 3 in E q . 3.7. The final l ibrary used incorporated sub-libraries generated with c 3 = 0,0.5,10,100 concatenated as il lustrated for the tem- perature solutions in Figure 3.4. The portion of the l ibrary wi th c 3 = 100 showed approximately a 10% increase in the magnitude of the temperature source term over the port ion with c 3 = 0 for the pre-ignition solutions; the change was much more pronounced (exceeding of 100%) after ignit ion occurs in the flamelet calculation. Figure 3.4: Flamelet Temperature [K] Contours for Complete L ibrary Here, we see a practical benefit of not enforcing dependence of an on n. A t a given time, it is conceivable that if the ideal combination of basis functions included a large amount of solution #50, we can see from Figure 3.4 that it might also include some amount of solution #125, which is very similar in terms of its conditional mean of temperature. In fact, since 75 flamelets were included using each scalar dissipation history and the solutions were extracted at the same time interval, solutions 50, 125, 200, and 275 all represent the same flamelet t ime, 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 stabil ization 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 ignit ion delay time wi th temper- ature is shown in Figure 3.5 along wi th the experimental data from the U B C shock tube facility [78]. The results shown were produced using an identical procedure; only the ini t ial temperature field in the flow solver and the ini t ia l oxidizer temperature used to bui ld the l ibrary 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 D el ay [m s] • D el ay [m s] c • b Ip jn itk • > * a b Ip jn itk • 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. A s a matter of interest, the fit of the simulation predictions to the experimental data was compared to the fit of a linear regression based on the available data. A s the object of this simu- lation exercise was to capture the effect of temperature on ignit ion 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. Whi le 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 addit ion to predicting the time at which ignit ion 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 ignit ion site appeared. Table 3.2 summarizes the t ime and location at which ignit ion sites are predicted. The laminar ig- nit ion delay, which is extracted from the l ibrary of basis functions used for each case, is included here to facilitate later discussion. Th is laminar delay is obtained by averaging of the ignit ion delays for the various histories in the appropriate l ibrary of basis functions, using the condit ionally averaged temperature increase criterion. Returning to the location data, a scatter plot of the experimentally observed ignit ion location against the predicted ignit ion location is inconclusive. Whi ls t some of the points coincide either exactly or nearly so, the majority of the experimental points fall between the simulation prediction and twice the distance downstream. 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. Addit ional ly, interpretation of the experimental data is diffi- cult due to many simultaneous or nearly simultaneous ignit ion sites. Wh i le it does appear that many of the experimentally observed ignit ion 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 ignit ion in the shock tube do not allow for measurement of the jet penetration, so confirming or disproving this significant possibil i ty 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 ignit ion occurred closer than Simulation Linear Predict ion Regression SSE K2 0.288 0.223 0.52 0.63 (SSyy = 0.596, fff^ = 1.29 ) Table 3.1: Al locat ion 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 ignit ion 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: Locat ion and Mixture Fract ion 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 K i m et al . [68] using condit ional moment closure, and an earlier implementation of condit ional source-term estimation reported by Bla i r [60]. A n third previous simulation undertaken by K i m et al. [69] also attempted to predict Fraser et al.'s results using a flamelet model. Those results very nearly coincide wi th those in [68], and so are omitted from the figure for clarity. It would appear prima facie that both the experimental results used for val idation 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 wi th the strong resemblance between Ishimyama et al.'s data and U B C ' 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 ignit ion delay on that pressure was found CHAPTER 3. NON-PREMIXED METHANE AUTO-IGNITION 61 + Fraser et al. -o- C S E - L F D • Kim et al. • U B C Shock Tube * Ishiyama et al. Blair *% r+ + ++ +++ ++ + 4.48 00 1 . 6 5 ^ >. 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 Pr ior Works o *̂ c c n 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 startl ing discrepancy, we must consider once again the question of what defines ignit ion. In obtaining the experimental data, Eraser et al. [63] uti l ized a measure of the pressure rise in the vessel to identify ignit ion. The more recent investigation discussed by Ishiyama et al . [67] reconsidered the question of determining ignit ion 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 l iquid fuels, for gaseous fuels detection of ignit ion 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. Th is means that the actual start of significant reaction must have occurred some time earlier. Further, the rate of heat release during this period following ignit ion is not necessarily constant between different oxidizer temperatures; specifically the delay between the true onset of ignit ion and that detected by the pressure rise wi l l vary wi th the amount of fuel mixed before ignit ion. In particular, by inspection of Figure 4 in [63], for longer ignit ion delays, the method of anal- ysis used to detect ignit ion from the pressure rise in the vessel employed by Fraser et al . appears to overestimate the ignit ion delay by nearly 2 ms for long ignit ion delay times, and by perhaps 0.3 ms for short ignit ion delays times. The measure used by Ishiyama et al . is much more likely to coincide wi th the onset of ignit ion determined by optical measures at U B C , as it defines ignit ion 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 st i l l expected to be longer than the delay detected at U B C as light emissions can be expected to be observed before the pressure responds. The earlier C M C and flamelet simulations are in reasonable agreement wi th Fraser et al.'s results, but not wi th the more recent experiments. A s discussed above, the criteria used there for ignit ion is one which wi l l identify ignit ion after the system has been burning for some significant amount of t ime. 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 w i th the matching data does not necessarily guarantee that they would be able to identify the onset of ignit ion using the more stringent criteria used in the latter experiments. After ignit ion (by the new definition) has occurred, the rate of pressure rise wi l l 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 t ime at which ignit ion actually occurred. Whi le 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 ini t ial flame propagation following igni- t ion. Wi thout such effects, it is unrealistic to continue the current simulation beyond ignit ion to the time at which a criterion corresponding to the earlier work could be applied. Despite this l imitat ion, which applied to the basis functions used by B la i r [60] as well, B la i r used a criteria to define ignit ion which would also occur much later than ignit ion as defined here — that the mean control volume average temperature exceeded 2000K. Th is can only occur well after the igni- t ion prediction based on the current criteria. The mechanism by which B la i r obtained ignit ion delays shorter by nearly an order of magnitude based on a criteria that should have delayed the identification of ignit ion 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 , B la i r used only temperature to perform the inversion for ~a; discussion elsewhere presents the difficulty dis- t inguishing between flamelet basis functions on temperature alone. However, B la i r d id not appear to suffer from a uniqueness problem problem related to using only temperature; Figure 4.3 in [60] indicates that in the l ibrary of basis functions used, the conditional temperature appears to be changing substantially between flamelets. However, the lack of flamelets wi th similar temperatures highlights a larger issue: in the l ibrary used by Bla i r , there are at best two flamelets (one of which is the ini t ial condition) present which have a peak temperature less than 1800K. This suggests that the ignit ion process is almost entirely unrepresented in the l ibrary of basis functions. A s well, the lack of ' ignited' solutions wi th similar temperatures is the result of using only one scalar dissipation history, which was increasing in t ime as opposed to the expected decrease in time in the problem under consideration, to bui ld the library. Th is suggests that it is quite possible that if an addit ional scalar were included, the minimizat ion of the least squares problem necessary to de- termine ~a may have had a very large residual indeed wi th 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 Bla i r '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 condit ional averages of any scalar up to the validity of the closure. A key feature of the ignit ion process is that some finite amount of t ime elapses without a noticeable increase in temperature and a slow bui ld up in reaction rate, fol- lowed by a rapid increase in temperature and reaction rate [72]. Whi le early implementations of C S E + L F D in K I V A - 3 V 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 condit ional averages of the three scalars used for the inversion process are shown in F ig - ure 3.7 along wi th the conditional average of the creation rate of the CH radical at several times prior to ignition. For clarity, only the condit ional average from the plane showing the largest increase over the ini t ia l condit ion is shown. Here we see that there is relatively minor change in the condit ional aver- age of temperature prior to ignit ion, followed by the development of a local increase in temperature near the stoichiometric mixture fraction. Following ignit ion, the increase in temperature spreads outward from the stoichiometric mixture fraction. The lack of change in the conditional average for temper- ature prior to ignit ion suggests, once again, that this is a poor scalar to base a determination of the state of the system on prior to ignit ion. The concentration of carbon monoxide also increases dramatically, but there is a more apparent increase in the concentration prior to ignit ion. The evolution of the condit ional average of the intermediate CH3OH is more interesting. The peak condit ional magnitude of this species is substantial well before ig- nit ion takes place, and changes location in mixture fraction space as ignit ion 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 ignit ion. Th is shift, as well as the non-negligible concentration of the scalar well prior to ignit ion, justifies the choice of this scalar to aid in selecting the appropri- ate combination of flamelets prior to ignit ion. Final ly, 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 U B C shock tube experiments. In the absence of soot, the hypothesis is that the light emissions observed are primari ly the result of CH radicals decaying from an excited state, and that the quantity of excited molecules is proport ional to the creation rate. The behaviour is promising for the detec- t ion of ignit ion; as wi th temperature, there is a very dramatic increase just as ignit ion occurs. Setting the absolute level at which to declare ignit ion has occurred is not obvious; in the U B C 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 condit ional average of the var- ious scalars at any mixture fraction and on any plane is plotted against t ime for a selection of oxidizer temperatures. The various quantities are normal- ized by a threshold chosen for the quantity, such that ignit ion was declared based on the various quantities when the quantity plotted exceeds a value of one. A l though 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. W i t h the thresholds chosen, there is very l itt le difference in the ignit ion delay identified by consideration of any of the three scalars, as is apparent from Figure 3.9. Given this comparison, the ignit ion 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 wi th the optical criteria used in the experiments. The data presented earlier (Figures 3.6, 3.5) was generated using this criterion. In addit ion to consistency in the predicted ignit ion delay, the location along the centreline at which ignit ion 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 i n c w 0.083) than the temperature criterion (Z\(T\Qmax i n c — 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 T O I = 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: Evolut ion of Peak Increase in Condit ional Averages 0.5 r E %. 0 o Temperature Increase - * C H Production Rate C O 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 Cr i ter ia 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 ignit ion is not immediately obvious. There are competing effects; increased turbulence levels wi l l increase the rate of mixing. Whi le this means that the fuel and oxidizer wi l l mix sooner, it also means that any radicals produced or temperature increase wi l l be dissipated quickly. Further, in a turbulent flow, there wi l l be a distr ibution of strain rates; wi th in any given region, there wi l l be a variety of strain histories present [72]. The laminar and turbulent ignit ion 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 ignit ion delay. A t low tem- peratures, where the delay is longer, the effect of turbulence is to decrease the delay. Interestingly, the ratio of turbulent to laminar ignit ion delays is almost linear wi th 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 ignit ion delay is short, the effect of tur- bulence dissipating temperature and radical bui ld up is dominant; although the presence of turbulence also suggests that there wi l l be regions of low strain, the likelihood of one of these regions having been already sufficiently mixed for ignit ion to occur is low. A t longer ignit ion delays, the majori ty of the field is well mixed prior to ignition. Then, the distr ibution of strain rates suggests that there is a strong l ikelihood that one of the well mixed regions wi l l encounter low strain, resulting in ignit ion promotion. 3.4 Conclusion A simulation has been undertaken to predict the ignit ion 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 R A N S paradigm takes advantage of concepts from both the C M C and flamelet models. A separate condit ional average is permitted (and encouraged) on each plane of cells ax- ial ly 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 wi l l be necessary. The predicted ignit ion delay was in excellent agreement wi th 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 wi th an addit ional 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 ignit ion event, as well as the occurrence of multiple, nearly simultane- ous ignit ion events in the experiment, prevented a meaningful val idation of the predicted ignit ion location. However, the l imited comparison that was possible did not present any evidence that we should be skeptical of the sim- ulat ion results. In order to produce simulations valid beyond ignit ion, it is necessary to develop libraries of basis functions which encompass the appro- priate physics to address the flame propagation following ignit ion, which are not captured in the solutions to the laminar flamelet equations. REFERENCES 71 REFERENCES [60] C . Bla i r . Implementation of conditional source-term estimatoin for prediction of methane ignition. Master's thesis, University of Br i t ish Columbia, 2003. [61] W . K. Bushe and H. Steiner. Condit ional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, Ju ly 1999. [62] W . K. Bushe and H. Steiner. Laminar flamelet decomposition for condi- t ional source-term estimation. Physics of Fluids, 15(6): 1564-1575, June 2003. [63] R. A . Fraser, D. L. Siebers, and C. F. Edwards. Autoigni t ion of methane and natural gas in a simulated diesel environment. SAE Technical Paper Series, 1991. [64] A . C . Hindmarsh. O D E P A C K , A systematized collection of ode solvers. In R. Stepleman et al . , editors, Scientific Computing, pages 55-64, A m - sterdam, 1983. North-Hol land. [65] J . Huang, P. H i l l , W . Bushe, and S. Munshi . Shock-tube study of methane ignit ion under engine-relevant conditions: experiments and modeling. Combustion and Flame, 136(1-2):25-42, 2004. [66] J . Iaconis. A n investigation of methane autoignition behaviour under diesel engine-relevant conditions. Master's thesis, University of Br i t ish Columbia, 2003. [67] T . Ishiyama, M . Shioj i , et al. Characteristics of spontaneous igni- t ion and combustion in unsteady high-speed gaseous fuel jets. In 2003 JSAE/SAE International Spring Fuels & Libricants Meeting, number J S A E 20030225 / S A E 2003-01-1922, Yokohama, Japan, May 19-22 2003. [68] S. H. K i m , K. Y . Huh , and R. A . Fraser. Model ing autoignit ion 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 . K i m , Y . Y u , et al. Numerical investigation of the autoignit ion of turbulent gaseous jets in a high-pressure environment using the multiple- rif model. Fuel, 83:375-386, 2004. [70] A . Kl imenko and R. Bilger. Condit ional 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- pat ion rate on ignit ion of non-premixed hydrogen/air mixtures in coun- terflow. Procedings of the Combustion Institute, 29:1629-1636, 2003. [72] E . Mastorakos, T . Bar i taud, and T . Poinsot. Numerical simulations of autoignit ion in turbulent mixing flows. Combustion and Flame, 109: 198-223, 1997. [73] J . D. Naber, D. L. Siebers, S. S. Di jul io, and C. K. Westbrook. Effects of natural-gas composition on ignit ion 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 . L inan. 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 wi th conditional source-term estimation. Physics of Fluids, 13(3): 754-769, March 2001. [77] G . Sull ivan. Personal communication. 2004. [78] G . Sull ivan. 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 appl ica- t ion of the method to a reacting flow wi th a steady ensemble averaged solu- t ion. Whi le 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- toignit ion, 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 wi th condit ional source-term estimation using laminar flamelet decomposition, the effects of arbitrar i ly complex chemistry can be accounted for wi th min imal computa- t ional cost. This is because the majority of the computational effort associ- ated wi th the chemistry is encountered during the l ibrary 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 l ibrary are entirely independent, the calculations for l ibrary genera- t ion 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 wi th condit ional moment closure using a single conditioning variable, such as extinction and ignit ion, can be captured in an unsteady simulation, a l - though there are yet to be resolved complications associated wi th capturing unsteady effects in a steady solver. In general, the model has broad applica- bi l i ty and should be useful in any circumstance where the necessary physics can be captured in the l ibrary of basis functions. In early formulations of conditional source-term estimation where the condit ional 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 condit ional 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 condit ional 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 wi l l 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 condit ional averages to be concerned with. The first is the approximately Gaussian variat ion about the conditional mean wi th 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. A l l first moment closures must assume that the former variation is small; second moment closures have been developed to address larger variations [83, 84]. The latter form of variat ion is more problematic, and causes difficulty for second moment C M C methods as well. To see why the latter form is so detrimental to first moment closures, suppose that within a given region, conditionally averaging individual points obtained from a fully resolved sample produces a conditional mean of tem- perature coinciding wi th the dashed line on the left of Figure 4.1. Further suppose that approximately half the points are clustered about the bot tom curve (solution A ) , and the remainder are clustered near the top curve (so- lut ion C ) . Th is 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 'Solut ion 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. A l though 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 . 0 2 — 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 condit ional 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 C M C method would compute the conditional reaction rates based on the conditional means indicated by the dashed lines, whilst a strictly val id closure would compute the reaction rate for each T, Yco pair of points and then conditionally average the result. For this situation, the first order C M C closure, as well as second order closures which presume a Gaussian distr ibution about the condit ional mean would be highly inaccurate. In C S E + L F D , faced wi th this situation, we have to choose the coefficients aA,<*B,ctc- Choosing = 0 , Q B = 1.0, ac = 0 would agree well wi th 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 C M C closure. However, using carbon monoxide also for inversion, the method would find a poor fit to the data for Yco- If the scalars chosen are sufficient to discriminate effectively between the basis solutions offered, the only way to obtain simultaneous agreement with both the temperature and carbon monoxide fields is to choose = 0.5, O B = 0.0, ac = 0.5. Then, the source term used wi l l be an average of the source term calculated based on (T\QA , ( ^ I O A and (T\QC , ( F j | 0 C , which is much closer to the ideal behaviour. W i t h this type of closure, the use of a second condit ioning variable is unnecessary; the physics which would normally be captured by using a second conditioning variable are implici t ly 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 w i th 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 al l possible physical conditions in the l ibrary used for de- composit ion, 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 fi ltering the un- derlying function in sufficiently varied ways that we choose between basis functions appropriately to reflect the behaviour of the condit ional 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- al ly 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 E q . 1.38, and can be wri t ten as: m in{ | | / 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 writ ten as a product of the appropriate coefficient, and the corresponding scalar value for that posit ion in basis function space; we are then try ing 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 E q . 4.2, it is apparent that the solution of E q . 4.1 for a(n*) is actually a functional dependent on (T\() (ri*), (Yi|C) (n*), • • •, (Yi\0 (ri*). 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- t ion 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 l imited 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 l ibrary of basis functions has been assembled. The presence of al l relevant physics in the l ibrary 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 inabi l i ty 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 F lame ' D ' , as well as continuing the autoignition simulations beyond ignit ion include edge flames [85], which can never be captured in a l ibrary built from simple solutions to the laminar flamelet equations. Th is highlights that although a reasonable place to start, the use of laminar flamelets as basis functions has l itt le future. Heat transfer between the flow and external objects is also of relevance in many practical flows, as are buoyancy effects — these also wi l l necessitate more comprehensive solutions for use as basis functions. Fortunately, the formulation is independent of the method used to construct the basis functions — they could be as well built from solutions constructed from C M C , or D N S for that matter. In general, C S E wi th decomposition can be viewed as a way to combine easily obtainable elementary solutions to bui ld a more realistic simulation. As more complex physics are added, it is l ikely that the number of scalars necessary for inversion may increase substantially to discriminate between the addit ional physics, as it is necessary to use a scalar for inversion that can discriminate between basis functions across al l phenomenological transitions. Once the appropriate combination of basis functions has been chosen, the mean reaction rate is obtained by virtue of an expression similar to E q . 1.29; that is: (4.3) 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 C S E methods, and it would be tr iv ia l to adjust the formulation to util ize more complex representations of the pdf. Compar ing E q . 4.3 (written for temperature) wi th the corresponding clo- sure using double-condition C M C (using scalar dissipation, wi th sample space X1 as the second conditioning variable): w « / I < w | C , X > ( C , x W x ' « 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 condit ion of being independent of Z, so that it captures variations in the flame structure not correlated wi th 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 l ibrary must be sufficiently fine grained such that a numerical integration over a' is not error-prone. Th is could be a problem where the flame structure varies rapidly wi th 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') s e e n m the integral of E q . 4.3 must be highly non-linear, so the constraint may be arduous. In Chapter 2, it was mentioned that when flamelets wi th high scalar dissi- pation were included in the l ibrary 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 possibil ity 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 t ime intervals, and the structure was changing more rapidly in flamelet t ime near extinction, the effective discretization became much coarser when the CHAPTER 4. DISCUSSION AND CONCLUSION 80 extinguishing flamelets were included. Al though 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 l ibrary 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: T A f - i E ( « n ( X o ) n ) _ n = l It is not yet clear if this addit ional constraint would produce an improve- ment — hopefully, given the assumptions of the model discussed so far, it should be impl ic i t ly satisfied based on the scalars already employed for the inversion. It may be useful as a check to evaluate E q . 4.6 and test it 's validity; if it is violated, then likely there is a violation of the earlier assumptions. Final ly, 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. A l though 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 wi l l always incorporate the a pri- ori assumption to some extent. In searching for a solution which is steady in the Euler ian sense, our choice seems logical. It wi l l provide stabil i ty for the inversion as the solution evolves, l imit the rate at which it can change X> / % ( C ) Jo / ( C R (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 t ime 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- t ion in the Lagrangian sense would have to produce a solution equal to the past solution in the Euler ian sense for a problem wi th an Euler ian steady solution. A n 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 wi th lami- nar flamelet decomposition has been presented and used successfully in the simulation of turbulent non-premixed methane ignit ion. The framework has also been used to develop a simulation of the Sandia Flame ' D ' . Th is for- mulat ion 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 l ibrary 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 l ibrary 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 ignit ion de- lay for an ignit ing methane jet, as well as obtaining a steady solution for the Sandia Flame ' D ' in the R A N S paradigm effectively demonstrates the potential of the model, and provides a framework for future development. REFERENCES 82 REFERENCES [80] C . Bla i r . Implementation of conditional source-term estimatoin for prediction of methane ignit ion. Master's thesis, University of Br i t i sh Columbia, 2003. [81] W . K. Bushe and H. Steiner. Condit ional moment closure for large eddy simulation of nonpremixed turbulent reacting flows. Physics of Fluids, 11(7):1896-1906, Ju ly 1999. [82] W . K. Bushe and H. Steiner. Laminar flamelet decomposition for condi- t ional source-term estimation. Physics of Fluids, 15(6):1564-1575, June 2003. [83] S. H. K i m , K. Y . Huh, and R. W . Bilger. Second-order condit ional 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 autoignit ion of turbulent flows. Physics of Fluids, 10(6), June 1998. Letter. [85] G . Ruetsch, L. Vervisch, and A . L inan. Effects of heat release on triple flames. Physics of Fluids, 7(6): 1447-1454, June 1995. RANS FIELDS APPENDIX A. RANS FIELDS 84 In Figure A . l , the R A N S fields showing contours of the uncondit ional mean temperature and CH creation rate are given at the time of ignit ion. 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. R A N S Temperature Field, 1150K Oxidizer R A N S C H Creation Field, 1150K Oxidizer 0.04 0.05 400 500 600 700 800 900 1000 1100 R A N S Temperature Field. 12O0K Oxidizer 0.04 0.05 400 500 600 700 8O0 900 1000 1100 1200 R A N S Temperature FieW, 1300K Oxidizer 400 500 600 700 800 900 1000 1100 1200 1300 R A N S Temperature Field, 1400K Oxidizer 0.015 _ 0.01 *0.005 0 0 0 5 1 1.5 2 2 5 3 3JS 4 4.5 R A N S C H Creation Field, 1200K Oxidizer 2 3 4 5 R A N S C H Creation Field, 1400K Oxidizer 400 600 1000 12O0 1400 Figure A . l : R A N S Temperature and C H Creation Contours at td
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
France | 5 | 2 |
Republic of Korea | 3 | 0 |
United States | 3 | 16 |
United Kingdom | 3 | 3 |
Germany | 2 | 1 |
China | 1 | 1 |
City | Views | Downloads |
---|---|---|
Unknown | 11 | 3 |
Oxford | 3 | 3 |
Beijing | 1 | 1 |
Mountain View | 1 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: