Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Combustion modeling using conditional source-term estimation with flamelet decomposition and low-dimensional… Wang, Mei 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2007-268193.pdf [ 8.1MB ]
Metadata
JSON: 831-1.0080430.json
JSON-LD: 831-1.0080430-ld.json
RDF/XML (Pretty): 831-1.0080430-rdf.xml
RDF/JSON: 831-1.0080430-rdf.json
Turtle: 831-1.0080430-turtle.txt
N-Triples: 831-1.0080430-rdf-ntriples.txt
Original Record: 831-1.0080430-source.json
Full Text
831-1.0080430-fulltext.txt
Citation
831-1.0080430.ris

Full Text

C O M B U S T I O N M O D E L I N G U S I N G C O N D I T I O N A L S O U R C E - T E R M E S T I M A T I O N W I T H F L A M E L E T D E C O M P O S I T I O N A N D L O W - D I M E N S I O N A L M A N I F O L D S b y M E I W A N G A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y i n T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Mathematics) T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A November 2006 © M e i Wang, 2006 Abstract C o m b u s t i o n m o d e l i n g is performed w i t h C o n d i t i o n a l Source-term Es-t imation (CSE) u s ing both Lamina r Flamelet Decompos i t ion (LFD) and L o w - d i m e n s i o n a l Man i fo lds . C S E w i t h Lamina r Flamelet Decompos i -t ion (LFD) is used i n the Large E d d y S imula t ion (LES) context to s tudy the non-premixed Sandia D-flame. The results show that the flame temperature and major species are w e l l predicted w i t h both steady and unsteady flamelet l ibraries. A m i x e d l ibrary composed of steady and unsteady flamelet solutions is needed to get a good predic t ion of N O . That the L F D m o d e l a l lows for tun ing of the results is found to be s ig-nificant d rawback to this approach. C S E is also used w i t h a Trajectory Generated L o w - d i m e n s i o n a l M a -n i fo ld ( T G L D M ) to simulate the Sandia D-flame. B o t h G R I - M e c h 3.0 and G R I - M e c h 2.11 are found to be able to predict the temperature and major species w e l l . However , on ly G R I - M e c h 2.11 gives a good predic-tion of N O . That G R I - M e c h 3.0 failed to give a good pred ic t ion of N O is i n agreement w i t h the findings of others i n the literature. The Stochastic Particle M o d e l (SPM) is used to extend the T G L D M to l o w temperatures where the or ig ina l con t inuum T G L D M failed. A new method for generating a trajectory for the T G L D M b y averaging different realizations together is proposed. The n e w T G L D M is used i n s imulat ions of a p remixed laminar flame and a perfectly stirred reac-tor. The results show that the new T G L D M significantly improves the predict ion. Final ly, a t ime filter is app l i ed to i n d i v i d u a l S P M realizations to i i eliminate the small time scales. These filtered realizations are tabulated into T G L D M which are then used to predict the autoignition delay time of a turbulent methane/air jet in RANS using CSE. The results are compared with shock tube experimental data. The TGLDMs incorpo-rating SPM results are able to predict a certain degree of fluctuations in the autoignition delay time, but the magnitude is smaller than is seen in the experiments. This suggests that fluctuations in the ignition delay are at least in part due to turbulent fluctuations, which might be better predicted with LES. iii Table of Contents Abstract i i Table of Contents iv List of Tables vii List of Figures viii List of Nomenclature xi Acknowledgements xiii Chapter 1. Introduction and Thesis Outline 1 1.1 Introduct ion 1 1.2 Out l ine 3 Chapter 2. Literature survey 5 2.1 Introduct ion 5 2.2 Turbulent combust ion 6 2.3 Compu ta t iona l F l u i d Dynamics 11 2.3.1 Direct N u m e r i c a l S imula t ion 11 2.3.2 Reynolds -Averaged Navier-Stokes 13 2.3.3 Large E d d y S imula t ion 15 2.3.4 H y b r i d C F D methods 17 2.4 Turbulent combust ion mode l i ng 18 2.4.1 L a m i n a r flamelet m o d e l 20 2.4.2 C o n d i t i o n a l M o m e n t Closure 23 2.4.3 C o n d i t i o n a l Source-term Es t imat ion 27 2.5 S impl i f ica t ion of reaction mechanisms 31 2.5.1 Q S S A and part ial equ i l i b r i um 31 2.5.2 Computa t iona l singular perturbat ion 35 2.5.3 Intrinsic L o w - D i m e n s i o n a l M a n i f o l d 36 2.5.4 Trajectory Generated L o w - D i m e n s i o n a l M a n i f o l d . . 40 2.5.5 I C E - P I C 43 2.6 Probabil is t ic descript ion of combust ion 44 2.6.1 Probabi l i ty Densi ty Func t ion (PDF) m o d e l 45 2.6.2 Stochastic Particle M o d e l 48 2.7 S u m m a r y 56 i v Chapter 3. Conditional source-term estimation with laminar flame-let decomposition in large eddy simulation of a turbulent non-premixed flame 57 3.1 Introduction 57 3.2 Formulation 57 3.2.1 Large Eddy Simulation 57 3.2.2 Laminar Flamelet Decomposition 58 ^3.2.3 Flamelet libraries 60 3.3 Simulation results and discussion 61 3.3.1 Steady Flamelet Library Results 61 3.3.2 Unsteady Flamelet Library Results 71 3.3.3 Mixed Library 76 3.4 Conclusion 78 Chapter 4. Simulation of a turbulent non-premixed flame using conditional source-term estimation with trajectory gen-erated low-dimensional manifold 80 4.1 Introduction 80 4.2 Models : 80 4.2.1 CSE-TGLDM 80 4.3 Implementation 82 4.3.1 Large Eddy Simulation ' 82 4.3.2 T G L D M 82 4.3.3 Experimental Setup 84 4.4 Results and discussion 85 4.4.1 Computational cost 85 4.4.2 Temperature and major species 87 4.4.3 Minor species 92 4.4.4 NO prediction 94 4.5 Conclusion 95 Chapter 5. Extending the trajectory generated low-dimensional manifolds to low temperature with the stochastic particle model 99 5.1 Introduction 99 5.2 Calculating trajectories with the SPM 101 5.2.1 Average trajectories 104 5.3 Validation 107 5.3.1 Laminar premixed flame 107 5.3.2 Perfectly stirred reactor calculations I l l 5.4 Conclusions 114 v Chapter 6. Prediction of autoignition fluctuation using trajectory generated low-dimensional manifolds generated using the stochastic particle model 115 6.1 Introduction 115 6.2 Results and discussion 116 6.2.1 Filtering individual realizations 116 6.2.2 Predicting fluctuations in autoignition 121 6.3 Conclusions • 135 Chapter 7. Conclusions and Future Work 137 7.1 Conclusions 137 7.2 Future Work 139 Bibliography 141 Appendix A . Inverting an integral equation in L F D 155 vi List of Tables 4.1 Break-down of computational cost in simulation 85 5.1 Computation cost of for averaging temperature T = 1200K .. 105 6.1 Experimental conditions of shock tube 124 v i i List of Figures 2.1 Borgh i d i ag ram for non-premixed combust ion 9 2.2 Regimes i n non-premixed turbulent combust ion 10 2.3 Exper iment measurement of scattered and condi t iona l aver-age temperature of Sandia D-flame at x/d = 30 24 2.4 Time scales i n reacting f low 37 2.5 C o m p a r i n g one d imens iona l mani fo ld generated b y I L D M and T G L D M 42 2.6 So lu t ion of a Hypo thes i s s imple system 51 3.1 Three flamelets solutions of temperature, C O and tempera-ture source term 59 3.2 Center l ine profile of temperature and mixture fraction 62 3.3 R a d i a l profiles of mixture fraction and temperature calculat-ed w i t h S F L 63 3.4 R a d i a l profiles of H 2 0 mass fraction calculated w i t h S F L . . . . 65 3.5 C o n d i t i o n a l averages of H 2 0 calculated w i t h S F L 66 3.6 C o n d i t i o n a l averages of C O and CO2 mass fractions calculat-ed w i t h S F L 67 3.7 C o n d i t i o n a l averages of C H 4 and O2 mass fractions calculat-ed w i t h S F L 68 3.8 C o n d i t i o n a l averages of O H and H 2 mass fractions calculat-ed w i t h S F L 69 3.9 C o n d i t i o n a l averages of N O calculated w i t h S F L 70 3.10 C o n d i t i o n a l averages of H 2 0 calculated w i t h U F L 71 3.11 C o n d i t i o n a l averages of C O and C 0 2 mass fractions calculat-ed w i t h U F L 72 3.12 C o n d i t i o n a l averages of C H 4 and 0 2 mass fractions calculat-ed w i t h U F L 73 3.13 C o n d i t i o n a l averages of O H and H 2 mass fractions calculat-ed w i t h U F L 74 3.14 C o n d i t i o n a l averages of N O calculated w i t h U F L 75 3.15 C o n d i t i o n a l averages of N O mass fraction calculated w i t h m i x -ed flamelet l ib ra ry 77 4.1 T G L D M l ibrary at stoichiometric mixture fraction 84 4.2 Instantaneous temperature field 86 4.3 Center l ine profile of temperature and mixture fraction 87 4.4 Temperature and mixture fraction rad ia l profile 88 v i i i 4.5 C o n d i t i o n a l averaged CO2 and H 2 0 mass fraction at differ-ent locations downs t ream 89 4.6 Center l ine profile of mass fraction of CO2 and H 2 0 90 4.7 Mass fraction of C 0 2 and H 2 0 radia l profile at different loca-tions downs t ream 91 4.8 C o n d i t i o n a l species mass fractions f rom interpolat ing T G L D M l ib ra ry 92 4.9 C o n d i t i o n a l species mass fractions from interpolat ing T G L D M l ibrary 93 4.10 T G L D M l ibrary of N O mass fraction at stoichiometric m i x -ture fraction 95 4.11 C o n d i t i o n a l averaged mass fraction of N O at different loca-tions downs t ream us ing G R I - M e c h 3.0 96 4.12 C o n d i t i o n a l averaged mass fraction of N O at different loca-tions downs t ream us ing G R I - M e c h 2.11 97 5.1 C o m p a r i s o n of c o n t i n u u m and S P M w i t h an in i t i a l tempera-ture of 1200K 102 5.2 C o m p a r i s o n of con t inuum and S P M w i t h an in i t i a l tempera-ture of 298K 103 5.3 Piecewise least squares fit of S P M w i t h an in i t i a l tempera-ture 1200K 106 5.4 Piecewise least squares fit w i t h l o w temperature (T=298K) i n the C 0 2 and H 2 0 phase space 107 5.5 Scalar profiles through the laminar p remixed flame w i t h i n i -t ia l reactant temperature of 298K 109 5.6 Scalar profiles through the laminar p remixed flame w i t h i n i -t ial reactant temperature of 298K 110 5.7 P S R calculat ion results: temperature and mass fractions of C 0 2 a n d O H 113 6.1 The wa i t t ime (ms) for every chemical reaction event as a function of the event number for a single real izat ion of the S P M w i t h an in i t i a l temperature of 298K 117 6.2 The wai t t ime (ms) for every chemical reaction event as a function of the C 0 2 mass fraction for a single real izat ion of the S P M w i t h an in i t i a l temperature of 298K 117 6.3 React ion rates filtered w i t h different filter size 120 6.4 Ten realizations of S P M w i t h an in i t ia l temperature of 1200K compared to the con t inuum trajectory 122 6.5 Profile of the jet at t = 1.30ms 126 6.6 Profile of the jet at t = 1.56ms 127 6.7 Profile of the jet at t = 1.73ms (with real izat ion 11) 128 ix 6.8 Profile of the jet at t = 1.73ms (wi th real izat ion 12) 129 6.9 Var ia t ion i n ign i t ion delay times of a transient methane jet injected into air at an in i t i a l temperature of T = 1150K as calculated w i t h R A N S us ing the C S E / T G L D M method . . . . 130 6.10 The chemical source term of the mass fraction of C 0 2 as a function of the mass fraction of O H d u r i n g autoigni t ion of a stoichiometric mixture of methane /a i r w i t h an in i t i a l tem-perature of 1150K at 30-bar , . . . . 132 6.11 C o m p a r i s o n of the ign i t ion delay t ime of methane jets issu-i n g into quiescent air calculated w i t h filtered S P M generated trajectories w i t h experiment data for an in i t ia l air tempera-ture T ~ 1150K '. . . . 1 3 4 x List of Nomenclature Symbol Meaning Units D Diffus iv i ty 2 —1 m • s Da D a m k o h l e r number 2 —1 m • s E A c t i v a t i o n Energy c a l - m o l - 1 F B o d y force h Specific En tha lpy I Identi ty M a t r i x J Jacobian M a t r i x 3h Entha lpy diffusivi ty J / (m 2 -s ) 3k Species molecular diffusivi ty k g / ( m 2 - s ) k React ion Rate S " 1 Ka K a r l o v i t z N u m b e r N o n - p r e m i x e d reference length scale m V Pressure bar P Probabi l i ty Densi ty Func t ion P Projection M a t r i x Pr Prand t l number Re Reynolds number Sc Schmidt number tc C h e m i c a l t ime scale s */ Integral t ime scale s K o l m o g o r o v time scale s Uref N o n - p r e m i x e d reference veloci ty m / s V veloci ty vector m / s V Eigen(Schur) Basis M a t r i x S Sensi t ivi ty T Temperature K W React ion rate s - 1 X Cartes ian coordinates Y Mass Fract ion Z M i x t u r e Fract ion (Y\Z = () C o n d i t i o n a l Average of Y w i t h C o n d i -t ion Z = ( x i Symbol Meaning Units Z"2 Variance of Mix tu re Fract ion C M C C o n d i t i o n a l M o m e n t Closure C S E C o n d i t i o n a l Source-term Est imat ion D N S Direct N u m e r i c a l S imula t ion F C T F l u x Corrected Transport I L D M Intr insical ly L o w Dimens iona l Man i f -o l d L E S Large E d d y S imula t ion L F D L a m i n a r Flamelet Decompos i t ion P S R Perfectly Stirred Reactor Q S S A Quasi-steady State A p p r o x i m a t i o n R A N S Reynolds averaged Navier-Stokes S F L Steady Flamelet L ib ra ry S P M Stochastic Particle M o d e l T G L D M Trajectory-Generated L o w Dimens io -na l M a n i f o l d U F L Uns teady Flamelet L ib ra ry A D iagona l M a t r i x of Eigenvalues s - 1 P Dens i ty k g - m ~ 3 T Igni t ion De l ay or Time Dura t i on /is,ms or s to Phys ica l t ime scale for P S R S " 1 LO C h e m i c a l Source Term s - 1 X Scalar Diss ipa t ion Rate s - 1 C Sample space for mixture fraction x i i Acknowledgements D u r i n g four years of m y P h D study, I met m a n y friends w h o have g iven me good memories . I w o u l d l ike to take this oppor tun i ty to thank a l l w h o he lped me through m y P h D study. I w o u l d l ike to thank the financial support f rom the U n i v e r s i t y of Br i t i sh C o l u m b i a a n d the N a t i o n a l Research C o u n c i l of Canada . I w o u l d l ike to thank Drs . Steve Rogak, C a r l F. O l l i v i e r - G o o c h and M i c h a e l J. W a r d for agreeing to.be members of m y research committee. I w o u l d l ike to thank a l l m y colleagues and friends i n Rus ty H u t 101 and Green Col lege . Thanks to Dr . Jian H u a n g for the discussions and valuable suggestions. Thanks to G o r d McTagga r t -Cowan for proof-reading part of the thesis. Thanks Dr . A n d r e a Frisque, Heather Jones and Bei J ing for precious comradery and help i n a l l aspects. I w o u l d l ike to thank m y supervisor, Dr . K e n d a l Bushe, for his f i rm support i n m a n y aspects of m y P h D study. Wi thou t h i m , this thesis w o u l d not be possible. I w o u l d l ike to express m y special gratitude to m y fami ly members - m y parents, sisters, sister-in-laws and nephews, for their f i rm sup-port on m y w a y of p u r s u i n g higher degrees. Final ly, I w o u l d l ike to thank m y boyfr iend, Feng X i a , whose deep unders tanding and support is greatly appreciated. x i i i C h a p t e r 1 Introduction and Thesis Outline 1.1 Introduction We are now facing serious environment problems and a severe energy shortage. Since 2000, oil prices have increased dramatically, first as a re-sult of strong demand growth due to economies booming in countries such as China and as a result of disruption i n the M i d d l e East which is the largest oil supply region. A record price of $78.40 per barrel was reached on July 13, 2006. Large effort have been put into searching for alternative energy resources. However, as the energy consumption detailed in US department of energy 2006 overview show, Petroleum products (38.96%), natural gas (23.04%) and coal (22.38%) remain the major sources [2]; the technologies used to generate power from these energy sources virtually all make use of combustion. Pollutants and green-house gases generated from combustion pro-cesses threaten the environment. Combustion produces pollutants such as oxides of nitrogen (NO x ) , soot, and unburnt hydrocarbons (HC); these pollutants can have significant impact on our health. In addition, the Environmental Protection Agency (EPA) has declared that carbon dioxide is a pollutant which contributes to global warming. More strin-gent regulations are being proposed by EPA for future years which 1 seek to protect the environment; these force manufacturers of energy convers ion devices to look for ways to reduce the pollutants . Scientists and researchers continue to search for w a y s to improve the efficiency and reduce pol lutant emissions f rom the b u r n i n g of these energy resources. The design of efficient and l o w emi t t ing equipment needs a thorough unders tanding of the combust ion process. W h i l e ex-periments are useful tools for testing the final product of the design process and the emergence of new technologies such as non-intrusive diagnostic techniques a l l ow researchers to explore phenomena not ac-cessible b y experiment before,,these are relat ively expensive and l i m -i ted as useful design tools. M o d e l i n g is more controllable for gener-at ing exact condi t ions b y changing parameters. Therefore it is n o w w i d e l y used as a des ign tool i n add i t ion to experiments. However , combus t ion mode l i ng is a compl ica ted a n d interdisci-p l ina ry subject. M o s t combus t ion takes place w i t h i n a turbulent rather than a laminar f low field. Turbulence itself is ve ry difficult to m o d e l : due to the large variety of scales, it is usual ly not feasible to solve a l l the scales (as one w o u l d do i n direct numer ica l s imulat ion) w i t h cur-rent computat ional capabilities. Instead, it is usua l ly necessary to solve for an averaged or filtered f low field, i n w h i c h case models need to be p rov ided for unclosed terms. Descr ib ing the chemical kinet ic mecha-n i s m is another challenge; detailed chemistry is an ongo ing research field to extend available chemical mechanisms f rom certain s imple fu-els to other more complex fuels b y matching available experiment data. A detailed chemical kinetic mechanism can easily inc lude tens or hun-dreds of species and hundreds or thousands of reactions. These reac-2 t i ons o c c u r w i t h a l a r g e r a n g e o f d i f fe ren t t i m e sca les w h i c h l e a d s to a st iff O D E s y s t e m tha t i s d i f f i c u l t to s o l v e . F u r t h e r m o r e , i n cases w h e r e d i r ec t n u m e r i c a l s i m u l a t i o n is n o t p o s s i b l e , the r e s u l t i n g t u r b u l e n c e a n d c h e m i s t r y i n t e r a c t i o n s h a v e to be m o d e l e d . C l a s s i c a l c o m b u s t i o n m o d e l s a s s u m e that the r e a c t i n g f l u i d c a n b e t rea ted as a c o n t i n u u m a n d so c a n n o t p r e d i c t m i c r o s c o p i c p h e n o m e n o n . H o w e v e r , the u n d e r l y i n g f l u c t u a t i o n s c a n be v e r y i m p o r t a n t to c o n -s i d e r d u r i n g the d e s i g n p r o c e s s for d e v i c e s tha t are s e n s i t i v e for i g n i -tion d e l a y t i m e f l u c t u a t i o n s tha t are w i d e l y o b s e r v e d i n e x p e r i m e n t s . M i c r o s c o p i c s i m u l a t i o n is e x t r e m e l y c o m p u t a t i o n a l l y e x p e n s i v e ; e v e n m e s o s c o p i c p a r t i c l e m o d e l i n g i s n o t feas ib le for a p p l i c a t i o n s c l o s e to p r a c t i c a l u s a g e (e.g. t he s i m u l a t i o n o f a jet). T h e r e f o r e , s i m u l a t i o n o f the s m a l l f l u c t u a t i o n s i n c o m b u s t i o n i s a b i g c h a l l e n g e . T h i s thes i s focuses o n the s t u d y o f c o m b u s t i o n m o d e l i n g w i t h af-f o r d a b l e cost a n d h i g h e r accuracy . It add re s se s s e v e r a l p r o b l e m s o f c o m b u s t i o n m o d e l i n g m e n t i o n e d a b o v e ; c l o s u r e for the m e a n r e a c t i o n rate , r e d u c e d c h e m i c a l m e c h a n i s m s a n d a p r o b a b i l i s t i c d e s c r i p t i o n o f the f l u i d are a l l e x p l o r e d . M o s t o f the s i m u l a t i o n cases i n v o l v e n o n -p r e m i x e d c o m b u s t i o n . A l l s t u d y cases use m e t h a n e - w h i c h i s the m a -jo r c o m p o n e n t o f n a t u r a l gas as fue l ; n eve r the l e s s , these m o d e l s c a n be e x t e n d e d to a n y g a s e o u s f u e l w i t h r e l a t i v e l y l o w effort . 1.2 Outline C h a p t e r 2 p r o v i d e s a r e v i e w of the l i t e r a tu re r e l e v a n t to the w o r k d e -s c r i b e d l a t e r i n t h i s thes i s . It f i rs t r e v i e w s d i f fe ren t C o m p u t a t i o n a l F l u i d D y n a m i c s ( C F D ) m e t h o d s . T h e n d i f fe ren t c o m b u s t i o n m o d e l s are re-3 v iewed . Me thods for reducing chemical reaction mechanisms are dis-cussed. The chapter concludes w i t h a survey of probabil is t ic descrip-t ion of combust ion. In Chapter 3 results are presented us ing C o n d i t i o n a l Source-term Est imat ion (CSE) w i t h L a m i n a r Flamelet Decompos i t ion (LFD) i n the Large E d d y S imula t ion (LES) of the Sandia D-flame. The w o r k i n this chapter has been submit ted for publ ica t ion i n Physics of F l u i d s as " C o n -di t iona l Source-term Es t imat ion w i t h Lamina r Flamelet Decompos i -tion i n Large E d d y S imula t ion of a Turbulent N o n - p r e m i x e d F lame" by M . W a n g and W. K . Bushe. In Chapter 4, C S E is used w i t h a Trajectory Generated L o w - D i m e n -sional M a n i f o l d ( T G L D M ) to simulate the same flame. A paper based on this w o r k has been pub l i shed as "S imula t ion of a Turbulent N o n -premixed F lame us ing Cond i t i ona l Source-term Es t imat ion w i t h Tra-jectory Generated L o w - d i m e n s i o n a l M a n i f o l d " b y M . W a n g , J. H u a n g and W. K . Bushe i n the 2006 Proceedings of C o m b u s t i o n Institute. In Chapter 5, the Stochastic Particle M e t h o d (SPM) is used to ex-tend the T G L D M to l o w temperatures where the O D E system is too stiff and the or ig ina l T G L D M method failed. In Chapter 6, the S P M is used i n combina t ion w i t h the T G L D M method to predict the fluctua-tions i n the of autoigni t ion delay time. The w o r k i n Chapters 5 and 6 has been combined and submit ted for publ ica t ion as "Trajectory G e n -erated L o w - d i m e n s i o n a l Man i fo lds generated us ing the Stochastic Par-ticle M o d e l " b y M . Wang , A . Frisque, J. H u a n g and W. K . Bushe i n the Combus t i on Theory and M o d e l l i n g . 4 Chapter 2 Literature survey 2.1 Introduction N u m e r i c a l s imula t ion of reacting flows is an impor tant research area because it contributes to our unders tanding of combust ion. It is also w i d e l y used i n indus t ry as a w a y to design more efficiently and to re-duce the need for expensive experiments. M o s t reacting f lows of practi-cal interest are turbulent. Turbulent combust ion is difficult to s imulate because it is an in terdisc ip l inary combina t ion of turbulent f l u i d dy-namics and chemistry. Turbulence itself poses a significant challenge: the existence of a w i d e variety of scales makes the direct numer ica l s imula t ion of turbulent f low fields w i t h h i g h Reyno lds number com-putat ional ly expensive. Chemis t ry is another major challenge: not on ly can there be a large number of species and reactions i n v o l v e d , but these lead to a large range of timescales, result ing i n a complex, stiff system that is computa t ional ly difficult to solve. Enormous effort has been put into reducing the complexi ty of chemical kinetics. W h e n consider ing the s imula t ion of turbulent combust ion, unless one uses direct numer-ical s imula t ion - w h i c h is extremely costly and can on ly be app l i ed to relat ively s imple f lows - there is an unclosed term i n the transport equation w h i c h is the average reaction rate w h i c h needs to be m o d -5 e l e d . B e c a u s e the r e a c t i o n rate i s a h i g h l y n o n l i n e a r f u n c t i o n o f m a s s f r a c t i o n o f spec i e s a n d t e m p e r a t u r e , s u b s t i t u t i n g the m e a n m a s s frac-t ions a n d t e m p e r a t u r e i n t o the A r r h e n i u s r e a c t i o n rate e x p r e s s i o n w i l l cause a n u n a c c e p t a b l y l a r g e error . T h e c l a s s i c d e t e r m i n i s t i c o r c o n t i n u u m a p p r o a c h s o l v e s the react-i n g s y s t e m o n m a c r o s c o p i c scales . F r o m the m i c r o s c o p i c sca le p o i n t o f v i e w , m o l e c u l e s i n a c h e m i c a l s y s t e m r a n d o m l y c o l l i d e a n d react w i t h c e r t a i n p r o b a b i l i t i e s . I d e a l l y , o n e w o u l d w a n t to s i m u l a t e the d y n a m -ics o f e v e r y m o l e c u l e , as o n e does i n M o l e c u l a r D y n a m i c s ( M D ) [3]; u n f o r t u n a t e l y , M D is n o t y e t c a p a b l e o f h a n d l i n g the c o m p l e x i t y o f the c h e m i s t r y n e e d e d for c o m b u s t i o n p r o b l e m s i n a n e n g i n e e r i n g con tex t . A s a c o m p r o m i s e , the S tochas t i c P a r t i c l e M o d e l ( S P M ) i s a m e s o s c o p i c a p p r o a c h w h e r e the s y s t e m i s d i s c r e t i z e d i n t o f in i t e v o l u m e s i n w h i c h the n u m b e r o f m o l e c u l e s o f e a c h spec ies a re c a l c u l a t e d . T h e s tochas -tic a p p r o a c h c a n a c c o u n t for the i n h e r e n t f l u c t u a t i o n s a n d c o r r e l a t i o n s tha t are i g n o r e d i n the d e t e r m i n i s t i c a p p r o a c h . T h i s c h a p t e r w i l l f i rs t r e v i e w d i f fe ren t C F D m e t h o d s . T h e n c o m -b u s t i o n m o d e l s i n c l u d i n g the l a m i n a r f l a m e l e t m o d e l , C o n d i t i o n a l M o -m e n t C l o s u r e a n d C o n d i t i o n S o u r c e - t e r m E s t i m a t i o n w i l l b e r e v i e w e d . M e t h o d s for r e d u c i n g the c h e m i c a l r e a c t i o n m e c h a n i s m w i l l b e d i s -c u s s e d . T h e c h a p t e r w i l l f i n i s h w i t h a s u r v e y o f t he p r o b a b i l i s t i c m o d -e l i n g a p p r o a c h e s . 2.2 T u r b u l e n t c o m b u s t i o n M o s t r e a c t i n g f l o w s o f p r a c t i c a l in te res t are t u r b u l e n t . T u r b u l e n t c o m -b u s t i o n i n v o l v e s v a r i o u s l e n g t h , v e l o c i t y a n d t i m e sca les . T h e s e p a -6 r a m e t e r s c a n b e u s e d to p l o t a d i a g r a m w h i c h n o t o n l y g i v e s a n i n -s i g h t i n t o the f l o w b u t a l s o p r o v i d e s a g u i d e l i n e for the f i tness o f a c e r t a i n m o d e l . D i f f e r e n t d i a g r a m s h a v e b e e n d e v e l o p e d for p r e m i x e d a n d n o n - p r e m i x e d f l a m e s . A s a l l t u r b u l e n t f l o w s s t u d i e d i n t h i s thes is are n o n - p r e m i x e d ones , o n l y n o n - p r e m i x e d d i a g r a m s are r e v i e w e d i n th i s l i t e r a tu r e s u r v e y . V a r i o u s c h a r a c t e r i s t i c scales h a v e b e e n d e f i n e d for the d i f f e ren t d i -a g r a m s d e v e l o p e d . S o m e u s e a R e y n o l d s n u m b e r to c h a r a c t e r i z e the t u r b u l e n t f l o w a n d a D a m k o h l e r n u m b e r to c h a r a c t e r i z e the r e a c t i o n z o n e [4]. A n o t h e r o p t i o n i s to u se a p a i r o f r a t i o s to d i v i d e the d i a g r a m : the r a t i o o f the r e l a t i v e f l u c t u a t i o n s at a m e a n s t o i c h i o m e t r i c m i x t u r e f r a c t i o n to the d i f f u s i v e t h i c k n e s s i n m i x t u r e f r a c t i o n space , a n d the time sca le r a t i o o f e x t i n c t i o n sca la r d i s s i p a t i o n ra te to the c o n d i t i o n a l F a v r e m e a n sca l a r d i s s i p a t i o n ra te [5]. A l t e r n a t i v e l y , the B o r g h i d i a -g r a m de f ines a v e l o c i t y r a t i o a n d a l e n g t h r a t i o to d e l i n e a t e b e t w e e n r e g i m e s [6]. T h e B o r g h i d i a g r a m w a s f i rs t u s e d fo r p r e m i x e d c o m b u s -t i o n , u s i n g p a r a m e t e r s c o m p o s e d o f a n o n c o m i n g t u r b u l e n c e v e l o c i t y n o r m a l i z e d b y the l a m i n a r b u r n i n g v e l o c i t y a n d the l a rges t t u r b u l e n t e d d y l e n g t h sca le n o r m a l i z e d b y the l a m i n a r f l a m e t h i c k n e s s . H o w -ever , a d i f f u s i o n f l a m e does n o t p r o p a g a t e ; therefore i t h a s n o c h a r a c -te r i s t ic f l a m e s p e e d . It i s a l s o d i f f i c u l t to d e f i n e f i x e d l e n g t h sca les for a d i f f u s i o n f l a m e . L e n t i n i [7] s u g g e s t e d r e v i s e d d e f i n i t i o n s to d e l i n e a t e a B o r g h i d i a -g r a m for a d i f f u s i o n f l a m e . H e d e f i n e d the reference l e n g t h sca le as: Iref = \/vtc ( 2 . 1 ) 7 a n d reference v e l o c i t y as: U r e f = i - (2.2) w h e r e tc d e n o t e s the c h e m i c a l t i m e sca le a n d v d e n o t e s the k i n e m a t i c v i s c o s i t y . A n o n - p r e m i x e d reference t i m e sca le c a n t h e n b e d e f i n e d as Uef = T h e d i a g r a m i s c r ea t ed u s i n g log(^~) ( a n d logiy1--) as the axes. S e v e r a l cons t an t s are d e f i n e d to d i v i d e the d i a g r a m i n t o d i f fe ren t r e g i o n s , as s h o w n i n F i g . 2 .1 . W h e n the t u r b u l e n t R e y n o l d s n u m b e r is s m a l l , a l a m i n a r f l a m e i s o b s e r v e d . T h e d o m a i n o f t u r b u l e n t c o m b u s -t i o n i s fu r the r d i v i d e d i n t o three z o n e s . A t u r b u l e n t K a r l o v i t z n u m b e r , Ka, i s d e f i n e d as the r a t i o b e t w e e n the c h e m i c a l t i m e sca le o f the l a m i -n a r f l a m e to the K o l m o g o r o v t i m e sca le (tk): Ka = ^ . (2.3) W h e n the s m a l l e s t sca le o f t u r b u l e n c e , the K o l m o g o r o v l e n g t h , i s l a r g e r t h a n the t h i c k n e s s o f the r e a c t i o n z o n e , o r w h e n the c h e m i c a l t i m e sca le is s m a l l e r t h a n the K o l m o g o r o v time sca le , Ka < 1, the s y s t e m is t h e n d e s c r i b e d as a l o c a l l y l a m i n a r f l a m e e m b e d d e d i n a t u r b u l e n t f l o w . T h i s z o n e i s s h o w n as " i s l a n d f o r m a t i o n a n d w r i n k l e d l a m i n a r f l ame , f l a m e l e t r e g i m e " o n F i g . 2.1. T h e l a m i n a r f l a m e l e t m o d e l , w h i c h w i l l be d e s c r i b e d later , c a n be u s e d for f l a m e s i n t h i s r e g i m e . T h e tu r -b u l e n t D a m k o h l e r n u m b e r (Da) is d e f i n e d as the r a t i o o f the i n t e g r a l f l o w t i m e sca le to the t i m e sca le o f the c h e m i c a l r e a c t i o n , Da=£- (2.4) w h e r e tj d e n o t e s the i n t e g r a l t i m e scale . W h e n Da < 1, w h i c h i s the r e g i o n a b o v e the l i n e Da = 1 o n F i g . 2 .1 , s m a l l sca le e d d i e s p e n e t r a t e 8 Figure 2.1: Borgh i d iagram for non-premixed combus t ion and disrupt the chemical reaction zone and the flamelet concept is not applicable. This is cal led the "thickened flame regime". Between the thickened flame regime and the flamelet zone is the " torn flame front" region. In this region, the turbulence interacts w i t h the flame but the f low does not completely o v e r w h e l m the flame. Lent in i ' s Borgh i d i ag ram for non-premixed combus t ion does not consider the effect of m i x i n g . Peters [5] defined the non-premixed tur-bulent combust ion regimes b y compar ing the ratio of Z'st/(AZ)p and the t ime scale ratio Xg/Xst i n a log- log p lo t s h o w n i n F i g . 2.2. Here , Z'st is the fluctuation at the mean stoichiometric mixture , (AZ)p is the diffusion thickness i n mixture fraction space defined as: (AZ)F = \VZ\st£D, (2.5) 9 0.1 0.01 para ted (1 ante lets n ^ flame extinction \ .... connected flame zones connected reaction zones : i H&Z)K= V 1 10 KM) Xt// Xsi Figure 2.2: Regimes i n non-premixed turbulent combus t ion [5] where £Q is diffusion thickness i n phys ica l space defined by: to = ( ^ ) 1 / 2 , (2-6) and D S T is the diffusion coefficient, a is the strain rate; Xq is the ex-t inction scalar d iss ipa t ion rate, Xst is the condi t ional Favre mean scalar diss ipat ion rate. Several lines are defined o n the d i a g r a m to separate different zones. The first l ine ' i s at Z'ST/(AZ)p = 1. For large mixture fraction fluctuations, where Z'ST > (AZ)p, w h i c h is the zone above the Z'ST/(AZ)p = 1 l ine , fluctuations i n mixture fraction space extend to sufficiently lean a n d r i ch mixtures that the diffusion layers a round the reaction zone are separated. W h e n the mixture fraction fluctuations are smal l (Z'ST < (AZ)F/ w h i c h is be low the Z'J(AZ)F = 1 l ine), the reac-t ion zones are connected. Ano the r l ine defined b y Z'ST/(AZ)R = 1 is also plotted, where (AZ)R = 1 is the reaction zone thickness. This l ine 10 corresponds to mixture fraction fluctuations that are equal to the thick-ness of the reaction zone. O n Fig . 2.2 it is a l ine w i t h slope of —1/4 [5]. Be low this l ine the mixture fraction fluctuations are smal ler than the reaction zone thickness and the reaction zones are connected, w h i c h means the mixture fraction field is essentially homogeneous. Ano the r l ine is plot ted at Xq = Xst w h i c h is the vert ical l ine i n F ig . 2.2. W h e n the scalar diss ipat ion is larger than the extinct ion scalar dissipat ion, w h i c h is the zone left of the vert ical l ine, the flame is ex-t inguished. A diffusion flame is d r a w n on F i g . 2.2 to illustrate the dif-ferent zones of a typical diffusion jet. 2.3 Computational Fluid Dynamics Flows are governed b y equations describing the conservation of mass and energy and the Navier-Stokes equations, w h i c h describe the con-servation of m o m e n t u m . It is very difficult to f ind exact analyt ical solu-tions for these equations i n f low fields of interest. Instead, a numer ica l solut ion is usua l ly used. D u e to the large variety of scales i n turbu-lence, it is computa t ional ly difficult to resolve i n the smallest scales; there are var ious methods that seek to p rov ide efficient and accept-ably accurate results. There are three basic ways to s imulate turbu-lent f lows: Direct N u m e r i c a l S imula t ion (DNS) , Reyno lds -Averaged Navier-Stokes ( R A N S ) methods and Large E d d y S imula t ion (LES). 2.3.1 Direct Numerical Simulation A l l f lows are fu l ly described b y the basic set of balance equations i n -c l u d i n g the continuity, m o m e n t u m , species and energy transport equa-tions [8]: 11 • Mass conservation (continuity): % + V • ( P v ) = 0 (2-7) where p and v denote density and veloci ty respectively. • M o m e n t u m : dt + V • (pvv) = - V p + V r + F (2.8) where r denotes the viscous stress tensor and F denotes the b o d y force. • Species transport equation: dpYk dt + V - ( p v y f c ) = - V Jf c + ^ (2.9) where Yk denotes species mass fraction, J f c denotes species diffu-sive f lux a n d iok is species chemical source-term. • Energy equation: d-~ + V • (pvh) = ~ + V • (Jh + v r ) + v F (2.10) where h denotes total enthalpy and Jh denotes enthalpy diffu-sion. For N e w t o n i a n f luids , r is defined b y Newton ' s l a w : T = / i ( V v + ( V v ) r - ^ ( V - v ) I ) (2.11) where the molecular viscosity, / i , is a property of the f lu id . The flux Jk can be approximated w i t h F ick ' s l aw: Jk = -pDkS7Yk (2.12) 12 w h e r e Dk i s t he m o l e c u l a r d i f f u s i o n coef f ic ien t o f the spec i e s k. T h e e n t h a l p y d i f f u s i o n , 3h, c a n b e a p p r o x i m a t e d w i t h F o u r i e r ' s l a w : k=l k w h e r e Pr i s t he P r a n d t l n u m b e r a n d Sck i s the S c h m i d t n u m b e r o f spec ies k. T h e So re t effect a n d m o l e c u l a r t r a n s p o r t d u e to p r e s s u r e g r a d i e n t s are t y p i c a l l y n e g l e c t e d , a n d a u n i t y L e w i s n u m b e r i s o f ten a s s u m e d to s i m p l i f y m o d e l i n g . D N S s o l v e s e q u a t i o n s E q . 2.7 to E q . 2.10 d i r e c t l y w i t h g r i d s su f f i -c i e n t l y f ine as to s o l v e the s m a l l e s t e d d y w i t h o u t the n e c e s s i t y o f a n y ex t r a m o d e l s . H o w e v e r , to r e s o l v e a l l the scales i n the f l o w , e s p e c i a l l y i n p r a c t i c a l h i g h R e y n o l d s n u m b e r t u r b u l e n t f l o w s , t he m e s h h a s to b e v e r y f ine . A c c o r d i n g to K o l m o g o r o v ' s l a w , the s i z e o f p e r s i s t i n g e d d i e s i s 0 ( i ? e ~ 3 / 4 ) i n three d i m e n s i o n s . F o r a f l o w w i t h Re « 1 0 6 , the m e s h s i ze h a s to be o n the o r d e r o f 10~2 to r e s o l v e a l l the sca les . T h e r e f o r e , the to ta l n u m b e r s o f m e s h ce l l s i n a h e x a h e d r o n u n i t n e e d to b e o n the o r d e r o f 1 0 1 3 [9]. W i t h p r a c t i c a l a p p l i c a t i o n s h a v i n g e v e n h i g h e r R e y n o l d s n u m b e r s , the c o m p u t a t i o n a l cost far e x c e e d s c u r r e n t c o m -p u t a t i o n a l c a p a b i l i t i e s . A n o t h e r d i s a d v a n t a g e o f D N S i s tha t t he i n i -tial a n d b o u n d a r y c o n d i t i o n s are o f t en d i f f i c u l t to d e f i n e fo r p r a c t i c a l f l o w s . T h e r e f o r e , D N S is m o s t l y u s e d for r e l a t i v e l y s i m p l e cases , s u c h as h o m o g e n e o u s i s o t r o p i c t u r b u l e n t f l o w , t u r b u l e n t b o u n d a r y l a y e r s a n d c h a n n e l f l o w [10-12] . 2.3.2 Reynolds-Averaged Navier-Stokes R e y n o l d s - a v e r a g e d N a v i e r - S t o k e s ( R A N S ) m o d e l s o l v e s for the t i m e -o r e n s e m b l e - a v e r a g e d v a l u e s . E a c h q u a n t i t y Q i s d e c o m p o s e d i n t o a 13 mean and a devia t ion: Q = Q + Q' (2.14) where Q = 0. To a v o i d explici t mode l i ng of the veloci ty and density correlation, a Favre average Q = & can be in t roduced such that: Q = Q + Q" (2.15) w i t h Q" = 0. A p p l y i n g the Favre average to the basic balance equations yields: ^ + V • (pv) = 0 (2.16) ^ + V • (pvv) = - V ( p V V ' ) - V p + V r + F (2.17) ^ + V • (pvYk) = -V{-p^YJ!) - V ^ + uTk (2.18) ^ + V • (pvh) = - V • ( p v ^ f ) + ~ + ( J ^ + v r ) + v F (2.19) a t a t The laminar diffusive fluxes J f c and Jh are usua l ly smal l compared w i t h turbulent transport at sufficiently h igh turbulent Reyno lds n u m -bers. The unclosed terms generated from averaging are the Reynolds stresses v " v " , species turbulent fluxes v"Y f e " , temperature fluxes v " T " and mean chemical source-term u>k- R A N S models are used to solve the first three of these unclosed terms. Closure for the mean chemical source-term is discussed i n section 2.4. The most w i d e l y used R A N S m o d e l is the k — e m o d e l , w h i c h is based on a gradient transport as-sumpt ion . In this m o d e l , two transport equations for the turbulent k i -netic energy k and its diss ipat ion rate e are also so lved. A l t h o u g h the k — e m o d e l is one of the most c o m m o n l y used models , it has m a n y uncertainties. It is based on the assumption of homogeneous isotropic turbulence. Furthermore, there are several empi r i ca l constants whose 14 v a l u e s are t u n e d b y m a t c h i n g to s i m p l e f l o w s . A l s o , i t i s n o t y e t c l ea r h o w the exac t b o u n d a r y c o n d i t i o n for the t w o a d d i t i o n a l e q u a t i o n s s h o u l d be d e f i n e d . W i l c o x [13] g i v e s a c o m p r e h e n s i v e r e v i e w a b o u t th is m o d e l . U n s t e a d y R e y n o l d s - a v e r a g e d N a v i e r - S t o k e s ( U R A N S ) o r T r a n s i e n t R e y n o l d s - a v e r a g e d N a v i e r - S t o k e s ( T R A N S ) s i m u l a t i o n s are u s e d for t r ans i en t f l o w s [14]. I n U R A N S , the u s u a l R e y n o l d s d e c o m p o s i t i o n i n E q . 2.15 i s e m p l o y e d b u t the t r ans i en t t e r m s are r e t a i n e d i n E q s . 2.16 t h r o u g h 2.19. T h e d e p e n d e n t v a r i a b l e s are n o t o n l y a f u n c t i o n o f the space , b u t a l s o a f u n c t i o n o f t i m e . I n U R A N S , o n e n e e d s to b e c a u t i o u s a b o u t u s i n g the s t a n d a r d k — e m o d e l w h i c h c a n to b e t o o d i s s i p a -t i v e [14]. 2.3.3 Large Eddy Simulation B e c a u s e R A N S o n l y s o l v e s for the a v e r a g e d sca la r s i n the f l o w f i e l d , i t o n l y p r o v i d e s i n f o r m a t i o n o n m e a n q u a n t i t i e s ; t h u s , i t i s n o t a l w a y s su f f i c ien t to d e s c r i b e t u r b u l e n t c o m b u s t i o n p r o b l e m s . F u r t h e r m o r e , R A N S m o d e l s t e n d to b e h e u r i s t i c a n d e m p i r i c a l i n n a t u r e . A n a l t e r n a t i v e is to u se L a r g e E d d y S i m u l a t i o n ( L E S ) w h i c h i s a p r o m i s i n g t e c h n i q u e be-cause i t i s c o m p u t a t i o n a l l y a f fo rdab le a n d g i v e s a m o r e p h y s i c a l . m o d e l for the R e y n o l d s stress. T h e b a s i c i d e a o f L E S i s to s o l v e d i r e c t l y for the l a rge sca les w h i c h c a r r y the m a j o r i t y o f the t u r b u l e n t k i n e t i c e n e r g y a n d m o d e l the s m a l l sca les w h i c h t e n d to b e m o r e i s o t r o p i c a n d eas ier to m o d e l . T h i s m e a n s the m e s h n e e d s to be f ine e n o u g h to c a t c h the l a rge e d d i e s ; n e v e r t h e l e s s , the m e s h i s s t i l l o f t en coa r se e n o u g h to be c o m p u t a t i o n a l l y a f fo rdab l e . I n L E S , a v a r i a b l e $ is f i l t e r e d i n the spec -t r a l s p a c e o r i n the p h y s i c a l space b y w e i g h t e d a v e r a g i n g i n a g i v e n 15 volume. The spatial filter operator on a g iven variable $ is defined as: $ = / G , (x ,x ' )$ (x , ) r fx ' (2.20) JD where D denotes the integration d o m a i n and G ( x , x ') is the filter ker-nel . Standard filters kernels inc lude the top hat filter, the box filter and the Gauss ian filter. A top hat filter is defined as: f 1, if | x - x ' | < f ; G ( x , x ' ) = ^ ' ~ 2 (2.21) I 0, otherwise. where A is the filter size. A n y variable $ can be decomposed into a filtered component $ and a fluctuating component $ ' : $ = $ + $ ' . (2.22) Here can be non-zero depending on the filter chosen w h i c h is differ-ent from R A N S . A p p l y i n g Favre average fil tering for the instantaneous balance equat ion gives: ff + V • (pv) = 0 (2.23) + V • (pvv) = - V ( p ( v v - vv ) ) — V p + V r + F (2.24) d p Y k + V • (pvYk) = - V ( p ( v y f c - v7 f c )) + ZTk (2.25) dt ~ + V • (pvh) = - V ( p ( v 7 J - wht)) + ^ + V • (J* + v r ) + v F . (2.26) The unclosed terms i n L E S are the subgrid-scale (SGS) Reyno lds stress term (vv — v v ) , unresolved species fluxes (vYk — v Y f c ) , enthalpy fluxes (vht—vht), filtered laminar diffusion fluxes J k and J h , and mean chem-ical source-term cbk. The most c o m m o n l y used SGS m o d e l is the Smagor in-sky m o d e l [15] w h i c h assumes the SGS stress fol lows a gradient-diffusion 16 process . T h e m a i n l i m i t a t i o n o f the S m a g o r i n s k y m o d e l i s tha t the m o -d e l c o n s t a n t m u s t b e t u n e d for d i f fe ren t f l o w s . G e r m a n o et a l . [16] p r o p o s e d a d y n a m i c s u b g r i d - s c a l e e d d y v i s c o s i t y m o d e l w h i c h deter -m i n e s the coef f i c ien t as a f u n c t i o n o f space a n d t i m e . M o i n et a l . [17] fu r the r e x t e n d e d th i s t e c h n i q u e to e d d y d i f f u s i v i t y . A d y n a m i c p r o -c e d u r e for d e t e r m i n i n g the f i l te r w i d t h r a t i o , a m o d e l p a r a m e t e r tha t a p p e a r s i n the d y n a m i c S m a g o r i n s k y m o d e l a n d tha t i s t r a d i t i o n a l l y a s s u m e d to be cons t an t , h a s a l so b e e n d e v e l o p e d [18]. L E S has b e e n u s e d w i t h success for a l a rge v a r i e t y o f p r o b l e m s . H o w e v e r , i t h a s s o m e r e m a i n i n g q u e s t i o n s a n d i t s p o t e n t i a l i s n o t f u l l y d e v e l o p e d . L E S i s a l s o a p p e a l i n g for t u r b u l e n t c o m b u s t i o n . A t h i g h R e y n o l d s a n d D a m k o h l e r n u m b e r s , the r a t e - c o n t r o l l i n g p roces ses o f m o l e c u l a r m i x i n g a n d c h e m i c a l r eac t ions o c c u r at the s m a l l e s t sca les w h i c h are n o t r e s o l v e d b y L E S ; t h u s these n e e d to m o d e l e d j u s t as i n R A N S [19]. H o w e v e r , L E S g i v e s be t ter p r e d i c t i o n s o f the s c a l a r m i x -i n g p r o c e s s a n d d i s s i p a t i o n rates c o m p a r i n g w i t h R A N S . P i t s c h [20] r e p o r t e d tha t the f i l t e r e d sca la r d i s s i p a t i o n rate p r e d i c t e d b y L E S l e a d s to a s u b s t a n t i a l l y i m p r o v e d p r e d i c t i o n , e s p e c i a l l y for m i n o r spec ies . P i t s c h [21] h a s a l s o p r o v i d e d a c o m p r e h e n s i v e r e v i e w o f u s i n g L E S i n t u r b u l e n t c o m b u s t i o n . 2.3.4 Hybrid C F D methods L E S p r o v i d e s a m o r e p h y s i c a l m o d e l t h a n R A N S b u t t he c o m p u t a -t i o n cost i s s t i l l q u i t e h i g h for p r a c t i c a l i n d u s t r y a p p l i c a t i o n s . A s w e l l , L E S i s d i f f i c u l t to a p p l y for n e a r - w a l l r e g i o n s w h i c h r e q u i r e f ine g r i d s , as m e n t i o n e d b y C h a p m a n [22]. S o m e n e w C F D m e t h o d s h a v e b e e n p r o p o s e d b e s i d e s D N S / R A N S / L E S , s u c h as D e t a c h e d E d d y S i m u l a -17 t ion (DES) w h i c h is based on a combinat ion of L E S and R A N S . D E S was first p roposed b y Spalart et al . [23,24]. It was deve loped as a technique for realistic s imulat ions of h igh Reynolds number turbu-lent f lows w i t h massive separation. It ut i l izes the R A N S i n the attached boundary layers and L E S i n the shear layers and separated f low re-gions. D E S predict ions of three-dimensional and t ime-dependent fea-tures i n mass ive ly separated flows are superior to R A N S [23,24]. In general, these D E S models use a transfer function to effect t ransi t ion from the standard R A N S turbulence m o d e l to the L E S sub-gr id type mode l . 2.4 Turbulent combustion modeling In both R A N S and L E S of a non-premixed reacting flow, closure for the mean chemical source-term of the reacting f low is cha l lenging because the reaction rate is a nonlinear function of temperature a n d species mass fractions. Subst i tut ing averaged values direct ly into an Ar rhen ius -type equation that describes reaction rate usua l ly does not lead to a good approx imat ion of the mean of the reaction rate. Several models have been proposed to resolve this, i n c l u d i n g conserved scalar m o d -els, the probabi l i ty densi ty function (PDF) m o d e l (which is discussed i n section 2.6.1) and the linear eddy mode l . Fast chemistry, the laminar flamelet m o d e l and C M C models are a l l conserved scalar models . The l inear eddy m o d e l reduces a l l the scalars onto a loca l ly one-d imens iona l descr ipt ion [25]. The significant feature of this m o d e l is that it expl ic i t ly models turbulent st irr ing, molecular diffusion and chemical reaction at a l l scales i n a one-dimensional equation. This mo-18 del has been app l i ed to L E S b y M c M u r t r y [26], w h o in t roduced an i n -d i v i d u a l l inear e d d y i n each computat ion cell to describe the subgr id m i x i n g , diffusion and reaction progress. The convective process was mode led by a " s p l i c i n g " of events by w h i c h l inear e d d y elements are copied to and f rom neighbor ing g r id cells. The disadvantage of this m o d e l is that it adds another d imens ion w h i c h can make the compu-tation quite expensive, par t icular ly because the g r i d must be ve ry fine i n the extra d imens ion . Conse rved scalar models are usual ly based o n the mixture fraction Z w h i c h describes the state of m i x i n g between fuel and oxidizer . M i x -ture fraction can be defined based on a combinat ion of element mass fractions Z = ^ ~ (2.27) where £ is the mass fraction of element i, and £ l ( o and £itp are the mass fractions of the element i n the oxidizer and fuel streams respectively. W i t h the assumpt ion of un i ty L e w i s number, the mixture fraction is independent of the choice of the element. A s E q . 2.27 shows, Z = 0 i n the oxid izer stream and Z = 1 i n the fuel stream. The fast chemistry m o d e l [27] makes the assumpt ion that the chem-ical reaction rate is inf ini te ly fast; therefore, the species mass fractions and temperature are determined b y the m i x i n g state, that is " m i x e d is burned" . A l l the reacting scalars can then be related b y the mixture fraction. A n inf ini te ly fast single step chemical reaction [28] was first used but it m a y be replaced b y chemical equ i l i b r i um condi t ion to han-dle multi-step chemistry [29]. However , this m o d e l lacks abi l i ty to pre-dict any finite rate chemistry effect. 19 2.4.1 Laminar flamelet model I n fast c h e m i s t r y , the r e a c t i o n t i m e scales are a s s u m e d to be i n f i n i t e l y s m a l l c o m p a r e d w i t h the t i m e sca les o f t r a n s p o r t . T h e r e a c t i o n z o n e i s a n i n f i n i t e l y t h i n in te r face . T h e f l a m e l e t m o d e l r e l axes the fast che-m i s t r y a s s u m p t i o n s o m e w h a t b y i n t r o d u c i n g a sca l a r d i s s i p a t i o n rate as a p a r a m e t e r to m e a s u r e the d e p a r t u r e f r o m the e q u i l i b r i u m state. H o w e v e r , t he f l a m e l e t m o d e l s t i l l h a s the u n d e r l y i n g a s s u m p t i o n tha t the t i m e sca les o f the c h e m i c a l k i n e t i c s are m u c h s h o r t e r t h a n the time scales o f c o n v e c t i o n a n d d i f f u s i o n ; therefore the r e a c t i n g in te r face c a n be v i e w e d as a n e n s e m b l e o f l a m i n a r d i f f u s i o n f l ame le t s . It i s v a l i d w i t h i n the " f l a m e l e t r e g i o n " as s h o w n i n F i g . 2 .1 . Pe t e r s [30] g i v e s a n e x t e n s i v e r e v i e w o f the f l a m e l e t m o d e l . A n e w c o o r d i n a t e Z w h i c h i s l o c a l l y n o r m a l to the sur face o f the s t o i c h i o m e t r i c m i x t u r e i s i n t r o -d u c e d . W e t h e n le t Z 2 = X2, Z3 = x 3, r = t, a n d a p p l y the f o l l o w i n g t r a n s f o r m a t i o n : d dZ d dx\ dxi dZ _d_ _ _d_ 0Z_Jl_ ( k _ 2 3 ) dxk dZk dxk dZ' A s s u m i n g tha t d i f f e r e n t i a l d i f f u s i o n effects c a n b e n e g l e c t e d ( that the L e w i s n u m b e r Let w h i c h r ep resen t s the r a t i o o f t h e r m a l to m a s s d i f -f u s i v i t y i s u n i t y for a l l spec ies ) , the s p a t i a l t r a n s p o r t e q u a t i o n c a n b e t r a n s f o r m e d to the m i x t u r e f r a c t i o n space : " f = f 0 + * - ^ ) < i 2 9 > dT PXd2T 1 r8p , , x 20 w h e r e p i s d e n s i t y , Yl i s the m a s s f r a c t i o n o f spec i e s i, cjj i s c h e m i c a l r e a c t i o n rate , T i s t e m p e r a t u r e , 1^ i s the e n t h a l p y o f spec i e s i, cp i s the spec i f i c hea t c a p a c i t y at cons t an t p r e s s u r e , a n d qR i s r a d i a t i o n ; x de -notes the sca l a r d i s s i p a t i o n ra te w h i c h r ep resen t s d i f f u s i v e t r a n s p o r t i n the f l a m e l e t e q u a t i o n . T h e r a d i a t i o n a n d p r e s s u r e f l u c t u a t i o n t e r m i n the e n e r g y e q u a t i o n are u s u a l l y o m i t t e d . T h e R(Yi) a n d R{T) t e r m s c o n t a i n d e r i v a t i v e s w i t h r e spec t to Z2 a n d Z 3 w h i c h a re o f a l o w e r o r -d e r c o m p a r e d w i t h the f i rs t t e r m o n the R H S o f E q s . 2.29 a n d 2.30 a n d c a n be n e g l e c t e d . T h e f i n a l f l a m e l e t e q u a t i o n s are: p a 7 = T a z ^ ( 2 3 1 ) dT P X d2T >A hi . ^ = Y d z ^ - ^ 7 p u - ( Z 3 2 ) i=i r There fo re , the o r i g i n a l t r a n s p o r t e q u a t i o n s w i t h f o u r d i m e n s i o n a l co -o r d i n a t e s {x\,xi,xz,t) a re n o w r e d u c e d to t w o d i m e n s i o n s o f (Z,T). T h i s p r o c e d u r e g r e a t l y r e d u c e s the c o m p u t a t i o n a l cost , s u c h that s o l v i n g E q s . 2.31 a n d 2.32 w i t h d e t a i l e d c h e m i s t r y i s p o s s i b l e . H o w e v e r , th i s r e d u c t i o n d e p e n d s o n the a s s u m p t i o n that the r e a c t i o n z o n e i s a o n e -d i m e n s i o n a l s t r u c t u r e n o r m a l to the sur face o f the s t o i c h i o m e t r i c m i x -ture. T h e i n s t a n t a n e o u s sca la r d i s s i p a t i o n rate x i s d e f i n e d b y : X = 2DVZ • V Z , (2.33) w h e r e D d e n o t e s m i x t u r e f r a c t i o n d i f f u s i v i t y . T h e sca l a r d i s s i p a t i o n i n c o r p o r a t e s the i n f l u e n c e o f c o n v e c t i o n a n d d i f f u s i o n n o r m a l to the sur face o f the f l a m e . It i s a p a r a m e t e r tha t re la tes the t w o c o o r d i n a t e sy s t e ms . T h e s c a l a r d i s s i p a t i o n i s t y p i c a l l y a s s u m e d to t ake the f o r m X = XoX*(Z)> w h e r e x*(Z) i s s o m e p r e s u m e d f u n c t i o n . P e t e r s [30] o b -21 tains an analyt ical solut ion for x* f rom an analyt ical so lu t ion of a one-d imens iona l l aminar m i x i n g layer: exp{-2[erfc- 1 (2Z)] 2 } X 1 ; exp{-2[er fc - 1 (2Z 0 ) ] 2 } 1 j Some other analyt ical solutions of x have also been proposed [31]. Pitsch obtained the condi t ional scalar diss ipat ion rate f rom the instan-taneous uncondi t iona l scalar diss ipat ion rate f ield b y inver t ing an i n -tegral equation [20]. P i t sch [31] also extended the flamelet equation to inc lude differ-ential diffusion effects w i t h m u c h more compl icated expressions than Eqs. 2.28 a n d 2.29. H e implemented the Lagrang ian Flamelet m o d e l for a steady turbulent diffusion C H 4 / H 2 / N 2 - a i r flame [32] and found that differential diffusion can cause a strong history effect and therefore can not be neglected i n that flame. However , he found that the differential diffusion effect is negl igible i n the Sandia D flame [33]. B y assuming that the flame is local ly steady, the t ime dependent term can be omit ted, w h i c h gives the Steady L a m i n a r Flamelet M o -del ( S L F M ) . A l ib ra ry of steady state flamelets can then be generated for arbitrary complex chemistry and tabulated i n terms of the scalar diss ipat ion rate xo and Z. A p resumed shape of the joint P D F of scalar diss ipat ion and mixture fraction w i t h the assumpt ion of statistical i n -dependence of Z and xo can be used to obtain the mean value of a species mass fraction. Yz= Yl(Z;X)P(Z)P(x)dZdX. (2.35) Jo Jo The S L F M relaxes the infinite fast chemistry m o d e l considerably and has been used w i t h some success. However , the S L F M is insuffi-cient w h e n dea l ing w i t h unsteadiness, extinction and re-ignit ion. N o r 22 c a n i t a c c o u n t for the h i s t o r y effect o f the s ca l a r d i s s i p a t i o n rate , i .e. , that the f l a m e i s i n f l u e n c e d n o t j u s t b y the c u r r e n t s c a l a r d i s s i p a t i o n b u t a l s o b y i ts h i s t o r y [34-36] . I n the U n s t e a d y L a m i n a r F l a m e l e t M o d e l ( U L F M ) , the u n s t e a d y t e r m is r e t a i n e d a n d a L a g r a n g i a n t i m e i s i n t r o d u c e d to a c c o u n t for h i s t o r y effects i n the f l a m e l e t s t r uc tu r e [37]. T h e L a g r a n g i a n t i m e c a n be d e t e r m i n e d b y i n t r o d u c i n g n o t i o n a l p a r t i c l e s at the i n l e t o f f l o w a n d t r a n s p o r t d o w n s t r e a m . P i t s c h [38] u s e d the u n s t e a d y f l a m e l e t m o d e l w i t h L a g r a n g i a n t i m e to cor re la te w i t h the E u l e r i a n t i m e c o o r d i n a t e s y s t e m . A d i f f u s i o n t i m e tx w a s u s e d to d e s c r i b e the t i m e n e e d e d to e x c h a n g e m a s s a n d e n e r g y o v e r the f l a m e t h i c k n e s s i n m i x t u r e f r a c t i o n space . If tx i s s h o r t c o m p a r e d w i t h the L a g r a n g i a n time, the u n s t e a d y t e r m i n the f l a m e l e t e q u a t i o n s c a n b e n e g l e c t e d . R e p r e s e n t a t i v e In ter -a c t i v e F l a m e l e t (RTF) [39] a n d a s i m i l a r a p p r o a c h c a l l e d E u l e r i a n P a r -t ic le F l a m e l e t M o d e l [40] h a v e b e e n d e v e l o p e d b a s e d o n the u n s t e a d y f l amele t . T h e R I F m o d e l h a s b e e n a p p l i e d to i n t e r n a l c o m b u s t i o n e n -g i n e s w i t h s o m e success [39]. 2.4.2 Conditional Moment Closure B o t h fast c h e m i s t r y a n d the l a m i n a r f l a m e l e t m o d e l s u g g e s t tha t c o n d i -t i o n a l a v e r a g e s w i t h i n the m i x t u r e f r a c t i o n s p a c e are u s e f u l for s t u d y -i n g t u r b u l e n t c o m b u s t i o n . C o n d i t i o n a l m o m e n t c l o s u r e ( C M C ) i s a l s o a c o n s e r v e d sca l a r m o d e l b a s e d o n c o n d i t i o n a l a v e r a g i n g tha t w a s first p r o p o s e d b y B i l g e r [41] a n d K l i m e n k o [42] i n d e p e n d e n t l y . A n ex t en -s i v e d i s c u s s i o n o f C M C is g i v e n b y K l i m e n k o a n d B i l g e r i n [43]. T h e c o n d i t i o n a l a v e r a g e p r o v i d e s a be t te r d e s c r i p t i o n o f the a c t u a l c o n d i -t i o n t h a n a g l o b a l a v e r a g e , r e s u l t i n g i n s m a l l e r f l u c t u a t i o n s a r o u n d the 23 F i g u r e 2.3: E x p e r i m e n t m e a s u r e m e n t o f s ca t t e red a n d c o n d i t i o n a l a v -erage t e m p e r a t u r e o f S a n d i a D - f l a m e at x/d = 30 (g rey do t s : sca t te red da ta ; S o l i d b l a c k l i n e : c o n d i t i o n a l a v e r a g e t e m p e r a t u r e ) . a v e r a g e as s h o w n i n F i g . 2.3. B i l g e r [41,43] u s e d a d e c o m p o s i t i o n a p p r o a c h to d e r i v e the c o n d i -t i o n a l a v e r a g e t r a n s p o r t e q u a t i o n . T h e b a s i c d e c o m p o s i t i o n h e u s e d is: F ( x , t) = Q(Z{x, t),x, t) + Y"(x, t) (2.36) w h e r e Q =< Y(x, t ) | Z = £ > is the c o n d i t i o n a l e x p e c t a t i o n o f Y w h e n the c o n d i t i o n i n g v a r i a b l e Z equa l s the v a l u e o f C; Y" i s the f l u c t u a t i o n w i t h respec t to the c o n d i t i o n a l m e a n . T h e d e r i v a t i o n i n v o l v e s d i f f e r e n -t i a t i n g E q . 2.36 as: dY _dQ ^ dQdZ { dY" dt dt d( dt dt { ' } 24 or: V F = V Q + | | VZ + VY" (2.38) T h e m o l e c u l a r d i f f u s i o n t e r m i s : V • (pDVY) = V • {pDVQ) + • {pDVZ) . (2.39) + V - (pDVY"). S u b s t i t u t i n g the a b o v e t r a n s f o r m a t i o n s i n t o the t r a n s p o r t e q u a t i o n a n d t a k i n g the c o n d i t i o n a l a v e r a g e g i v e s the u n c l o s e d f o r m o f the C M C e q u a t i o n : dO d20 < p\( >~^+< P\C >< v|C > - V Q - < p\C >< N\C > ~ ^ (2.40) = < p|C > < u\C > +eQ + ey . w h e r e < iV|£ > = D(yZ)2 deno te s the sca la r d i s s i p a t i o n ra te . T h i s h a s b e e n m o d e l e d i n v a r i o u s w a y [44-47] . T h e u n c l o s e d t e r m eQ d e n o t e s the d i f f u s i v e effects: eQ=< V{pDVQ)+pDVZ-V^\C> . (2.41) It i s s m a l l a n d therefore c a n b e n e g l e c t e d w h e n R e i s l a r g e . A n o t h e r u n c l o s e d t e r m , e y , i s d e f i n e d as: eY = -< p — + p u • V Y " - V • ( D p V Y " ) | C > . (2.42) at T h e p r o p o s e d c l o s u r e for ey i s [43]: v - ( < pic > < u " Y " K > P(Q) v - ( < p 1 C > ( - A V Q ) p ( 0 ) ^ ( 0 < p I C > ~ P(()<P\(> (2.43) 25 The condi t ional averaged value < Y^ > describes the characteristics of the reactive scalar fields m u c h better than the uncondi t iona l values < Yi >. In cases where the condi t ional fluctuations Y" can be neglected, the condi t ional averaged reaction rate can be closed by: C J I C ^ C T I C X I C P I C ) ( 2 - 4 4 ) where £ is the sample space for mixture fraction, ui is the chemical reaction rate, T\( is the condi t ional ly filtered temperature, Yi\( is the condi t ional ly filtered mass fraction vector and p\( is the condi t iona l ly filtered density. The enthalpy C M C equat ion is s imi lar to E q . 2 . 4 3 . K l i m e n k o ' s ap-proach to de r iv ing the C M C transport equation used the P D F transport equation to obtain a very s imi lar final result as Bi lger ' s [ 4 2 ] . W i t h i n a homogeneous field, the C M C equation reduces to: d < Y\C > d2 < Y\C > at' - < W > dC2 = < * K > ( 2 - 4 5 ) w h i c h is s imi lar to the unsteady flamelet equation; however their un -de r ly ing phys ica l assumptions differ significantly. C M C o n l y makes the assumpt ion that the var ia t ion about the condi t iona l averages is smal l . W h i l e the C M C m o d e l us ing first moment closure has p roven to be successful i n m a n y cases [ 4 4 , 4 6 , 4 8 - 5 1 ] , it is insufficient w h e n the fluctuation about the condi t ional mean become significant, such as i n flames w i t h extinct ion and reignit ion [ 5 , 4 3 ] . To account for the fluc-tuations w h i c h are not correlated w i t h fluctuations i n Z, different ap-proaches have been proposed; second-order closure of the condi t iona l chemical source-term is one opt ion [ 5 2 , 5 3 ] . K i m [ 4 7 ] app l i ed second-order condi t iona l moment closure mode l i ng to turbulent non-premixed 2 6 h y d r o c a r b o n flames w i t h s i g n i f i c a n t l o c a l e x t i n c t i o n a n d r e i g n i t i o n . H i s r e su l t s s h o w tha t the s e c o n d - o r d e r C M C p r e d i c t i o n h a s be t te r agree-m e n t t h a n the f i r s t -o rde r C M C for the i n t e r m e d i a t e s a n d r e a c t i o n o f fue l . H o w e v e r , i t r e q u i r e s m o d e l i n g s e v e r a l u n c l o s e d t e r m a n d the v a r i -a t i o n b a s e d o n a s i n g l y c o n d i t i o n a l a v e r a g e s t i l l m i g h t b e l a rge . A sec-o n d c o n d i t i o n i n g v a r i a b l e c a n b e i n t r o d u c e d for s u c h cases. C h a [54] uses the s c a l a r d i s s i p a t i o n ra te as the s e c o n d c o n d i t i o n i n g v a r i a b l e ; the resu l t s s h o w tha t i t c a n d e s c r i b e the e x t i n c t i o n s e e n o n a v e r a g e b u t p r e d i c t s the o n s e t o f r e i g n i t i o n too ear ly . 2.4.3 Conditional Source-term Estimation B u s h e a n d S t e i n e r [55] p r o p o s e d C S E to o b t a i n the c o n d i t i o n a l ave r -a g e d v a l u e s w i t h o u t s o l v i n g t r a n s p o r t e q u a t i o n s . T h e y b e g i n b y n o t i n g that, i n the i d e n t i t y 7 ( x , t ) « /171CP(C;x,t)dC, (2.46) Jo w h e r e / is t a k e n to b e d e n s i t y , t e m p e r a t u r e o r a s p e c i e s m a s s f r a c t i o n , o n e i s , i n s o m e w a y , s o l v i n g a t r a n s p o r t e q u a t i o n for the left h a n d s i d e , s u c h tha t the q u a n t i t y i s k n o w n . A l s o , o n e i s l i k e l y to b e s o l v i n g for the m e a n a n d v a r i a n c e o f m i x t u r e f r ac t ion , f r o m w h i c h i t i s p o s s i b l e to r e c o n s t r u c t a n a p p r o x i m a t i o n to P(C; x , t), t he f i l t e r e d d e n s i t y f u n c -t i o n ( F D F ) , u s i n g a /3-funct ion. I n m a n y a p p l i c a t i o n s o f C M C m e t h o d s , / | C i s a s s u m e d to h a v e a s m a l l s p a t i a l g r a d i e n t as t h i s h a s b e e n s h o w n to b e t rue i n e x p e r i m e n t s a n d i n D N S . I n m a n y f l o w s , there i s a s p a t i a l h o m o g e n e i t y i n f\(, s u c h as i s p r e s e n t i n a n a x i s y m m e t r i c jet. A s s u m -i n g the c o n d i t i o n a l a v e r a g e d v a l u e s for s o m e e n s e m b l e A o f M p o i n t s ( w h i c h w o u l d b e c o n t r o l v o l u m e ave rages i n a f in i te v o l u m e a p p r o a c h ) 27 is homogeneous, then E q . 2.46 can be wri t ten for each poin t m w i t h i n the ensemble: where <> denotes the ensemble average w i t h i n the ensemble A. This is an integral equat ion - a F r edho lm equation of the first k i n d - and , w h e n the integral is replaced b y a numer ica l quadrature of TV points, one obtains a l inear system of M equations for TV variables, w h i c h one can then solve. In a free jet, a reasonable choice for the ensemble is a plane at con-stant distance x downs t ream of the nozzle since the variances of condi -tional temperature and species are found to be smal l on that plane. The mean mixture fraction Z(x, t) is obtained b y so lv ing its transport equa-tion; the SGS favre variance of mixture fraction Z " 2 ( x , t) can be m o d -eled i n L E S us ing dynamic procedure [17]. A s s u m i n g the /5-function approximates the F D F of mixture fraction sufficiently w e l l , < f\( > can be obtained b y inver t ing E q . 2.47 w i t h l inear regular iza t ion and substituted into the C M C hypothesis E q . 2.44 to obtain the condi t ional reaction rate. The mean chemical source term Co is then obtained b y integrating E q . 2.47 for each point i n the ensemble and the transport equations for a l l reactive scalars can thus be closed. W h i l e this approach has s h o w n promise i n prev ious L E S work , it has some significant drawbacks. Firstly, so lv ing an inverse p rob lem such as that posed b y the invers ion of E q . 2.47 requires some k i n d of regular izat ion i n order to obtain a stable solut ion for the condi t iona l average. One approach to this regular izat ion is to s imul taneous ly re-duce the error i n the approximat ion of E q . 2.47 a n d a der ivat ive of (2.47) 28 the c o n d i t i o n a l a v e r a g e w i t h respec t to m i x t u r e f r a c t i o n [56]. T h i s a p -p r o a c h g e n e r a l l y l e a d s to a s o l u t i o n that i s too s m o o t h , s u c h that i t has the a p p e a r a n c e o f b e i n g o v e r l y d i f f u s i v e . T h i s i s p a r t i c u l a r l y t rue w h e n a p p l i e d to i n v e r s i o n for the m a s s f r a c t i o n o f a m i n o r spec i e s . A n o t h e r a p p r o a c h i s to a s s u m e tha t the c o n d i t i o n a l a v e r a g e c h a n g e s s l o w l y i n t i m e , s u c h tha t o n e s i m u l t a n e o u s l y m i n i m i z e s the e r r o r i n the a p p r o x -i m a t i o n o f E q . 2.47 a n d the c h a n g e i n the c o n d i t i o n a l a v e r a g e o v e r the las t t i m e s tep . A n o t h e r d r a w b a c k i s tha t o n e h a s to s o l v e a t r a n s p o r t e q u a t i o n for e v e r y spec i e s i n the c h e m i c a l k i n e t i c m e c h a n i s m . T h u s , w h e n i n v e r t i n g E q . 2.47 d i r e c t l y for the c o n d i t i o n a l l y f i l t e r e d m a s s frac-t ions , o n e i s l a r g e l y l i m i t e d to r e d u c e d c h e m i c a l k i n e t i c m e c h a n i s m s . L a m i n a r F l a m e l e t D e c o m p o s i t i o n is a d i f fe ren t a p p r o a c h to s o l v -i n g E q . 2.47. R a t h e r t h a n s o l v i n g d i r e c t l y for the c o n d i t i o n a l ' a v e r a g e , i n L F D , o n e d e c o m p o s e s the c o n d i t i o n a l a v e r a g e i n t o a l i n e a r c o m b i -n a t i o n o f s o l u t i o n s to the l a m i n a r f l a m e l e t e q u a t i o n [30 ,57] . T h u s , the i n v e r s i o n o f E q . 2.47 b e c o m e s a ma t t e r o f i n v e r t i n g for a b a s i s f u n c t i o n coef f ic ien t m a t r i x at / r o ( x , t)= Yl a«e*(0 P - ( C ; x> *) dC- (2-48) Jo t=l w h e r e 0 , i s o n e s o l u t i o n o f the l a m i n a r f l a m e l e t e q u a t i o n w i t h i n a l i -b r a r y o f Nf s o l u t i o n s . H e r e a g a i n , / c a n be a spec i e s m a s s f r a c t i o n , t e m p e r a t u r e o r dens i t y . T h e c o n d i t i o n a l a v e r a g e o f a s ca l a r i s t h e n f\Z = C = 5^^0,(0- (2-49) i=r T h i s m e t h o d h a s s o m e i m p o r t a n t a d v a n t a g e s . F i r s t l y , b e c a u s e a l l spec i e s c a n b e s a v e d i n the l i b r a r y i n a d v a n c e , L F D c a n s o l v e fo r m i n o r 2 9 species wi thou t so lv ing transport equations for them, resul t ing i n s ig-nificant savings i n computat ional effort. Secondly, it is possible to store the reaction rates associated w i t h the flamelet solutions i n the l ibrary and look up the rates i n the l ibrary to obtain the condi t iona l average of the reaction rates us ing E q . 2.49. This represents another significant savings i n computa t ional effort. For details of h o w to invert the inte-gral equation, please see appendix A . The flamelet decomposi t ion method was used i n an a priori test and p roved to be able to p rov ide acceptable closure for R A N S b y compar-i n g the results w i t h D N S data [58]. Grout et al . [59] have found that this flamelet decompos i t ion can prov ide good predict ions of the autoigni-t ion of turbulent methane /a i r jets i n R A N S . I n that w o r k , they used a more general T ikhonov Regular iza t ion w i t h the assumpt ion that the condi t ional average of the current timestep shou ld not be significantly different f rom that of the previous one (or the previous i teration i n a steady-state calculation). This form of regular izat ion has a better phys -ical interpretation than assumptions about smal l gradients i n mixture fraction. U s i n g flamelet solutions as a l ibrary leaves some remain ing ques-tions though. The first question is the va l id i ty of the assumptions made i n L F D - specifically, an assumption about the flame be ing i n the flame-let regime m a y be i m p l i e d . The second question pertains to h o w the flamelet l ib ra ry shou ld be generated. Ideally one w o u l d l ike to include al l of the possible flamelets i n this l ibrary w h i c h is clearly imposs i -ble. There has been no rigorous attempt made to establish h o w m a n y flamelets are needed i n order to get acceptable accuracy. A l s o , i n a 3 0 L E S c a l c u l a t i o n , o n e w o u l d s u r e l y f a v o r a n u n s t e a d y f o r m u l a t i o n to a c c o u n t for the effects o f u n s t e a d i n e s s o n the f l a m e . T h e n , the s ca l a r d i s s i p a t i o n as a f u n c t i o n o f time s h o u l d m a t c h the f l o w f i e l d u n d e r c o n s i d e r a t i o n . H u a n g a n d B u s h e [60] u s e d C S E w i t h T G L D M i n c o r p o r a t i n g the m o r e g e n e r a l T i k h o n o v r e g u l a r i z a t i o n o f G r o u t et a l . [61] i n p r e d i c t i n g the a u t o i g n i t i o n d e l a y time i n t u r b u l e n t jets w i t h a R A N S m e t h o d w i t h s o m e success . T h e m e a n i g n i t i o n d e l a y a n d N O e m i s s i o n s w e r e r e a s o n -a b l y w e l l - p r e d i c t e d , b u t , b ecause a R A N S m e t h o d w a s u s e d , f l u c t u a -t ions i n the i g n i t i o n d e l a y s c o u l d n o t be p r e d i c t e d . 2.5 Simplification of reaction mechanisms C h e m i c a l r e a c t i o n m e c h a n i s m s t y p i c a l l y i n v o l v e a l a r g e n u m b e r o f spec ies a n d r eac t ions . F o r e x a m p l e , G R I M e c h 3.0 [62], w h i c h d e s c r i b e s the c h e m i s t r y o f n a t u r a l gas r e a c t i n g w i t h air , i n c l u d e s 53 spec i e s a n d 325 e l e m e n t a r y r e a c t i o n s for m e t h a n e / a i r c o m b u s t i o n . T h e a u t o i g n i -t i o n o f D i e s e l f u e l , h a s t h o u s a n d s o f e l e m e n t a r y r eac t i ons . B e c a u s e o f the l a rge v a r i e t y o f time scales i n the c h e m i c a l r e a c t i o n m e c h a n i s m , the s y s t e m is o f ten st i ff a n d c o m p u t a t i o n a l l y e x p e n s i v e to s o l v e . A m o n g a l l the e l e m e n t a r y r eac t ions , o n l y a f e w d e t e r m i n e the ra te o f the o v e r -a l l p roces s ; these are r e f e r r e d as r a t e - l i m i t i n g r eac t i ons . B y f o c u s i n g o n these r eac t i ons , a r e d u c e d c h e m i c a l m e c h a n i s m c a n b e d e v e l o p e d to save c o m p u t a t i o n a l cost . 2 . 5 . 1 Q S S A a n d p a r t i a l e q u i l i b r i u m T h e t w o m o s t c o m m o n m e a n s to s i m p l i f y c h e m i c a l r e a c t i o n m e c h a -n i s m s are the Q u a s i - s t e a d y State A p p r o x i m a t i o n ( Q S S A ) a n d the p a r -31 t ial e q u i l i b r i u m method. The Q S S A was first a p p l i e d to chemical k i -netic schemes b y Bodenstein [63]. B y assuming the rates of consump-tion and format ion are the same for some species, one can obtain an approximate analyt ical solut ion for the kinetic differential equations. For example, for a s imple case w i t h two elementary, irreversible steps: r< &12 c 2^3 0 Oi >02 ^03 the rate l a w gives: d[S; - = -kniS!] (2.50) dt ^ = * i 2 [S i ] - k23[S2] . (2.51) d-^ = k23[S2}. (2.52) A s s u m i n g S2 is a h igh ly reactive species and the rate of consumpt ion is approximate ly equal to the rate of formation, then ^ ~ = k12[S1]-k23[S2}^0. (2.53) A n analyt ical so lu t ion for Si and 5*3 can then be deduced f rom E q . 2.50 to E q . 2.53. [Si] = [S1]0-exp(-k12-t). (2.54) [S3} = {S1]0-[l-exp(-k12-t)}. (2.55) C o m p a r i n g E q . 2.54 and E q . 2.55 w i t h the exact analyt ical solut ion shows that, where the assumpt ion is v a l i d , Q S S A can be a good ap-proximat ion , and that the solut ion f rom Q S S A approaches the exact solut ion i n a short t ime. App l i ca t ions of Q S S A have been used to convert the stiff differen-tial equations into non-stiff systems. Peters and Kee [64] propose a sys-32 t e m a t i c w a y to u s e Q S S A for c o n s t r u c t i n g r e d u c e d c h e m i c a l m e c h a -n i s m s . B y l o o k i n g at the m a g n i t u d e o f i n t e r m e d i a t e spec i e s , t h e y i d e n -t i fy c e r t a i n spec i e s w h i c h sa t i s fy the Q S S A . I n p a r t i a l e q u i l i b r i u m , the a s s u m p t i o n i s n o t m a d e for spec ies , b u t for e l e m e n t a r y r eac t ions . Pe te r s et a l . [64] i n t r o d u c e t h i s m e t h o d to i d e n t i f y fast p roces ses . T h e n e t r e a c t i o n rate o f the fast r e ac t i ons is s m a l l w h e r e a s the f o r w a r d a n d b a c k w a r d r e a c t i o n ra tes are l a rge . F o r e x a m p l e , i f a p a r t i a l e q u i l i b r i u m a s s u m p t i o n i s m a d e fo r r e a c t i o n : A + B ^ C + D t h e n kl[A][B] = k2[C}[D}. (2.56) B y u s i n g the p a r t i a l e q u i l i b r i u m a s s u m p t i o n , the c o n c e n t r a t i o n o f i n -t e r m e d i a t e s c a n b e e x p r e s s e d i n t e r m s o f the c o n c e n t r a t i o n o f the s tab le spec ies . T h e Q S S A a n d the p a r t i a l e q u i l i b r i u m a s s u m p t i o n s are u s e d i n m a n y a p p l i c a t i o n s [65,66] . Pe te r s et a l . [64] u s e d p a r t i a l e q u i l i b r i u m for s o m e r e a c t i o n s to d e v e l o p a four - s tep r e d u c e d m e c h a n i s m . T h i s r e d u c e d m e c h a n i s m i s ab l e to p r e d i c t the i m p o r t a n t p r o p e r t i e s o f m e t -h a n e / a i r d i f f u s i o n f l a m e s , s u c h as the o v e r a l l s t r u c t u r e a n d e x t i n c t i o n l i m i t s . W a r n a t z [67] a p p l i e d th i s m e t h o d to a p r e m i x e d s t o i c h i o m e t r i c h y d r o g e n / a i r f l a m e at 1 b a r p r e s s u r e , u n b u r n e d m i x t u r e t e m p e r a t u r e 2 9 8 K w i t h the i n t e r m e d i a t e s c a l c u l a t e d u s i n g the p a r t i a l e q u i l i b r i u m a s s u m p t i o n a n d t h e n c o m p a r e d w i t h d e t a i l e d c h e m i s t r y . T h e resu l t s s h o w that the p a r t i a l e q u i l i b r i u m a s s u m p t i o n g i v e s a g o o d p r e d i c t i o n at t e m p e r a t u r e s h i g h e r t h a n 1 6 0 0 K b u t n o t for l o w e r t e m p e r a t u r e s . A n -o the r l i m i t a t i o n o f the m e t h o d i s tha t the r eac t ions w h i c h are i m p o r t a n t c h a n g e i n t i m e a n d s p a c e d u r i n g the c o m b u s t i o n p r o c e s s ; there fore , i t 3 3 is necessary to have specific reduced schemes for specific applications. G o o d knowledge of the chemical kinetics and the flame are essential for a p p l y i n g either the Q S S A or the part ial e q u i l i b r i u m assumptions appropriately. Several analysis methods are available to identify fast reactions. Sensit ivity analysis identifies the rate l i m i t i n g reaction steps b y analy-sis of the dependence of the species concentration solutions or veloci ty on the system parameters, such as the rate coefficients of the elemen-tary reactions [68,69]. If the reactions are found to have h i g h sensi t ivi-ties, they have a h i g h impact on the rate of the system and need to be retained i n the reduced chemistry. Reactions that have l o w sensitivities are less important and can be e l iminated f rom the detai led mechanism. In this way, a s impl i f i ed chemical mechanism is generated. The d raw-back of sensit ivi ty analysis is that it cannot determine the m a i n reaction path the system takes to evolve. Integral reaction f low analysis and local reaction f low analysis can be done for f lows d u r i n g the numer ica l s imula t ion of the combust ion process. The relative contributions of different reactions to the forma-tion and consumpt ion of the chemical species are then analysed [8]. This method can be used to determine the characteristic reaction path. Eigenvalue analysis considers the Jacobian matr ix of the reaction system to d is t inguish the t ime scales of reactions. M e t h o d s depend-i n g on eigenvalue analysis include Compute r S ingular Perturbat ion (CSP) [70,71], the Intrinsic L o w - D i m e n s i o n a l M a n i f o l d s ( I L D M ) , the Trajectory Generated L o w - D i m e n s i o n a l M a n i f o l d ( T G L D M ) and the Invariant Cons t ra ined-equi l ib r ium Edge Pre-Image C u r v e (ICE-PIC) 3 4 method [72,73]. 2.5.2 Computational singular perturbation The Computa t iona l Singular Perturbation (CSP) method was proposed b y L a m and Gouss is [70,71]. U n l i k e Q S S A and par t ia l e q u i l i b r i u m , this method does not depend on in tu i t ion and experience. C S P resolves the stiffness p rob lem b y cont inuously moni to r ing the contr ibut ion to the change of the system state by each of its components (species mass fractions and temperature). A t any g iven time, the system moni tored is G = % = SrFr, (2.57) at where y = (yi,y2, ••.,yN-i,T)T is the array of mass fractions of N - l species a n d temperature; S r denotes the stoichiometric vector; Fr de-notes the reaction rate of the r t h reaction, r = 1...M (where M is the total number of reactions). The system can be decomposed into fast and s low subspaces: . G = Gfast + Gsiow = a;! - 1 + Gsiou, (2.58) where a t is an iV x m matr ix whose co lumns are the basis vectors that span the fast subspace, m denotes the fast subspace d imens ion , and f is an m x m matr ix representing the magni tude of the projected G vector on the fast subspace. The C S P me thod starts b y integrating E q . 2.57 u s ing sufficiently smal l t ime steps. Then the Jacobian matrix is computed and decom-posed to f ind the large negative eigenvalues w h i c h can be decoupled from the system, thus de termining the d imens ion of the fast subspace. A differential equation for f is solved us ing sui tably sma l l t ime steps. 35 T h i s s tep i s i t e r a t e d s e v e r a l t i m e s u n t i l G / a s t i s i n s i g n i f i c a n t . T h e reac-t ions r e l a t e d to the fast m o d e are b y p a s s e d a n d E q . 2 .57 i s i n t e g r a t e d u s i n g l a r g e t i m e s teps . A t e a c h t i m e s t e p , the J a c o b i a n m a t r i x i s c o m -p u t e d . W h e n e v e r the d i m e n s i o n o f fast s u b s p a c e c h a n g e s , i t i s nec -es sa ry to r e t u r n to the i n i t i a l p o i n t a n d r e c a l c u l a t e a l l the p a r a m e t e r s a g a i n . C S P h a s b e e n u s e d to c o n s t r u c t the g l o b a l r e d u c e d m e c h a n i s m for a l a m i n a r p r e m i x e d C H 4 / a i r f l a m e f r o m a c o m p l e x d e t a i l e d k i n e t i c m e c h a n i s m c o n s i s t i n g o f 279 reac t ions a n d 49 spec i e s [74]. T h e s e v e n -s tep r e d u c e d m e c h a n i s m is ab l e to r e p r o d u c e the spec i e s p r o f i l e s a n d l a m i n a r b u r n i n g v e l o c i t y a c c u r a t e l y o v e r a w i d e r a n g e o f m i x t u r e frac-t ions a n d t e m p e r a t u r e s w h i l e a l so s i g n i f i c a n t l y r e d u c i n g the c o m p u t a -t i o n a l cost . L u et a l . e m p l o y e d C S P to genera te a 4-s tep a n d a 10-step r e d u c e d m e c h a n i s m for h i g h - t e m p e r a t u r e H 2 / a i r a n d C H 4 / a i r o x i d a -t i o n [75]. W h e n u s e d to c o m p u t e p e r f e c t l y s t i r r e d reac tors a n d the o n e -d i m e n s i o n a l p l a n a r p r o p a g a t i n g p r e m i x e d f l a m e s the m e t h o d s h o w e d g o o d a g r e e m e n t w i t h the f l a m e s p e e d , f l a m e t e m p e r a t u r e a n d f l a m e s t ruc tu re . C S P c a n a l s o b e u s e d for the a n a l y s i s o f the c o m p l e x b e h a v i o r a s s o c i a t e d w i t h t w o - s t a g e i g n i t i o n o f l a rge h y d r o c a r b o n m o l e c u l e s [76]. C S P i s u s e f u l for i n v e s t i g a t i n g a n d r e d u c i n g m e c h a n i s m s . H o w e v e r , i t d o e s n o t p r o v i d e a c o m p u t a t i o n a l l y ef f ic ient t e c h n i q u e for c o m b u s -t i o n m o d e l i n g [77]. 2.5.3 Intrinsic Low-Dimensional Manifold I n a r e a c t i n g f l o w , the c h e m i c a l r e a c t i o n p r o c e s s c o v e r s a l a r g e v a r i e t y o f t i m e sca les , t y p i c a l l y r a n g i n g f r o m 1 0 ~ 1 0 s (e.g. r a d i c a l r eac t ions ) to m o r e t h a n 1 s ( N O f o r m a t i o n ) . T h i s l a rge r a n g e o f t i m e sca les i s i n d i c a -36 Chemical time scales H Physical time scales F i g u r e 2.4: T i m e scales i n r e a c t i n g flow t i ve o f the h i g h st iffness o f the s y s t e m w h i c h r e q u i r e s a st i ff s y s t e m s o l v e r a n d s m a l l t i m e steps. H o w e v e r , the p h y s i c a l t i m e sca les o f tu r -b u l e n t c o m b u s t i o n flows u s u a l l y c o v e r a s m a l l e r r a n g e f r o m 1 0 ~ 4 to 1 0 - 2 as s h o w n i n F i g . 2.4. Spec i e s i n v o l v e d i n r eac t i ons faster t h a n the p h y s i c a l p roces ses c a n be e l i m i n a t e d . M a a s a n d P o p e [78] f irst p r o p o s e d the I n t r i n s i c L o w - D i m e n s i o n a l M a n i f o l d ( I L D M ) b a s e d o n the a n a l y s i s o f l o c a l t i m e sca les o f t he react-i n g s y s t e m . I L D M n e e d s as i n p u t o n l y the f u l l k i n e t i c m e c h a n i s m a n d the d e s i r e d n u m b e r o f degrees o f f r e e d o m of the r e d u c e d s y s t e m . F o l -l o w i n g M a a s a n d P o p e [78], e x p a n d i n g ns spec i e s w i t h e n t h a l p y a n d p re s su re , a n n = ns + 2 s y s t e m is g i v e n b y : ^ = F ( ^ ) + S ( ^ , V ^ , A ^ ) (2.59) w h e r e 'tp = (h,p, ip\, ip2, i>na)T', F ( ^ ) is the c h e m i c a l s o u r c e - t e r m a n d H deno te s p h y s i c a l p rocesses s u c h as c o n v e c t i o n , d i f f u s i o n , etc. T h e i r a n a l y s i s focuses o n a s i m p l e h o m o g e n e o u s , a d i a b a t i c , i s o b a r i c s y s t e m | = F ( ^ (2.60) T h e e i g e n v a l u e o f the J a c o b i a n (fj = i d e n t i f i e s d i f fe ren t t i m e scales o f the c h e m i c a l r e a c t i o n s y s t e m . T h e y use S c h u r d e c o m p o s i t i o n 37 for the Jacobian matrix: J = V A V (2.61) where V is the n x n r ight Schur vector matr ix of J , A is an upper-triangular matr ix whose diagonal entries are the eigenvalues of J , and V is the inverse matr ix of V ; V and V are sorted i n descending order according to the corresponding eigenvalues A . The inverse of the eigenvalues are the characteristic t ime scales: r = \ - \ (2.62) W h e n A has a large magni tude, it represents a relat ively fast process and can be decoupled f rom the detailed system. A reaction space nr = ns—ne and fast time-scale numbers nj = nr—nc, where nc is the number of s low time-scales, are also defined. The l o w d imens iona l man i fo ld is defined by: V ; %T{w) = 0. (2.63) where V / is nf x n matr ix taken from the last cor responding rows of V . This equation can be so lved w i t h the system constraints a n d progress variables. The final I L D M equation is wri t ten as: / v>(?/o\ G ( V , r ) = = 0 , (2.64) where V / is the matr ix of Schur vectors corresponding to the large neg-ative eigenvalues and P(if), r ) are 2+ne+nc addi t ional parametric equa-tions, w h i c h inc lude the adiabaticity and constant pressure condit ions, n e element conservation equations and n c extra parameter equations. In compar ing the I L D M results of a C O / H 2 / a i r system w i t h conven-tional methods, M a a s and Pope [79] report that the I L D M provides 38 accura t e r e s u l t s n o t o n l y i n areas c l o s e to e q u i l i b r i u m b u t a l s o at l o w t e m p e r a t u r e s w h e r e c o n v e n t i o n a l m e t h o d s f a i l . M a a s a n d P o p e [ 7 9 ] w e n t o n to p r o p o s e a g e n e r a l p r o c e d u r e to c o u -p l e the r e d u c e d c h e m i c a l m e c h a n i s m w i t h p h y s i c a l p r o c e s s e s s u c h as f l o w a n d m o l e c u l a r t r anspo r t . T h e y f o r m u l a t e a p r o j e c t i o n o p e r a t o r Pr w h i c h pro jec ts a p h y s i c a l p e r t u r b a t i o n E(ip) o n t o a p e r t u r b a t i o n T(6) w i t h i n the m a n i f o l d PrE = T . T h e n the m a n i f o l d i s set b y : ^ = S ( 0 ) + r r > ( 0 ) , vm, WW) ( 2 - 6 5 ) w h e r e 9 = (0i,62, . . . ,0jv) i s the p a r a m e t r i z a t i o n o f the m a n i f o l d , t y p i -c a l l y w i t h N < ns. T h e y d e m o n s t r a t e tha t t h i s p r o c e d u r e c a n b e suc -c e s s f u l l y u s e d i n a p e r f e c t l y s t i r r e d f l o w reac to r c a l c u l a t i o n . A s the s e p a r a t i o n b e t w e e n c h e m i s t r y a n d d i f f u s i o n p r o c e s s e s be -c o m e s less s i g n i f i c a n t , the d i m e n s i o n o f the m a n i f o l d h a s to b e i n -c r e a s e d to get a c c e p t a b l e accu racy , for e x a m p l e , at l o w t e m p e r a t u r e s . S i m p l e l i n e a r p r o l o n g a t i o n i n t h i s area i s s t r a i g h t f o r w a r d a n d e a s y to i m p l e m e n t b u t s o m e t i m e s i t c a n n o t p r o v i d e a n a c c u r a t e s o l u t i o n . G i c -q u e l [ 8 0 ] p r o p o s e s c o n s t r u c t i n g the m a n i f o l d b y u s i n g the s o l u t i o n o f l a m i n a r p r e m i x e d free f l a m e s for l o w t e m p e r a t u r e s to e x t e n d the m a n -i f o l d s g e n e r a t e d w i t h I L D M . T h i s p r o c e d u r e i s c a l l e d F l a m e P r o l o n g a -t i o n o f I L D M (FPI ) . G i c q u e l i m p l e m e n t e d F P I s u c c e s s f u l l y i n p r e d i c t -i n g the r e s p o n s e o f a p r e m i x e d f l a m e to s t r a i n i n g u s i n g the d o u b l e -p r e m i x e d c o u n t e r f l o w f l a m e c o n f i g u r a t i o n . A s i m i l a r i d e a b y V a n O i -j e n a n d d e G o e y [ 8 1 ] i s c a l l e d the f l a m e l e t - g e n e r a t e d m a n i f o l d m e t -h o d ( F G M ) . T h e s e m e t h o d s d o n o t u s e m a t h e m a t i c a l a n a l y s i s to p r o v e h o w m a n y p r o g r e s s v a r i a b l e s are suff ic ient . P h a s e - s p a c e I L D M (PS-I L D M ) [ 8 2 ] c o m b i n e s the i d e a o f I L D M a n d F G M w h i c h i n c l u d e s the 3 9 d i f f e r e n t i a l d i f f u s i o n effects a n d h a s a s o u n d m a t h e m a t i c a l a n a l y s i s . I L D M h a s the a d v a n t a g e tha t i t is l o c a l l y d e t e r m i n e d , w h i c h m e a n s that g i v e n a p a r t i c u l a r c o m p o s i t i o n , the c o r r e s p o n d i n g p o i n t o n the m a n i f o l d c a n be d e t e r m i n e d w i t h o u t c o n s t r u c t i n g the w h o l e m a n i f o l d . L o c a l l y d e t e r m i n e d m e t h o d s c a n b e e a s i l y a p p l i e d to h i g h e r d i m e n -s i o n s a n d in situ a d a p t i v e t a b u l a t i o n [83]. H o w e v e r , E q . 2.63 i s d i f f i c u l t to s o l v e a n d I L D M c a n n o t gua ran t ee the ex i s t ence a n d u n i q u e n e s s o f the m a n i f o l d . A n o t h e r d i s a d v a n t a g e o f I L D M is tha t i t i s n o t a n i n e r -t i a l m a n i f o l d u n l e s s the time sca le s e p a r a t i o n a p p r o a c h e s i n f i n i t y . T h i s m e a n s tha t the r e a c t i o n v e c t o r s a l o n g the I L D M are m o s t l y n o t i n the t angen t p l a n e o f the m a n i f o l d . A m a n i f o l d is i n e r t i a l i f the r e a c t i o n t ra -j e c t o r y f r o m a n y p o i n t o n the m a n i f o l d r e m a i n s o n the m a n i f o l d . F o r a n i n e r t i a l m a n i f o l d , the ra te o f c h a n g e o f the r e d u c e d v a r i a b l e s i s d e -t e r m i n e d b y the f u l l c h e m i c a l k i n e t i c s w i t h o u t a p p r o x i m a t i o n . F o r a n o n - i n e r t i a l m a n i f o l d , the ra te o f c h a n g e o f the r e d u c e d v a r i a b l e s de -p e n d s s t r o n g l y o n the c h o i c e o f p a r a m e t e r s u s e d to p ro jec t the f u l l ra te of c h a n g e v e c t o r o n the m a n i f o l d . S i n g h et a l . [84] s h o w e d tha t I L D M i s o n l y a n a p p r o x i m a t i o n o f the m o r e f u n d a m e n t a l s l o w i n v a r i a n t m a -n i f o l d ( S I M ) . D a v i s a n d S k o d j e p r o p o s e d u s i n g a m o d i f i e d F r a s e r a l g o -r i t h m to i m p r o v e the I L D M [85]. 2.5.4 Trajectory Generated Low-Dimensional Manifold P o p e a n d M a a s [86] a l s o d e v e l o p e d the Tra jec to ry G e n e r a t e d L o w -D i m e n s i o n a l M a n i f o l d ( T G L D M ) m e t h o d to genera te a l o w - d i m e n s i o n -a l m a n i f o l d . I n T G L D M , the m a n i f o l d i s d e f i n e d u s i n g t ra jector ies , w h i c h are the p a t h t h r o u g h c o m p o s i t i o n s p a c e t a k e n b y the s y s t e m f r o m the b o u n d a r y o f the r e a l i z a b l e c o m p o s i t i o n s p a c e to c h e m i c a l 4 0 equ i l ib r ium. Each trajectory is evo lved by so lv ing E q . 2.60 u s ing a stiff O D E solver; therefore a converged solut ion is guaranteed. The manif-o l d is then tabulated i n a lower d imensional parameter space. If an at-tractive man i fo ld does exist, any trajectory whose in i t i a l state is on the mani fo ld stays on the mani fo ld . Therefore the T G L D M method gener-ates an iner t ia l man i fo ld equivalent to the S I M . Pope and M a a s [86] define the in i t i a l points of the trajectory us ing the extreme-value-of-major-species method. The extreme bounda ry v a l -ues are so lved b y the constraint of element conservation. For a two d i -mensional man i fo ld , ne + 2 species are chosen according to the global reaction of the system. Then a linear equation system that satisfies the element mass conservation is constructed. ne+2 Yc{j, i)Y(i) = Ye(j) j = 1 , n e , (2.66) i=i where n e denotes the chosen number of elements, Yc denotes the mass fraction of element j i n species i, Y(i) denotes the mass fraction of i, and Ye(j) is the mass fraction of element j. It is n o w possible to solve the extreme value of the mass fractions of bounda ry u s ing linear p r o g r a m m i n g techniques [87]. H u a n g [60] proposed a different means w h i c h is to use constrained equ i l i b r ium to obtain the in i t i a l points. F ig . 2.5 show the compar i son of I L D M w i t h T G L D M [1]. This one d imens iona l man i fo ld demonstrates good agreement between the I L D M and the T G L D M . H u a n g [1] also compared the performance of I L D M and T G L D M b y analys ing a hypothetical reaction system. The analysis showed that T G L D M w i t h op t imized in i t ia l condi t ions gives a better approximat ion of detai led chemistry than the I L D M under moderate perturbation. H u a n g constructed a T G L D M from a detai led reaction 41 } Manifold with MP method ; Manifold with current method j \ ••••••• Reaction trajectories ^ 4 o r o^mol/kg] 6 7 8 9 Figure 2.5: C o m p a r i n g one d imensional mani fo ld generated b y I L D M and T G L D M [1] mechanism [88] w h i c h contains 55 species and 278 elementary reac-tions. F ive major species ( C H 4 / 0 2 / C O , C 0 2 and H 2 0 ) were selected to construct a two d imens iona l mani fo ld . C 0 2 and H 2 0 were chosen as progress variables because of their relat ively l o n g and s imi la r t ime scales w h i c h avoids s k e w i n g the projected mani fo ld . E lement conser-va t ion for C , H a n d O is used i n E q . 2.66 to solve for the in i t i a l values. U n l i k e Pope a n d mass ' w o r k [86], w h i c h used transformed coordinates i n tabulating the mani fo ld , H u a n g proposed us ing Yco2 and YH2o as the table entry for mass fraction and reaction vectors directly. The mass fraction and reaction rate of any species and temperature can be i n -terpolated from the table b y us ing De launay t r iangulat ion w h i c h per-forms a poin t search and surface interpolat ion on a two-d imens iona l unstructured g r i d [89]. The density of the unstructured mesh varies 42 natural ly w i t h the change of stiffness i n the reaction system. H u a n g [1] app l i ed the constructed T G L D M to three different me-thane/a i r reaction systems: an unstrained laminar p r e m i x e d flame, a laminar non-premixed flame and a perfectly stirred reactor. The results are generally i n good agreement w i t h detailed chemistry, a l though as the phys ica l t ime scales approach the chemical time scales i n the fast subspace, relat ively large errors occur. H u a n g [60] also app l i ed T G L D M w i t h C S E to s imulate exper imen-tal methane jets under engine-relevant condit ions. It was found to be able to predict the experimental igni t ion delay and ign i t ion kernel lo -cations w e l l . The total N O ^ predic t ion was found to be p romis ing . 2 . 5 . 5 I C E - P I C Recently, R e n et al . p roposed an Invariant Cons t ra ined-equ i l ib r ium Edge Pre-Image C u r v e ( ICE-PIC) me thod [72,73] to reduce the detai led chemical mechan i sm w h i c h has the advantages of bo th the I L D M and the T G L D M . A s ment ioned before, the T G L D M generates an inert ial mani fo ld compar ing w i t h the I L D M ; however, it lacks the advantage the I L D M has: the I L D M is local ly determined. This proper ty is espe-cial ly important for constructing higher d imens ion manifo lds because an IS A T me thod [83] can then be used to save storage. The I C E - P I C me thod introduces a trajectory based on ideas from T G L D M and Rate-controlled Constra ined E q u i l i b r i u m ( R C C E ) [90-92]. Starting from the constrained chemical equ i l i b r i um calculated as one does i n R C C E , a man i fo ld is constructed b y so lv ing the O D E system as i n the T G L D M method. Then a Pre-Image C u r v e s me thod [93] is used to per form species reconstruction local ly f o l l o w i n g the reaction 43 t ra jec tory a l o n g the I C E m a n i f o l d . I C E - P I C m e t h o d h a s b e e n t e s ted i n a n i d e a l i z e d h y d r o g e n / a i r sy s -t e m w i t h s i x c h e m i c a l spec i e s w i t h the resu l t s c o m p a r e d to three o the r d i m e n s i o n - r e d u c t i o n m e t h o d s [72]. T h e resu l t s s h o w tha t the e r ro r i n -c u r r e d b y the I C E - P I C m e t h o d i s s m a l l ac ross the w h o l e f l a m e , e v e n i n the l o w t e m p e r a t u r e r e g i o n . It w a s a l so t e s ted i n a s t o i c h i o m e t -r i c m e t h a n e / a i r a u t o i g n i t i o n a n d the resu l t s are f a v o r a b l e i n c o m p a r -i s o n w i t h Q S S A a n d R C C E [73]. I C E - P I C has a l s o b e e n t e s ted i n a o n e - d i m e n s i o n a l p r e m i x e d l a m i n a r f l ame , w h e r e a " c l o s e - p a r a l l e l " as-s u m p t i o n i s i n t r o d u c e d w h i c h l e a d s to a " t r a n s p o r t c o u p l i n g " t e r m . T h e r e su l t s s h o w tha t th i s p r o c e d u r e s i g n i f i c a n t l y r e d u c e s the d i m e n s i o n -r e d u c t i o n e r ro r s [73]. W h i l e the I C E - P I C m e t h o d has s h o w n grea t p r o m i s e , i t h a s n o t y e t b e e n u s e d i n the s i m u l a t i o n o f a t u r b u l e n t f l a m e . 2.6 Probabilistic description of combustion P r o b a b i l i s t i c d e s c r i p t i o n o f c o m b u s t i o n , w h i c h u s u a l l y i n v o l v e s u s i n g a M o n t e C a r l o a p p r o a c h to s o l v e the c o m b u s t i o n s y s t e m , are a n at-t r ac t ive o p t i o n , e s p e c i a l l y for t u r b u l e n t c o m b u s t i o n . It c a n be u s e f u l for i n v e s t i g a t i n g s y s t e m s w i t h m a n y degrees o f f r e e d o m (e.g., P r o b a -b i l i t y D e n s i t y F u n c t i o n m o d e l ) , w h i c h are d i f f i c u l t to s o l v e w i t h f in i t e v o l u m e o r f in i t e d i f f e rence m e t h o d . A n o t h e r u s e f u l fea ture i s i t s n a t u -r a l g e n e r a t i o n o f f l u c t u a t i o n s ( i r r e g u l a r p a r t i c l e t ra jector ies) , w h i c h are d i f f i c u l t to o b t a i n u s i n g d e t e r m i n i s t i c m e t h o d s . I n th is s e c t i o n , t w o m o d e l s are d e s c r i b e d : the P r o b a b i l i t y D e n s i t y F u n c t i o n ( P D F ) m o d e l a n d the S tochas t i c P a r t i c l e M o d e l ( S P M ) , b o t h 4 4 of w h i c h ful ly u t i l i ze the advantages of the M o n t e Ca r lo method. There w i l l be greater emphasis on the S P M m o d e l since it is used extensively i n this thesis. 2.6.1 Probability Density Function (PDF) model The probabi l i ty densi ty function (PDF) method is a statistical approach to close the mean reaction rate problem. If the joint probabi l i ty densi ty function of a l l scalars is k n o w n , the mean reaction rate can be deter-m i n e d b y a s imple integration. For example, a joint p robabi l i ty density function P ( v & i , \ I / n ; x , t ) , where is the value of scalar i (density, temperature and mass fraction), is in t roduced. The mean reaction rate is: £ * . ( x , t ) = y , . . . y ' c J * i ( * 1 , . . . , * n ) P ( * 1 ) . . . ) * n ; x , t ) d * 1 . . . d * n . (2.67) N o w the quest ion is h o w to obtain P . One w a y is to assume the P D F shape as a g iven function f rom available mean quantities [4]. For ex-ample, P cou ld be assumed to be a c l ipped Gauss ian or /3-function de-termined b y the mean and variance of a variable. The two moments -mean and variance can be solved from a transport equation. Howeve r , the p resumed P D F usua l ly oversimplifies the actual P D F w h i c h causes large errors. Ano the r more elegant w a y is to solve a transport equat ion for the P D F [94-97] whose shape is a consequence of the f lu id m i x i n g and re-action. Pope [95,97] de r ived a conservation equation for the joint P D F of velocities a n d scalars. The one-point joint P D F P ( u , vi/; x , t)d\xd^ is the probabi l i ty that at t ime t and one spatial po in t x , the f l u i d has ve-loci ty components i n the range u — du/2 < u < u + ciu/2 a n d the values 45 of scalars i n the range * - d * / 2 < * < + The transport equa-tion is [95] + v • (PuP) + ( P g - vp) • v u p + ]T A [ W i P ] = i=l 1 n Pi V u • [< - V • r + V p ' | u , * > P] - ^ — [ < V • ( p D V ^ ) | u , * > P], (2.68) where g is the gravitat ional force, is the source-term for the scalar z, r is the stress tensor, the symbol V u is the divergence operator w i t h re-spect to the three components of the velocity, and the angular brackets denote condi t ional averages, condi t ioned w i t h respect to fixed values of u and * . A l l terms on the L H S are closed, representing the change of the P D F w i t h time, convect ion, transport i n the veloci ty space due to gravita-tion and pressure gradients and the chemical source-term, respectively. That the chemical source-term is closed is a most attractive feature of this mode l . The terms on the R H S , however, are not c losed and repre-sent respectively viscous stress and pressure fluctuations and the fluc-tuating part of the mic ro -mix ing . These unclosed terms are generated by assuming the one-point P D F is a sufficient descr ipt ion of the f low (neglecting spatial correlations). The molecular m i x i n g term is the most difficult one to m o d e l and several models have been proposed, i n c l u d -i n g the Interaction b y Exchange w i t h the M e a n ( IEM) or L inea r -Mean Square Es t imat ion models ( L M S E ) [6,98], Coalescence-Dispersion (C-D) m o d e l [99, 100] (also cal led C u r l ' s m i x i n g model) , m a p p i n g clo-sure [101] and Euc l idean M i n i m u m Spanning Trees (EMST) [102]. One characteristic of E q . 2.68 is its h igh dimensional i ty : the veloci ty 46 components and a l l scalars enter the P D F evolu t ion equat ion as inde-pendent variables. The finite vo lume and finite difference methods are usual ly not efficient i n so lv ing the P D F equation because the size of memory increases exponent ial ly as a function of d imens ion . Instead, most numer ica l solutions of the transport equation of P D F invo lve us-i n g M o n t e - C a r l o methods [95], where the m e m o r y requirement is l i n -early dependent o n the number of d imensions . In this approach, the P D F is represented b y a large number of stochastic particles. These particles evolve sequential ly i n time according to convect ion, chemical reaction, molecular transport and b o d y forces, and thus can simulate the evolu t ion of the P D F . In pract ical applicat ions, several k inds of P D F s are used w h e n con-s ider ing different sets of variables. In the simplest form, the joint P D F P (u , \ P ;x , £)• is reduced to a P D F for composi t ion on ly ( P ( * ; x , t ) ) to treat chemical reaction exactly; the veloci ty field is computed b y a stan-dard turbulence m o d e l (e.g., the k — e model) so lv ing the averaged Navier-Stokes equations [103-105]. This is also cal led the compos i t ion P D F method. Ano the r P D F method is based on the joint P D F of compo-sitions a n d veloci ty [106]. Alternat ively, a complete P D F me thod based on the joint P D F of velocity, compos i t ion and turbulence frequency has recently been proposed [107,108]. Cao et al . [107] used the joint P D F of velocity, compos i t ion and turbulence w i t h different chemical mecha-nisms to simulate the Sandia D , E and F flames; a S impl i f i ed L a n g e v i n M o d e l (SLM) for the evolu t ion of the particle veloci ty was used, and the stochastic frequency m o d e l of V a n Slooten et a l . [109] was used for the turbulence frequency of particles; the E M S T [102] m i x i n g m o d e l for 47 molecular diffusion on the composi t ion, w i t h an in situ adapt ive tabu-lat ion (ISAT) a lgor i thm [83] for implement ing the chemistry g iven b y the different mechanisms. H i s s imula t ion results have good agreement w i t h experiments w h i c h suggests that the joint P D F m o d e l can rep-resent the turbulence-chemistry effects i n these flames. H o w e v e r , the biggest disadvantage of joint P D F models is that the computa t ional cost is quite h igh . 2.6.2 Stochastic Particle Model The determinist ic approach to combust ion s imula t ion regards the t ime evolut ion as a continuous and w h o l l y predictable process w h i c h is governed b y a set of coupled , ordinary differential equations. It solves the reaction system as an in i t ia l value problem: x = Vp(x), x(t0) = x0. (2.69) where x(t) denotes a unique deterministic pa th based o n the in i t ia l condit ions and chemical kinetics, V denotes the stoichiometric vector and p(x) is the reaction rate. A solut ion for x can be so lved b y integrat-i n g E q . 2.69. The stochastic approach regards the time evolu t ion as a r andom-w a l k process w h i c h is governed b y a single differential equation, cal led the chemical master equation. Gi l lespie ' s stochastic me thod [110] starts from a chosen in i t ia l state. Then the probabi l i ty P that the next reaction w i l l occur after r seconds is computed. The step length r is chosen ac-cording to this probabi l i ty dis tr ibut ion. The next step is to determine what event w i l l happen next. Subsequently, this event occurs and the system is upda ted to a new state. This i teration continues un t i l cer-48 tain convergence criteria are met. The above procedure is also cal led a M a r k o v process. Bunker et al . [ I l l ] and Gi l lespie [110] proposed us ing stochastic methods to calculate homogeneous reaction systems separately. G i l l e -spie [112] conc luded that stochastic s imula t ion has significant advan-tages over the determinist ic method, i n that it can capture microscopic fluctuations. Howeve r , it also has several disadvantages; the execu-t ion time for the stochastic s imula t ion is usua l ly m u c h longer and i n -creases propor t iona l to the number of reaction channels and molecules of species to be s imulated. However , this fact is somewhat offset b y the smal l amount of memory it requires. The a lgor i thm requires a " re l i -able" uni t - in terval u n i f o r m r a n d o m number generator. H o w e v e r , there is no specific cr i ter ion to measure this reliability. The stochastic ap-proach also requires a large number of runs to estimate an ensemble average, w h i c h means it is computat ional ly expensive. Di rec t ly so lv ing the master equat ion is difficult; therefore Gi l lesp ie p roposed a stochas-tic s imula t ion approach. Later Gi l lesp ie gave a r igorous der iva t ion of the chemical master equation [113]. This approach is a mesoscopic de-scr ipt ion of the chemical reaction. The difference between the deterministic, or c o n t i n u u m approach and Gi l lesp ie ' s me thod is s h o w n i n the fo l lowing hypothet ical reaction system. A ^ B fc2 where k\ and k2 are the forward and backward reaction rates w h i c h are set to 10 and 1 respectively. Let x\ and x2 be the concentrations of A and B respectively. The in i t ia l condi t ion is x\(t = 0) = M and x2(t = 0) = 0. 4 9 The classic determinist ic method solves the O D E system: - h k2 Xi h -k2 - X 2 . (2.70) where xx + x2 = M. Subst i tut ing x2 b y x2 = M - xx generates a one d imens iona l O D E for xx: (2.71) (2.72) x\ = —kxxi + k2(M — x\). This has a single asymptot ical ly stable steady state solut ion: k2M X l fci + k2' Let P(X]_, t) denote the probabi l i ty of the concentration of A to be Xx at time t. A s s u m i n g the state of the system is Xx at t ime t and is X at time t + dt, one the fo l l owing must have occured d u r i n g t ime dt: • X = Xx + 1 and reaction A —> B happened • X = Xi — 1 and reaction B —> A happened • X = . X i , the concentration does not change, no reaction hap-pened Then the temporal evolu t ion of the probabi l i ty d is t r ibut ion for this s imple system can be approximated by: P(Xi,t + dt) =P(Xi + l,t)ki(Xi + l)dt (2.73) + P(Xi - 1, t)k2(M -Xi + l)dt + P(XU t)(l - kxXdt - k2(M - Xi)dt) Let t ing dt —> 0 results i n the chemical master equation: d P ( ^ b t ) =P(Xi + l^hiXi + 1) - P(Xut)kiXi (2.74) + P{Xi - I, t)k2{M -Xi + 1)- P(XU t)k2(M - Xx). 50 1 0 ' ' • 1 0 0 . 5 1 1 .5 Time F i g u r e 2.6: S o l u t i o n o f a H y p o t h e s i s s i m p l e s y s t e m ( T h i c k l i n e : de te r -m i n i s t i c m o d e l ; t h i n l i n e : G i l l e s p i e ' s m o d e l ) . T h e state shi f t v e c t o r s are vi — —1 a n d v2 — \; the p r o p e n s i t y f u n c t i o n s are ax[X) = kxXx a n d a2(Xl) = k2(M - Xx). A c o m p a r i s o n o f the s o l u t i o n s o f the c o n t i n u u m a n d G i l l e s p i e m e t h -o d s is s h o w n i n F i g . 2.6. W h i l e the c o n t i n u u m s o l u t i o n i s a s m o o t h c u r v e , the S P M c a n p r e d i c t the f l u c t u a t i o n s w h i c h are a r o u n d the c o n -t i n u u m s o l u t i o n . A s K r a f t ' s a n a l y s i s i n d i c a t e s , the d e t e r m i n i s t i c e q u a -t i o n i s d e r i v e d f r o m the s tochas t i c l i m i t w h e n the n u m b e r o f p a r t i c l e s goes to i n f i n i t y [114]. H e r e a n e x t e n s i o n o f G i l l e s p i e ' s a l g o r i t h m , i n c l u d i n g t r a c k i n g t e m -p e r a t u r e a n d e n e r g y as c a n b e f o u n d i n [115], i s g i v e n b e l o w . T h e d i s -51 crete m a s t e r e q u a t i o n d e s c r i b i n g a d i sc re te M a r k o v p r o c e s s is g i v e n ^ = ]C[^n,n'P„<(*) - Wn',nPn(t)] (2.75) n' w h e r e p n ( t ) d e n o t e s the p r o b a b i l i t y d e n s i t y o f a s tate n a t t i m e t. I n the state v e c t o r n = (XX X2 ... XM T), X i s the m o l e c u l e n u m b e r o f spec ies j o u t o f M t o t a l spec ies a n d T is t e m p e r a t u r e ; n' i s the s y s t e m state after e v e n t I o c c u r s c h a n g i n g the state f r o m state n; Wn^n i s the ra te o f c h a n g e w h i c h c a n be i n t e r p r e t e d as the f r e q u e n c y o f a n e v e n t I o c c u r r i n g w h i l e the s y s t e m is i n state n. A c c o r d i n g to G i l l e s p i e ' s a l -g o r i t h m the t w o i m p o r t a n t q u e s t i o n s w h i c h n e e d to b e a n s w e r e d are: w h e n w i l l the n e x t e v e n t h a p p e n a n d w h i c h e v e n t w i l l h a p p e n g i v e n the s y s t e m state n? T h e t i m e s t e p c a n b e d e t e r m i n e d f r o m the m a s t e r e q u a t i o n E q . 2.75 for the d e c a y o f a state n b y the s y s t e m l e a v i n g i t t o w a r d n' after e v e n t / ^ J T = " £  Wn',nPn(t) = - £ Wl>npn(t). (2.76) n' i T h e d i s t r i b u t i o n o f r e s i d e n c e t i m e r i s the e x p o n e n t i a l f u n c t i o n : Pn0(r) = Wnoe-w»°\ (2.77) w h e r e Wno i s the t o t a l f r e q u e n c y to l e a v e the state n 0 a n d Wno=YWl>»o- (2-78) / H e r e , Wi>no i s the ra te o f c h a n g e o f r e a c t i o n I at state n 0 . T h e a v e r a g e w a i t t i m e i s < r > = (2.79) T h i s e x p o n e n t i a l d i s t r i b u t i o n c a n b e g e n e r a t e d b y : r = - ^ - l n ( l - n ) , (2.80) 52 where 0 < r\ < 1 is a uniform random number. To determine which event will occur, the simplest method is to add up all relative frequencies and then multiply by a uniformly dis-tributed random number r 2 . The event that will be selected nseiected r a t e is then chosen by: "^selected rate Umax "^selected rate~r~l n = l n = l n = l where n m a x is the total number of chemical reactions in the system. This simple method works well for relatively small systems with few rates. The rate for chemical reactions is then calculated by the Arrhenius-type equation wr = K- c M , (2.82) where kr is the reaction constant, c; is the concentration of species i, and vr>i is the reaction coefficient of species i in reaction r. After nseiected rate is found, the system is updated from state n to n' by adjusting the molecule numbers: NiReactants * ^iReactants ^r,ifieactants (2.83) NiProducts * Niproducts -h VTtiReactants • (2.84) For the homogeneous system, the temperature change is also tracked by calculating through the change in the internal energy, the approxi-mate temperature rise, given by: A T r = - ^ - ( Y, vr,iHi(Toid) - Yl KilMToid) (2.85) v,old i=iProducts i—iReactants - RTold{NP The algorithm for a homogeneous reaction system with tracking of temperature can be summarized by the following steps: 53 1. A n in i t ia l state n is chosen; the chemical reaction rate of each pos-sible state change is then calculated. 2. A wa i t time r is determined b y E q . 2.80. 3. One of the reaction events is chosen b y E q . 2.81. 4. The chosen reaction is carr ied out, the system state is upda ted b y E q . 2.83 and the temperature rise is calculated b y E q . 2.85. This a lgor i thm is carr ied out repeatedly un t i l certain convergence cr i -teria are met. The S P M can be v i e w e d as a m u c h higher resolut ion s i m -ula t ion of a chemical system than one obtains us ing a c o n t i n u u m ap-proach, since it tracks every i n d i v i d u a l react ion. .However , it shou ld be noted that every t ime the a lgor i thm is i n v o k e d (for each realization), the system evolves through a different trajectory, since, every time it is i nvoked , one expects'different reactions to occur w i t h different w a i t i n g times. To improve the time performance of Gi l lesp ie a lgori thms, whose computa t ion cost increases propor t ional ly to the number of molecules of species and reaction channels, G i b s o n and Bruck proposed the "next reaction m e t h o d " [116]. G i b s o n identif ied the calcula t ion w h i c h is re-peated i n every i teration of the computat ion, such as i n step 1, where every chemical reaction rate is recalculated. H e then in t roduced sev-eral changes to improve the efficiency. For example, recalculat ing the chemical reaction rate on ly if it changes, s toring r and reus ing it where appropriate and us ing a dependency graph, etc. H o w e v e r , the depen-dency graph used to store the relationship between changes i n reaction rate to a g iven channel and indexed pr ior i ty queues is a trade-off i n this 54 algor i thm. Gi l l e sp ie proposed the "Tau-Leap m e t h o d " to save c o m p u -tation cost w i t h an acceptable loss i n accuracy [117]. The o r ig ina l G i l l e -spie a lgor i thm solves the chemical master equation exactly and obtains the exact t ime r and determines the exact event to happen. Howeve r , it is sometimes unnecessary to obtain such accuracy. Instead of f ind-ing out w h i c h reaction happens at w h i c h time step, Gi l l e sp ie proposed calculat ing h o w m a n y of each reaction are occurr ing i n a certain time interval . If the t ime interval chosen is large, significant computa t ion costs can be saved. However , this is a trade-off of detai l and compu-tation cost; usually, r needs to be relatively smal l to get acceptable re-sults. Gi l l esp ie further proposed an i m p r o v e d procedure for determin-i n g the m a x i m u m leap size for a specified degree of accuracy [118]. To offset the h i g h computat ional demand of the stochastic particle mode l , the a lgor i thm can be para l le l ized to u t i l ize cluster / g r i d tech-nology [119]. A h y b r i d a lgor i thm has also been proposed [120]; it par-titions the system into fast and s low reaction sub-systems, and then uses determinist ic me thod to handle the fast dynamics and Gi l l e sp ie algori thms to solve the s low dynamics . The appl ica t ion of stochastic algorithms based on [110] have spread into m a n y areas. Reaction-diffusion problems have been s tudied us ing stochastic algori thms i n conjunction w i t h components account ing for the diffusion process b y D a v i d [121] and Breuer et al . [122]. The forma-tion of soot us ing a coagulat ion reaction m o d e l has been s tudied [123]. Surface processes have been s tudied by us ing stochastic a lgori thms i n Frenklach's w o r k [124]. 55 2.7 Summary C o m b u s t i o n s i m u l a t i o n i s a m u l t i d i s c i p l i n a r y subjec t o f f l u i d d y n a m -ics a n d c h e m i s t r y . T h i s c h a p t e r h a s p r o v i d e d a r e v i e w o f the l i t e r a tu r e w h i c h i s r e l a t e d to the s t u d y i n th i s thes is . S e v e r a l o p p o r t u n i t i e s for n e w re sea rch c a n b e i d e n t i f i e d : • W h i l e C S E u s i n g L F D h a s b e e n a p p l i e d to R A N S , i t h a s n o t b e e n a p p l i e d to L E S . • C S E u s i n g a l o w d i m e n s i o n a l m a n i f o l d i n L E S h a s a l s o n o t b e e n a t t e m p t e d before . • T h e S P M offers a n e w w a y to genera te t ra jector ies i n a T G L D M w h i c h c o u l d offer s i g n i f i c a n t n e w p o s s i b i l i t i e s . I n the res t o f t h i s thes i s , these o p p o r t u n i t i e s for r e s e a r c h are e x p l o r e d . 56 r Chapter 3 Conditional source-term estimation with laminar flamelet decomposition in large eddy simulation of a turbulent non-premixed flame 3.1 Introduction This chapter focuses on the appl icat ion of C o n d i t i o n a l Source-term Est imat ion (CSE) w i t h Lamina r Flamelet Decompos i t ion (LFD) to the Large E d d y S imula t ion (LES) of the Sandia D-flame, a p i lo ted jet flame that has been s tudied extensively i n the literature. The purpose is to test this combinat ion of models against h i g h qual i ty exper imental data. 3.2 Formulation 3.2.1 Large Eddy Simulation For filter operat ion and balance equation of L E S , please refer to C h a p -ter 2 Section 2.3.3. The force tensor r is mode led by: r = A ( ( V v ) + ( V v f ) - H/iV • vS (3.1) where 8 is the uni t tensor. The molecular viscosi ty fj, is m o d e l l e d to be 57 a function of the filtered temperature fo l lowing a power l aw: fj. ~ f °-7. (3.2) The unclosed sub-gr id scale (SGS) fluxes are mode l ed u s ing e d d y vis-cosity and eddy diffusivi ty models: p ( w - vv) - - / w t ( V v + ( V v ) r - | ( V • v)I) (3.3) p{vZ - \Z) = -pDtVZ (3.4) where vt and Dt are the sub-gr id eddy viscosi ty a n d diffusivi ty re-spectively, bo th of w h i c h can be calculated by us ing the d y n a m i c mo-del [16,17]. The L E S balance equations are solved us ing the same code as was used i n the w o r k of Steiner and Bushe [56], itself a der ivat ive of a D N S code or ig ina l ly coded b y Boersma [125]. The spatial discret izat ion uses a second-order finite v o l u m e method w i t h a staggered structured g r id i n spherical coordinates. The g r id is 192 x 84 x 48 i n the rad ia l , tangen-tial and az imutha l direct ion respectively. The integrat ion i n t ime uses a second-order predictor-corrector-projector me thod [126]. The in f low boundary condi t ion is set according to experimental data for fuel, p i -lot and coflow quantities. The outf low boundary condi t ion is a convec-tive boundary w h i c h convects vort ical structures out of the computa-t ional doma in . A traction-free condi t ion is used for the lateral b o u n d -ary [125]. The code ut i l izes the Message-Passing Interface (MPI) [127] for para l le l iza t ion and was r u n on eight nodes of a P4 -Xeon B e o w u l f L i n u x cluster. 3.2.2 Laminar Flamelet Decomposition For details of C S E w i t h L F D , please refer to Chapter 2 Section 2.4.3. 58 2500i 6| 0 0 O1 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 z Figure 3.1: Three flamelets solutions of temperature, C O and tempera-ture source term It is possible to invert E q . 2.48 w i t h several different scalar fields w i t h i n each control v o l u m e for a c o m m o n a;. Grou t [59] found that the difference i n the temperature source term can be quite significant w i t h very litt le difference i n the condi t ional temperature. Therefore, i n E q . 2.48, / is a combina t ion of temperature and the C O mass fraction: F i g 3.1 shows the temperature, mass fraction of C O and tempera-ture source term of three flamelets from so lv ing the flamelet equation. The var ia t ion of temperature is smal l , but the var ia t ion of the corre-spond ing temperature source term is b ig . W i t h the higher var ia t ion of C O , invers ion w i t h both temperature and C O w i l l improve the accu-racy of basis function coefficient [59]. Altogether, at each time-step, the fo l lowing operations are under-taken: 1. E q . 2.49 is inver ted for a{ u s ing the no rma l i zed temperature and C O mass fraction. (3.5) 59 2. di is used to solve for the condi t ional mean source-terms of tem-perature and C O mass fraction. E q . 2.48 is used to obtain the un-condi t ional source-terms, w h i c h are then used to close the trans-port equations for temperature and C O mass fraction. 3. Transport equations for the filtered means of the mixture fraction, temperature and C O mass fraction are then integrated i n time to the next time-step. 3.2.3 Flamelet libraries A n unsteady flamelet l ibrary (UFL) was bui l t b y so lv ing the unsteady flamelet equations [32] s h o w n i n E q . 2.29 and 2.30 as was done i n p rev i -ous R A N S s imula t ion of the Sandia D-flame [59]. Differential diffusion effects was neglected since they have been found to have little in f lu -, ence on the results [57]. Grou t [59] estimated the mean scalar diss ipa-t ion rate his tory b y so lv ing the non-reacting f low field and the mean scalar d iss ipat ion rate at the stoichiometric interface as a function of distance f rom the nozzle . A mean scalar diss ipat ion rate i n Lagrang ian time was obtained b y assuming the stoichiometric interface was prop-agated downs t ream at the mean axis veloci ty of the stoichiometric i n -terface. The mean scalar diss ipat ion rate his tory for the Sandia D-flame can also be found i n Pitsch's w o r k [32]. To b u i l d a l ibrary, the reacting f low field was in i t i a l i zed w i t h a m i x i n g solut ion w i t h bounda ry con-ditions at Z = 0 match ing the coflow, at Z = 1 match ing the fuel f ield and at the p i lo t mixture fraction matching the pi lo t f ield. The bounda ry . condit ions were set to match both the fuel and ox id ize r streams. The flamelet was so lved un t i l it reached steady state at the far field mean 60 sca la r d i s s i p a t i o n ra te , t h e n the s ca l a r d i s s i p a t i o n ra te w a s i n c r e a s e d l i n e a r l y to gene ra t e a s o l u t i o n c o r r e s p o n d i n g to h i g h s t r a i n tha t l e a d s to e x t i n c t i o n . T h e s o l u t i o n s w e r e s a m p l e d at r e g u l a r i n t e r v a l s a n d w e r e s i m p l y c o n c a t e n a t e d toge the r i n t o the o n e u n s t e a d y f l a m e l e t l i b r a r y . A s t e a d y l i b r a r y - c o n s t r u c t e d w i t h s o l u t i o n s to E q s . 2.29 a n d 2.30 w i t h o u t the t i m e d e r i v a t i v e o n the L H S - w a s a l s o b u i l t . W h i l e a s t e a d y f l a m e l e t l i b r a r y ( S F L ) is n o t e x p e c t e d to p r e d i c t t r a n s i e n t p h e n o m e n a w e l l , i t s h o u l d s t i l l be ab l e to p r e d i c t the t e m p e r a t u r e a n d m a j o r spec ies s ince the t ime- sca l e s a s s o c i a t e d w i t h t he i r f o r m a t i o n are u s u a l l y m u c h s m a l l e r t h a n the f l o w t ime-sca le s . T h e sca la r d i s s i p a t i o n ra te x w a s v a r -i e d to get s o l u t i o n s f r o m e q u i l i b r i u m to e x t i n c t i o n w h i c h w e r e b u i l t i n t o a s t e a d y l i b r a r y . T h e c h e m i c a l k i n e t i c m e c h a n i s m u s e d i n b u i l d i n g the f l a m e l e t l i b r a r i e s w a s G R I - M e c h 3.0 [62]. 3.3 Simulation results and discussion T h e f u e l jet i n the S a n d i a D - f l a m e i s c o m p o s e d o f th ree p a r t s a i r a n d o n e p a r t C H 4 b y v o l u m e w i t h a n a x i a l v e l o c i t y o f 49.6 m / s . T h e d i a m e -ter o f the f u e l jet n o z z l e i s 7.2 m m a n d i t i s e n c l o s e d w i t h a n a i r c o f l o w h a v i n g a v e l o c i t y 0.9 m / s . A p i l o t i s p l a c e d i n the v i c i n i t y o f the jet to s t a b i l i z e the f l a m e . T h e jet has a R e y n o l d s n u m b e r o f 2 2 , 400 . A d e -t a i l e d d e s c r i p t i o n o f t h i s f l a m e a n d the e x p e r i m e n t a l d a t a a re d e s c r i b e d b y B a r l o w et a l . [128]. 3.3.1 Steady Flamelet Library Results F i g . 3.2 s h o w s the t i m e - a v e r a g e d c e n t e r l i n e p r o f i l e o f m i x t u r e frac-t i o n a n d t e m p e r a t u r e . A l l t e m p e r a t u r e p r o f i l e s are n o r m a l i z e d b y T0 = 291K. T h e s i m u l a t i o n r e su l t s for m i x t u r e f r a c t i o n agree w e l l w i t h the 61 N 0 . 5 0 ^ * A N A - - „ ' ^ \ A • / X / * ^ A \ / A \ At ' O s ' A . O 1 0 5t; 0 2 0 4 0 6 0 0 X /D Figure 3.2: Center l ine profile of temperature and mixture fraction (cir-cles: experimental mixture fraction; so l id l ine: calculated mixture frac-tion; triangles: n o r m a l i z e d experimental temperature; dashed line: cal-culated no rma l i zed temperature). experiment data. The temperature is somewhat over-predicted a long the centerline. F i g . 3.3 shows the temperature and mixture fraction predict ions at four downs t ream locations as functions of the radius. The mixture frac-tion is w e l l predicted throughout the f low field. The temperature over-predic t ion on the centerline extends to some other locations i n the flow. H o w e v e r it appears that the over-predict ion of temperature is worst at the centerline and, at some locations away from the centerline, the tem-perature is predic ted w e l l . F ig . 3.4 shows the Favre average of H 2 0 mass fraction at the four 62 Figure 3 . 3 : R a d i a l profiles of mixture fraction and temperature at four downst ream locations (circles: experimental mixture fraction; so l id l ine: calculated mixture fraction; triangles: no rma l i zed experimental temperature; dashed line: no rma l i zed temperature calculated w i t h SFL) . 63 downst ream locations w h i c h are generally w e l l predicted. There is a slight overpredic t ion of the peak average H 2 0 mass fraction and hence an overpredict ion of the H 2 0 mass fraction on the centerline at a l l four stations. F i g . 3.5 shows the condi t ional average of H 2 0 mass fraction. The location of the peak mass fraction is shifted s l ight ly to the r ich side. However , the magni tude of the peak is, for the most part, w e l l predicted. A s found i n earlier w o r k [56], i f the condi t iona l average of a scalar is wel l -predic ted, then its uncondi t iona l average is w e l l -predicted as s h o w n i n F i g . 3.4. This is an ind ica t ion that the mixture fraction mean a n d variance are reasonably wel l -predic ted b y the L E S s imula t ion a n d that the /3-function based o n these quantities is an ac-ceptable approx imat ion to the F D F . Figs. 3.6 through 3.9 show the condi t ional averages of var ious mass fractions at four different locations downs t ream obtained f rom the l i -brary us ing the coefficient basis function and integrat ing E q . 2.48. The predic t ion of C O is quite good , a l though this is the one species for w h i c h a transport equation was solved, so it w o u l d be expected that this species w o u l d be predicted w e l l . The predict ions of H 2 0 , 0 2 and C H 4 are also reasonably good, a l though the C 0 2 mass fraction is under-+-predicted at a l l stations. The O H mass fraction is reasonably w e l l predicted; the H 2 is somewhat over-predicted, a l though this m a y be attributed to i gno r ing the effects of differential diffusion. The N O mass fraction is s ignif icantly under-predicted at a l l stations, w h i c h is to be expected f rom the S F L results, since N O chemistry is s low and w i l l not s be w e l l represented w i t h a steady state calculation. 64 x/D=7.5 x/D=15 0.12i • . . . . 1 0.12, . . Figure 3.4: R a d i a l profiles of H 2 0 mass fraction at four downs t ream locations (circles: experimental Favre average of H 2 0 mass fraction; so l id l ine: n o r m a l i z e d Favre average of H 2 0 mass fraction calculated w i t h SFL) . 65 Figure 3.5: C o n d i t i o n a l averages of H 2 0 mass fraction at four stations downst ream of the nozz le (circles: experimental H 2 0 ; so l id l ine: H 2 0 calculated w i t h SFL) . 66 Figure 3.6: C o n d i t i o n a l averages of C O and C 0 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental C O ; so l id l ine: C O calculated w i t h S F L ; triangles: experimental C 0 2 ; dashed l ine: C 0 2 calculated w i t h SFL) . 67 x/D=30 0.25 0.15 x/D=45 0.05 0.4 0.6 Mixture fraction 0.4 0.6 Mixture fraction Figure 3 . 7 : C o n d i t i o n a l averages of C H 4 and 0 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental C H 4 ; so l id line: C H 4 calculated w i t h S F L ; triangles: experimental 0 2 ; dashed line: 0 2 calculated w i t h SFL) . 6 8 0.4 0.6 Mixture fraction 0.4 0.6 Mixture fraction . x 10 x/D=30 x10 ' x/D=45 0.4 0.6 Mixture fraction 0.4 0.6 Mixture fraction Figure 3.8: C o n d i t i o n a l averages of O H and H 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental H 2 ; so l id l ine: H 2 calculated w i t h S F L ; triangles: experimental O H ; dashed l ine: O H calculated w i t h SFL) . 6 9 Figure 3 . 9 : C o n d i t i o n a l averages of N O mass fraction at four stations downst ream of the nozz le (circles: experimental N O ; so l i d l ine: N O calculated w i t h SFL) . 7 0 0.14, - • '—, , _ , _ x/D=15 0.14. • . Figure 3.10: C o n d i t i o n a l averages of H 2 0 mass fraction at four stations downst ream of the nozz le (circles: experimental H 2 0 ; so l i d l ine: H 2 0 calculated w i t h U F L ) . 3.3.2 Unsteady Flamelet Library Results The s imula t ion results us ing the unsteady flamelet l ib ra ry ( U F L ) are s h o w n i n Figs. 3.10 through 3.14. For the most part, the predict ions are significantly i m p r o v e d b y us ing the U F L . The one notable excep-t ion is the pred ic t ion of N O . U s i n g the U F L , there is a significant over-predic t ion of N O . A t first glance, the N O results are somewhat d isappoin t ing . H o w -ever, the N O x chemistry i n G R I - M e c h 3.0 is k n o w n to have problems i n 71 Figure 3.11: C o n d i t i o n a l averages of C O and C 0 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental C O ; so l i d l ine: C O calculated w i t h U F L ; triangles: experimental C 0 2 ; dashed line: C 0 2 calculated w i t h U F L ) . 72 0.25, x/D=7.5 0.25j x/D=15 A Figure 3.12: C o n d i t i o n a l averages of C H 4 and 0 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental C H 4 ; so l id l ine: C H 4 calculated w i t h U F L ; triangles: experimental 0 2 ; dashed l ine: 0 2 calculated w i t h U F L ) . 73 Figure 3.13: C o n d i t i o n a l averages of O H and H 2 mass fractions at four stations downs t ream of the nozz le (circles: experimental H 2 ; so l id l ine: H 2 calculated w i t h U F L ; triangles: experimental O H ; dashed l ine: O H calculated w i t h U F L ) . 7 4 Figure 3 . 1 4 : C o n d i t i o n a l averages of N O mass fraction at four stations downst ream of the nozz le (circles: experimental N O ; so l id l ine: N O calculated w i t h U F L ) . 75 predic t ing N O [129] for this type of flame; indeed, i n other s imulat ions us ing C S E w i t h a l o w d imens iona l mani fo ld method to reduce the che-mist ry (which w i l l be discussed i n Chapter 4), s imi la r over-predictions of N O have been found us ing this chemistry [130]. In that context, the N O results s h o w n i n F i g . 3.14, be ing consistent w i t h predict ions ob-tained us ing var ious other methods, are actually encouraging. Thus, altogether, the results obtained w i t h the U F L are very good and are clearly superior to those obtained w i t h the S F L . 3 . 3 . 3 M i x e d L i b r a r y L o o k i n g at the predict ions of the N O field f rom the two different l i -braries, an average of these two fields w o u l d l i ke ly come ve ry close to the experimental results. Thus a mixed l ibrary was assembled b y con-catenating the two previous libraries together to fo rm a n e w library. One possible rationale for pu t t ing together such a l ib rary (aside from that it seems l ike that combinat ion might improve the N O predictions) is this: it is possible that, for a large fraction of the flow, the chemis-try is l i ke ly to be i n a near-steady-state, since the chemistry w i l l l i ke ly be fast relative to the local chemistry; however, it is also l i ke ly that some smaller fraction of the f low w i l l have shorter time-scales, so that flamelets f rom the unsteady l ibrary w o u l d be more appropriate. The C S E method w i t h L F D should , i n pr inciple , be able to select the right balance between these two scenarios; this m i x e d l ib rary provides an oppor tuni ty to-test the effect of the makeup of the flamelet l ib ra ry on the predict ions u s ing C S E w i t h L F D . The predict ions of temperature and major species remain good us-i n g the m i x e d library, be ing essentially identical to those obtained us-7 6 0 0.2 0.4 0.6 0.8 1 Mixture fraction 0 0.2 0.4 0.6 0.8 1 Mixture fraction 0 0.2 0.4 0.6 0.8 Mixture fraction 0 0.2 0.4 0.6 0.8 1 Mixture fraction Figure 3.15: C o n d i t i o n a l averages of N O mass fraction at four stations downst ream of the nozz le (circles: experimental N O ; so l i d l ine: N O calculated w i t h m i x e d flamelet l ibrary). i n g the U F L . The results s h o w n i n F i g . 3.15 indicate that i n c l u d i n g both steady and unsteady solutions does improve the pred ic t ion of mass fraction of N O significantly, par t icular ly at the two stations furthest from the nozzle . That a good predic t ion of N O has been obtained here b y s i m p l y m i x i n g two libraries together is l ike ly fortuitous, rather than indicat ive of some universa l ly v a l i d approach to b u i l d i n g an appropriate library. Indeed, as ment ioned above, the N O x chemistry i n G R I - M e c h 3.0 is k n o w n to have problems i n predic t ing N O . This br ings into quest ion 77 the va l id i t y of these results. B y combin ing the two libraries i n the arbi-trary w a y done here, the flamelet l ibrary m a y s i m p l y have been tuned to fit the experimental results. Thus, despite the results of the m i x e d l ibrary calculat ion p r o v i d i n g an excellent pred ic t ion of N O , this test has identif ied a potential shor tcoming i n the L F D approach: the results of the s imula t ion can be influenced b y wha t flamelets are i nc luded i n the library. One important a i m of L E S combust ion models shou ld be to avo id , wherever possible, the need to tune aspects of the m o d e l to match experimental results and, i n this respect, C S E w i t h L F D appears to fail . Altogether, it w o u l d seem that the or ig ina l formula t ion of C S E , us-ing , for example, a l o w d imens iona l mani fo ld to reduce the chemistry, is a superior approach to L F D . This other k i n d of approach does not require (or, i n fact a l low) tun ing of any parameters. 3.4 Conclusion C S E w i t h L F D was used i n L E S to obtain three different predict ions of the Sandia D-f lame. It is found that both the steady a n d unsteady libraries can predict the temperature and most species w e l l . The steady l ibrary l ed to an under-predic t ion of N O and the unsteady one led to an over-predict ion of N O . B y us ing a l ibrary compr ised of both steady and unsteady flamelets, the predict ion of N O is greatly i m p r o v e d . However , L F D has been s h o w n to have several disadvantages. D u e to l imitat ions of the flamelet mode l , it cannot deal w i t h phys ica l pro-cesses such as ign i t ion , extinction and re-ignit ion. P r io r informat ion is needed regarding the l i ke ly time-history of scalar d iss ipa t ion i n the 7 8 f low to b u i l d an unsteady flamelet library. Final ly , there appears to be a significant p r o b l e m w i t h L F D , i n that it is possible to influence the predict ions of the method b y altering the composi t ion of the flamelet l ibrary: a good predic t ion of N O was obtained w i t h a kinet ic mech-an i sm that is k n o w n to have problems predic t ing N O x chemistry b y m o d i f y i n g the compos i t ion of the flamelet library. Despite these problems w i t h the L F D method, it has been s h o w n here that excellent predict ions of a non-premixed flame are possible us ing the combina t ion of C S E and L F D . The combinat ion shows some promise that it c o u l d be used to simulate other non-premixed turbulent flames if used w i t h appropriate caution. 7 9 Chapter 4 Simulation of a turbulent non-premixed flame using conditional source-term estimation with trajectory generated low-dimensional manifold 4 . 1 I n t r o d u c t i o n In the last chapter, a significant shor tcoming i n the L F D approach to CSE was identif ied. In this chapter, Cond i t i ona l Source-term Es t ima-t ion (CSE) w i t h a Trajectory Generated L o w - d i m e n s i o n a l M a n i f o l d ( T G -LDM) is used i n Large E d d y Simula t ion of the Sandia D-f lame w i t h the hope that d i spens ing of the flamelet l ibrary w i l l e l iminate the short-coming w h i l e preserving the h i g h qual i ty of the predict ions. 4 . 2 M o d e l s 4 . 2 . 1 C S E - T G L D M In this study, the Trajectory Generated L o w - d i m e n s i o n a l M a n i f o l d ( T G -LDM) me thod [1,86] has been used to generate the table for C S E w h i c h shares the same basic idea as the Intrinsic L o w - d i m e n s i o n a l Man i f -old ( I L D M ) method [78]. Chemica l reaction processes cover a large 80 range of t ime scales, typica l ly ranging from 10~ 1 0 s to more than 1 s w i t h m a n y reactions and species. O n l y a smal l fraction of these scales can interact w i t h turbulence. Thus, the fast scales can be decoupled from the system. C o m p a r e d w i t h early efforts i n r educ ing the chemi-cal reacting system, such as the quasi-steady state and par t ia l equi l ib-r i u m assumptions, the I L D M and T G L D M methods on ly need the ful l kinetic mechanism^and the desired number of degrees of freedom of the reduced system as input . However , the I L D M me thod has diff icul-ties w i t h convergence and other problems as po in ted out b y Nafe and Maas [131]. In the T G L D M [1] method, the starting point of the trajectory is decided by some problem-specific op t imiza t ion criteria subject to the conservation of element mass fraction. Once the starting po in t is de-c ided, the trajectory is automatical ly evo lved i n the compos i t ion space b y so lv ing the chemical reaction system us ing the detai led chemistry w i t h a stiff o rd inary differential equation (ODE) system solver [132]. The trajectory w i l l eventual ly end at the e q u i l i b r i u m point . C o m p a r e d to the I L D M method , the T G L D M method has the advantage that it guarantees convergence a n d the reaction vector is a lways tangent to the trajectory. The T G L D M table is bu i l t i n advance of the reacting f low c o m p u -tation. Transport equations are then solved for continuity, m o m e n t u m , temperature, Y b c v YH2O and any other scalars of par t icular interest. A t each time-step i n so lv ing these transport equations, 1. The condi t iona l average of temperature a n d the two progress variables used for T G L D M lookup , < Yco2\i > a n d < Y # 2 o | £ >, 81 are obtained us ing E q . 2.47. 2. The condi t iona l reaction source-terms are interpolated f rom the T G L D M b y us ing the condi t ional values of the two progress va r i -ables. To overcome stability problems caused b y unsmooth inter-po la t ion i n the table, the timestep is spl i t into smaller t imesteps us ing D V O D E [132] so that the interpolat ion process migh t be done m a n y times. 3. The uncondi t iona l mean source-term is obtained from an equa-t ion l ike E q . 2.47 and substituted into the transport equation for temperature and species. Because on ly those condi t iona l values on each plane need to be interpolated, the computa t ional cost is relat ively low. Further computa t ional savings can be real ized b y interpola t ing for the condi t ional average of the mass fractions of species direct ly f rom the T G L D M table rather than so lv ing a transport equation. 4 . 3 I m p l e m e n t a t i o n 4.3.1 Large Eddy Simulation The code of Boersma [133], w h i c h was used i n an earlier imp lemen-tation of C S E i n L E S [56] and i n the w o r k described i n Chapter 3, has been further modi f ied . The discretization, m o d e l i n g approaches and boundary condit ions are the same as described i n Section 3.2.1. 4.3.2 T G L D M The T G L D M table was generated us ing < Yco2\€ > a n d < YH2O\€ > as progress variables because they have s imi lar and relat ively l ong 82 formation t ime scales, w h i c h makes the projection a n d interpolat ion more accurate [1]. M i x t u r e fraction space was discret ized w i t h forty-one points on a stretched g r id such that most points were clustered around the stoichiometric mixture fraction Zst = 0.351. The in i t ia l con-dit ions were defined to match the unburned mixtures i n the Sandia D-flame coflow and fuel jet. F ive major species C H 4 / 0 2 , C O , C 0 2 and H 2 0 were chosen to construct the in i t ia l starting point . D e p e n d i n g on the bounda ry condi t ion , each mixture fraction has be tween 38 a n d 57 trajectories. E a c h trajectory has 200 points evenly d is t r ibuted i n order to get a smooth interpolat ion f rom the table. i One difficulty i n generating a T G L D M for the current configurat ion occurs w h e n the in i t i a l temperature is low. In this zone, the separation among the var ious chemical reactions becomes so large that even the dedicated stiff O D E solver breaks d o w n . The general ly used method i n I L D M to solve this p rob l em is l inear prolongat ion w h i c h is not appro-priate i n this case because the r ight in i t ia l starting points are needed. Instead of l inear prolongat ion, the in i t ia l point is l inear ly extended b y a d d i n g a sma l l fraction of the e q u i l i b r i u m so lu t ion to raise the tem-perature un t i l the system can evolve automatical ly again. The T G L D M and the l o w temperature zone are i l lustrated i n F i g . 4.1. The advantage of this me thod is that it ensures the conservation of enthalpy. This l o w temperature zone is not important for most of the flame under consid-eration, but it can affect the results close to the nozz le where the strain rate is h i g h and local ex t inc t ion/ re igni t ion occurs. Both G R I - M e c h 3.0 and G R I - M e c h 2.11 [62] have been used to b u i l d tables that are approximate ly 765MB a n d 600MB i n b ina ry format w i t h 83 (a) (b) Figure 4.1: T G L D M l ibrary at stoichiometric mixture fraction (a). The enlarged g raph (b) shows the linear extended at l o w tempera-tures (dashed l ine is the linear prolongated part) al l mass fractions and chemical source-terms of species stored. The i n -terpolation was performed b y d i v i d i n g the who le compos i t ion space into sma l l triangles connected b y adjacent points. A g iven poin t is then found to be i n one certain triangle and the value is interpolated us ing the values at the three corners. 4 . 3 . 3 E x p e r i m e n t a l S e t u p The fuel jet i n the Sandia D-flame is composed of three parts air a n d one part C H 4 b y v o l u m e w i t h a veloci ty of 49.6 m / s . The diameter of the fuel jet nozz le is 7.2 m m and it is enclosed w i t h an air coflow hav ing a veloci ty 0.9 m / s . A pi lo t is p laced i n the v i c in i t y of the jet to stabilize the flame. The jet has a Reynolds number of 22,400. A detai led descript ion of this flame and experimental data is g iven b y B a r l o w et al . [128]. 84 Table 4.1: B reak -down of computat ional cost i n s imula t ion Ca lcu la t ion Transport C S E Inversion Chemis t ry Average time/<5£ (s) 6.11 4.27 5.31 Peak t i m e ( s ) 7 6 1923 4 . 4 R e s u l t s a n d d i s c u s s i o n 4 . 4 . 1 C o m p u t a t i o n a l c o s t A l t h o u g h the f low being s imulated is essentially s teady L E S is inher-ently unsteady. To compare w i t h the experiments, it was necessary to run the unsteady calculat ion unt i l the mean results converged, w h i c h took several weeks. Table 4.1 breaks d o w n the computat ional cost per timestep of solv-ing the cont inui ty and Navier-Stokes equations and five transport equa-tions for Z, T, N O , C 0 2 and H 2 0 w i t h invers ion and interpolat ion us-ing 8 processors and 2000 timesteps. W h i l e the computa t ional costs of integrat ing the transport equations and the invers ion for C S E have low standard deviations, the cost of chemistry ( interpolating w i t h i n the T G L D M , then integrating the chemistry w i t h D V O D E ) is ve ry h i g h in some timesteps. The reason for this is that the results interpolated from the table are not a lways sufficiently'smooth. Thus the O D E solver might choose smaller sub-timesteps to converge the so lu t ion i n each timestep. The t ime consumed on the invers ion for C S E and the subse-quent chemistry calculat ion is almost double that of so lv ing the trans-port equations. 85 86 0< • • 1 1 1 1 i In 0 10 20 30 40 50 60 70 80 x/D Figure 4.3: Center l ine profile of temperature and mixture fraction (o: measured Z; so l id l ine: calculated Z; A : measured temperature; — : calculated temperature) 4.4.2 Temperature and major species The fo l lowing s imula t ion results were found us ing G R I - M e c h 3.0 w i t h G R I - M e c h 2.11 g i v i n g s imi lar results. F i g . 4.2 shows the instantaneous temperature f ield. F i g . 4.3 shows that the centerline predict ions of mix -ture fraction and temperature are i n reasonable agreement w i t h the ex-per imental data. The agreement is also good as a function of the radius at different locations downst ream (Fig. 4.4), the results be ing s imi lar to those of Steiner and Bushe [56]. Cond i t i ona l averages of CO2 and H 2 0 are s h o w n i n F i g . 4.5. The agreement is good on the lean side and around the stoichiometric mix -ture fraction i n general. There is a slight underpredic t ion on the r ich side of stoichiometric at a l l locations downstream, especial ly at x/D = 7.5. 8 7 x/D=7.5 x/D=15 Figure 4 . 4 : Temperature and mixture fraction radia l profile at differ-ent locations downs t ream (o: measured Z; so l id l ine: calculated Z; A : measured temperature; — : calculated temperature) 8 8 X/D=7.5 0.12i • . 0 0.2 0.4 0.6 0.8 1 mixture fraction x/D=30 0.14i 0 0.2 0.4 0.6 0.8 1 mixture fraction x/D=15 0 0.2 0.4 0.6 0.8 1 mixture fraction x/D=45 0 0.2 0.4 0.6 0.8 1 mixture fraction Figure 4 . 5 : C o n d i t i o n a l averaged mass fraction of C 0 2 and H 2 0 at dif-ferent locations downs t ream ( o : measured H 2 0 ; so l id l ine: calculated H 2 0 ; A : measured C 0 2 ; — : calculated C 0 2 ) 8 9 0 . 1 4 Figure 4.6: Center l ine profile of mass fraction of C 0 2 and H 2 0 (o: mea-sured H 2 0 ; so l i d l ine: calculated H 2 0 ; A : measured C 0 2 ; — : calcu-lated C 0 2 ) Figs. 4.6 and 4.7 show reasonably good agreement w i t h the experi-mental data for the centerline and radia l profiles of the mass fractions of C 0 2 and H 2 0 . It is clear that i f the condi t iona l averaged va lue of species mass fraction can be correctly interpolated f rom the T G L D M table d i rec t ly then the spatial mean of the mass fraction of that species can be found b y m a p p i n g the corresponding condi t iona l values into phys ica l space. This suggests that the predic t ion of the filtered density function is reasonably good. F i g . 4 . 8 shows the condi t ional species mass fractions of C H 4 and 0 2 obtained b y interpolat ing from the T G L D M table. The pred ic t ion of C H 4 shows good agreement w i t h the experiment at a l l locations d o w n -9 0 Figure 4 . 7 : Mass fraction of C 0 2 and H 2 0 rad ia l profile at different locations downs t ream (o: measured H 2 0 ; so l id l ine: calculated H 2 0 ; A : measured C 0 2 ; — : calculated C 0 2 ) 91 x/D=7.5 0.25i • • 1 mixture fract ion x/D=30 0.3i • • — 0.25 0 0.2 0.4 0.6 0.8 1 mixture f ract ion x/D=15 0.25i . . 0 0.2 0.4 0.6 0.8 1 mixture f ract ion x/D=45 0.3. . . . — , 0.25, 0 0.2 0.4 0.6 0.8 1 mixture f ract ion Figure 4.8: C o n d i t i o n a l species mass fractions f rom interpolat ing the T G L D M l ibrary compared w i t h experimental data (o: measured C H 4 ; — : calculated C H 4 ; A : measured 0 2 ; — : calculated 0 2 ) stream. The 0 2 is w e l l predicted i n general, except at x/D = 7.5, where the predic t ion is poor, r i ch of the stoichiometric mixture fraction i n a manner s imi lar to wha t was seen i n F ig . 4.5. 4 . 4 . 3 Minor species A s ment ioned earlier, it is possible to interpolate the temperature and species mass fraction from the T G L D M table direct ly rather than solv-i n g transport equations to save computat ional cost. F i g 4.9 shows the condi t ional species mass fractions of H 2 , C O a n d the O H radica l ob-tained b y interpolat ing f rom the table directly. The interpolat ion of the 92 F i g u r e 4 . 9 : C o n d i t i o n a l spec ies m a s s f r ac t ions f r o m i n t e r p o l a t i n g the T G L D M l i b r a r y c o m p a r e d w i t h e x p e r i m e n t a l d a t a ( o : m e a s u r e d H 2 x 10; : c a l c u l a t e d H 2 x 10; A : m e a s u r e d O H x 10; — : c a l c u l a t e d O H x 10; +: L I F m e a s u r e d C O ; : c a l c u l a t e d C O ) c o n d i t i o n a l a v e r a g e o f O H agrees w i t h the e x p e r i m e n t a l d a t a i n g e n -e r a l b u t g i v e s a s l i g h t o v e r p r e d i c t i o n a r o u n d the s t o i c h i o m e t r i c m i x -ture f r a c t i o n . W i t h C O , a g o o d a g r e e m e n t o n the l e a n s i d e a n d a r o u n d the s t o i c h i o m e t r i c m i x t u r e f r a c t i o n is f o u n d , b u t C O i s , for the m o s t p a r t , u n d e r p r e d i c t e d o n the r i c h s i d e o f s t o i c h i o m e t r i c . I n c o n s t r u c t i n g the T G L D M , i t h a s b e e n i m p l i c i t l y a s s u m e d tha t a l l spec ies d i f f u s i v i t i e s are e q u a l . T h i s a s s u m p t i o n a p p e a r s to a d v e r s e l y affect the p r e d i c t i o n o f H 2 , w h e r e the p r e d i c t i o n sugges t s a s t r o n g p e a k i n H 2 j u s t r i c h o f the s t o i c h i o m e t r i c m i x t u r e f r a c t i o n w h i c h i s n o t p r e s e n t 9 3 i n the experimental results. A t the three stations closest to the burner, the experimental results a l l have a more diffusive character than our predict ion, consistent w i t h a higher species diffusivity. Thus , the pre-dictions of H 2 are not very good, par t icular ly under r i ch condit ions. Some of the condi t ional species mass fractions obtained b y interpo-la t ing f rom the T G L D M l ibrary/especial ly H 2 at x/D — 45, are not very smooth i n mixture fraction space. The reason is that the condi t ional species mass fraction is interpolated from a T G L D M generated for a part icular mixture fraction. In b u i l d i n g the T G L D M at each mixture fraction, a compromise has to be made i n representing the man i fo ld i n regions w i t h h i g h d imens iona l topology (cusps, folds etc). Specifically, wherever the man i fo ld folds, a single trajectory is - somewhat arbi-trari ly - selected to represent the mani fo ld i n that region. Pope [86] has suggested that these folds and cusps can be neglected i n the manner done here. A t different values of mixture fraction, different compro-mises are made. Thus , at adjacent mixture fractions, it is possible that a fold i n the man i fo ld is represented b y rad ica l ly different trajectories leading to an unsmoothness i n the mani fo ld i n mixture fraction space. 4 . 4 . 4 N O p r e d i c t i o n Because of the s low time scales of YNO compared to Yco2 a n d Yn2o, it is difficult to interpolate for the mass fraction of N O us ing just two progress variables, as can be seen i n F i g . 4.10. The N O mass fraction re-mains frozen as Y c o 2 a n d ^ H 2 O develop toward the e q u i l i b r i u m point at early stages. A s the trajectories evolve toward the e q u i l i b r i u m point , the mass fractions of C 0 2 and H 2 0 stabilize, but the N O mass frac-t ion increases sharply. To overcome this, a transport equation for N O 94 x 10 2-0.2 Figure 4.10: T G L D M l ibrary of N O mass fraction at stoichiometric mix -ture fraction is so lved w i t h the source-term obtained b y interpolat ing it f rom the T G L D M table. F ig . 4.11 shows that the predict ion of < YNO\€ > at a l l locations us ing the G R I - M e c h 3.0 mechanism significantly overpredicts the ex-per imental data. S imi la r f indings that G R I - M e c h 3.0 signif icantly over-predicts the N O mass fraction were also discussed at the Fif th Inter-national Workshop on Measurement and C o m p u t a t i o n of Turbulent N o n - p r e m i x e d Flames (TNF5) [134]. F i g . 4.12 shows that the predic-t ion of N O us ing G R I - M e c h 2.11 is greatly improved . 4.5 C o n c l u s i o n C S E - T G L D M has been used i n a L E S of the Sandia D-flame. It has been found that C S E - T G L D M can be used w i t h L E S i n non-premixed tur-bulent combus t ion to p rov ide good predict ions of temperature a n d of 95 Figure 4.11: C o n d i t i o n a l averaged mass fraction of N O at different lo cations downs t ream us ing G R I - M e c h 3.0 (symbols: measured; l ine: cal culated) 9 6 Figure 4.12: C o n d i t i o n a l averaged mass fraction of N O at different lo-cations downs t ream us ing G R I - M e c h 2.11 (symbols: measured; l ine: calculated) 9 7 al l species for w h i c h transport equations are so lved. Significant com-putat ional savings can be real ized b y interpolat ing for m i n o r species condi t ional mass fractions direct ly from the T G L D M table rather than so lv ing transport equations; these can then be m a p p e d back into the phys ica l space, where agreement w i t h experimental results is s t i l l very good. G o o d predict ions of N O require so lv ing a transport equat ion and interpolat ing the condi t ional source-term from the T G L D M table; i n this case, G R I - M e c h 2.11 predicts N O m u c h better G R I - M e c h 3.0. 9 8 Chapter 5 Extending the trajectory generated low-dimensional manifolds to low temperature with the stochastic particle model 5.1 Introduction W h e n the in i t i a l temperature of the methane /a i r system is l ower than 400K, the ign i t i on de lay t ime becomes exceedingly l o n g w h i c h makes it extremely difficult for a stiff O D E solver to obtain a solut ion. This makes it impossible to incorporate the autoigni t ion trajectory into a T G L D M under these condit ions. One w a y to c i rcumvent this p rob lem, as was done i n Chapter 4, is to shift the in i t ia l po in t s l ight ly into the composi t ion space b y m i x i n g the unburned mixture w i t h a sma l l frac-tion of the e q u i l i b r i u m composi t ion. This has the effect of increasing the temperature and reducing the in i t ia l stiffness i n the system. F r o m this in i t ia l point , the O D E solver has no difficulty e v o l v i n g the system to equ i l i b r i um and generating a trajectory. Between the o r ig in and this new in i t ia l point , s imple l inear interpolat ion was used. The resul t ing trajectory is fit into the man i fo ld as was s h o w n i n F i g . 4.1. W h i l e the me thod of m i x i n g unburned and e q u i l i b r i u m compos i -99 t i ons toge the r d o e s a l l o w u s to e a s i l y c a l c u l a t e a t rajectory, i t i s a n u n -s a t i s f y i n g s o l u t i o n to the p r o b l e m . It i s c l e a r f r o m F i g . 4.1 tha t the t ra -j e c t o r y i n i t i a t i n g at the o r i g i n takes a n a b r u p t c h a n g e i n d i r e c t i o n w h e n u s i n g the O D E s o l v e r - a n a b r u p t c h a n g e i n d i r e c t i o n that , l o o k i n g at o the r t ra jector ies i n the m a n i f o l d , does n o t a p p e a r to b e n a t u r a l . H e n c e , i t is l i k e l y tha t t h i s o n e t ra jec tory i s n o t o n the t rue m a n i f o l d w e seek, p a r t i c u l a r l y a t l o w t e m p e r a t u r e s . T h e l o w - t e m p e r a t u r e r e g i o n o f the m a n i f o l d m a y n o t b e i m p o r t a n t for a p p l i c a t i o n s w h e r e the c h e m i s t r y i s a l w a y s c l o s e to e q u i l i b r i u m , s u c h as i n n o n - p r e m i x e d f l a m e s far f r o m e x t i n c t i o n . I n cases w h e r e the l o w t e m p e r a t u r e r e g i o n o f the m a n i f o l d is i m p o r t a n t , h o w e v e r - as i n the case o f a f ree ly p r o p a g a t i n g p r e m i x e d f l a m e , for e x a m p l e - t h i s t ra jec tory i s l i k e l y g o i n g to l e a d to s i g n i f i c a n t e r rors . F u r t h e r , i t i s p o s s i b l e tha t the e r rors i n the t ra jec tory c o u l d be s i g n i f i c a n t e v e n at r e l a t i v e l y h i g h t e m p e r a t u r e s , there fore , there i s a n e e d to f i n d a n a l t e r n a t i v e s t r a t egy for m a k i n g a t r a jec to ry tha t d o e s n o t r e q u i r e the u n p h y s i c a l m i x i n g o f e q u i l i b r i u m a n d u n b u r n e d c o m -p o s i t i o n s . T h e S tochas t i c P a r t i c l e M o d e l ( S P M ) [115] uses the c h e m i c a l M a s -ter e q u a t i o n [113] to d e s c r i b e the r e a c t i o n s y s t e m i n i t s p r o b a b i l i s t i c n a t u r e , as e x p l a i n e d i n S e c t i o n 2.6.2. B e c a u s e i t u se s a M o n t e - C a r l o t e c h n i q u e , i t i s i m m u n e to p r o b l e m s o f st iffness i n the r e a c t i o n m e c h a -n i s m ; h o w e v e r , the c a l c u l a t i o n t i m e to get the S P M to a p p r o a c h c h e m i -c a l e q u i l i b r i u m is m a n y t i m e s grea ter t h a n i s n e e d e d w h e n s o l v i n g the c o n t i n u u m e q u a t i o n s (for e x a m p l e , i n c a l c u l a t i n g a u t o i g n i t i o n ) . I n t h i s chap te r , t he S P M w i l l b e u s e d to e x t e n d the T G L D M to l o w t e m p e r a t u r e s b y u s i n g i t to c a l cu l a t e the a u t o i g n i t i o n trajectory. I n d o -100 i n g so , a n u m b e r o f q u e s t i o n s a b o u t the a u t o i g n i t i o n t ra jec to ry w i l l b e e x p l o r e d that a r i se w h e n o n e c o n s i d e r s w h a t m i g h t b e l e a d i n g to the d i f fe ren t i g n i t i o n d e l a y t i m e s b e t w e e n d i f fe ren t r e a l i z a t i o n s : • D o d i f fe ren t r e a l i z a t i o n s o f the i g n i t i o n t ra jec tory f o l l o w a d i f fe r -ent p a t h t h r o u g h c o m p o s i t i o n space? • If t h e y d o f o l l o w a d i f fe ren t p a t h , h o w d i f f e ren t i s the p a t h ? D o the p a t h s c o n v e r g e o n t o a s i n g l e m a n i f o l d be fo re t h e y r e a c h the e q u i l i b r i u m c o m p o s i t i o n o r d o t h e y u l t i m a t e l y o n l y c o n v e r g e at e q u i l i b r i u m ? • If a n a v e r a g e t ra jec tory i s to be c o m p u t e d b y a v e r a g i n g m a n y re-a l i z a t i o n s together , d o e s the r e s u l t c o n v e r g e w i t h a n a c c e p t a b l y s m a l l n u m b e r o f r e a l i z a t i o n s ? D o e s the a v e r a g e t ra jec tory t ake the s a m e p a t h as the t ra jec tory that r e su l t s f r o m s o l v i n g the c o n -t i n u u m e q u a t i o n s ? 5.2 Calculating trajectories with the S P M A l l S P M resu l t s s h o w n i n th is p a p e r are for s t o i c h i o m e t r i c m i x t u r e s o f m e t h a n e i n a i r w i t h a n i n i t i a l p r e s s u r e P = 1 bar . a n d v o l u m e V = • (0.1 / i m ) 3 . T h e c h a r a c t e r i s t i c l e n g t h o f th i s v o l u m e is l a r g e r t h a n the m e a n free p a t h o f the m i x t u r e , w h i c h i s a p p r o x i m a t e l y O ( 1 0 ~ 3 nm)3 for a m i x t u r e o f i n i t i a l t e m p e r a t u r e b e t w e e n 2 9 8 K a n d 1 3 0 0 K . T h e c h e m i -c a l m e c h a n i s m u s e d i s G R I 2.11 [135]. S i n c e there i s a s i g n i f i c a n t c o m -p u t a t i o n a l cos t fo r e a c h r e a l i z a t i o n c a l c u l a t e d w i t h t he S P M , o n l y the t ra jec tory tha t d e s c r i b e s a u t o i g n i t i o n w a s s t u d i e d . F i g . 5.1 s h o w s t e n S P M c a l c u l a t i o n s i n ( C Q 2 , H 2 Q ) c o m p o s i t i o n space 101 Figure 5.1: C o m p a r i s o n of con t inuum and S P M results w i t h an i n i -t ial temperature of 1200K (Thick l ine: con t inuum trajectory; T h i n lines: S P M realizations). compared to the trajectory calculated w i t h the c o n t i n u u m l i m i t w i t h an in i t ia l temperature of T = 1200K. There appear to be fair ly w i d e v a r i -ations between the different realizations throughout their trajectory. In part, what is be ing seen here are fluctuations w i t h i n each trajectory that are associated w i t h the smallest time-scales i n the chemistry. In general, the ten S P M results fo l low roughly the same pa th through the composi t ion space as the con t inuum solut ion w h i c h is a lways i n the mids t of the S P M results. F ig . 5.2 shows ten S P M calculations i n ( C 0 2 , H 2 0 ) compos i t ion space compared to the trajectory calculated w i t h the c o n t i n u u m l i m i t (as was s h o w n i n F i g . 4.1) w i t h an in i t ia l temperature of T = 298K. Here , the S P M results and the con t inuum trajectory do not overlap over a broad range. O n l y w h e n the temperature has r isen significantly does the con-t i n u u m trajectory overlap w i t h the S P M results. A s s u m i n g that the S P M results are t racing the true trajectory, F i g . 5.2 suggests that the 102 0.1 0.08 2™ 0.06 0.04 0.02 0 0.02 0.04 \ O 0 6 • Y C 0 2 F i g u r e 5.2: C o m p a r i s o n o f c o n t i n u u m a n d S P M re su l t s w i t h a n i n i t i a l t e m p e r a t u r e o f 2 9 8 K ( T h i c k l i n e : c o n t i n u u m t ra jec tory; T h i n l i n e s : S P M r e a l i z a t i o n s ) . m i x i n g o f u n b u r n e d a n d e q u i l i b r i u m c o m p o s i t i o n s to c i r c u m v e n t the p r o b l e m w i t h c a l c u l a t i n g a u t o i g n i t i o n for the l o w i n i t i a l t e m p e r a t u r e causes a s i g n i f i c a n t e r r o r i n the t ra jec tory o v e r a w i d e r a n g e . T h e S P M re su l t s i n d i c a t e tha t the a u t o i g n i t i o n d e l a y t i m e for t h i s s y s t e m w i t h a n i n i t i a l t e m p e r a t u r e o f 2 9 8 K i s o n the o r d e r o f 1 0 2 6 s, o r 1 0 1 8 y e a r s , m e a n i n g tha t a r e a c t i o n i s e x t r e m e l y u n l i k e l y . H o w e v e r , after th i s e x t r e m e l y l o n g i n d u c t i o n t i m e , the i g n i t i o n p r o c e s s ( r a i s i n g the t e m p e r a t u r e f r o m jus t o v e r 2 9 8 K to o v e r 1 0 0 0 K ) t a k e s p l a c e o v e r a v e r y s h o r t t i m e ( I O - 3 s) . T h i s i s w h y the st i ff O D E s o l v e r i s u n a b l e to i n t eg ra t e the s y s t e m at t h i s l o w t e m p e r a t u r e . I n t r y i n g to i n c o r p o r a t e the S P M r e su l t s i n t o a T G L D M , a s i g n i f i c a n t p r o b l e m i s : a s i n g l e r e a l i z a t i o n u s i n g the S P M r e q u i r e s o v e r 100,000 even t s ; w i t h o u t p o s t - p r o c e s s i n g , a l l o f these w o u l d h a v e to b e s t o r e d i n the T G L D M , w h i c h w o u l d v a s t l y inc rease the s t o r age r e q u i r e m e n t s a n d u l t i m a t e l y defea t the p u r p o s e o f the T G L D M m e t h o d . T h u s , to 103 generate a usable trajectory f rom the S P M , it is necessary to either use only one real izat ion and a p p l y a fi l tering operat ion (to cut d o w n o n the number of points needed to represent the trajectory i n the manifold) -w h i c h w o u l d result i n a single smooth trajectory for the real izat ion - or assemble a sufficient number of independent runs of the S P M together to obtain a statistically converged average trajectory. 5.2.1 Average trajectories Averag ing trajectories together is not entirely t r iv ia l : as can be seen i n Figs. 5.1 and 5.2 the different trajectories can take significantly differ-ent paths through the composi t ion space, sometimes even reversing their direction. The different trajectories have very different ign i t ion delay times. In this work , piecewise linear least squares fits are taken to m a n y realizations f rom the S P M . The compos i t ion space is d i v i d e d into 200 discrete zones (using C 0 2 mass fraction) w i t h a h i g h densi ty of points at lower temperatures. M a n y realizations f rom the S P M s i m -ula t ion are generated w i t h the same in i t ia l condi t ion and the results are sorted into the different zones. Realizat ions are generated un t i l the least squares fits for a l l mass fractions of interest i n each zone have converged. The resul t ing fits are then evaluated at the zone centers to form an "averaged" trajectory. This technique is app l i ed first to the h i g h in i t i a l temperature of 1200K. E n o u g h realizations are averaged together to meet the conver-gence criteria. F i g . 5.3 shows the result of the least squares fit w h e n the m a x i m u m residual (for any species) is smaller than 0.01. The agree-ment between the c o n t i n u u m trajectory and the least squares fit of S P M is very good i n both the C 0 2 - H 2 0 and C 0 2 - O H compos i t ion spaces. 104 Table 5.1: C o m p u t a t i o n cost of for averaging temperature T = 1200K M a x i m u m Res idua l 0.40 0.20 0.01 Requi red real izat ion 90 273 2270 Similar results are found for other species as w e l l . The number of re-alizations required to meet var ious different levels of convergence is l is ted i n Table 5.1. The number of realizations needed increases dra-mat ical ly as the m a x i m u m residual decreases. To achieve a m a x i m u m residual smaller than 0.01 requires several days of compu t ing on a sin-gle 2 . 4 G H z processor. That the S P M results agree ve ry w e l l w i t h the con t inuum solut ion is a strong va l ida t ion of the approach used here, but the extreme computat ional cost suggests that this me thod shou ld on ly be used i n circumstances where other methods are not available. The autoigni t ion trajectory for l o w in i t ia l temperatures cannot be cal-culated us ing the con t inuum approach; however, the c o n t i n u u m ap-proach does w o r k w h e n the temperature can be raised (albeit artifi-cially). Thus, it is sufficient to calculate the trajectory w i t h the S P M only i n the l o w temperature zone near the o r ig in i n the compos i t ion space: once the temperature has r isen enough, i t is possible to swi t ch to us ing the c o n t i n u u m solut ion approach, w i t h the end po in t of the S P M result as a pseudo- ini t ia l point. F ig . 5.4 shows the result of us-i n g this procedure. The averaged S P M trajectory is calculated on ly up to the poin t i n the graph indicated b y the large square symbol . A f -ter that point , the trajectory is calculated w i t h the c o n t i n u u m solut ion. The n e w trajectory matches w e l l w i t h the S P M realizations s h o w n and clearly represents a significant improvement over the c o n t i n u u m tra-jectory calculated w i t h the combined unburned and e q u i l i b r i u m in i t ia l 105 0 0 . 0 2 0 . 0 4 0 . 0 6 v T co 2 (a) x 1 0 (b) Figure 5.3: Piece wise least squares fit of S P M realizations w i t h in i t ia l temperature 1200K; (a) C 0 2 - H 2 0 ; (b) C 0 2 - O H (Line: c o n t i n u u m tra-jectory; Circles: piecewise least squares fit). 106 Figure 5.4: Piecewise least squares fit w i t h l o w temperature (T=298K) i n the C 0 2 and H 2 0 phase space (Thin lines: ten S P M realizations; Dashed l ine: trajectory starting from arbitrary mixture of unburned and e q u i l i b r i u m composit ions; So l id line: combined S P M / c o n t i n u u m trajectory; Square symbol : the poin t where the trajectory changes from S P M to c o n t i n u u m solution). composi t ion. 5 . 3 Validation Two cases, a p r emixed laminar flame and a perfectly st irred reactor, are chosen to test the performance of the SPM-extended T G L D M . 5 . 3 . 1 L a m i n a r p r e m i x e d f l a m e A freely propagat ing laminar p remixed flame is chosen to test the n e w T G L D M . The govern ing equations described i n Section 2.3.1 are so lved 107 i n one d imens ion us ing a second-order centered finite difference scheme The bounda ry condit ions for the unburned mixture are Yi = YUji, T = TU at x —• - c o . (5.1) where F denotes the mass fraction of species i, Yu<i denotes the unburn t mass fraction of species i, Tu denotes the unburnt temperature. The boundary condit ions for the burned end are dYi dT - 7 ^ = 0, — = 0 at x —> oo. (5.2) ox ox Both the T G L D M calculated us ing the mixture of u n b u r n e d a n d equi-l i b r i u m composi t ions i n the place of the autoigni t ion trajectory a n d the SPM-extended T G L D M tables have been used for m o d e l i n g Co{. The predic t ion of temperature a n d major species u s i n g bo th tables match w e l l w i t h the detai led chemistry a n d are quite s imi lar ; there-fore, on ly the predict ions obtained w i t h the SPM-ex t ended T G L D M are s h o w n i n F i g . 5.5. Some minor species such as H a n d C H are w e l l predicted b y bo th models , as s h o w n i n F i g 5.6a. The O H mass fraction is overpredicted w h i c h migh t be because of the lack of diffusion i n the T G L D M . The pred ic t ion of C H 3 0 is better w i t h the SPM-ex tended T G L D M than w i t h the T G L D M calculated u s i n g the mix ture of un -bu rned a n d e q u i l i b r i u m composit ions i n the place of the autoigni t ion trajectory, as s h o w n i n F i g . 5.6b, but s t i l l underpredic ted b y half w h e n compared w i t h the detai led chemistry result. The flame veloci ty w i t h the T G L D M calculated u s ing the mixture of unburned and e q u i l i b r i u m composit ions i n the place of the autoigni-tion trajectory is 20.33 c m / s , w h i c h is under-predicted compared w i t h detailed chemistry, w h i c h gives 36.83 c m / s . The SPM-ex tended T G -108 0.25 (a) 2 5 0 0 2 3 4 5 x [ m m ] (b) Figure 5.5: Scalar profiles through the laminar p remixed flame w i t h in i t ia l reactant temperature of 298K: (a) Major species; (b) Temperature (Squares: SPM-ex tended T G L D M ; So l id l ine: Deta i led chemistry). 109 Figure 5 . 6 : Scalar profiles through the laminar p remixed flame w i t h in i t i a l reactant temperature of 298K: (a) M i n o r species (Sol id l ine: D e -tailed chemistry; Symbols : SPM-extended T G L D M ) ; (b) C H 3 0 (Solid l ine: Deta i led chemistry; • : SPM-ex tended T G L D M ; x : L i n e a r l y inter-polated T G L D M ) . 110 L D M i m p r o v e s the p r e d i c t i o n to 22.27 c m / s . A t f i rs t g l a n c e , th i s m o d -est i m p r o v e m e n t i s s o m e w h a t d i s a p p o i n t i n g . H o w e v e r , i t is k n o w n that, at l o w t e m p e r a t u r e s , the d i m e n s i o n o f the I L D M i n b o t h s t r a i n e d a n d u n s t r a i n e d p r e m i x e d f l a m e s inc reases r a p i d l y b e c a u s e i t i s h a r d to d e c o u p l e fast sca les a t l o w t e m p e r a t u r e s [136]. S e v e r a l d i f f e r en t m e t h -o d s h a v e b e e n p r o p o s e d to a d d r e s s th i s p r o b l e m i n I L D M m e t h o d s , i n -c l u d i n g the F l a m e P r o l o n g a t i o n o f I L D M (FPI) m e t h o d [80], the F l a m e -l e t - G e n e r a t e d M a n i f o l d m e t h o d ( F G M ) [81] a n d P h a s e - S p a c e I L D M (PS-I L D M ) [82]. T h e s e m e t h o d s e a c h u s e the s o l u t i o n o f o t h e r f l a m e s w h i c h u s u a l l y h a v e d i f f u s i o n t e r m s a n d t h e n t abu la t e the r e s u l t s i n t o a m a -n i f o l d . T h e T G L D M m e t h o d , l i k e the o r i g i n a l I L D M m e t h o d , d o e s n o t i n c l u d e d i f f u s i o n . T h u s , i t i s l i k e l y tha t the i m p r o v e m e n t i n the p r e -d i c t i o n o f the p r e m i x e d f l a m e s p e e d i s l i m i t e d b e c a u s e the l o w d i -m e n s i o n a l T G L D M i s f u n d a m e n t a l l y i n s u f f i c i e n t to c a l c u l a t e the f l a m e v e l o c i t y a c c u r a t e l y at s u c h a l o w t e m p e r a t u r e . I n d e e d , as the i n i t i a l t e m p e r a t u r e decreases , the a c c u r a c y o f the p r e d i c t i o n o f v e l o c i t y d e -creases [1]. 5.3.2 Perfectly stirred reactor calculations T h e p e r f e c t l y s t i r r e d r eac to r ( P S R ) i s a n i d e a l i z e d c o n t i n u o u s l y - s t i r r e d f l o w t a n k reactor . I n f l o w c o n t i n u o u s l y f l o w s i n t o the r eac to r to m i x w i t h the f l u i d i n s i d e the r eac to r w h i l e r e a c t i o n p r o d u c t s a re e x h a u s t e d . T h e P S R i s a v a l u a b l e test case for r e d u c e d c h e m i s t r y a n d i t i s u s e d he re to test w h e t h e r o r n o t the a v e r a g e d S P M a u t o i g n i t i o n t ra jec tory i m p r o v e s the p e r f o r m a n c e o f the T G L D M at l o w t e m p e r a t u r e s . T h e 111 balance equat ion for P S R is g iven by: dY — = F ( Y ) - w ( Y - Y 0 , (5.3) where Y = [Yi Y 2 ••• YW] is the mass fraction vector and N is the total number of species; F is the reaction rate vector; to represents the rate at w h i c h hot combust ion gases f low into the reactor; Y ; is the mass fraction vector of the inflow. The in i t ia l state of the P S R is set to be a stoichiometric methane /a i r mixture w i t h a temperature of 298K; the system pressure is set to 1 bar. If the in f low is set to be the same as the in i t ia l state, the system w i l l not be able to ignite w i t h i n an acceptable time; therefore an in f low w i t h hot product gases - the equ i l i b r ium composi t ion and temperature of the mixture i n the in i t i a l state - was cont inuously fed into the system. Over t ime, the temperature rises because of the m i x i n g w i t h the hot in f low and eventual ly the system ignites. For this test case, ui = 200 s - 1 is chosen, w h i c h leads to igni t ion i n the P S R after approximate ly 5 ms. The results of three P S R calculations are s h o w n i n F i g . 5.7. The or ig-ina l T G L D M , where the autoigni t ion trajectory is calculated in i t i a l ly b y m i x i n g the unburned and equ i l i b r ium composi t ions together, over-predicts the ign i t ion time for the P S R - this despite the fact that this PSR test is very s imi la r i n its construction to the w a y this man i fo ld was constructed. The T G L D M incorporat ing the averaged S P M trajectory for au to igni t ion provides a significant improvement i n the P S R cal-culat ion, reproduc ing the results w i t h fu l l chemistry almost perfectly. This indicates that the n e w method for constructing the autoigni t ion trajectory i n the l o w temperature T G L D M is very successful. 112 3000 (b) Figure 5 . 7 : P S R calculat ion results: temperature a n d mass fractions of CO2 and O H (sol id l ine: detailed chemistry; dash dot l ine: averaged S P M T G L D M ; dashed line: T G L D M calculated u s ing arbitrary mixture of unburned and e q u i l i b r i u m composit ions i n place of au to igni t ion tra-jectory); (b) is an enlargement of (a). 113 5.4 Conclusions I n th i s s t u d y , the S tochas t i c P a r t i c l e M o d e l h a s b e e n u s e d to gener -ate t rajectories for Tra jec tory G e n e r a t e d L o w D i m e n s i o n a l M a n i f o l d s . W h e n m a n y t ra jector ies are a v e r a g e d together , t he r e s u l t c o n v e r g e s w i t h a n a c c e p t a b l y s m a l l n u m b e r o f r e a l i z a t i o n s a n d th i s a v e r a g e t ra -j e c t o r y t akes the s a m e p a t h as the t ra jec tory tha t r e su l t s f r o m s o l v i n g the c o n t i n u u m e q u a t i o n s . T h u s , there a p p e a r s to b e a n a l t e rna te m e -t h o d to s o l v i n g the c o n t i n u u m e q u a t i o n s for o b t a i n i n g the a u t o i g n i -t i o n t ra jec tory tha t i s n o t af fected b y st iffness i n the c h e m i c a l k i n e t i c m e c h a n i s m o f f e r e d b y the S P M . It h a s b e e n s h o w n tha t t he m e t h o d c o m b i n i n g u n b u r n e d a n d e q u i l i b r i u m c o m p o s i t i o n s to a v o i d st iffness i n the a u t o i g n i t i o n t ra jec tory at e x t r e m e l y l o w i n i t i a l t e m p e r a t u r e s i s c l e a r l y d e f i c i e n t a n d tha t the a v e r a g e d S P M t ra jec tory p r o v i d e s a v i -ab le a l t e r n a t i v e . T h e a v e r a g e d S P M trajectory, w h e n i n s e r t e d i n t o the T G L D M , l e a d s to a n e x c e l l e n t p r e d i c t i o n o f a P S R w h e n c o m p a r e d to the P S R c a l c u l a t e d w i t h f u l l c h e m i s t r y , w h e r e a s the p r e d i c t i o n u s i n g the T G L D M c o n s t r u c t e d u s i n g the m e t h o d c o m b i n i n g u n b u r n e d a n d e q u i l i b r i u m c o m p o s i t i o n s i n t r o d u c e d a s i g n i f i c a n t e r r o r i n the P S R c a l -c u l a t i o n s . It offers a l i m i t e d i m p r o v e m e n t i n the p r e d i c t i o n o f a l a m i -n a r p r e m i x e d f l a m e w h i c h m i g h t be a t t r i b u t e d to the i n e f f i c i e n c y o f the t w o - d i m e n s i o n a l m a n i f o l d u s e d here . 114 Chapter 6 Prediction of autoignition fluctuation using trajectory generated low-dimensional manifolds generated using the stochastic particle model 6.1 Introduction In Chapter 5, the Stochastic Particle Model (SPM) was used for extend-ing the Trajectory Generated Low-dimensional Manifold (TGLDM) me-thod to low temperatures by averaging many realizations together to obtain one trajectory. In this chapter, the potential for using individ-ual SPM realizations in a TGLDM will be explored. Here the questions being addressed are: • Is there some way that a single realization can be incorporated into a TGLDM? • If TGLDMs incorporating different realizations of the autoigni-tion trajectory are used in CFD calculations of turbulent autoigni-tion, do the results of the different calculations give different ig-nition delay times? • If so, do different realizations used in CFD lead to fluctuations of the same order as are seen in experiments? 115 6.2 Results and discussion 6.2.1 Filtering individual realizations Filtering realizations from the S P M offers the possibility to compare different realizations, to study the effects of chemical fluctuations on the manifold and possibly to include the effects of chemical fluctua-tions into C F D calculations by including different realizations from the S P M into different manifolds. The very premise of low dimensional manifolds is that, when the smallest time-scales are eliminated, the resulting system can be well approximated by a lower dimensional manifold. Fig. 6.1 shows the wait times for every chemical step as a function of the event number for a single realization of the S P M . Fig. 6.2 shows the wait times as a function of the C 0 2 mass fraction. These figures show very clearly how vast is the range of time-scales in the system during autoignition. They also shows how there is a clear separation between long and short time scales. There are three main zones in which the different wait times fall. There is a very long time-scale zone, initially on the order of 1026s (these points only occur very early in the realization so that the points are obscured by the left hand axis in the graph in Fig. 6.1). There is a somewhat shorter intermediate time-scale zone initially on the order of 102os. Finally, there is a very short time-scale zone with time-scales shorter than 10 _ 1s. Ideally, very short time-scales should be filtered away. A filter of a size within the intermediate time-scale zone would be preferred. U n -fortunately, there is no fixed time-scale that one could use that would be appropriate throughout the evolution of the trajectory: as the mass 116 0 1 2 3 4 5 Events . „ 4 x 10 Figure 6.1: The wa i t t ime (ms) for every chemical reaction event as a function of the event number for a single real izat ion of the S P M w i t h an in i t i a l temperature of 298K 0 2 4 v 6 8 10 Figure 6.2: The wa i t t ime (ms) for every chemical reaction event as a function of the C 0 2 mass fraction for a single real izat ion of the S P M w i t h an in i t i a l temperature of 298K 117 fraction of C 0 2 increases, the temperature rises and the reaction rates increase, w h i c h causes the two longer time-scale zones to merge and drop towards the sma l l time-scales. A t the time of the early events i n the realizat ion, filter at a time-scale of a round 10 1 9 s is expected; later that filter w o u l d be far too broad. A n alternative to direct t ime-filtering is to filter i n event space. A species mass fraction can be expressed as a function of the number of chemical events n (as the wai t time was plot ted i n F i g . 6.2) and the filter can be app l i ed i n n space: Y(n)=. Y{n')g(n,n')dn' (6.1) Jo where g(n,n') w o u l d be the filter kernel and N w o u l d be the total number of events available. F i l te r ing i n this w a y has a number of ad-vantages. Foremost among these is that the event space filter can be v i e w e d as a time-filter where the filter scale changes according to the local time-scales since the local wa i t t ime is calculated according to the local time-scales. It is also a filter that is s imple to a p p l y to the S P M results. The filter kernel w i l l be taken to be a centered top-hat function. The system tends to have long wai t times fo l lowed b y a large n u m -ber of very short wa i t times. In selecting a d imens ion Ne (the number of events o n either side of the top-hat function) for the filter, i t is found that it must be larger than the number of events between the events w i t h l ong wa i t times. Otherwise , the filter is ineffective at e l imina t ing the smal l time-scales. To meet this cri terion for the realizations s h o w n i n Figs. 5.1 a n d 5.2, Ne s h o u l d be larger than 54. A n alternate me thod for establishing the m i n i m u m value for Ne is to consider the d is t r ibut ion of wai t times as r a n d o m numbers that are 118 samples from an underlying distribution - indeed, they are: they fol-low the exponential distribution - and calculate the number of sam-ples needed to calculate the average of the distribution with some error bound and a degree of confidence. According to the central l imit theo-rem, the average value approaches a normal distribution if the sample size N is big enough. Assuming that Wn = 1 in Eq. 2.77, choosing a confidence level of 0.95 and absolute error \Xn - n\ of 0.1, Ne = 384. If Ne = 54, the error is 0.27. In order to establish the optimal value of Ne, different values of Ne are applied to the trajectories. The first test considers how well reac-tion rates are predicted (given that these must be well predicted for the T G L D M to be useful); this test is shown in Fig. 6.3. The original un-filtered reaction rates from this trajectory show very large fluctuations, especially i n the source-term for temperature. A p p l y i n g a filter with Ne — 54 (not shown) does have some smoothing effect, but the trajec-tory in (CGs, H2O) space still fluctuates excessively: the main goal in applying the filter is to smoothen the trajectory sufficiently that it can be described wi th a reasonably small number of points. Filtering with Ne = 1080 or 2700 meets this goal and the resulting trajectories match the reaction rates reasonably well . Interestingly, these two curves are difficult to distinguish from one another in Fig. 6.3. It is also interest-ing to see that these curves are somewhat similar to the continuum solution. When the filter dimension is raised further to Ne = 5400, it is seen that the prediction of the temperature source-term no longer fol-lows the trend in the unfiltered data. Clearly, there is a maximum value of i V e that should be used to avoid overly smoothing the trajectory. 119 (b) F i g u r e 6.3: R e a c t i o n ra tes f i l t e r e d w i t h d i f fe ren t f i l t e r s i z e for a n i n i t i a l t e m p e r a t u r e o f 1 2 0 0 K ( G r a y do ts : o n e u n f i l t e r e d r e a l i z a t i o n f r o m the S P M ; g r a y s o l i d l i n e : Ne = 1080; b l a c k s o l i d l i n e : JV e = 2700; d a s h - d o t l i n e : = 5400; d a s h l i n e : the c o r r e s p o n d i n g c o n t i n u u m trajectory;) . 120 F i g . 6.4 s h o w s the r e s u l t s o f f i l t e r i n g t e n d i f f e ren t S P M tra jector ies w i t h NE = 1080 a n d NE = 2700. W i t h NE = 2700, a l l t e n t ra jector ies d e p a r t s i g n i f i c a n t l y f r o m the c o n t i n u u m s o l u t i o n a r o u n d YQO2 = 0 . 0 1 5 . G i v e n th i s , i t i s c o n c l u d e d that - at leas t fo r t h i s s y s t e m - NE s h o u l d b e set to 1080. W i t h NE = 1080, a l l o f the t ra jector ies s t i l l f o l l o w the c o n t i n u u m t ra jec tory r e a s o n a b l y c lose ly , a l t h o u g h there are s i g n i f i c a n t v a r i a t i o n s b e t w e e n e a c h f i l t e r e d trajectory. 6.2.2 Predicting fluctuations in autoignition A s w a s m e n t i o n e d i n C h a p t e r 2, s i m u l a t i o n s u s i n g a R A N S a p p r o a c h w i t h C S E a n d T G L D M to o b t a i n c l o s u r e for the c h e m i c a l s o u r c e - t e r m s h a v e d e m o n s t r a t e d a r e a s o n a b l e p r e d i c t i o n o f the a v e r a g e i g n i t i o n d e -l a y t i m e o f t r a n s i e n t t u r b u l e n t jets; the c o m b i n a t i o n , b y i t s n a t u r e , w a s n o t ab l e to p r e d i c t a n y f l u c t u a t i o n s . E x p e r i m e n t a l d a t a i n d i c a t e that there are s i g n i f i c a n t f l u c t u a t i o n s i n the i g n i t i o n d e l a y t i m e . T h e r e are t w o p o s s i b l e sou rces o f the f l u c t u a t i o n s i n the i g n i t i o n d e l a y t i m e . 1. I g n i t i o n o f a n o n - p r e m i x e d s y s t e m t y p i c a l l y r e q u i r e s tha t the l o -c a l s ca l a r g r a d i e n t s are s u f f i c i e n t l y l o w for a s u f f i c i e n t l y l o n g t i m e - e x a c t l y h o w l o w a n d for h o w l o n g d e p e n d s o n the che -m i s t r y g o v e r n i n g the s y s t e m . I n a t u r b u l e n t f l o w , the l o c a l s ca l a r g r a d i e n t s f l uc tua t e w i t h the l o c a l v e l o c i t y f i e l d . T h e r e f o r e , tu r -b u l e n t f l u c t u a t i o n s c a n i m p e d e the i g n i t i o n p r o c e s s ; v a r i a t i o n s i n the t u r b u l e n t f l u c t u a t i o n s b e t w e e n d i f fe ren t r e a l i z a t i o n s o f the f l o w f i e l d c o u l d l e a d to f l u c t u a t i o n s i n the i g n i t i o n d e l a y t i m e . 2. O n the o t h e r h a n d , there are a l so s i g n i f i c a n t f l u c t u a t i o n s b e t w e e n 121 0 0 . 0 2 0 . 0 4 0 . 0 6 Y c o 2 (a) 0 . 1 -— 0 . 0 8 Mr 0 . 0 6 ii 0 . 0 4 1 1 0 . 0 2 f 0 " 0 0 . 0 2 0 . 0 4 0 . 0 6 Y c o 2 (b) Figure 6.4: Ten realizations of S P M w i t h an in i t ia l temperature of 1200K compared to the con t inuum trajectory: (a) iV e = 1080; (b) Ne — 2700 (Thick line: c o n t i n u u m trajectory; T h i n lines: S P M realizations). 1 2 2 di f fe ren t r e a l i z a t i o n s i n the i g n i t i o n d e l a y t i m e o f h o m o g e n e o u s s y s t e m s ; these f l u c t u a t i o n s are d u e to the r a n d o m n a t u r e o f the m o l e c u l a r c o l l i s i o n s tha t l e a d to c h e m i c a l r e a c t i o n s , p a r t i c u l a r l y at e a r l y t i m e s i n the i n d u c t i o n time w h e r e c o l l i s i o n s tha t l e a d to a c h a n g e i n c h e m i c a l c o m p o s i t i o n are ra re even t s . It i s p o s s i b l e tha t the f l u c t u a t i o n s i n the i g n i t i o n d e l a y are d u e to a c o m b i n a t i o n o f the t w o p o s s i b l e sources , w h e r e t u r b u l e n t f l u c t u a t i o n s w o u l d l e a d to v a r i a t i o n s i n the t i m e the l o c a l s ca l a r g r a d i e n t s are su f f i -c i e n t l y l o w a n d c h e m i c a l f l u c t u a t i o n s w o u l d l e a d to v a r i a t i o n s i n h o w l o n g the g r a d i e n t s w o u l d h a v e to b e l o w for i g n i t i o n to b e s u c c e s s f u l . H e r e , the h y p o t h e s i s tha t the f l u c t u a t i o n s i n the i g n i t i o n d e l a y t i m e are p r i m a r i l y d u e to c h e m i c a l f l u c t u a t i o n s a n d tha t these c a n b e w e l l -r e p r e s e n t e d i n a s i m u l a t i o n b y u s i n g d i f fe ren t i n d i v i d u a l ( f i l tered) re-a l i z a t i o n s o f t he S P M i n the T G L D M w i l l be tes ted . T h e s a m e a p p r o a c h ( a n d c o d e ) as w a s u s e d b y H u a n g a n d B u s h e [60] i s u s e d he re . T h e R e y n o l d s - a v e r a g e d c o n s e r v a t i o n e q u a t i o n s for a f u l l y c o m p r e s s i b l e tu r -b u l e n t r e a c t i v e f l o w are s o l v e d o n a t w o - d i m e n s i o n a l , a x i s y m m e t r i c m e s h i n c y l i n d r i c a l c o o r d i n a t e s . A s t a n d a r d k — e t u r b u l e n c e m o d e l [137] i s u s e d . T h e g o v e r n i n g p a r t i a l - d i f f e r e n t i a l e q u a t i o n s are s o l v e d u s i n g a F l u x - C o r r e c t e d T r a n s p o r t ( F C T ) a l g o r i t h m [138] tha t i s g e n -e r a l l y f o u r t h - o r d e r accu ra t e i n space w i t h a f i n i t e - v o l u m e r ep re sen t a -tion. A s e c o n d - o r d e r R u n g e - K u t t a time a d v a n c e s c h e m e is u s e d for the t e m p o r a l d i s c r e t i z a t i o n . A d e n s i t y - b a s e d s c h e m e is u s e d to a c h i e v e the p r e s s u r e / v e l o c i t y c o u p l i n g . T r a n s p o r t e q u a t i o n s are s o l v e d for the m a s s f r ac t ions o f CO2 a n d H2O; c l o s u r e for the c h e m i c a l s o u r c e - t e r m s i n these t r a n s p o r t e q u a t i o n s i s o b t a i n e d w i t h C S E u s i n g the T G L D M 123 Table 6.1: Exper imenta l condit ions of shock tube. N o . Runs ^injector [mm] p. " Pb 1 o •L i T e 28 1.1 75 30 1.5 300 1150-1400 "injection pressure [bar] &Back pressure[bar] (pressure behind reflected shock) ^injection duration[ms] rfFuel temperature Temperature behind reflected shock [K] method. Each r u n of this R A N S code requires approximate ly one day on a single 2 . 4 G H z processor. The experiments of S u l l i v a n et al . [139] are s imulated. In these shock-tube experiments, a reflected shock-wave was used to compress air to relat ively h i g h temperatures and a pressure of 30 bar; methane was i n -jected short ly after the shock reflection into the essentially quiescent air w i t h a pressure ratio of 2.5 for 1.5 ms. The detection of ign i t ion delay is captured through optical observation instead of pre- igni t ion mea-surements [140,141] or sudden change of pressure rise [142]. Table 6.1 summarizes the experimental condit ions and m a i n parameters i n these experiments. A n in i t i a l air temperature of 1150K is chosen for these s imulat ions because there were a number of experiments at this temperature to compare to. Fi l tered S P M realizations are used to substitute o n l y the autoigni-tion trajectory i n the T G L D M . In pr inc ip le , filtered S P M trajectories cou ld be used for the who le mani fo ld , but the ign i t ion de lay t ime is a l -most exclus ively determined b y the mani fo ld i n the reg ion nearest the or ig in . Ten different T G L D M s were generated each us ing a different 124 filtered S P M real izat ion for the autoigni t ion trajectory. Each T G L D M table has values for 50 different discrete values of mixture fraction. Two add i t iona l T G L D M s were generated b y sort ing through the 10 T G L D M s and p i c k i n g out the filtered S P M trajectory w i t h the shortest and longest ign i t ion delay time at each mixture fraction and assem-b l i n g a combina t ion us ing a l l of the shortest delay times and a l l of the longest delay times. Ano the r s imula t ion was per formed us ing a T G -L D M w h i c h is constructed b y substi tuting the autoigni t ion trajectory calculated b y averaging S P M realizations i n the manner described ear-lier. Final ly , for compar ison, a s imula t ion has been per formed us ing the or ig ina l T G L D M calculated us ing the c o n t i n u u m solu t ion for a l l trajectories. F i g . 6.5 to F i g . 6.7 are the profiles of temperature, mass fraction of O H , mixture fraction and mixture fraction variance f ield at different times calculated w i t h the T G L D M S P M predic t ing the m a x i m u m i g n i -t ion delay times (realization 11 i n F ig . 6.9). F i g . 6.8 shows the profiles at about the same times as F i g . 6.7 but w i t h the T G L D M S P M pred ic t ing m i n i m u m ign i t ion delay time predic t ion (realization 12 i n F i g . 6.9). Several w a y s to decide the igni t ion delay t ime from the s imula t ion results have been proposed [59,60,143]. H u a n g found the ign i t ion de-lay time b y p lo t t ing the m a x i m u m mean of a scalar (temperature, O H or C 2 H 2 mass fraction) as a function of the time, then extrapolat ing the m a x i m u m slope of a scalar history back to the level at the onset of i n -jection [60]. In this work , ign i t ion is defined to be w h e n the m a x i m u m mean temperature i n any control vo lume exceeds 2000K. H u a n g [60] also calculated the ign i t i on delay t ime based o n the mass fractions of 125 F i g u r e 6.5: P r o f i l e o f the jet c a l c u l a t e d u s i n g T G L D M w i t h m a x i m u m i g n i t i o n d e l a y t i m e s be fo re i n j ec t i on s tops (at t = 1 . 30ms) 126 Figure 6 . 6 : Profile of the jet calculated us ing T G L D M w i t h m a x i m u m igni t ion delay times short ly after injection stops (at t = 1 . 5 6 m s ) 127 5 10 x [cm] s a Mixture fraction 5 10 x [cm] 0 .3 0 .2 0.1 0 Mixture fraction variance 5 10 x [cm] Figure 6.7: Profile of the jet calculated us ing T G L D M w i t h m a x i m u m igni t ion de lay times short ly after igni ted (at t = 1.73ms) 128 T[K] -3 x 10 Mixture fraction variance Figure 6.8: Profile of the jet calculated us ing T G L D M w i t h m i n i m u m igni t ion delay times short ly after igni ted (at t = 1.73ms) 129 2 . 4 2 . 2 1 . 8 1 . 6 r—i " I—r 1.4 II—u—U—u—U—U—U—u—u—U—U—U—I 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 S P M r e a l i z a t i o n s Figure 6.9: Var ia t ion i n ign i t ion delay times of a transient methane jet injected into air at an in i t i a l temperature of T = 1150K as calculated w i t h R A N S us ing the C S E / T G L D M method. (Dashed l ine: T G L D M us ing only c o n t i n u u m solutions; so l id l ine: averaged S P M autoigni t ion trajectory i n T G L D M ; realizations 1-10: T G L D M s us ing filtered S P M autoigni t ion trajectories; real izat ion 11: T G L D M us ing filtered S P M au-toigni t ion trajectories w i t h m a x i m u m igni t ion delay times; real izat ion 12: T G L D M us ing filtered S P M autoigni t ion trajectories w i t h m i n i m u m igni t ion delay times.) O H and C2H2. Howeve r , the trends i n the ign i t ion de lay time are s imi -lar to the one calculated us ing temperature and the difference is smal l . Therefore, on ly results us ing temperature as the ign i t ion cri ter ion are s h o w n here. There are indeed significant fluctuations i n the ign i t ion delay time predicted us ing the different T G L D M s evident i n F i g . 6.9. The results obtained us ing the filtered S P M trajectories a l l seem to agree w e l l w i t h the result u s ing the averaged S P M trajectory. The m a x i m u m ign i t ion delay time a m o n g the S P M results comes from us ing the T G L D M con-130 structed us ing the realizations w i t h the m a x i m u m ign i t ion delay times; s imilarly, the m i n i m u m comes from the T G L D M constructed w i t h the realizations w i t h the m i n i m u m delay. These two results suggest a strat-egy for est imating the range of ign i t ion delays to be expected: use the m i n i m u m a n d m a x i m u m ign i t ion delays f rom a large n u m b e r of S P M runs to construct a " m i n i m u m T G L D M " and a " m a x i m u m T G L D M " . It is interesting that the ign i t ion delay time r predic ted b y us ing the T G L D M w i t h the averaged S P M is smaller than that predic ted us-i n g the c o n t i n u u m T G L D M . One possible cause of this discrepancy is the smal l number of molecules of radical species at the very begin-n ing of ign i t ion . In the S P M calculat ion, there are discrete steps i n the concentrations of a l l species, since the S P M calculates the number of molecules i n the control vo lume and this must be a non-negative inte-ger. In the c o n t i n u u m calculation, the concentrations can take on any non-negative real value. Some reaction rates can be non-l inear func-tions of the species concentrations, w h i c h c o u l d l ead to a discrepancy between the average of the reaction rate. A further invest igat ion of different radicals i n the system reveals the number of molecules of O H calculated us ing the S P M ranges from 0 to 4, w i t h 0 and 1 be ing b y far the most frequent number found. The chemical source-term for C 0 2 is a h igh ly nonl inear function of mass fraction of O H as s h o w n i n F i g . 6.10. This source-term is tabulated i n the T G L D M ; the ign i t ion delay time calculated b y the R A N S code is clearly h i g h l y sensitive to this source-term because it is one of on ly two rates used to change the chemical composi t ion i n the R A N S code. A s can be seen i n F i g . 6.10, us ing the con t inuum approx imat ion leads 131 " 0 0 . 2 0 . 4 0 . 6 0 . 8 1 M a s s f r a c t i o n o f O H ^ n - e F i g u r e 6.10: T h e c h e m i c a l s o u r c e t e r m o f the m a s s f r a c t i o n o f CO2 as a f u n c t i o n o f the m a s s f r a c t i o n o f O H d u r i n g a u t o i g n i t i o n o f a s t o i c h i o -m e t r i c m i x t u r e o f m e t h a n e / a i r w i t h a n i n i t i a l t e m p e r a t u r e o f 1 1 5 0 K at 30 bar. ( S o l i d l i n e : c o n t i n u u m trajectory; d a s h e d l i n e : l i n e a r in t e r -p o l a t i o n b e t w e e n the O H c o n c e n t r a t i o n s e q u i v a l e n t to z e r o a n d o n e m o l e c u l e b e i n g p r e s e n t i n the S P M . ) 132 to a m u c h lower source-term for CO2 than one gets u s ing the averaged S P M results, w h i c h , early i n the induc t ion per iod , w o u l d be the av-erage between the rates calculated w i t h zero and one molecule i n the S P M vo lume . Th i s appears to exp la in w h y the c o n t i n u u m a n d S P M results differ. It does not offer any guidance, however, as to w h i c h of the two re-sults is the correct one. In pr inc ip le , the result obtained w i t h the S P M shou ld be more accurate, since it is somewhat more phys ica l ly accu-rate. H o w e v e r , the S P M implementa t ion used here makes the some-what unrealistic assumpt ion that the one, tiny control v o l u m e is adia-batic. A more realistic S P M implementa t ion w o u l d calculate for m a n y control vo lumes and a l l ow for diffusion amongst these. The average of these control vo lumes might then be used i n place of the averaged S P M . Unfortunately, a d d i n g diffusion to the S P M appears to increase the computa t ional cost dramat ical ly (clearly, the computa t ional cost cannot scale on ly w i t h the number of control volumes) . The quest ion of w h i c h of the two predict ions of the igni t ion delay time is the better one remains unresolved and is left for future research. F ig . 6.11 compares the s imula t ion results to the exper imental data. In the experiments, it is impossible to precisely control the condit ions after the shock-reflection, so that there are some m i n o r variat ions i n the temperature between different runs; the error i n the est imation of the post-shock-reflection air temperature is approximate ly 15K. It shou ld be stressed that the error i n the experimental ign i t ion delay is a round 50 /is: the exper imental ly measured fluctuations i n the ign i t ion delay are real. There is an u n d e r l y i n g functional dependence of the ign i t ion 133 1140 1145 1150 1155 T [ K ] 1160 1165 Figure 6.11: C o m p a r i s o n of the ign i t ion delay t ime of methane jets is-su ing into quiescent air calculated w i t h filtered S P M generated trajec-tories w i t h experiment data for an in i t ia l air temperature T w 1150K (•: experimental results of S u l l i v a n et al . [139]; x: R A N S results us ing T G L D M s w i t h filtered S P M autoigni t ion trajectories; *: R A N S result us ing T G L D M w i t h averaged S P M autoigni t ion trajectory; o: R A N S result us ing T G L D M based on con t inuum approximat ion) . delay on temperature, w i t h the igni t ion delay d r o p p i n g w i t h increas-i n g temperature. However , it is evident from the figure that the fluctu-ations i n the ign i t ion delay for this system at this temperature range are about 700 /is. Inc lud ing chemical fluctuations i n the manner adopted here does lead to significant fluctuations i n the ign i t ion de lay time. However , the magni tude of the fluctuations is lower than seen i n the experiments by about a factor of two or three. It w o u l d appear that the chemical fluctuations do have a significant impact on the ign i t ion delay of the turbulent jet, but that to predict the fu l l fluctuations accu-rately w o u l d require the abi l i ty to include turbulent fluctuations. Large 134 r1 E d d y S i m u l a t i o n ( L E S ) w o u l d b e o n e s t r a t egy for d o i n g t h i s , a l t h o u g h L E S i s c o n s i d e r a b l y m o r e c o m p u t a t i o n a l l y e x p e n s i v e t h a n the R A N S t e c h n i q u e u s e d h e r e s i n c e i t r e q u i r e s a t h r e e - d i m e n s i o n a l g r i d . I n a n L E S con tex t , i t w o u l d l i k e l y b e p r e f e r a b l e to i n c l u d e s e v e r a l d i f f e ren t f i l t e r ed r e a l i z a t i o n s o f a u t o i g n i t i o n i n the o n e T G L D M , r a t h e r t h a n i n -c l u d i n g o n l y the o n e r e a l i z a t i o n i n the T G L D M , as w a s d o n e he re . O n e c o u l d t h e n c h o o s e a m o n g these d i f fe ren t t ra jector ies a t r a n d o m , w h i c h w o u l d b e a m o r e r ea l i s t i c r e p r e s e n t a t i o n o f the c h e m i c a l f l u c t u a t i o n s t h a n w a s d o n e i n the s i m p l e R A N S s i m u l a t i o n s p e r f o r m e d here . 6.3 C o n c l u s i o n s It has b e e n f o u n d that d i f fe ren t r e a l i z a t i o n s o f the a u t o i g n i t i o n t ra-j e c to ry f o l l o w a d i f fe ren t p a t h t h r o u g h c o m p o s i t i o n space . A f t e r f i l -t e r i n g o u t the s m a l l t ime- sca l e s i n the t ra jector ies , the p a t h s d o u l t i -m a t e l y c o n v e r g e o n t o a s i n g l e m a n i f o l d be fo re t h e y r e a c h the e q u i -l i b r i u m c o m p o s i t i o n , h o w e v e r t h e y take v e r y d i f fe ren t p a t h s t h r o u g h c o m p o s i t i o n space at l o w e r t e m p e r a t u r e s . B y f i l t e r i n g the t ra jec tory i n e v e n t space , i t i s p o s s i b l e to i n c o r p o r a t e a s i n g l e r e a l i z a t i o n f r o m the S P M i n t o a T G L D M . W h e n T G L D M s i n c o r p o r a t i n g d i f f e r en t r e a l i z a -t i ons o f the a u t o i g n i t i o n t ra jec tory are u s e d i n c a l c u l a t i o n s o f t u r b u l e n t a u t o i g n i t i o n , the r e su l t s o f the d i f fe ren t c a l c u l a t i o n s g i v e d i f f e r en t i g n i -t i o n d e l a y t i m e s . H o w e v e r , the T G L D M s i n c o r p o r a t i n g d i f f e ren t r e a l -i z a t i o n s o f the a u t o i g n i t i o n trajectory, w h e n u s e d i n the R A N S c o d e , d o n o t l e a d to f l u c t u a t i o n s o f the s a m e o r d e r as are s e e n i n e x p e r i m e n t s . T h i s i n d i c a t e s tha t t u r b u l e n t f l u c t u a t i o n s l i k e l y a l s o h a v e a n i m p o r t a n t effect; to a c c o u n t for the c o m b i n e d effects o f t u r b u l e n t a n d c h e m i c a l 135 fluctuations, the use of Direct N u m e r i c a l S imula t ion (DNS) or Large E d d y S imula t ion (LES), w h i c h is computat ional ly cheaper and more plausible, is indicated. 136 Chapter 7 Conclusions and Future Work 7 . 1 C o n c l u s i o n s This thesis has explored combust ion m o d e l i n g w i t h C o n d i t i o n a l Source-term Es t imat ion us ing both L amina r Flamelet Decompos i t i on and L o w -dimens iona l M an i fo ld s . First, C S E was used w i t h L F D i n L E S to obtain three different predict ions of the Sandia D-flame. It was found that both steady and unsteady libraries can predict the temperature and most species w e l l . A l ibrary comprised of both steady a n d unsteady flamelets was needed to get a good predic t ion of N O . L F D was s h o w n to have several disadvantages. The most impor tant of these is that it appears to be possible to influence the predict ions by al tering the compos i t ion of the flamelet library. A n impor tant goal i n L E S m o d e l development shou ld be to avo id the need (or even the ability) to tune the m o d e l to fit experimental results. In this respect, C S E w i t h L F D appears to have a significant short-coming. This (in part) mot ivated the development of C S E - T G L D M for L E S w h i c h was also app l i ed to simulate the Sandia D-f lame. It was found that C S E - T G L D M can be used w i t h L E S i n non-premixed tur-bulent combust ion to p rov ide good predict ions of temperature and of a l l species for w h i c h transport equations are so lved . B y interpo-137 la t ing for m i n o r species condi t ional mass fractions direct ly from the T G L D M rather than so lv ing transport equations, significant c o m p u -tational savings were real ized. W h e n these condi t ional averages are m a p p e d back into the phys ica l space, the agreement w i t h exper imental results is ve ry good. G o o d predictions of N O require so lv ing a trans-port equation for that species mass fraction and obta in ing the condi -t ional source-term from the T G L D M table. G R I - M e c h 2.11 predicts N O m u c h better G R I - M e c h 3.0 because of deficiencies that have been iden-tified i n the N O x chemistry i n G R I - M e c h 3.0. In order to address a shortcoming i n the T G L D M method at l o w temperatures, the Stochastic Particle M o d e l was used to generate tra-jectories. W h e n m a n y realizations of the autoigni t ion trajectory pre-dicted b y the S P M were averaged together, the average trajectory takes the same path as the trajectory that results from so lv ing the c o n t i n u u m equations. Thus , the S P M provides an alternate me thod to solve the con t inuum equations to obtain the autoigni t ion trajectory that is not affected by stiffness i n the chemical kinetic mechanism. It was s h o w n that the method c o m b i n i n g unburned and e q u i l i b r i u m composi t ions to a v o i d stiffness i n the autoigni t ion trajectory at extremely l o w i n i -tial temperatures is clearly deficient and that the averaged S P M trajec-tory provides a viable alternative to this approach. The averaged S P M trajectory, w h e n inserted into the T G L D M , leads to an excellent predic-t ion of a P S R w h e n compared to the P S R calculated w i t h fu l l chemistry, whereas the pred ic t ion us ing the T G L D M constructed u s ing the met-h o d combin ing unburned and equ i l i b r ium composi t ions in t roduced a significant error i n the P S R calculations. 138 It was also found that different realizations of the autoigni t ion tra-jectory fo l low a different path through compos i t ion space. After f i l -tering out the smal l time-scales i n the trajectories, the paths ul t imately converge onto a single man i fo ld w e l l before they reach the e q u i l i b r i u m composi t ion. H o w e v e r , they take very different paths through compo-si t ion space at l ower temperatures. B y fil tering the trajectory i n event space, it was s h o w n to be possible to incorporate a single real izat ion from the S P M into a T G L D M . W h e n T G L D M s incorpora t ing different realizations of the autoigni t ion trajectory were used i n calculations of turbulent autoigni t ion, the results of the different calculations gave dif-ferent ign i t ion delay times. However , w h e n they were used i n a R A N S code, the T G L D M s incorporat ing different realizations of the autoigni-tion trajectory d i d not lead to fluctuations of the same order as seen i n experiments. This indicates that turbulent fluctuations l i ke ly also have an important effect; to account for the combined effects of turbulent and chemical fluctuations, use of Large E d d y S imula t ion is indicated. 7.2 Future Work M o r e cases s h o u l d be r u n to test the homogenei ty assumpt ion i n C S E , such as a bluff b o d y flame. A more complete flamelet l ib ra ry for L F D can be bui l t , e.g. to inc lude re igni t ion or edge flame phenomenon to predict more severe cases such as the Sandia E-flame where extinction is more significant. C S E w i t h T G L D M might be more p romis ing because the T G L D M is more flame independent . However , it is important to i m p r o v e the per-formance of manifolds , especially i n regions where the phys ica l t ime 139 scales are comparable to the time scales i n the fast subspace. H i g h e r d imens ion manifolds m a y be necessary for such a case, but the storage and retrieval cost w o u l d have to be addressed. The S P M i n combina t ion w i t h T G L D M can predict fluctuations i n the autoigni t ion de lay time. However , u s ing R A N S , the magni tude was smaller than is found i n the experimental data. Future s tudy is needed to use a more advanced C F D method, such as D N S or L E S , to predict the effects of turbulent fluctuations w h i c h appear to be a major factor causing the variations observed i n the experiments. 140 Bibliography [1] J. H u a n g , N a t u r a l gas combust ion under engine-relevant condi -tions, P h . D . thesis, Un ive r s i ty of Br i t i sh C o l u m b i a (2006). [2] The A n n u a l Energy Out look 2006 w h i c h is available at: h t t p : / / www.e i a .doe .gov /o i a f / aeo / index .h tml . [3] D . C . Rapaport , The art of molecular dynamics s imula t ion , C a m -br idge U n i v e r s i t y Press, 1997. [4] P. A . Libby , F. A . W i l l i a m s , Turbulent reacting f lows, A c a d e m i c press L o n d o n , 1994. [5] N . Peters, Turbulent combust ion, 1st E d i t i o n , C a m b r i d g e U n i -versi ty Press, 2000. [6] R. Borgh i , Turbulent combust ion mode l l i ng , Progress i n Energy C o m b u s t i o n a n d Science 14 (1988) 245-292. [7] D . Len t in i , Assessment of the stretched laminar flamelet ap-proach for nonpremixed turbulent combust ion, C o m b u s t i o n Sci-ence a n d Technology 100 (1994) 95-122. [8] J. Warnatz , U . Maas , R. W. Dibb le , Combus t i on , phys ica l and chemical fundamentals, mode l i ng and s imula t ion , experiments, pol lutant formation, 3rd Ed i t ion , Springer, 2001. [9] P. M o i n , K . Mahesh , Direct numer ica l s imula t ion: a tool i n tur-bulence research, A n n u a l R e v i e w of F l u i d Mechan ics 30 (1998) 539-578. [10] V. Eswaran , S. B . Pope, A n examinat ion of forcing i n direct numer ica l s imulat ions of turbulence, C o m p u t a t i o n a l F l u i d s 16 (1988)257-278. [11] R. M . J. K i m , P. M o i n , Turbulence statistics i n fu l ly deve loped channel flow, Journal of F l u i d Mechanics 177 (1987) 133-166. [12] P. R. Spalart, Direct numer ica l s tudy of leading-edge contami-nation, F l u i d D y n a m i c s of Three-Dimensional Turbulent Shear F lows and Transit ion. 141 [13] D . C . Wi l cox , Turbulence mode l l i ng for C F D , D C W Industries, L a Canada , 1998. [14] S. H . Johansson, L . D a v i d s o n , E . Olsson, N u m e r i c a l s imula t ion of vortex shedd ing past triangular cyl inders at h i g h reynolds n u m -ber us ing a k-e turbulence mode l , International Journal for N u -merical Me thods i n F lu ids 16 (10) (1993) 859-878. [15] J. Smagorinsky, Genera l circulat ion experiments w i t h the p r i m i -tive equations, M o n t h l y Weather R e v i e w 91 (3) (1963) 99-164. [16] M . Germano , U . P i o m e l l i , P. M o i n , W. H . Cabot , A dynamic subgrid-scale e d d y viscosi ty mode l , Physics of F lu ids 3 (1991) 1760-1765. [17] P. M o i n , K . Sqires, W . Cabot , S. Lee, A d y n a m i c subgrid-scale m o d e l for compressible turbulence and scalar transport, Physics of F lu ids 3 (1991) 2746-2757. [18] A . E . Tejada-Martiz, K . E . Jansen, A d y n a m i c smagor insky m o d e l w i t h d y n a m i c determinat ion of the filter w i d t h ratio, Physics of F l u i d s 16 (2004) 2514-2528. [19] S. B. Pope, Ten questions concerning the large-eddy s imula t ion of turbulent f lows, N e w Journal of Physics 6 (2004) 35-59. [20] H . Pi tsch, Improved pollutant predict ions i n large-eddy s imu-lations of turbulent non-premixed combust ion b y cons ider ing scalar d iss ipat ion rate fluctuations, Proceedings of the C o m b u s -t ion Institute 29 (2002) 1971-1978. [21] H . Pitsch, Large-eddy s imula t ion of turbulent combust ion, A n n u . Rev. F l u i d M e c h . 38 (2006) 453-482. [22] D . R. C h a p m a n , Computa t iona l aerodynamics development and outlook, A I A A J. 17 (1979) 1293-1294. [23] P. R. Spalart, W. H . Jou, M . Stretlets, S. R. A l l m a r a s , Commen t s on the feasibili ty of L E S for wings and on the h y b r i d R A N S / L E S approach, Proceedings of the First A F O S R International Confer-ence on D N S / L E S . 142 [24] P. R. Spalart, Strategies for turbulence m o d e l i n g and s imulat ions, International Journal of Heat and F l u i d F l o w 21 (2000) 252-263. [25] A . R. Kers te in , A l inear-eddy mode l of turbulent scalar transport and m i x i n g , C o m b u s t i o n Science and Technology 60 (1988) 391. [26] P. A . M c M u r t r y , S. M e n o n , A . R. Kers te in , L inea r e d d y sub-gr id m o d e l for turbulent reacting flows: appl ica t ion to hydrogen-air combust ion, Proceedings of the C o m b u s t i o n Institute (1992) 271 -278. [27] R. W. Bilger, Turbulent flows w i t h non-premixed reactants, Top-ics A p p l i e d Physics 44 (1980) 65-113. [28] S. Burke , T. Schumann, Diffusion flames, Industr ia l and E n g i -neer ing Chemis t ry 20 (1928) 998-1004. [29] J. C o u p l a n d , C . H . P r i d d i n , M o d e l l i n g of f low and combus t ion i n a p roduc t ion gas turbine combustor, Turbulent Shear F l o w s 5 (1987) 310-323. [30] N . Peters, L a m i n a r diffusion flamelet models i n non-premixed turbulent combust ion, Progress i n Energy C o m b u s t i o n Science 10 (3)(1984)319-339. [31] H . Pi tsch , N . Peters, A consistent flamelet fo rmula t ion for non-premixed combust ion consider ing differential diffusion effects, C o m b u s t i o n a n d Flame 114 (1998) 26-40. [32] H . Pi tsch, Uns teady flamelet mode l i ng of differential diffusion i n turbulent jet diffusion flames, C o m b u s t i o n a n d F lame 123 (3) (2000)358-374. [33] H . Pi tsch, H . Steiner, Large-eddy s imula t ion of a turbulent p i -lo ted methane /a i r diffusion flame (Sandia flame D ) , Physics of F lu ids 12 (2000) 2541-2554. [34] N . Peters, C o m m e n t and reply on the "Assessment of combus-t ion and submodels for turbulent nonpremixed hydrocarbon flames", C o m b u s t i o n and F lame 116 (1999) 675-676. 143 [35] N . Swamina than , R. Bilger, Assessment of combus t ion and sub-models for turbulent nonpremixed hydrocarbon flames, C o m -bus t ion a n d F lame 116 (1999) 519-545. [36] N . Swamina than , Flamelet regime i n non-premixed combust ion, C o m b u s t i o n and Flame 129 (2002) 217-219. [37] F. M a u s s , D . Kel ler , N . Peters, A lagrangian s imula t ion of flamelet extinction and re-ignit ion i n turbulent jet diffusion flames, Proceedings of the Combus t i on Institute 23 (1990) 693-698. [38] H . Pi tsch, M . C h e n , N . Peters, Uns teady flamelet m o d e l i n g of turbulent hydrogen-air diffusion flames, Proceedings of the C o m b u s t i o n Institute 27 (1998) 1057-1064. [39] H . Barths, C . Hasse, S imula t ion of combust ion i n direct injection diesel engines us ing an eulerian particle flamelet m o d e l , Pro-ceedings of the Combus t i on Institute 28 (2000) 1161-1168. [40] P. Coe lho , N . Peters, N u m e r i c a l s imula t ion of a m i l d combus t ion burner, C o m b u s t i o n and Flame 124 (2001) 503-518. [41] R. W . Bilger, C o n d i t i o n a l moment closure for turbulent reacting f lows, Physics of F lu ids 5 (1993) 436-444. [42] A . Y. K l i m e n k o , Mul t i componen t diffusion of var ious admix-tures i n turbulent flow, F l u i d Dynamics 25 (3) (1990) 327-334. [43] A . K l i m e n k o , R. Bilger, Cond i t i ona l moment closure for tur-bulent combust ion, Progress i n energy and combus t ion science 25 (6) (1999) 595-687. [44] S. K i m , K . Y. H u h , Use of the condi t ional moment closure m o d e l to predict no formation i n a turbulent C H 4 / H 2 flame over a bluff body, C o m b u s t i o n and Flame 130 (2002) 94-111. [45] S. H . K i m , K . Y. H u h , R. A . Fraser, M o d e l i n g autoigni t ion of a turbulent methane jet b y the condi t ional moment closure mode l , Proceedings of the Combus t i on Institute 28 (2000) 185-191. 144 [46] C . D e v a u d , K . N . C . Bray, Assessment of the appl icabi l i ty of con-di t iona l moment closure m o d e l to a lifted turbulent flame: first order m o d e l , C o m b u s t i o n and Flame 132 (2003) 102-114. [47] S. H . K i m , K . Y. H u h , R. W. Bilger, Second-order condi t ional moment closure m o d e l i n g of local ext inct ion and re igni t ion i n turbulent non-premixed hydrocarbon flames, Proceedings of the C o m b u s t i o n Institute 29 (2002) 2131-2137. [48] M . R o o m i n a , R. Bilger, Cond i t i ona l moment closure m o d e l i n g of turbulent methanol jet flames, Combus t i on Theory and M o d e l -i n g 3 (1999) 698-708. [49] M . R o o m i n a , R. W . Bilger, Cond i t iona l moment closure ( C M C ) predic t ion of a turbulent methane-air flames, C o m b u s t i o n and flame 125 (2001) 1176-1195. [50] M . Fairweather, R. M . Wooley, First-order condi t iona l moment closure m o d e l i n g of turbulent, non-premixed methane flames, C o m b u s t i o n and Flame 138 (2004) 3-19. [51] S. K i m , K . Y. H u h , L . Tao, A p p l i c a t i o n of the el l ipt ic condi t iona l moment closure m o d e l to a two-dimens iona l nonpremixed methanol bluff-body flame, Combus t i on and F lame 75-90. [52] N . Swamina than , R. W. Bilger, Cond i t i ona l variance equation and its analysis, Proceedings of the C o m b u s t i o n Institute (1997) 1191-1198. [53] E . Mastorakos , R. W. Bilger, Second-order condi t iona l moment closure for the autoigni t ion of turbulent f lows, Physics of F lu ids 10 (1998)1246-1248. [54] C . M . C h a , G . Kosaly , H . Pi tsch, M o d e l i n g ext inct ion and reig-n i t ion i n turbulent nonpremixed combust ion u s ing a doub ly -condi t ional moment closure approach, Physics of F l u i d s 13 (12) (2001)3824-3834. [55] W. K . Bushe, H . Steiner, Cond i t iona l moment closure for large eddy s imula t ion of nonpremixed turbulent reacting flows, Physics of F l u i d s 11 (1999) 1896-1906. 145 [56] H . Steiner, W. K . Bushe, Large eddy s imula t ion of a turbulent reacting jet w i t h condi t ional source-term estimation, Physics of F lu ids 13 (2001) 754-769. [57] H . Pi tsch, Improved pol lutant predict ions i n large-eddy s imu-lations of turbulent non-premixed combust ion b y consider ing scalar d iss ipat ion rate fluctuations, Proceedings of the C o m b u s -t ion Institute 29 (2003) 1971-1978. [58] W. K . Bushe, H . Steiner, Lamina r flamelet decompos i t ion for con-di t ional source-term estimation, Physics of F lu ids 15 (2003) 1564-1575. [59] R. Grout , C o m b u s t i o n mode l l i ng w i t h condi t ional source-term estimation and laminar flamelet decomposi t ion, Mas te r ' s thesis, Un ive r s i ty of Br i t i sh C o l u m b i a (2004). [60] J. H u a n g , W . K . Bushe, S imula t ion of transient turbulent met-hane jet ign i t ion and combust ion under engine-relevant condi -tions us ing condi t ional source-term estimation w i t h detai led ch-emistry, Submit ted to Combus t i on Theory and M o d e l i n g . [61] R. Grout , W. K . Bushe, C . Blair , Lamina r flamelet decomposi t ion i n autoigni t ion predic t ion, Submit ted to C o m b u s t i o n Theory and M o d e l l i n g . [62] G . P. Smi th , D . M . Go lden , M . Frenklach, N . W . M o r i -arty, B . Eiteneer, M . Goldenberg , C . T. B o w m a n , R. K . H a n s o n , S. Song, W. C . Gardiner, V. V. L i s s i ansk i , Z . Q . h t t p : / / w w w . m e . b e r k e l e y . e d u / g r L m e c h / . [63] M . Z . Bodenstein, E ine theorie der photochemischen reaktions-geschwindigkei ten , Physics of C h e m 85 (1913) 329-397. [64] N . Peters, R. Kee , The computa t ion of stretched laminar methane-air diffusion flames us ing a reduced four-step mech-anism, C o m b u s t i o n and Flame 68 (1) (1987) 17-30. [65] M . D . Smooke, Reduced kinetic mechanisms and asymptot ic ap-proximat ions for methane-air flames, Springer-Verlag, 1991. [66] N . Peters, B. Rogg , Reduced kinetic mechanisms for applicat ions i n combust ion systems, V o l . 15, Springer-Verlag, 1991. 146 [67] J. Warnatz , Concentrat ion, pressure, and temperature-dependence of the flame veloci ty i n hydrogen-oxygen-ni t rogen mixtures, C o m b u s t i o n Science and Technology 26 (5-6) (1981) 203-213. [68] U . N o w a k , J. Warnatz , Sensit ivity analysis i n al iphat ic hydrocar-b o n combust ion, Progress i n Astronaut ics and Aeronaut ics 113 (1988) 87-103. [69] M . Nehse , J. Warnatz , C . Chevalier , Kine t ic m o d e l i n g of the ox i -dat ion of large al iphatic hydrocarbons, Proceedings of the C o m -bust ion Institute (1996) 432-459. [70] S. H . L a m , D . A . Goussis , Unders tand ing complex chemical k i -netics w i t h computat ional s ingular perturbation, Proceedings of the C o m b u s t i o n Institute (1988) 931-941. [71] S. H . L a m , D . A . Goussis , The C S P method for s i m p l i f y i n g kinet-ics, International journal of chemical kinetics 26 (1994) 461-486. [72] Z . Ren , S. B . Pope, A . V l a d i m i r s k y , J. M . Guckenheimer , A p p l i c a -t ion of the I C E - P I C method for the d imens ion reduct ion of chem-ical kinetics, Journal of C h e m . Physics of 124 (2006) 114111. [73] Z . Ren , S. B. Pope, A . V lad imi r sky , J. M . Guckenheimer , A p p l i c a -t ion of the I C E - P I C method for the d imens ion reduct ion of chem-ical kinetics coupled w i t h transport, Proceedings of the C o m b u s -t ion Institute. [74] A . Massiasa , D . Diamant isa , E . Mastorakos, D . A . Goussisa , A n a lgor i thm for the construction of global reduced mechanisms w i t h C S P data, Combus t i on and Flame 117 (1999) 685-708. [75] T. L u , Y. Ju, C . K . L a w , C o m p l e x C S P for chemistry reduct ion a n d analysis, C o m b u s t i o n and Flame 126 (2001) 1445-1455. [76] A . Kazakov , M . Chaos , Z . Zhao , F. L . Dryer , C o m p u t a t i o n a l s in-gular per turbat ion analysis of two-stage ign i t ion of large hydro -carbons, Journal of Physics of Chemis t ry 110 (21) (2006) 7003 -7009. 147 [77] A . S. Toml in , T. Turanyi , M . J. P i l l i n g , Low-temperature combus-t ion and autoigni t ion, Mathemat ica l tools for the construction i n -vestigation and reduct ion of combust ion mechanisms 57 (1997) 1531-1556. [78] U . Maas , S. B . Pope, S impl i fy ing chemical kinetics: Intrinsic l o w -d imens iona l manifolds i n composi t ion space, C o m b u s t i o n and Flame 88 (1992) 239-246. [79] U . Maas , S. B . Pope, Implementat ion of s impl i f i ed chem-ical kinetics based on intrinsic low-d imens iona l manifolds , Int. S y m p . on combust ion 24 (1992) 103-112. [80] O . G icque l , N . Darabiha , D . Thevenin , L a m i n a r p remixed hydro-g e n / air counterf low flame simulat ions us ing flame prolongat ion of I L D M w i t h differential diffusion, Proceedings of the C o m b u s -tion F lame 28 (2000) 1901-1908. [81] J. A . v a n Oijen, L . P. H . de Goey, M o d e l i n g of p remixed laminar flame us ing flame-generated low-d imens iona l manifolds , Seven-teenth International C o l l o q u i u m on the D y n a m i c s of Explos ions and Reactive Systems (1999) 25-30. [82] H . Bongers, J. A . V. Oijen, L . P. H . de Goey, Intrinsic l o w -dimens iona l man i fo ld method extended w i t h diffusion, Pro-ceedings of the C o m b u s t i o n Institute 29 (2002) 1371-1378. [83] S. Pope, Computa t iona l ly efficient implementa t ion of com-bust ion chemistry us ing i n situ adaptive tabulat ion, C o m b u s -tion Theory M o d e l l i n g 1 (1997) 41-63. [84] S. S ingh , J. Powers , S. Paolucci , O n s low manifolds of chemical ly reactive systems, Journal of C h e m . Physics of 117 (2002) 1482-1496. [85] M . D a v i s , R. Skodje, Geometr ic invest igat ion of l ow-d imens iona l manifolds i n systems approaching equ i l i b r ium, Journal of C h e m . Physics of 111 (1999) 859-874. [86] S. B . Pope, U . Maas , S impl i fy ing chemical kinetics: Trajec-tory generated low-d imens iona l manifolds , Mechan ica l and Aerospace Engineer ing Report: F D A 93-11, C o r n e l l Univers i ty . 148 [87] J. P. Igniz io , T. M . Caval ier , L inear p rog ramming , Prentice H a l l International Series i n Industr ial and Systems Engineer ing , Pren-tice H a l l , E n g l e w o o d Cliffs , N J , 1994. [88] J. H u a n g , P. G . H i l l , W. K . Bushe, S. R. M u n s h i , Shock-tube s tudy of methane ign i t ion under engine-relevant condit ions: experi-ments a n d mode l ing , Combus t i on and Flame 136 (2004) 25-42. [89] R. Renka , Tripack: A constrained two-d imens iona l de launay tr i-angulat ion package, A C M Transactions on Mathemat ica l soft-ware 22 (1996) 1-8. [90] J. C . Keck , D . Gi l l esp ie , Rate-controlled pa r t i a l -equ i l ib r ium me-thod for treating reacting gas mixtures, C o m b u s t i o n a n d F lame 17(1971)237-241. [91] J. C . Keck , Rate-controlled const ra ined-equi l ibr ium theory of chemical reactions i n complex systems, Prog . Energy Combus t . 16(1990)125-154. [92] Q . Tang, S. B . Pope, Implementat ion of combust ion chemistry b y i n s i tu adapt ive tabulation of rate-controlled constrained equi-l i b r i u m manifolds , Proceedings of the C o m b u s t i o n Institute 29 (2002)1411-1417. [93] Z . Ren , S. B. Pope, Species reconstruction u s ing pre-image curves, Proceedings of the Combus t i on Institute 30 (2005) 1293-1300. [94] C . D o p a z o , E . O ' B r i e n , A n approach to the descr ip t ion of a tur-bulent mixture , A c t a A s t r o n 1 (1974) 1239. [95] S. B. Pope, P D F methods for turbulent reactive f lows, Progress i n Energy C o m b u s t i o n Science 11 (1985) 119-192. [96] J. Y. C h e n , W . K o l l m a n n , R. W . Dibb le , P D F m o d e l i n g of turbu-lent nonpremixed methane jet flames, C o m b u s t i o n Science and Technology 64 (1989) 315-346. [97] S. B . Pope, Computa t ions of turbulent combust ion: Progress and challenges, Proceedings of the C o m b u s t i o n Institute 23 (1990) 591-612. 149 [98] C . D o p a z o , Probabi l i ty density function approach for heated tur-bulent axisymmetr ic jet centerline evolut ion , Physics of F lu ids 18 (1975)397-404. [99] R. L . C u r l , D i spe r sed phase mix ing : 1. theory a n d effects i n s im-ple reactors, A l C H E 9 (1963) 175-181. [100] J. Janicka, W . K o l b e , W . K o l l m a n n , Closure of the transport equa-t ion for the probabi l i ty density function of turbulent scalar fields, Journal of N o n - E q u i l i b r i u m Thermodynamics 4 (1979) 47-66. [101] H . C h e n , S. C h e n , R. H . Kra ichnan , Probabi l i ty d is t r ibut ion of a stochastically advected scalar f ield, Physics R e v i e w Lett. 63 (1989)2657-2660. [102] S. Subramaniam, S. B. Pope, A m i x i n g m o d e l for turbulent re-active f lows based on euc l id ian m i n i m u m spann ing trees, C o m -bust ion and Flame 115 (1998) 487-514. [103] R. Cao , S. Pope, N u m e r i c a l integration of stochastic differen-tial equations: weak second-order mid-po in t scheme for appl ica-t ion i n the compos i t ion P D F method, Journal of C o m p u t a t i o n a l Physics 185 (2003) 194-212. [104] J. C h e n , W . C h a n g , M o d e l i n g differential diffusion effects i n tur-bulent nonreact ing/react ing jets w i t h stochastic m i x i n g models , C o m b u s t i o n Science and Technology 133 (1998) 343-375. [105] A . M a s r i , R. Cao , S. Pope, G . G o l d i n , P D F calculations of turbu-lent l ifted flames of h 2 / n 2 fuel i s su ing into a v i t ia ted co-flow, C o m b u s t i o n Theory M o d e l l i n g 8 (2004) 1-22. [106] V. Raman , R. Fox, A . Harvey, H y b r i d f in i te -volume/ t ranspor ted P D F s imulat ions of a par t ia l ly p remixed methane-air flame, C o m b u s t i o n and Flame 136 (2004) 327-350. [107] R. Cao , S. B . Pope , A . R. M a s r i , Turbulent l if ted flames i n a v i t i -ated coflow investigated us ing joint P D F calculations, C o m b u s -tion and F lame 142 (2005) 438-453. [108] K . L i u , S. Pope, D . Caughey, Calculat ions of bluff-body s tabi l ized flames us ing a joint P D F m o d e l w i t h detai led chemistry, C o m -bus t ion and Flame 142 (2005) 89-117. 150 [109] P. V. Slooten, Jayesh, S. Pope, Advances i n P D F m o d e l i n g for i n -homogeneous turbulent f lows, Physics of F l u i d s 10 (1998) 246-265. [110] D . T. Gi l l esp ie , A general method for numer ica l ly s imula t ing the stochastic t ime evolu t ion of coupled chemical reactions, Journal of Compu ta t iona l Physics 22 (1976) 403-434. [ I l l ] D . L . Bunker , B . Garrett, T. Kle indiens t , G . S. L . I l l , Discrete s im-ula t ion methods i n combust ion kinetics, C o m b u s t i o n and Flame 23 (1974) 373-379. [112] D . T. Gi l l esp ie , Exact stochastic s imula t ion of coup led chemical reactions, Journal of Physics of Chemis t ry 81 (1977) 2340-2361. [113] D . T. Gi l l esp ie , A rigorous der ivat ion of the chemical master equation, Phys ica A 188 (1992) 404-425. [114] M . Kraft , M . Wanger, N u m e r i c a l s tudy of a stochastic particle method for homogeneous gas-phase reactions, Compute r s and M a t h , w i t h App l i ca t i ons 45 (2003) 329-349. [115] A . Frisque, J. Schnakenberg, J. H u a n g , W. K . Bushe, Stochastic s imula t ion of the variations i n the autoigni t ion de lay t ime of pre-m i x e d methane and air, Combus t ion Theory and M o d e l i n g . [116] M . A . G i b s o n , J. Bruck, Efficient exact stochastic s imula t ion of chemical systems w i t h m a n y species and m a n y channels, Journal of C h e m i c a l Physics 104 (2000) 1876-1889. [117] D . T. Gi l l esp ie , A p p r o x i m a t e accelerated stochastic s imula t ion of chemical ly reacting systems, Journal of C h e m i c a l Physics 115 (2001) 1716-1733. [118] D . T. Gi l l esp ie , L . R. Pe tzold , Improved leap-size selection for ac-celerated stochastic s imula t ion , Journal of C h e m i c a l Physics 119 (2004) 8229-8234. [119] M . Schwehm, Para l le l stochastic s imula t ion of whole-ce l l m o d -els, Proceedings of 2nd Int. C o n f Systems B io logy (ICSB2001) (1996)333-341. 151 [120] E . L . Hasel t ine , J. B . Rawl ings , A p p r o x i m a t e s imula t ion of cou-p l e d fast and s low reactions for stochastic chemical kinetics, Journal of C h e m i c a l Physics 117 (2002) 6959-6969. [121] D . D a v i d , J.-P. Boon , Y . - X . L i , Lattice-gas automata for coupled reaction diffusion equations, Physics of Rev. Lett. 66 (19) (1991) 2535-2538. [122] H . P. Breuer, W. Huber , F. Petruccione, F luctuat ion effects on wave propagat ion i n a reaction-diffusion process, Phys ica D 73 (1994) 259-273. [123] H . P. Breuer, J. H o n e r k a m p , E Petruccione, A stochastic approach to complex chemical reactions, C h e m . Physics of Lett. 190 (3) (1992)199-201. [124] M . Frenklach, S imula t ion of surface reactions, Pure A p p l . C h e m . 70 (1998)477-484. [125] B. J. Boersma, Direct s imula t ion of a jet diffusion flame, Center for Turbulence Research ,Annual research briefs (1998) 47-56. [126] H . N . N a j m , P. S. Wyckoff , O . M . K n i o , A semi- impl ic i t numer ica l scheme for reacting flow, Journal of Compu ta t iona l Physics 143 (1998)381-402. [127] Information available at: h t tp : / / w w w - u n i x . m c s . a n l . g o v / m p i / . [128] Information available at: http: / / w w w . c a . s a n d i a . g o v / T N F / D a t a A r c h / F l a m e D . h t m l . [129] Proceedings of the T h i r d International Workshop o n Measure-ment and C o m p u t a t i o n of Turbulent N o n p r e m i x e d Flames, Sec-t ion 2 - Sandia / Darmstadt P i lo ted C H 4 / A i r F lame D . [130] M . Wang , J. H u a n g , W. K . Bushe, S imula t ion of a turbulent non-p remixed flame us ing condi t ional source-term est imation w i t h trajectory generated low-d imens iona l man i fo ld , Proceedings of the C o m b u s t i o n Institute. [131] J. Nafe , U . Maas , I l d m paper, Proceedings of the C o m b u s t i o n In-stitute 29(2002)1379-1385. 152 [132] A . C . H i n d m a r s h , O d e p a c k , a s y s t e m a t i z e d c o l l e c t i o n o f o d e s o l v e r s , S c i e n t i f i c C o m p u t i n g 1 (1983) 5 5 - 6 4 . [133] B . J . B o e r s m a , G . B r e t h o u w e r , F. T. M . N i e u w s t a d t , A n u m e r i c a l i n v e s t i g a t i o n o n the effect o f the i n f l o w c o n d i t i o n s o n the self-s i m i l a r r e g i o n o f a r o u n d jet, P h y s i c s o f F l u i d s 10 (1998) 8 9 9 - 9 0 9 . [134] I n f o r m a t i o n a v a i l a b l e at h t t p : / / w w w . ca . s a n d i a . g o v / T N F / 5 t h W o r k s h o p / T N F 5 . h t m l . [135] I n f o r m a t i o n a v a i l a b l e at: h t t p : / / d i e s e l . m e . b e r k e l e y . e d u / ~ g r i _ m e c h / n e w 2 1 / v e r s i o n 2 1 . [136] D . S c h m i d t , T. B l a s e n b r e y , U . M a a s , I n t r i n s i c l o w - d i m e n s i o n a l m a n i f o l d s o f s t r a i n e d a n d u n s t r a i n e d f l a m e s , C o m b u s t i o n T h e -o r y M o d e l l i n g 2 (1998) 135 -152 . [137] W . P. Jones , B . E . L a u n d e r , T h e p r e d i c t i o n o f l a m i n a r i z a t i o n w i t h a t w o - e q u a t i o n m o d e l o f t u r b u l e n c e , I n t e r n a t i o n a l J o u r n a l o f H e a t a n d M a s s Trans fe r 15 (1972) 301 -314 . [138] J. P. B o r i s , D . L . B o o k , F l u x - c o r r e c t e d t r a n s p o r t , i . shas ta , a f l u i d t r a n s p o r t a l g o r i t h m tha t w o r k s , J o u r n a l o f C o m p u t a t i o n a l P h y s i c s 11 (1) (1973) 3 8 - 6 9 . [139] G . D . S u l l i v a n , J. H u a n g , T. X . W a n g , W . K . B u s h e , S. N . R o g a k , E m i s s i o n s v a r i a b i l i t y i n . ga seous f u e l d i r e c t i n j e c t i o n c o m p r e s -s i o n i g n i t i o n c o m b u s t i o n , S o c i e t y o f A u t o m o t i v e E n g i n e e r s 114 (2005) 7 8 0 - 7 8 9 . [140] R . A . F rase r , D . L . S iebe r s , C . F. E d w a r d s , A u t o i g n i t i o n o f m e -t h a n e a n d n a t u r a l gas i n a s i m u l a t e d d i e s e l e n v i r o n m e n t , S A E T e c h n i c a l P a p e r Ser ies . [141] J . D . N a b e r , D . L . S i ebe r s , S. S. D i j u l i o , C . K . W e s t b r o o k , Effec ts o f n a t u r a l - g a s c o m p o s i t i o n o n i g n i t i o n d e l a y u n d e r d i e s e l c o n d i -t i o n s , C o m b u s t i o n a n d F l a m e 99 (1994) 192-200 . [142] T. I s h i y a m a , M . S h i o j i , et a l , C h a r a c t e r i s t i c s o f s p o n t a n e o u s i g -n i t i o n a n d c o m b u s t i o n i n u n s t e a d y h i g h - s p e e d g a s e o u s f u e l jets, J S A E / S A E I n t e r n a t i o n a l S p r i n g F u e l s & L u b r i c a n t s M e e t i n g . 153 [143] C . B l a i r , I m p l e m e n t a t i o n o f c o n d i t i o n a l s o u r c e - t e r m e s t i m a t i o n for p r e d i c t i o n o f m e t h a n e i g n i t i o n . , M a s t e r ' s thes i s , U n i v e r s i t y o f B r i t i s h C o l u m b i a (2003). [144] W . H . P re s s , S. A . T e u k o l s k y , W . T. V e t t e r l i n g , B . P. F l a n n e r y , N u -m e r i c a l R e c i p e s i n F o r t r a n , C a m b r i d g e U n i v e r s i t y P re s s , C a m -b r i d g e , U K , 1992. 154 Appendix A Inverting an integral equation in LFD I n L a m i n a r F l a m e l e t D e c o m p o s i t i o n ( L F D ) , the c o n d i t i o n a l a v e r a g e s are c a l c u l a t e d f r o m a c o m b i n a t i o n o f s o l u t i o n s i n the l i b r a r y as s h o w n i n E q . 2.48. T h e i n v e r s i o n is d o n e for the coef f ic ien t m a t r i x a; b y i n v e r t -i n g ~ r1 N f /= / Y^a&iiQPdC. ( A . l ) J o i=i T h e a l g o r i t h m d e s c r i b e d he re c a n b e f o u n d i n [59 ,144] . A p p l y i n g a q u a d r a t u r e r u l e to E q . A . l : ' l P(09iQdC * £ P ( a M C * H « i P ( C ' X ] p ( a ) ( A . 2 ) '° fc=l fc=l J<k K fc=i w h e r e K d e n o t e s the to t a l n u m b e r o f d i s c r e t i z e d p o i n t s i n £ space a n d wk d e n o t e s the w e i g h t o f the q u a d r a t u r e r u l e a n d 0 = 5^0^(0. ( A . 3 ) i=i 155 E q . A . 2 c a n be w r i t t e n for J c o m p u t a t i o n a l ce l l s o n a su r f ace u p o n w h i c h a i s a s s u m e d to be h o m o g e n e o u s : ( < /|C >*=!,»=! ••• </|C>l,JV \ V <T\C>K,I . . . <T\C>K,N ) I n br ie f : / « B a ( A . 4 ) ( A . 5 ) w h e r e B i s d e f i n e d a c c o r d i n g to E q . A . 5 , w h i c h c a n b e e x t e n d e d to a c o m b i n a t i o n o f m u l t i p l y sca la rs as d e f i n e d i n E q . 2.49. R w M a ( A . 6 ) w h e r e M is i n a f o r m s i m i l a r to B . T h e s o l u t i o n o f E q . A . 6 b y u s i n g E q . 3.5 i n a leas t s q u a r e s sense i s g i v e n b y s o l v i n g the l i n e a r s y s t e m o f e q u a t i o n s : ( M * M + A I ) a n = M * ^ + A a n _ i , ( A . 7 ) w h e r e I i s a n i d e n t i t y m a t r i x a n d A is a w e i g h t i n g coef f ic ien t . 156 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080430/manifest

Comment

Related Items