' .7, n’ "7" IS TflESl S 2:33 . LIBRARY 548/ was: Michigs " 'jjtate Umversnty This is to certify that the thesis entitled DIRECT NUMERICAL SIMULATION OF PREMIXED TURBULENT METHANE FLAMES MODELED WITH VARIOUS REDUCED MECHANISMS presented by NEGIN ANGOSHTARI has been accepted towards fulfillment of the requirements for the MS . degree In MECHANICAL ENGINEERING /fl/7~/ Major Professor's Signature 570 3/200 3 Date _ ... _V-._.-----.--a----.-u-.-o----o-.-.-o-.-.a--.-.---------o-.-.-----o-.-.-.-.-.-.-.-.---n- MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN Box to remove this checkout from your record. TO AVOID FINE; return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 6/01 cJCIRC/DateDue.p65—p.15 DIRECT NUMERICAL SIMULATION OF PREMIXED TURBULENT METHANE FLAMES MODELED WITH VARIOUS REDUCED MECHANISMS By NEGIN ANGOSHTARI A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE MECHANICAL ENGINEERING 2003 ABSTRACT DIRECT NUMERICAL SIMULATION OF PREMIXED TURBULENT METHANE FLAMES MODELED WITH VARIOUS REDUCED MECHANISMS By NEGIN ANGOSHTARI The mutual effects of turbulence and combustion are studied via analysis of computational data. generated by direct. numerical simulation (DNS) of a premixed methane—air flame in a shearless turbulence mixing layer. The chemical reaction is modeled with a 1-step. a 5-step and a 12-step reduced kinetics scheme. The 12- step mechanism is that of Sung et. al (1998) and is shown to accurately predict the characteristics of laminar premixed flames such as the flame thickness and the flame propagation speed. The flame behavior in various turbulent environments is investi- gated by changing the peak of the initial energy spectrum. Flame surface densities, flame strain and curvature statistics, the local structure of the flamelets. the local extinction and re—ignition, the turbulence generation and destruction by combustion and some other important. phenomena in turbulent combustion are investigated. The importance of the chemical kinetics models in representing the flow/flame features and the accuracy of “simple” (l-step and 5-step) reduced schemes are also studied by comparing the results obtained the more detailed kinetics mechanisms. For Kian: his love. guidance. support, and patience. iii Acknowledgcments I would like to thank the National Science Foundation and the American Chemical Society for the financial support of my graduate education. I would also like to thank my principal professor Dr. FA. J aberi for being a. good adviser and for his timely advices. I would like to thank my committee members Dr. Mei Zhuang. Dr. Tom Shih(.) and specially Dr. Indrek \Vichman for their helpful guidance through my research and also for their enlightening insights during my graduate studies. I am forever in the debt of Dr. Renate Snider, who showed me the way of technical writing. It was mainly with her excellent guidance that I gained the insight and courage of writing what. I have found in my research. I would like to thank RT. Edwards. for the permission to print some figures and a table. I would also thank Mr. Craig Gunn. for helping me brush up my writing and also for reminding me -with his pleasant attitude - to always seek the pleasant. parts in everything I do. I would like to thank all my professors and teachers. I still remember the first day of kintergarten, and I will always remember an immense number of sweet. memories from my school days. I would like to thank my family: my mom , my dad, and my three sisters for always believing in me and for their patience during my long school years. I can not even start to describe how much I appreciate their encouragement and sacrifices. It. was because of their encouragement, that I could made the long journey to United States to pursue my graduate studies. It was their belief in me that gave me the strength to beat the sad feeling of l‘iomesickness and be successful in pusuing my goals. And finally, I would like to thank my dear husband. Kian I\Iehravaran for being iv the best friend, partner, and soul—mate I could wish for. It. is because of what I have learned from him that I could do what I did. It. is in his way of life that. we can achieve anything if we are persistant enough. It is a sweet moment when one sees any once-thought-impossilfle task coming to reality. That is my best. lesson in life-so far! Table of Contents ABSTRACT ii LIST OF TABLES viii LIST OF FIGURES ix 1 Introduction 1 2 Governing Equations 7 2.1 Introduction ................................ 7 2.2 Conservation Equations ......................... 7 2.2.1 Conservation of mass ....................... 7 2.2.2 Conservation of momentum ................... 8 2.2.3 Conservation of energy ...................... 8 2.2.4 Conservation of species ...................... 9 2.2.5 Gas mixtures and ideal gas relations .............. 9 2.2.6 Stress Tensor ........................... 11 2.2.7 Conduction of heat and molecular diffusion of species ..... 11 2.2.8 Chemical kinetics ......................... 15 2.2.9 Stoichiometry and equivalence ratio ............... 18 3 Numerical Solution 20 vi 3.1 Improvements ............................... 21 4 Laminar l-D Results 23 4.1 Introduction ................................ 23 4.2 Conservation equations .......................... 24 4.3 Steady-state results ............................ 25 4.4 Results from transient solver ....................... 29 4.5 Premixed flame thickness ......................... 30 5 Turbulent Flames 33 5.1 Introduction ................................ 33 5.2 Initialization ................................ 34 5.3 Results ................................... 36 5.3.1 Flame contours .......................... 37 5.3.2 Statistics of the strain rate and the flame ............ 50 5.3.3 Local flame structure ....................... 53 6 Conclusion 63 APPENDIX 65 LIST OF REFERENCES 67 vii 1.1 3.1 3.2 Q"! i.. 5.2 5.3 List of Tables Classical regimes of turbulent premixed combustion. Source: ref. [1] Courtesy of RT. Edwards Publishers ................... Percent of time spent in each subroutine in the original formulation and after all the approximations/ optimizations. ............ Detailed distribution of CPU time inside the stress subroutine, before and after the modifications. ....................... Velocity statistics of the developed 2-D turbulence field ........ Summary of performed simulations ................... Calculated correlations between a) pubs" pmn ............. viii 22 36 53 1.1 1.2 1.3 4.1 4.2 4.3 4.4 4.5 4.6 4.7 5.1 List of Figures Flame wrinkling by turbulence. Source: ref. [1] Courtesy of RT. Edwards Publishers. ........................... Variation of the turbulent flame speed with RMS turbulent speed. SL is the laminar flame speed for the corresponding fuel. The flame. speed increases as the turbulent rms velocity increases, up to the middle of the bending zone. It levels off in the quenching zone, which eventually will lead to quenching. Source: ref. [1] Courtesy of RT. Edwards Publishers .................................. Classical turbulent. combustion diagram. Source: ref. [1]. Courtesy of RT. Edwards Publishers. ........................ Profiles of temperature. velocity. and density obtained by premix using the 12-step kinetics ............................ Temperature and mole fractions of CH4. C02 and CO obtained by premix using the 12-step kinetics .................... Temperature and mole fractions of C H 3, (71120 and (72114 obtained by premix using the 12-step kinetics .................... Temperature and mole fractions of H20. H 02 and OH obtained by premix using the 12—step kinetics .................... Temperature and mole fractions of 02. ['12 and H obtained by premix using the 12-step kinetics ......................... Evolution of temperature profiles in the 5-step mechanism.The solution moves to the cold side due to volumetric expansion of the flame and converges to a steady profile. ...................... The reaction rates for H20. CO. (.702. and CH4 calculated using the 12-step mechanism (ref. [9]) and our unsteady DNS code. The results are consistent with the calculations using GRI 2.11 as in ref. [23] . . . Schematic of the geometry of the problem ............... ix 26 27 27 28 28 3O 37 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 Distortion of the flame surface by turbulence as a. time sequence Instantaneous contour plots of a)vorticity b) temperature c) CH4 mass fraction (1) C 02 reaction rate for case 1 at time t/TO = 1 ....... Instantaneous contour plots of a)vorticity b) ten'iperature c) CH4 mass fraction (1) C ()2 reaction rate for case 2 at. time {/70 = 1 ....... Instantaneous contour plots of a)vorticity b) temperature c) CH4 mass fraction (1) C02 reaction rate for case 3 at time t / To = 1. The reaction is more present in the hot side compared to the sin'nilations with smaller 6/1 ratios. ................................. Instantaneous contour plots of a) vorticity b) temperature c) CO re- action rate (I) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 4 at time t/TO = 1 ................... Instantaneous contour plots of a) vorticity b) temperature c) CO re- action rate d) OH mass fraction e) CO mass fraction f) C H4 mass fraction for case 5 at time {/70 = ................... Instantaneous contour plots of a) vorticity b) temperature c) CO re- action rate (I) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 6 at time t/TO = 1 ................... Instantaneous contour plots of a) vorticity b) temperature c) CO re- action rate d) OH mass fraction e) C 0 mass fraction f) CH4 mass fraction for case 7 at time t/TO = 1 ................... Instantaneous contour plots of a) vorticity 1)) temperature c) CO re- action rate (1) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 8 at. time t / T0 = 1 ................... Instantaneous contour plots of a.) vorticity 1)) temperature c) C 0 re- action rate (1) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 9 at time t/To = 1 ................... The y-averaged temperature is plotted vs. x. along with the corre- sponding laminar profile for the 1-step simulations ........... The y-averaged temperature is plotted vs. x, along with the corre- sponding laminar profile for the 5-step simulations ........... 38 41 43 44 45 46 47 48 49 51 5.14 5.15 5.16 5.17 5.18 5.19 5.20 5.22 5.23 5.24 5.25 5.26 The y-averaged temperature is plotted vs. x, along with the corre- sponding laminar profile for the 12-step simulations Consumption speed ratio 3,, versus normalized tangential strain rate for case 1. Linear fit shows the trend. . . Consum tion 5 )eed ratio 5 versus normalized tangential strain rate P I n o for case 2. Linear fit shows the trend. . . Consumption speed ratio 3,, versus normalized tangential strain rate for case 3. Linear fit shows the trend. . . Consumption speed ratio 3n versus normalized tangential strain rate for case 4. Linear fit shows the trend. . . Consumption speed ratio 5,, versus normalized tangential strain rate for case 5. Linear fit shows the trend. . . Consumption speed ratio 3,, versus normalized tangential strain rate for case 6. Linear fit shows the trend. . . Consumption speed ratio 3,, versus normalized tangential strain rate for case 7. Linear fit shows the trend. . . Consumption speed ratio 3,, versus normalized tangential strain rate for case 8. Linear fit shows the trend. . . Consumption speed ratio sn versus normalized tangential strain rate for case 9. Linear fit shows the trend. . . Local profiles of OH reaction rate for a) case 5 b) case 6 c) case 8 d) case 9. Shown by 1, 2, 3, 4 are sections that. will later be examined in Figures 5.25 and 5.26 . . . . Local profiles of OH reaction rates for a) case 5 b) case 6 c) case 8 (1) case 9 Local profiles of C H4 reaction rate for a) case 5 b) case 6 0) case 8 d) case 9 xi 54 55 57 58 60 61 62 Chapter 1: Introduction Turbulence is present in nearly all practical flows. hence a great deal of effort has been expended on defining the physical basis of turbulent flows (ref. [2]) . Despite all efforts to date. turbulence is still not a fully known phenomenon. Direct numerical simulations (DNS) of very low Reynolds number flows have been successfully performed using best of the computational power currently availabe (ref. [3]). However, simulation of higher Reynolds number flows. s1.)ecifically those with more complex geometries or chemistry. requires nmch more computational power than what is presently available. The complete set of reactions even for the simplest hydrocarbon fuel (methane) contains more than 300 steps(ref). \Vitli the present computational power it is nearly impossible to even simulate a three dimensional box of turbulent flame. The reduced mechanisms are a good approximation that bring these simulations within reach. A premixed flame is defined as a flame in which the fuel and oxidizer are mixed homogenously before combustion (ref. [4]). A laminar flame can be visualized as the propagation of a flame front into a fresh gaseous reactant mixture. The propagation speed is usually between 20-100 cm / s for most hydrocarbon fuels at atmospheric pressure. The flame front thickness is of the order 0.1 111111. When turbulence is introduced into the flame field, the flame front will interact with the turbulent eddies. wrinkling and distorting the flame. As a result of this interaction, both the fuel consumption rate and the overall flame thickness may be increased. but this might not always be the case. Poinsot and Veynante (ref. [1]) have made a comprehensive review of turbulent premixed flames. An abbreviated review is discussed here. Damkohler was the first person to systematically investigate the turbulent combustion. He was also the first to introduce the idea of flame wrinkling and its Volume V Turhulem flame propagating against the mean flow at speed s T Flame front Cross section A: premlxed flow entering at speed ’T Fresh gas Burnt gas Figure 1.1: Flame wrinkling by turbulence. Source: ref. [1], Courtesy of RT. Ed- wards Publishers. eflect on controlling turbulent flames. He considered wrinkling as the main mechanism for controlling turbulent flames. He described the turbulent flame speed (57) as the speed required in the flame control volume to keep the flame stationary (Fig. 1.1). The turbulent flame speed and the overall reaction rate in the control volume are directly proportional. A number of correlations between turbulent flame speed and turbulent root mean square(RMS) velocity have been documented (refs. [5]-[6]). All of them express a similar trend: turbulent flame speed increases linearly with turbulent RMS velocity in the low turbulent zone. However, if the turbulent RMS velocity is increased, the growth of the turbulent flame speed will level off and eventually the flame will be quenched (Fig. 1.2). Damkohler explains that every point on the flame surface is moving with the laminar flame speed, thus he considers the increased turbulent flame speed as a consequence of the increased flame surface caused by flame wrinkling. By closely examining the available data so far a universal analytical expression for the dependence of ST on u’ has not been obtained (ref. [7]). Several other A Turbulent flame speed (ST‘ / . s 0 . j _ . . .s - L i . TLildeetil RMS) \GlOt fly I ; m the fresh gas m; ‘ * > Low turbulence zone Bendu g zone Quenching srzau' llmll Figure 1.2: Variation of the turbulent flame speed with RMS turbulent speed. SL is the laminar flame speed for the corresponding fuel. The flame speed increases as the turbulent rms velocity ii'icreases, up to the middle of the bending zone. It levels off in the quenching zone, which eventually will lead to quenching. Source: ref. [1], Courtesy of RT. Edwards Publishers. factors like boundary conditions, turbulent spectrum. initial conditions, and fuel type are also strongly involved (ref. [8]). A highly wrinkled flame sheet separates the reactants from the products. This flame sheet is being convected and bent. The convection speed depends on local conditions such as curvature of the flame sheet and strain rate. The flame sheet, which is considered as a reacting front of very small but finite thickness, acts as a volume source that affects the turbulence field. (ref. [4]) C Gunter-gradient. diffusion is the result of turbulent velocity fluctuations that force the reactants to move in a direction opposite to the mean reactant gradient. thr , as a result of density difference accross the flame, the velocity variance (turbulent energy) is increased by a factor of twenty. (ref. [4]) Turbulent premixed combustion can be described as the interaction of a flame front (with thickness 6 and speed 82) and an ensemble of eddies representing turbulence. These eddies range from the Kolmogorov scale (1”,) to the integral scale (It). The Danikohler number is defined as the ratio of the turbulent time scale to the flame time scale and is defined for the largest eddies. A high Damkiihler number indicates a. small chemical time scale compared to the eddy time scale, thus turbulence is not able to affect the flames inner structure, which still resembles the structure of the laminar flame. This condition is known as the flamelet limit. In this case the mean burning rate can be estimated from the burning rate of a laminar flame multiplied by the overall flame surface area. On the other hand, a small Damkohler number shows the flame will be highly affected by turbulence because velocity fluctuations are mixing the reactants and products, leading to the well stirred reactor limit. The Karlovitz number K0, which ('lescril‘)es the smallest eddies, is the ratio of the chemical time scale to the Kolmogorov time scale. W'ith Ka < 1, the chemical time scale is smaller than the turbulent time scale and the flame thickness is smaller than the Kolmogorov scale. As mentioned earlier, the flames inner structure is then similar to a laminar premixed flame, only it is now wrinkled by the turbulence. Based on the ratio of 111/82 this region can be divided into two categories: the wrinkled and the corrugated flamelet regimes. If u’ < U2, then the speed of turbulent motion is too low to wrinkle the flame up to the flame interactions. This is the wrinkled flamelet regime. However, when u.’ > SE, then the speed of turbulent motion is enough to wrinkle the flame front. up to the flame front interactions, leading to the formation of separate fresh and burnt gas pockets. When Ka > 1 and Do > 1, the turbulent time scale is still larger than the chemical time scale, but the Kolmogorov scales are smaller than flame thickness, and are capable of modifying the inner flame structure. This regime is called the thickened flame regime or the distributed reaction zones model. In this regime the stretch caused by the Kolmogorov scales is larger than the critical flame stretch. This can cause local flame quenching. A summary of all flame regimes is described Table 1.1: Classical regimes of turbulent premixed combustion. Source: ref. [1], Courtesy of RT. Edwards Publishers. ] K0 <1 (1),, >1) [ A}, >1 (1)., >1) 1),, <<1 7 Flamelets Thickened flames \Vell stirred reactor Flame is thinner Small turbulent scales All turbulent time scales than all turbulent scales may enter the flame front are smaller than the chemical time scale in Table 1.1 and Figure 1.3. The computational run-time in these DNS simulations can be reduced by using a reduced mechanism, which retains the most possible/ important reaction characteristics of the original reaction set. In the present work. DNS is used to solve the flame propagation problem. The only approximation here is the use of reduced kinetics, which is necessary because of today’s computational limitations. For this study we chose to use the 12-step Augmented Reduced Mechanism(AR.\I). (ref. [9]) This reduced mechanism is developed from GRI-AIECH 1.2, It has 16 SI)(‘(’1('.‘S(][2, I], ()2, 0”, I120, [[02, [1202, (Vilg, C'H4, C0, C02 . (11/20, (72112, 721/4, (lg/16, N2) and 1‘2 steps. It contains more non—steady-state intermediates than the conventional 4 0r 5—step mechanisms, therefore it predicts a wide range of combustion phenomena, including the global response of auto—ignition, laminar flame propagation and counter-flowing non-premixed systems under extensive thermodynamical parametric variations (ref. [10]). Validation comparisons show that the 10‘2—step mechanism with 22 species (ref. [11]) is the minimum size for a skeletal mechanism to be as comprehensive as the 12-step ARM. In order to assess the behavior of 12—step mechanism, some other simulations will be performed using both the stoichiometric 1-step (ref. [12]) methane-air mechanism and a 5-step reduced mechanism (ref. [13]). RMS velocity Well stirred reactor / flame speed ., o (U/SL) Gorrugated flamelets) ( Wrinkled flamelets D 1 Integral length scale / flame thickness (I /6) t Laminar combustion Figure 1.3: Classical turbulent. combustion diagram. Source: ref. [1], Courtesy of RT. Edwards Publishers. In this study, we will investigate the effect of an exothermic chemical reaction (in this instance a stoichiometric premixed methane-air mixture) on a developed homogeneous turbulent. field. First, We will optimize the code because even with a laminar 1-D configuration without optimization, the simulation takes almost 10 days to run on a Sun Enterprise 450. we will optimize the CPU time by examining a laminar premixed flame. The objective is to obtain the correct laminar propagation speed and determine whether the optimization affects the accuracy of the results. Next, we will introduce different. turbulent fields to study how the turbulence and chemical reaction fields interact. The turbulent fields will have different scales and intensities, which will enable us to investigate the effect of large and small scales of turbulent flow on the flame behavior. Chapter 2: Governing Equations 2. 1 Introduction In this section we. will describe the partial ('lifferential equations used in our simulations. The derivation of the equations can be found in standard texts (ref. [14]) in many different forms. We will concentrate on the final form of the equations used in this project. This involves choices for transport properties, chemical reaction mechanisms, inclusion or exclusion of radiation, bouyancy, and several other physical phenomena. The N avier-Stokes equations (mass and momentum transport) are used along with equations which represent the contribution of individual species to the total mass distribution. The total number of variables sums to N + 6. Having a large number of species (e.g. 30 species) has a large impact on the effort required to solve the governing equations. The total number of variables required to represent. a 3-D compressible multi-species flow are : density p, 3—D velocity field at, Energy 6, pressure p, and mass fractions Yk of the N species. As it will be shown in next section, pressure is obtained from an algebraic relation (the ideal gas law, eq. 2.6), which means we don’t need an extra differential equation in order to solve for pressure, thus N + 5 differential equations must be solved. 2.2 Conservation Equations 2.2.1 Conservation of mass The total mass conservation equation does not change in reacting flows: 0p apt: , ._ = 2. 1 2.2.2 Conservation of momentum The equation for conservation of momentum is the same in reacting and non-reacting flow, . , N7 IN! ('9 0 0]) 0T”- ‘ ()0! . ‘ — -+..—. =———+——+ Y..-=.’+ Eu".-. 2.2 01/1th 0.1310ou 01,] 0.” I); Afk.) 01.1 pk_1 Afr.) ( ) where kaj is the volume force acting on species L: in direction j. Even though this equation does not include explicit reaction terms. the flow is modified by combustion: the dynamic viscosity )1 changes because temperature varies in a ratio about 1 : 8. Density also changes in the same ratio and dilatation through the flame front increases. As a consequence, the local Reynolds number varies much more than a non reacting flow. Even though the momentum equations are the same with and without combustion, the behavior of reacting flows is very different. A typical example is found in jets: turbulent non reacting jets may become laminar once they are ignited. 2.2.3 Conservation of energy The energy conservation equation requires the greatest attention because multiple forms exist. The form used in this study is the total energy equation [)81 Oct 0 0Q: 0 , N I —=_ —" 1‘ :‘*“' _"‘ 1": ll .,',~, V.“ 2.3 p [)t at, + 81:1(pu 8‘) 01.1 + 817}- (Uju ) +Q +P; Lfk. (u + A.) f l where Q is the heat source term (caused, for example, by an electric spark or a radiation flux), not to be confused with the combustion heat release. Here, N p Z kakj (u,- + V“) 18 the power produced by the volume forces fk acting on k=1 species It. The complete form of the energy equation is not always needed and simplified forms are often utilized in combustion. These simplifications include the constant pressure (low Mach number) approximation, which leads to some terms being neglected; equal heat. capacites for all species; and constant heat capacity for the mixture only. The complete form of the total energy equation (2.3) is used without. any of the above a11)proximations. 2.2.4 Conservation of species Mass conservation for species A? is (9/)YK 0 —- ' ~ /.' '. = .’. 2. 0t +0.11, (p(u'+h")yk) M” ( 4) where Vim: is the 1” component of the diffusion velocity Vk of species A: and Q, is the reaction rate of species k. By definition: N N Z Yka. :0 and go, = 0. k=1 k=1 2.2.5 Gas mixtures and ideal gas relations The mass fraction is a suitable variable regularly used in reactive flow simulations, and is defined by Yk : _, (2.5) where mk is the mass of species It: and m is the total mass of the gas. The ideal gas equation of state is conveniently used in combustion calculations which does not introduce significant errors due to the conditions usually considered (low pressure and high temperatures): R p = pWT, (2.6) where R = 8.314 J . is the perfect gas constant, p is density ,and W is the mean molek molecular weight of the mixture: IV 1 ‘ Y,,. = , 2.7 W A2231”, < > where Wk is the molecular weight of species k. For a reacting flow, there are multiple possible variables to represent. energy or enthalpy. In this project, the total energy (sum of internal energy plus kinetic energy) will be used in the conservation equation for energy. The internal energy is related to the enthalpy by e = h — p/p, (2.8) , where N h = 2 mg, (2.9) k=1 and is the sum of two parts: T 11,, = cpde+ A129,]. , (2.10) To v chemical senatble where Ahg‘k is the enthalpy of formation of species Is: at a reference temperature To which is usually 298.15K, and cpk is the specific heat of species k at constant pressure, which usually increases with temperature. In this project, enthalpy and 10 specific heat. data were also taken from the CHEMKIN thermodynamic database. 2.2.6 Stress Tensor The viscous stress tensor TU- is 2 Oak - 011,- 011- ,.=__ _o,. ._ _——J— . 2.11 T" 3fl (9.17,, J + H (01')- + 01",) ( ) where the velocity components are u, i. = 1 - 3, ,u is the dynamic viscosity, and (51-1 is the Kroenecker delta function. The viscous stress tensor and pressure are often combined into one tensor 0,,- defined by: - - 2 011;, , 8:1» 811.- Uij : TU — p00 = —I)()U — glamou‘ + [1(0—11 + (7:) , (2.12) ., . ..J .4 where p is the static pressure (ref. [1]). In this thesis, the pure species viscosity is otained from (ref. [15]): r 7rm A. A: B F 5 ' ‘ = ———-.—-— 2.13 ml 16 ”052(22): ( ) In order to get the mixture averaged value, we use "— 2k: Xi,“ (2 14) __ k e . i=1 21:1 XJ‘I’kj 2.2.7 Conduction of heat and molecular diffusion of species The heat flux in direction i (q,- is) N 0T ‘ qt : —Aj()—;- + th'kf/kd's (2.15) " k=1 11 where (A) is the heat diffusion coefficient. The first term on the right hand side is the Fourier law of heat, conduction. The second term on the right. hand side is associated with the differential diffusion of species with different enthalpies. This is the Dufour effect which is the heat flux due to species mass fraction gradients and is often neglected. In this project, pure species thermal conductivities were taken from the TRANSPORT subroutine of the CHEMKIN package (ref. [15]). The mixture—averaged thermal conductivity is given as (ref: eq. 50 from TRANSPORT) k 1 1 kzl A'zl Ak In order to reduce the computational time without sacrifycing any accuracy one solution is to use mixture average conductivity coefficient. A polynomial is fitted using data from CHEMKIN to obtain pure species conductivities that are then employed to obtain the mixture averaged conductivity. The diffusion velocities are obtained from (ref. [14]) N X-X W p N vszz J k(vk—vj)+(is—Xj)—p—+;th(f.—fi-). (2.17) k=1 where Djk is the binary diffusion coefficient of species j into species k and X k is the mole fraction of species k7, i.e., X k :2 YkW/l/Vk. The Soret. effect (the diffusion of mass due to temperature gradients) is neglected. The system (2.17) is a system of linear equations with an 3N2 unknowns which would be very expensive to solve, considering that it has to be repeated for all points in the domain for all time steps. An approximation often employed in turbulent reacting flow codes is to use the Fick’s law through ka, the mixture averaged diffusion coefficient of species k, 12 (“)Y. a»; = -D........—*. (2.18) ()J',‘ where V,”- is the diffusion velocity of species A: in direction 1'. In most turbulent codes( since the turbulent mixinig dominates the molecular diffusion process), ka is obtained by assuming a constant Lewis number for each species k, A = pcpl)mk. (219) Lf’k The Lewis number Lek is interpreted as the ratio of molecular heat. diffusion to molecular species diffusion. In this study, the diffusion velocities are also obtained from mixture averaged diffusion coefficients, but. not by assuming a constant Lewis number for each species. This procedure gives results very close to the exact formulation here (2.17). V, = v, + m. + v... (2.20) where V), is the ordinary diffusion velocity (ref. [16]) (2.21) and where X k is the mole fraction, ka is the mixture averaged diffusion coefficient in terms of binary diffusion coefficient ka, which is defined as 1 — Y. ijk DA} The Soret effect (W'k) is non—zero only for low molecular weight species such as H, H2,H€Z 13 [2(1an 1 (IT 11".: ———.—.—, * X, T an (2.23) where 6;, is thermal diffusion ratio whose sign causes lower molecular weight species to diffuse from regions of low to higher temperature. In this study Wk is assumed negligible for all species. Finally, VC is the correction velocity and is included to ensure that mass fraction sum is unity, k 2 m1. = 0 (2.24) 1:21 In order to calculate the binary coefficients a 4”” degree polynomial was fitted to the data obtained from CHEMKIN for all the species involved in the simulation in the temperature range from 250K to 5000K. It can be shown (ref. [1]) that using Fick‘s law does not conserve the global mass. The conventional approach to solve. for N — 1 species and absorbing the difference in the last species as YN = 1 — [3:11 Yk can result in not conserving the total mass. In certain cases, counter-gradient transport for YN may be induced. This error is negligible, however, when the last species YN is a high concentration diluent such as N2 in air flames. Also, when all diffusion coefficients are equal or nearly equal, there is no significant error. This restriction is often used in turbulent combustion codes. Thus D), = D essentially causes the error in Eq. (2.4) to disappear. We calculated the correction velocity for species k as N Vch = —YA~ Z Vagkw (225) k=1 where V”, is the diffusion velocity of species k in the if” direction. The use of a correction velocity ensures that the species and mass conservation equations are 14 compatible. If we still solve for N — 1 species, the correct YN can be obtained. 2.2.8 Chemical kinetics Consider a system of N species reacting through M reactions: N N ZVIIUAII‘ : 214,121!" (2.26) k=1 k=l for " = 1,111, where Mk is a symbol for s )ecies k, V’, and V’.’ - are the molar . . L] k} stoichiometric coefficients of species A: in reaction j. Mass conservation enforces the condition: N N Z V1.1”; = 2 V1,]!- l'l'k (2.27) k=1 k=1 or N Z w..- W. = 0 (2.28) L21 forj=1...]l[. 14.] = Vii] — Vi.) (2.29) For sin‘iplicity, only mass reaction rates are used. For species 1:, this rate (12),. is the sum of rates oi, produced by all M reactions: where Qj is the rate of progress of reaction j. Summing all reaction rates to... one obtains: M N N 2% = E Q2 2 ”"sz = 0 (2.30) k=l k=1 1:1 showing that total mass is conserved. The progress rate (21- of reaction j is written as 15 n Vk) _ , ~ pli- . N M- . Qj—Afjl—I ff: —l\rJ-Hk=1 ff: . (2.31) k=1 where K!) and K”- are the forward and reverse rates of reaction j. and ka/W’. = lel), is the molar concentratitm of species 1:. These rate constants constitute a central problem of combustitfm modelling. They are usually modelled using the empirical Arrhenius law . ,. ' E ' 721‘ Ifo- = ADI”) exp -RT) = ADTHJ exp — ,1’). (2.32) Expressing the individual progress rates Q J for each reaction means providing data for the pre—exponential constants A“. the temperature exponents .3] and the activation temperatures T0}- for each of the j separate reactants. Before giving these constants, even identifying which species and which reactions should or should not be kept. in a given scheme is the challenge of chemical kinetics. The backwards rates Kr}- are computed from the forward rates through the equilibrium constants as follows K;- A'rj : \v J 9 (2:33) Pa Z V, , as“) All" ' (H’I‘) k l JCXI)( R - R'I‘ where Pa = 1 bar. The A symbols refer to changes occuring when passing from reactants to products in the 1” reaction. The quantities AH}Q and AS? are respectively enthalpy and entropy changes for reaction j. It is important to mention that data on Q]- correspond to models : except for certain reactions, these data are obtained experimentally. The values of kinetic parameters are often disputed in the kinetics community. For many flames, very different schemes are found in the literature with comparable accuracy. The choice of a scheme is thus a difficult and controversial task. Unftn‘tunately for numerical 16 combustion the values of the parameters used to calculated Q]- and the stiffness associated with the determination of Q )- create a central difficulty for numerical combustion; the space and time scales corresponding to the Q 1- terms are usually very small and require meshes and time steps which can be orders of magnitude smaller than in non-reacting flows. Particular attention must be paid to the activation energy Ej : the exponential dependence of rates to E 1- leads to considerable difficulties when E] takes large values (typically larger than 60 Area] / mole). A given combustion code for laminar flames may work perfectly with a given mesh when E j is small but will require many more points to resolve the flame structure when Ej is increased even by a small amount. In this study, we have chosen the following three mechanisms: the 1-, 5-, and 12-step schemes of [12], [13] and [9], respectively. The 1—step mechanism is assumed to be an irreversible global reaction CH4 + 20-; -—-> (720 + 21-1-20. (2.34) The chemical reaction rate is written as where Ta = 14906K, which yields an activation energy of E = 29.618kcal / mole and the frequency factor is BO, 2 5.2 x 1013 (ref. [12]). The 5-step chemical reaction mechanism is (ref. [13]): 311-; + 02 + CO =‘— 31-120 + CO 1‘12 + OH :1‘ 2H20 3H; + CO 2 H20 + CH4 ['12 + CO; 2 H20 + CO 3H2 + C0 + 2N0 = 3H.2() + CO + N2 17 '1 Tim, mac. (TIME; (2.35) The 12-step chemical reaction mechanism ([9]) consists of the following steps: O; + 2CO .——* 2COg H + 02 + CO = OH + C02 H; + ()2 + CO = H + OH + CO; HO; + CO = OH + CO; 02 +112()2 + CO = OH + 1102 + C02 02 + CgHg = H + C()2 02 + (111;, + CO + (72114 = CH4 + CO; + C1120 + C2H2 02 + 2le1;, = H; + CH4 + C02 02 + 2CH3 + CO = (II-1.4+ COQ + CHQO ()2 + CH3 + CO = H + C02 + CHgO ()2 + CO + CQHG = CH4 + C02 + CHQO H+OH .——: 1120 2.2.9 Stoichiometry and equivalence ratio The most important species present. in a flame are the fuel and the oxidizer. Thus, the definition of stoichiometric ratio and equivalence ratio are based on their mass fractions, although the definition of equivalence ratio depends also on the configuration of the burner (premixed or non premixed). 14F + V'OO —> Products, (2.36) 18 H A u‘ (if where V’F and MO are the the number of moles of fuel and oxidizer when all of the reactants are burnt to products. Thus, (Efi =J$22=3, man yp st VP” p where s is called the mass stoichiometric ratio. The equivalence ratio is defined as yrrp‘ (YI‘1)/()/F‘) (15 = 8 , = -:— , 2.38 3'0 y 0 } O st ( ) ¢=ST£, 93% 77.2.0 01' where mp and mo are the mass flow rates of fuel and oxidizer, respectively. The equivalence ratio is an important parameter. Rich combustion has to > 1 (the fuel is in excess) while lean combustion has gb < 1 (the oxidizer is in excess). In premixed hydrocarbon/ air flames, fuel is mixed with 02 and N2 with typically 3.76 moles of Nitrogen per 1 mole. of Oxygen. Since the sum of mass fractions must be unity, the fuel mass fraction is given by: 1 2}: (2m) 1 +5, (1+ 3.76%) 19 Chapter 3: Numerical Solution Compact finite differencing [17] was used to approximate the derivatives. Equations (2.1)-(2.4) were integrated using the extended MacCormack method [18]. A brief summary of this method can be found in the Appendix. Non—reflecting boundary conditions were used for the inflow [19] and the outflow [20] boundaries. One goal of this study was to reduce the CPU time in order to make the turbulent 2-D and 3—D calculations affordable. First, the original transport and thermal properties are reviewed, then the aptn‘oximat ions are explained and justified in the following section. As seen in Table 3.1 the ”stress" subroutine was consuming 86% of all the CPU time. out of which 87% was used to calculate the diffusion coefficients. The auxiliary relations required to solve the conservation equations are described fully in the previous chapter, both in the exact and in the. approximated form. An acceptable simplification for viscosity (ref. [15]) is to use k n— — Z—————— X” ~ZiInI. (3.1) k=lzj=JlX®kJ 2 Before this approximation the calculation CPU time for viscosity coefficient was high due to the two summations (one in the nominator and the other in the denominator). The value of 7),? is obtained from the polynomial fit to CHEMKIN data from 250 to 5000 K. In order to make the code fully independent of calling CHEMKIN, tables were constructed for specific heat (Cp) and enthalpy (h) values. 20 Table 3.1: Percent of time spent in each subroutine in the original fornmlation and after all the approximations/ optimizations. No. Subroutine Function percent CPU percent CPU name time before time after 1 Continuity Solves continuity 0.020 0.081 2 Y~Moment um Solves Y momentum 0.013 0.084 3 Energy Solves enery equation 0.015 0.067 4 Species Eq. Solves Species 0225 1.183 5 C alc. Interm. Updates intermediate 0.056 0.230 values 6 Avg. Mol. Wt. Calculates average 0.033 0.140 molecular weight 7 Cale. Temp. Calculates 0.174 0.972 temperature 8 Pressure Calculates pressure 0.002 0.010 9 B.Values Imposes B.C.’s 0.150 0.705 10 Stress Sub. Calculates thermal and 86.024 8.036 transport propertires 11 Chem. Kinetics Calculates species 13.287 88.492 reaction rates 3. 1 Improvements Table 3.1 shows the percentage of total time spent in various parts of the program. For example, evaluating thermodynamic properties and the stress tensor (stress subroutine took 86% of the total ro ram run time before the modifications P g were implemented. After the modifications, this subroutine took only 8% of the total run time. The overall effect of the optimizations was to reduce the total ro ram run time by about 70% without affectin the accuracv of results. P g . g . Table 3.2 shows the details of the transport properties as CPU time inside the stress subroutine. Most. of the CPU time was used to compute diffusion velocities. 21 Table 3.2: Detailed distribution of CPU time inside the stress subroutine, before and after the modifications. No. Subroutine Function percent CPU percent CPU time before time after 20 Tempfsz Vel. Der. Calculates velocity and 0.005 0.257 temperature derivatives 21 Mass F rac. Der. Calculation of mass 0.126 10.195 fraction derivatives 22 Enthalpy Calculates enthalpy 0.074 4.161 23 Viscosity Calculates Viscosity 12.192 3.114 24 Diffusion Vel. Calculates diffusion 87.193 74.008 velocity 25 Cond. Heat. flux Calculates conductive 0.408 8.265 heat flux 22 Chapter 4: Laminar 1-D Results 4. 1 Introduction The first step in the turbulent premixed simulations is having correct and steady solutions for the laminar non-strained premixed flame. Specialized computer codes like PREMIX have been developed for this purpose (ref. [21]). The CHEMKIN thermodynamic database (ref. [15]) and the corresponding subroutines (ref. [22]) are used in PREMIX to accurately obtain the thermodynamic and transport properties. The concentration and temperature profiles obtained are used to initialize the turbulent. flame calculations. Understanding the physics and finding a valid mathematical model to the one-dimensional laminar flame propagating into a premixed gas (methane-air for example), has been the subject of research for many years. There are four ways to obtain a laminar premixed flame solution. 1. Analytical solution, which is obtained by solving the 1-D equations in the flame frame of reference. This is only possible when the chemistry and transport models are relatively simple. Despite the fact that they work only for simplified situations, significant. insight can be achieved. 2. Steady-state numerical solution. This is a boundary—value problem in which the flame speed has to be obtained as the eigenvalue. This type of solution is used in computer codes like PREMIX (ref. [21]). 3. Transient solutions that are the more general form of the above method, as they have the transient term. When allowed to progress sufficiently in time, it can be used to obtain the steady-state solution. 4. Experiment; measuring the distribution of concentrations in a stable 23 configuration is a valuable source of data for constructing the detailed mechanisms. After a brief review of the cmiservation equations required for 1-D premixed flames, a discussion of their solution method will be presented, followed by results from the numerical simulation of 1—D premixed flames using different kinetics. 4.2 Conservation equations For laminar one-dimensit)nal premixed flames, the ccmservation equations are simplified as (ref. [1]): o Massconservation 0p Bpu 2% 0r ( ) 0 Species conservation. For k = 1 to N - 1: (9ka 0 l / _ , at. + a: (001- + lk) 1k) — Wk (4-2) 0 Energy conservation 0T er a arr an: N C —— w— =0." — A,— — —_—- C Y)" 4.3 p p(0t+u(2.r) T+Br( 0.1:) 0.1: (2:; Wk k ( ) with . I IV } . w' = - lkcdk T zit-=1 24 The momentum equation is not used unless it is necessary to calculate the pressure using the velocity field. These equations describe a wave that propagates from the burnt gas side to the fresh gas side. The propagation speed reaches a constant value S'L when transients are ignored. If the equations are. written in the reference frame of the flame. the flame will be fixed and the reactants will move towards the flame with speed SL. pu = constant = plSL (4.4) 0 , , . .—- (MU + H) 3k) = wk ('45) 01‘ 0T 3 0T 0T N C _— 2 LL95 +— A— — —+— C, .Y.V. 4.6 p path I + 0.1? < 8.1?) 6.1? pg "k A k ( ) Using chemical reaction and molecular diffusion models (such as Fick‘s law), the above set, of equations can be closed and solved. provided that the proper boundary conditions are utilized. There is a mathematical obstacle regarding the cold boundary difficulty that is discussed in ([14]). The solution of this system in an eigenvalue problem by nature where the flame speed SL is the eigenvalue. 4.3 Steady-state results In the last twenty years, many mathematical tools (refs) have been developed to solve Eqs.(4.4) to (4.6). When proper boundary conditions are set up and the problem is discretized on a finite difference grid. the resulting system is a strongly nonlinear boundary value problem. An excellent tool specifically developed for solving this nonlinear boundary I' I Y 3000 am“ ........................... at v (cm/s) ,. --_-‘- P (g/cm3)| 'I’ s 2500' 0 o 0.05 0.1 0.15 0.2 x (cm) Figure 4.1: Profiles of temperature, velocity, and density obtained by premix using the 12—step kinetics value problem is PREMIX ([21]). Figure 4.1 shows the profiles of temperature, density, and velocity obtained by PREMIX using the 12-step methane/ air mechanism ([9]). The flame speed obtained is 42c-m/s. The results for most of the 16 species involved in the 12—step mechanism are shown in Figs 4.2-4.5. The Temperature profile is also demonstrated in all of these figures as a reference to the location in the flame. PREMIX uses an adaptive grid to resolve the high-gradient locations by implementing closer grids. The adaptive process and Newton iteration are repeated until the desired accuracy is reached. The above process does not always converge to a solution, but once a solution is obtained for the prescribed conditions. the user can continue the process for other conditions (such as a different stoichiometric ratios) with fewer iterations. 26 mole fraction G I" 11-"" . 5“ o 0.05 0.1 0.15 0.2 x (cm) Figure 4.2: Temperature and mole fractions of C H4, C02 and CO obtained by premix using the 12-step kinetics 0.1 0.09 0.08 ‘ 0.07 C 30.06 8 “20.05 0 750.04 0.03 0.02 0.01 0 0.05 0.1 0.15 0.2 ka) Figure 4.3: Temperature and mole fractions of CH3, CH20 and C2H4 obtained by premix using the 12-step kinetics 27 0.2 .0 _a U" mole fraction 0 0.05 Figure 4.4: Temperature and mole fractions of 11-20, H 02 and OH obtained by premix using the 12-step kinetics 0.25 V V Y .0 3 mole fraction .0 0.05 O 0.08 0. I 0. | 5 0.2 Figure 4.5: Temperature and mole fractions of 02. H2 and H obtained by premix using the 12-step kinetics 28 4.4 Results from transient solver PREMIX solves the steady-state eigenvalue problem. Initial conditions and estimated profiles should be close enough to the final solution in order to obtain a converged solution. Specifying initial conditions can become frustrating for a detailed mechanism, since it involves a trial and error process. Also, in cases of non-elementary reactions, a converged solution is not always obtained. It is possible to use a transient solver (the same solver used in turlg)ulence calculations) in one-dimension to obtain the laminar premixed solution, provided that the diffusion velocities are calculated accurately and enough time is allowed for ignition of the flame. While not an efficient procedure for this purpose, it guarantees a solution without the necessity of a close initial estimate. In this study, PREMIX is used to obtain a solution for the 12-step mechanism. The solution is then fed to the transient solver to adjust it with the approximations used for the diffusion velocities. It is the latter solution that is used in the turbulent calculations. In the case of the global one—step reaction and the 5-step mechanism, the solution is obtained by assuming a hyperbolic tangent profile for the temperature and species distribution based on stoichiometric values of the reactants and products. For example, in the case of temperature, the following distribution is assumed: T0) = g (T. + T2) - g (T. — T2) tanh (d x y), (4.7) Where y varies between {—éLy. $151,] to assure that the flame is positioned at the origin (zero coordinates). The temperature is T1 at. the cold side and T2 at the hot side, which is lower than the adiabatic flame temperature. The quantity (i is a constant that determines the thickness over which the profiles vary. Similar distributions are used for species concentrations that change between stoichiometric 29 AMA “UV I I I I I I I 2000- 1800' 1600- 1400- T (K) 1200- time (intervals of 0.8 ms) 1000' 800- 600* 400+ “10.2 —0.15 —0.1 —0.05 0 0.05 0.1 0.15 0.2 x(cm) Figure 4.6: Evolution of temperature profiles in the 5-step mechanism.The solution moves to the cold side due to volumetric expansion of the flame and converges to a steady profile. parameter values. The ignition of the flame in the 5-step mechanism is shown in Fig. (4.6). The solution moves to the cold side due to volumetric expansion of the flame and also as a result of the flame propagation itself, and then converges to a steady profile. The CH4 reaction rate using the 12-step mechanism (ref [9]) and our unsteady DNS code (Fig. 4.7), gives the flame speed as 41.5 cm/s. The calculated value of laminar speed is 14.88 cm/s and 29 cm / s for the 1- and the 5-step mechanisms, respectively. 4.5 Premixed flame thickness To estimate the required mesh size in a direct numerical simulation, it is useful to define a flame thickness since the flame structure must be resolved adequately. 30 .ul U1 0 U1 0 I U1 I O ( 0.05 0.1 0.15 0.2 x (cm) molar production rate (10'3 gmol/cm3.s) 0 Figure 4.7: The reaction rates for I120, C0, C02, and CH4 calculated using the 12-step mechanism (ref. [9]) and our unsteady DN S code. The results are consistent with the calculations using GRI 2.11 as in ref. [23] The other usage of such a definition is the availability of a reference length in order to compare with the other scales of the flow, such as turbulence the scales. Different definitions exist for premixed flames thickness. Some of them are based on concentration or temperature profiles that obviously require a solution, while others are based only on scaling laws and can be used as first estimates. A flame thickness 6 based on scaling can be shown, /\1 5 = —, pICpSL (4.8) which merely gives the order of flame thickness for the simplified reactions. A temperature-based thickness may be constructed by defining the distance over which the reduced temperature varies from 0.01 to 0.99. The slow reactions taking place in the burnt gases will lead to bigger thickness values, which can be 31 misleading if used in computations. A more useful thickness based on temperatre can be defined, for example, by - T2 — T1 60 : ——.'\—-. 4.9 L Imagg) ( ) It is also possible to define the thickness based on the species or the radical concentratirms. However, these values would be very different. from values obtained based on the previously mentioned definitoins. because the radicals may occupy a much smaller region than the temperattire-variation region defined above. This issue is evident in detailed chemistry problems, due to the presence of radicals. 32 Chapter 5: Turbulent Flames 5. 1 Introduction Turbulent premixed flames behave in a very different way than laminar premixed flames. Providing that we move with the flame. all the laminar flame charecteristics for the flames discussed in the previous chapter (that is, freely propagating and fixed flames) such as concentraion or reaction rate distributions (and as a result flame speed) will remain the same. However, in turbulent flames there is a range of time and length scales that leads to flame distortion. This distortion along with the various curvatures present, will in turn cause strain on the flame sheet. The curvature variation also influences the flame behavior. On the other hand, flames affect the turbulence in their own way. The heat release caused by the flame will cause a volumetric expansion that smoothes out. the vorticity filed. Moreover, the heat release will lower the viscosity, causing the reduction of the effective Reynolds number, which leads to faster dissipation. Thus, the small scales will be diminished more rapidly. \Nith the presently available computational power, simulationg a three dimensional turbulent flame with a complex chemistry mechanism is not easily within reach. For economical reasons we are using 2-D turbulence in this work, despite the differences between 2—D and 3—D turbulence. In a 2-D turbulence there is no vortex stretching mechanism, thus the energy cascade that supplies the energy to the smaller size eddies is not complete. Therefore, the small scale vortices are not as influential on the flow and the flame as they are in 3—D turbulence. In addition, because of the broken energy cascade, the small scales do not follow the Kolmogorov law. However, further investigation in 3-D premixed flames reveals that the flame propagation topology tends to be locally two dimensional, especially in the highest 33 curvature regions. (refs. [24], [25], [26]). This fact can be used as a partial justification for utilizing 2-D turbulence (ref. [27]). 5.2 Initialization The simulation of the 1-D laminar premixed flame was discussed in the previous chapter. The same solution is used to initialize the turbulent 2-D premixed flame simulations. In order to achieve that, a well developed 2-D turbulence field is imposed on the presently laminar premixed flame. The resulting turbulent premixed flame will be allowed to evolve according to the corresponding the turbulence field‘s eddy turn-over time. Various phenomena might be observed after the mentioned time is reached, such as stretching and bending of the flame surface, and occasionally tearing of the flame surface and the formation of isolated packets of unreacted or quenched gas. The turbulence field is initialized as a random solenoidal velocity field with a zero mean and a Gaussian spectral density funtion 2 EU.) 2 a4 exp {—2 (5) }. A70 where c is a constant calculated based on the initial turbulent energy, I: is the wave number, k0 is the wave number containing the maximum energy (in this case k0 = 5), and E is the energy density (ref. [28]). Since this field is not compatible with the Navier-Stokes equations, it is allowed to evolve with a 2-D spectral code, until the enstrophy dissipation is maximized. The enstrophy dissipation is defined n=QVXaq, 34 Table 5.1: Velocity statistics of the developed 2-D turbulence field a, o n I l 1 L Rm] HeL H 0.359 8.94 2.75[0.06[0.427 9711230321] where w is the vorticity vector. The developed field has the charachteristics described in Table 5.1. Here, 0,, = %u,~u, is the incompressible turbulent energy or velocity variance, 9 = éwgwg is the enstrophy, and l = (MIMI/2 and L = 03/2n1/3 are the flow micro and macro scales, respectively, while Re, and ReL are their corresponsing Reynolds numbers. This field is initially specified as a non-dimensional field. The advantage of doing so lies in the fact that. the scaling techniques will enable us to use exactly the same field, although the actual physical size varies. This procedure provides the opportunity to experiment with various flame thickness to turbulent scale ratios, while maintaining the same turbulent field. Three different sizes are considered, and simulations are performed for each size with the three different reaction mechanisms, resulting in total of nine cases (Table 5.2). In all nine cases the rms of velocity (u’) is made equal to 60000m/s using the appropriate scaling factor. The 1-D laminar solution (from the previous chapter), is placed in the middle of the field (where :r = 0). Initially there is no variation in the y-direction. Then the turbulence field is imposed on the grids. Since the prepared turbulence field has 512by512 grid points, which is consistent with the flame field, it is possible establish a correspondence between each point in the turbulent field to one point in the flame field. Hence, the ratio of the turbulence length scale to the flame thickness can be varied without changing the turbulence field. The same unsteady code used in the previous chapter to simulate the laminar flame is also used in order to simulate the 2-D turbulent premixed flame. The flame 35 Table 5.2: Summary of performed simulations No. Steps L Int. scale [A scale ReL He, 7;: Tt(.S‘(3C) 6 (5 / L (cm) L(cm) 1(cm) (cm) 1 1 2.9204 0.2 0.0284 7762 1102 400 3.3e-5 0.052 1.83 2 1 1.4602 0.1 0.0142 3881 551 400 1.67e—5 0.052 3.66 3 1 0.7301 0.05 0.0071 1941 276 400 8.3e-6 0.052 7.32 4 5 2.9204 0.2 0.0284 7762 1102 194 3.3e—5 0.044 1.55 5 5 1.4602 0.1 0.0142 3881 551 194 1.67e-5 0.044 3.1 6 5 0.7301 0.05 0.0071 1941 276 194 8.3e-6 0.044 6.2 7 12 2.9204 0.2 0.0284 7762 1102 143 3.3e—5 0.038 1.34 8 12 1.4602 0.1 0.0142 3881 551 143 1.67e-5 0.038 2.68 9 12 0.7301 0.05 0.0071 1941 276 143 8.3e-6 0.038 5.32 propatation is mainly in the x-direction. The boundary condition in the x—direction is non-reflecting (ref. [19]), hence it allows for expansion in the x—direction. The boundary condition in the y-direction is periodic as there is not considerable flame propagation in the y-direction. A schematic of the problem geometry is shown in Figure 5.1. The MacCormack method is used to integrate the equations in time along with 4th order compact method for differencing. Further explanations can be found in Appendix A. As we further progress in time, the flame is gradually distorted, stretched and bent. by the turbulence. A sequence of the flame contours can be seen in Figure 5.2, which can also demonstrate tearing off of the flame and production of isolated flame packets. 5.3 Results In this section, we discuss the results for the turbulent combustion simulations. We group the results into three categories. First, the flame contours will provide the 36 l (_SiOlCl110melI ic [ mixture of methane and .m . [ l uompuoo Mepunoq ogpoued Periodic Non-reflecting boundary condition Figure 5.1: Schematic of the geometry of the problem opportunity to investigate the general view of the whole flame structure, without the necessity of looking into details. Second, in the statistical study we will demonstrate the trends associated with dependence of the flame speed on the tangential strain rate. Last, the flame structure is investigated through turbulent flame profiles in selected sections and is compared with the laminar distribution. Images in this thesis are presented in color. 5.3.1 Flame contours The results of the turbulent flame simulations are presented here in t / To = 1 and as flame contours for several parameters in Figures 5.3 to 5.11. Each set of figures is associated with a reaction mechanism and is plotted for vorticity, temperature, CH4 mass fraction for all mechanism. However, for the l-step mechanism the 002 reaction rate is plotted, while for the other two mechanisms (5- and 12—steps) the CO reaction rate and, the CO and OH mass fractions are plotted. The vorticity distribution does not vary dramatically with various mechanisms or 6/l ratios. However, as the 6/l increases, the vorticity dissipation rate will increase. This fact can be demonstrated by calculating the normalized vorticity 37 f) Figure 5.2: Distortion of the flame surface by turbulence as a time sequence 6) 38 (fl/{20) for the three different box sizes. These values are: 0.93, 0.83, and 0.77 for box sizes of 2.93mi. 1.46cm. and 0.736771, respectively. This happens as a result. of the drop in the effective Reynolds number as the box size decreases (Figures 5.3(a) to 5.11(a)). The flame thickness broadens as displayed in the temperature contours (Figures 5.3(b) to 5.11(b)). The same concept is better viewed by looking at the reaction rate contours (Figures 5.3(c) to 5.11(c)). In the l—step temperature contours((Figures 5.3(b) to 5.5(b))), the temperature in the burnt section is nearly uniform and no noticeable vartiation is observed. However, in temperature contours for 5— and 12-step mechanisms on the burnt side, there are isolated packets whose temperatures are lower that the adiabatic flame temperature. The proportionate area of these packets increases as the box size decreases. This observation is consistent with the fact that a increase in 6/1 will not provide the sufficient length for the flame to reach the adiabatic flame temperature. However, the laminar flame thickness of the 1-step mechanism is larger than that of 5- and 12-step mechanisms. Thus the usual expectation would be that the temperature would reach the adiabatic flame temperature in a longer length, but in reality this will not happen. The fact that l-step model leads to less temperature variations in the burnt section, demonstrates the smaller effect of turbulence on the flames modeled with the l-step mechanism. This can be further investigated by examining Figs. 5.15 to 5.23. In the reaction rate contours (Figures 5.3(c) to 5.11(c)). the thickening of the flame is observed. In these concours (e.g. Fig. 5.11(c)), it is very well demonstrated that the flame thickness varies based on the local curvature. In the region were that flame curvature is concave towards the burnt gas side, flame stretching and the near tearing off of the flame is observed. However, when in the. regions with a curvature is convex towards the burnt gases, the flame has thickened considerably. These discussions can also be applied to the other species mass fraction 39 contours. In the larger 6/l values, there are streaks of radicals present in the burnt gases, due to the fact that the flame did not have enough time to reach the same stage as the flame in the bigger boxes. This happens because all of these simulations are performed up to the same nondimensionalized time scale, but the actual flame time in the three different boxes is also different. The overall picture shows that the three mechanisms are not behaving the same. Moreover, the peak of the species concentration for each species is different in each of the three mechanisms. This also was seen for the laminar flames in the previous chapter. To compare the average turbulent flame thickness with laminar flame thickness, calculated using the three different mechanisms, the turbulent temperature average in y-direction is calculated for the three different 6/l ratios. It is then plotted along with the corresponding laminar temperature profile for each mechanism Figs 5.12-5.14. All are averaged at t / To = 1. The average temperature profiles for all of the mechanisms have the same trend. Although, as the U6 increases, the temperature average becomes smoother, which shows the smaller influence of turbulence on the flame. 40 8) Y'LL c) VAL 0° 0.5 I x/L 0.5 X/ L Figure 5.3: Instantaneous contour plots of a)vorticity b) temperature C) CH, mass fraction d) CO; reaction rate for case 1 at time t/ro = 41 3) V'LL 0) V'LL 0.5 x/ L Figure 5.4: Instantaneous contour plots of a)vorticity b) temperature C) CH, mass fraction d) CO; reaction rate for case 2 at time t/To = 1 42 a) C) o o 0.5 x/L 0.5 x/L Figure 5.5: Instantaneous contour plots of a)vorticity b) temperature c) CH4 mass fraction (1) ()0; reaction rate for case 3 at time t / To = 1. The reaction is more present in the hot side compared to the simulations with smaller 6/1 ratios. 43 a) b) 4 2 _J .J :0. o :0. -2 .4 4 6 oo 0.5 I "o x/ L C) in. >~ 0.5 x/ L e) do. > 0.5 x/ L 0.5 x/L Figure 5.6: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate d) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 4 at time t/To = 1 a) b) '4'- «asso'oeo co 9 v! .- ° C) e) 0.5 0.5 x/L x/L Figure 5.7: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate d) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 5 at time t/To = 1 45 a) 10) VAL bégkonao vy- .. O C) VAL 6) VAL 0 0 0.5 l o x/ L 0.5 x/L Figure 5.8: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate d) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 6 at time t/‘ro = 1 46 a) b) vy- 3 h o u A 0 Y4}- 50 VI o. x/L c) d) {4L VAL o.5 x/L 6) f) Y'LL as 05 x/L x/L Figure 5.9: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate d) 0H mass fraction e) CO mass fraction f) CH4 mass fraction for case 7 at time t/ro = 1 47 8) {LL 0) VAL . 0.5 x/L x/L 6) f) VAL 0.5 x/L Figure 5.10: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate (1) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 8 at time t / To = 1 48 a) b) VAL os x/L c) d) I 25 02 . . a15 Q1 .5 ‘ one o 1 ms ' 0.1 o o as l x/L VAL VAL xiL 6) f) VAL 05 05 x/L x/L Figure 5.11: Instantaneous contour plots of a) vorticity b) temperature c) CO reaction rate (1) OH mass fraction e) CO mass fraction f) CH4 mass fraction for case 9 at time t/To = 1 49 2200. — 8/I=1.83 , —- 5/l=3.66 200°" — 5n=7.32 1800- — laminar 1600 Q] 400> :1 zool 1000- 800- 600 400- A l 1 0.] 0.2 0.3 0.4 0.5 0.6 x (cm) Figure 5.12: The y-averaged temperature is plotted vs. x, along with the correspond- ing laminar profile for the 1—step simulations 5.3.2 Statistics of the strain rate and the flame In this section the statistic of consumption speed ratio is investigated. The consumption ratio is defined as ff; LDC”, dn .x) . (f—oo w’C114d72)0 Sn: 9 where n is the direction normal to the flame surface and 0 denotes the intial (laminar) state. The unit normal vector is defined as 77 _ va' ' (Veb-V¢)1/2’ where (,b is the variable used for defining the flame surface. This quantity is plotted versus normalized tangential strain rate ((1,) that is defined as 50 4 — 5/|=1.55 200° , — 5/l=3.1 13004 — 6/I=6.2 I i — laminar i 1600-“ A1400- 25,200 |.— 1000- 800' . 600' i 400- < 0.] 0.2 0.3 0.4 0.5 0.6 x (cm) Figure 5.13: The y-averaged temperature is plotted vs. x, along with the correspond- ing laminar profile for the 5—step simulations 2200.. — 5/I;1.34 T 1 i 4 ._ as: 1800-1 — laminar = l600 ’ 4 $1400 ' ‘ t-l 200 ' ‘ lOOO *- ‘ 800 r ‘ 600 * 400L 0.] 0:2 0:3 0:4 0:5 0.6 x (cm) Figure 5.14: The y-averaged temperature is plotted vs. x, along with the correspond- ing laminar profile for the lQ—step simulations 51 (1,217:V17, where F is the unit tangent vector on the flame surface. For all the simulations, scatter plots of 3,, versus normalized strain rate am) are given in Figures 5.15 to 5.23. In this study, the location of flame surface is defined where the temperature is equal to 1800K. In case 1 (l-step. smallest 6/1), a large degree of scatter is observed. However, an overall trend showing a decrease in the consumption speed as the tangential strain increases is also observed (ref [20]). For the larger 6/1 ratios, the scatter decreases and there is also a better correlation between 3,, and (1,70 (Figs. 5.16 and 5.17). In the 5-stcp mechanism and the smallest 6/1 (case 4, Fig. 5.18), much less scatter is observed compared to the its l-step counterpart. This scatter is further reduced as the 6/1 is increased (case 6, Fig 5.20), which is very close to a line. Vt’hen comparing to cases 4 and 5, it can be seen that. the slope of that line decreases, which means that. the effect of turbulence on the flame is reduced to a mere deformation. In the 12—step mechanism and the smallest 6/1 (case 7, Fig. 5.21), the correlation between 5,, and (1,70 is not as good as the 5-step case (case 4), although it is less scattered compared to case 1. In case 8 scatter amount is low (Fir 5.22), but again in case 9 the scatter is increased. This demonstrates the sensitivity of the 12-step mechanism to the flame thickness and the turbulence scale ratio. Again, the three mechanisms show different. trends. Depending on the location in the flame, the flame speed might. be affected by the flame front curvature. In order the to study this effect, it. is necessary to first. 52 Table 5.3: Calculated correlations between a) pawn pm" case no. 1 2 3 4 5 6 7 8 9 I pat," -().-11 -0.42 -0.57 —0.17 —0.(i3 -0.($(i —0.25 -0.:32 -0.55 pm," 0.05 -0.0()4 0.021 -0.0(i 0.03 0.083 0.003 0.09 0.017 calculate the curvature. The curvature is defined as 1. 0.92 ()s2 where (.r. y) are the coordinates of two points on the. flame surface and s is the arc length on the surface subject to the condition 8(1' = 0) = 0. To quantify this dependence, the correlation coefficient between the flame speed and the curvature and also between the flame speed and the tangential strain rate is calculated (Table. 5.3). For the simulations, no correlation is observed between the flame speed and the curvature. However, the correlation between tangential strain rate and flame speed is apparent. The strongest correlation is seen in the cases that use 5—step mechanism. Moreover, this correlation coefficient. increases as the 6/1 increases. which is consistent with the. previous observation of the scatter reduction as a result of increased 6/1. 011 the other hand, a smaller correlation indicates that parameters other than the tangential strain rate, such as the turli)ulent eddy size are it'ifiuential in the flame speed variation. 5.3.3 Local flame structure In order to study the local flame structures, the OH radical is a suitable choice, because it is a suitable scale to demonstrate the flame sensitivity to turbulence. This is due to the fact that the OH radical exists in a thin region in the flame sheet. Figure 5.15: Consumption speed ratio 5,, versus normalized tangential strain rate for case 1. Linear fit. shows the trend. Figure 5.16: Consumption speed ratio 3,, versus normalized tangential strain rate for case 2. Linear fit shows the trend. 54 4 eat To Figure 5.17: Consumption speed ratio 3,, versus normalized tangential strain rate for case 3. Linear fit shows the trend. 2 1 T I I f -\ Fl ,3 1.8~ , (Q o 1.6. V Cr [7‘6 C. - 1.4- - 1.2- - m“ 1. 4 0.3h i 0.6’ 1 0,4. ‘ 0.2» 4 o I l l l 1 -6 -4 -2 4 6 a I t: 0 Figure 5.18: Consumption speed ratio 3,, versus normalized tangential strain rate for case 4. Linear fit shows the trend. 55 \‘g’ ,A. 3’ f1- )0 0 £01“; '7' 1 8" g] '1‘- " .' {Jr 15‘ A, 1.6" F 5 C) 1'4. ‘ 1.2- ' U): 1' 0.8' 0.6“ d 0.4” d 0.2' ' o l g I, I L -6 - -2 0 4 6 2 a T t 0 Figure 5.19: Consumption speed ratio 3,, versus normalized tangential strain rate for case 5. Linear fit shows the trend. 1.8? 20° ‘ 1.4- 1_2. in“ 1' 0.8" 0.6' J 0.4" f 0.2' 4 96 -; _. 0 ‘ l 6 4 at to Figure 5.20: Consumption speed ratio 3,, versus normalized tangential strain rate for case 6. Linear fit shows the trend. 56 Figure 5.21: Consumption speed ratio 3,, versus normalized tangential strain rate for case 7. Linear fit shows the trend. Figure 5.22: Consumption speed ratio 3,, versus normalized tangential strain rate for case 8. Linear fit shows the trend. 57 l O i- i- . p t- Figure 5.23: Consumption speed ratio 3,, versus normalized tangential strain rate. for case 9. Linear fit shows the trend. and it is very much dependent on the distribution of the other radicals, and it. is present in both the 5- and the 12-step mechanisms. In other words, the OH radical tracks the flame front. The contours plotted in Figures 5.24 (a) to 5.24 ((1) show the OH reaction rate contours for the amount of —0.01 and +0.01. For the cases modeled with the 12-step mechanism the line representing the positive value of OH production is no longer continuous due to the fact that the OH radical production rate is not high enough. Therefore, it is necessary to use the existing radicals, rather than using the ones that are normally produced. This demonstrates the effect of turbulence in the flame. The l2-step mechanism demonstrates great sensitivity to the OH production. This sensitivity does not exist if the 5-step mechanism is used. On the other hand, turbulence with smaller length scales does not seem to have a substantial effect on the OH reaction rate. In Figures 5.24 (a) to 5.24 ((1), four different sections are chosen in the direction normal to the flame surface to show the flame properties distribution in that direction (Figs 5.25 and 5.26). It is also possible to compare these profiles with the laminar flame profiles. Compared to the laminar flame distributions. Fig. 5.25 shows that the profiles of the 12-step mechanism are changing much more than the 5-step flame. Furthermore. the peak amounts of the two mechanisms are different. The same behavior is also observed in ('11.; reaction rate profiles. The 5-step profiles are very much like the laminar case. which confirms the smaller influence of the turbulence on the 5-step flaine(cases 6 and 9 . Figs 5.26 (b) and ((1)). 59 a) 5-step (5/1 = 3.1 b) 5-step 6/1 = 6.2 \YL ‘. 7’. ’il’ Y w v . f I fi F " ‘10,, r' ‘u-S < ' {I ’3'.’ \" 4 Y 3‘ “\t/ I . -. . \ . . \ 3*. ' " . ' l -_ _ ‘ > \ ' .' 1 .1 \ ‘l i) ‘3' ’ \ v' PL? _\ ll / \‘j‘ 4 -‘ I", \ / , . ‘. 7 '. /' l. ‘L. P . .‘ ,’2' 4- d ‘ /'b: / .\ V‘.\\ .‘\ p- ‘ I. [it i! s l ,\\ \\' //‘~. “ . y“. ’ ///,\\ ’ ‘ .. /I- \ I .\ '\.- 1 \\.\.' .,.: L l A L A A \; 1' 1A A}; 4 A L c) 1'2-step 6/1 = 3.66 (l) 1‘2-step (5/1 = 7.32 T Yf TY/rfi v T 77 Y H V‘ ‘ Tj-V‘ ,"r -‘ \v' ,' l _ q . .‘. :l‘--"‘-. h f\ "" \ t . ,. ~ - . . , . / \_~ . ' .3 . . ‘ .1 , l) I _- ‘: l l l v . ‘x l l L ‘ \ 5.}, ' l 4 l )l 4 l \ If ,’ j . . \ ,- . f 4 . / '1‘1‘ , I .. \ { i.” L _‘ /’J ~—" 4.. 1 - '~ ’ 's . \‘ / I ‘ \‘~. J/ _¥",,-— / ‘ .‘ \ , , _ ‘_-, \ - ~\__ -’:,,\. Vi.“ \ \ __. . P». \ :1 . 25“ k‘ui—x ‘ ‘ .. R .H- 1 l" ‘ _ “ \ \ 4:12.. \‘ 1 “—\ \-2 \ ‘\~\ I \\,‘ l \h,‘ l ., l ‘ l ’ I , I . 4 > I, f' / . r ’7 4;“ , _ I ‘ III l 7 b 4 f V ‘- x‘ i‘ it \l ‘ l ' ' 1I L {J t _. ll l If"‘ \ I." ’ \ i l ‘ 11‘ \ ‘~ 4 if ‘ ‘ ‘ i \\ \‘ l - ‘\ . ‘ \_ / l ( ,,' 1‘ . A T . ‘ 4' i ‘ A i L, 4 A .. ‘1.._L._._.| .LJLL__'¢ as, ,L _ 7 _L _ Figure 5.24: Local profiles of OH reaction rate for a) case 5 1)) case 6 c) case 8 d) case 9. Shown by 1, 2, 3. 4 are sections that will later be examined in Figures 5.25 and 5.26 a.) 5—step (5/1 = 3.1 0.] I I‘Jax.’ Section 1 . . Section 2 « _ ' i —- Section 3 l J g 77* Section4 -0-l' ' ----~ Laminar ‘ -0.02 -o.01 0 0.071 ' _ ' 0.02 c) 1‘2-step (VI = 3.66 i / 0.03% / l 0.02. i. i -0.01l l--- Section] ' -* Section 2 l . — Section 3 i i Laminar -0.02} -o.03- —0.04 -0.03 o .<——.—~—A—__ .— .—L. .. -0.02 -o.01 0 0.01 0.02 o. b) 5-step (5/1 = 6.2 I ' 1 0.1 l "f .~'?rs.. . 0.05 07—“. ~\ r—-—— —-- -—-— ~———- ~‘:\ i —— Section] \\ ' Section 2 . : —— Section 3 * .\ . -0.0S‘ - l —--~ Laminar -0.i * 4| " 0.00s — 0.6" 73.01 S —0.005 0‘— d) 1°2-step 6/l = 7.32 -0.0| S -0.01 0.03 ,I l l A 0.02L 0.01} 0m§§ ' ‘\T.\ -0.01 I \. i Section 2 ; — Section 3 ' ---~ Laminar -0.02;» 003: \1 " 'i —— Section I— l l -0.03”"10.02 —0.01 o 0.01 0.02 0.03 Figure 5.25: Local profiles of OH reaction rates for a) case 5 b) case 6 c) case 8 d) case 9 61 a) 5-step (5/1 = 3.1 b) 5-step 6/l = 6.2 I ____ _.-._‘____.._ -0.1' s t — Section] ._.__._,__ ,3/ - Section 2 : £32; .. . Section 3 ' -0-l 5! — Section 3 I/ 1 Section 4 i ‘ Section 4 \‘v/ m---” i . __ 4 fl"? 1“"“Lammar _ l__..___._-_.-__- ___-__ A A ._ A ... -00. 15 -0.01 -0.005 0 0.005 0.0] 0.015 | | . 90.203? —0.01 0033 ‘ o 000?— o.oTi'/"o.ois c) 12—step 6/1 = 3.66 [__ 0 “\ . \i‘ -o.osi | «0.1 . ‘ — _ __ l ' ' 9 , ---- Section] i i — Section] l 1 . I; . - Section 2 0 ]Si _.’ ' Section 2 -o,15; \ I i —- Section 3 _ ' i 1‘ j" — Section 3 i V» me Lama: fl . L--- Laminar_J 1 ‘ , g, . . -o.03“ -o.02"_——0'.0—i o 0.01 0.02 ‘0-03 ‘0-02 '00] 0 0-0] 0-02 Figure 5.26: Local profiles of C H4 reaction rate for a) case 5 b) case 6 0) case 8 d) case 9 62 Chapter 6: Conclusion DN S of turbulent premixed flames were performed with three different chemical reaction mechanisms and also with different turbulent length scales. The thermal and transport formulae where adjusted for the maximum computational efficiency, while keeping the accuracy in a very acceptable level. The laminar flames with three different mechanisms were solved as an unsteady problem. The solution was continued until a steady state was reached. The species distributions were consistent with both experiments (ref [23]) and the laminar flame calculations of other investigators (refs [0] [10]). This laminar flame solution was used as the flame initialization for the turbluent field that has already evolved using a spectral code. The turbulent. field was scaled into three different physical box sizes. This was done to maintain the exact same characteristics of the turbulence field, while changing the size. Therefore, it was possible to experiment with different i‘nechanisms and different ratios of flame thickness to turbulent length scale. As expected, the l—step mechanism (lid not perform well, and sometimes it did not. even show the same trends as the other two mechanisms. There was much less variation in 1-step contour plots, for example, than in the other mechanisms’ contour plots. The l-step mechanism was specifically unable to show the small flame-like activities that were still present in the burnt side of flames modeled with 5- and 12—step mechanisms. In the contour plots, the flame stretching effect was clearly seen. It. "as also observed that when the flame curvature is convex to the hot side, the flame is thicker, and when it is convex to the cold (unburnt) side, the flame is thinner. Overall, the flame speed demonstrated a. noticeable dependency on the strain rate, especially for the smaller box sizes. This fact is stronger in the 5—step 63 mechanism. For the range of variable chosen in this work. the curvature did not show any correlation with the flame speed. In conclusion. the 1‘2—step i’r’iechanism has shown much a bigger sensitivity to turbulence. For example for the 1‘2-step contours the OH procution is very low, which raises the expectation that this mechanism might perform better at predicting the tear—offs and pocket formations whose existence has been showwn by ex1_)e.riments. Nevertheless, the final verification of which mechanism really performs better, and under what condtions. remains until simulations are performed with the full mechanism (GRI-Mech 3.0). 64 APPENDIX 65 The MacCormack scheme was extended by using forward-backward compact (.lifferencing (ref. [30]) Cf (B$A)f.’_ +(Bi/l)f’ I I I i z 1 1+1 o._ +a,+o, =—:l: , f’ l f f“ Ar Ar where o = :1, A = g B 2 :13, G = —V—:‘3. For the one-dimensional model equation, the solution for the new time step n + 1 is obtained from the solution in time step 71, through U; = U," — AA+F" U," = U,’ — AA‘F‘ 1 Lfer—l = _2_ [07:1 + [fin] 66 [ll [6] l7] l8] [9] [11] [12} [13] [141 [16] Bibliography T. Poinsot. and D. Veynante, editors, Theoretical and Numerical Combustion, R.T. Edwards, Philadelphia, PA, 2001. S. B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, UK, 2000. P. Mom and K. Mahesh, Direct Numerical Simulation .' A tool in turbulence research, Arum. Rev. Fluid Mech. 30, 539 (1998). S. Pope, Turbulent Premixed Flames, Ann. Rev. Fluid Mech. 19, 237 (1987). B. Hakberg and A. Gosman, Analytical determination of turbulent flame speed from combustion models, 20th (Intl.) Symposium on Combustion , The. Combustion Institute , 225 (1984). N. Peters, The turbulent burning velocity for large-scale and small-scale turbulence, J. Fluid Mech. 384, 107 (1999). K. Kuo, Principles of Combustion, John \K’iley, 1986. R. Cheng and I. Shepherd, The influence of burner geometry on prcmircd turbulent flame propagation, Combust.. Flame 85, 7 (1991). C. Sung, C. Law, and J. Chen. An Augmented Reduced Mechanism for Methane Oridation with comprehensive Global Parametric Validation, 27th Symposium on Combustion/ The Combustion Institute , 295 (1998). C. Sung, C. Law, and J. Chen, Further Validation of an Augmented Reduced Mechanism for Methane Oxidation .' Comparison of Global Parameters and Detailed Structure, Combust. Sci. and Tech. 156, 201 (2000). A. Kazakov and M. Frenklach, Reduced reaction sets based on GRI—Mcch 1.2. http: //www.ME.Berkeleyedu/drm . M. Bui-Pham, Studies in Structures of Laminar Hydrocarbon Flames, Ph.D Thesis, UCSD (1992). H. Mallampalli, Evaluation of CH4/NOJ: Global Mechanisms Used for Modeling Lean Premixed Turbulent Combustion of Natural Gas, MS. Thesis, Chemical Engineering Department, Brigham Young University, Provo, UT (1996). F. Williams. Combustion Theory. Perseus Publishing, 1985. R. J. Kee, F. M. Rupley, and J. A. Miller, The chemkin thermodynamic data base, Technical Report SAND87-8215B, Sandia National Laboratories. Livermore, CA, 1991. C. F. C urtiss and J. O. Hirschfelder. Journal of Chemical Physics 17, 550 (1949) 67 [17] [18] [19] [23] [‘24] [20'] [27] [28] [29] [30] S. K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comp. Phys. (1992). C. Kennedy and M. Carpenter, Comparison of Several numerical methods for simulation of compressible shear layers, NASA Technical paper . T. Poinsot and S. Lele, Boundary Conditions for Direct simulations of Compressible viscous flows, Journal of Computational physics (1992). D. Rudy and J. Strikwerda, Boundary conditions for subsonic compressible Navier-Stohes calculations, Computers and Fluids (1981). R. J. Kee, J. F. Grear, M. D. Smooke, and J. A. Miller, A fortran program for modeling steady laminar one—dimensional premixed flames, Technical Report SAN D85-8240, Sandia National Laboratories, Livermore, CA, 1991. R. J. Kee et al., A fortran computer code package for the evaluation of gas-phase multicomponent transport. properties, Technical Report SAND86—8246, Sandia National Laboratories, Livermore, CA, 1990. S. R. Turns, An introduction to combustion: concepts and applications, McGraw Hill, 2000. W. Ashurst, Geometry of premixed flames in three dimensional turbulence, Technical Report Proc. Summer Program, Center for Turbulence Research, Stanford University, CA, 1990. R. S. Cant, C. J. Rutland, and A. Trouve, Statistics for laminar flamelet modeling, Technical Report Proc. Summer Program, Center for Turbulence Research, Stanford University, CA, 1990. S. S. Girimaji and S. B. Pope, Propagating sufaces in isotropic turbulence, Journal of Fluid Mech. 234, 247 (1992). M. Baum, T. Poinsot, D. Haworth, and N. Darabiha, Direct Numerical Simulation of Hg/OQ/NQ flames with Complex Chemistry in Two—dimensional Turbulent Flows, J. Fluid Mech. 281, 1 (1994). F. A. Jaberi and S. James, Eflects of chemical reaction on two-dimensional turbulence, Journal of Scientific Computing 14, 31 (1999). M. Baum, T. Poinsot, and D. Thevenin, Accurate Boundary Conditions for Multicomponent Reactive Flows, J. Comp. Physics 116, 247 (1994). C. Kennedy and M. Carpenter, Several new numerical methods for compressible shear-layer simulations, Applied Numerical Methods . 68