. 255“ :4 K“ : .4. Huh: .. w a. In a. 1|va . (IL . it}? {fit . V mix } a; .7 a .1. wfixgfi. at: . ‘ ,s _‘ _ , Y but. . 2: .1 , x . .. a 5... .373 . SW?“ . , '5‘ 1.4 ‘ SQ.“ J4 . ,. 9... .. an: .. . fiw ,ufiu? ‘ can“? I . :1. . . 011:4“ tn. .1 EM. , .93. hasub 3 .1. a. £1 42mg. ... 1 2.3.9.7.... .mk A .u... 3...... . a I hag/mus... i .2 '1 ..§...r... . $.51. 7 3»? g. ..... 15.1.15! a .1: x... a. hr- 0.»!b‘haxl kg L a? _. f i éfiégggé 9% ; LIDHAHY Michigan State University This is to certify that the dissertation entitled MATHEMATICAL MODELING AND NUMERICAL SIMULATIONS OF CHEMICALLY REACTING TURBULENT JETS presented by KIAN MEHRAVARAN has been accepted towards fulfillment of the requirements for the PhD. degree in Mechanical Engineering V (Major “nature MON 2 G . 20 o 5 a / Date MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN Box to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE MATHEMATICAL MODELLING AND NUMERICAL SIMULATIONS OF CHEMICALLY REACTING TURBULENT JETS by KIAN MEHRAVARAN A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY MECHANICAL ENGINEERING 2005 ABSTRACT MATHEMATICAL MODELLING AND NUMERICAL SIMULATIONS OF CHEMICALLY REACTING TURBULENT JETS BY KIAN MEHRAVARAN Chemically reacting turbulent flows are investigated using numerical techniques. The main objectives of this work could be classified into two categories: 1. To understand the impact of gravity on transitional and turbulent jet flames. 2. To develop models and algorithms for accurate and efficient simulations of high- speed turbulent flows with hydrocarbon combustion. In (1), both Direct Numerical Simulations (DNS), and Large Eddy Simulations (LES) of turbulent jet flames have been performed under different gravity conditions. Simpli- fied chemistry without using a model was employed in DNS, while in LES the Filtered Mass Density Function (FMDF) method was used as a closure for the composition and reaction of methane/ air combustion. Both the DN S, and the LES results are consistent with previous findings and in- dicate that in the absence of gravity, combustion damps the flow instability; hence reduces ”turbulence production” and jet growth. However, in the ”finite-gravity” conditions, combustion generated density variations may promote turbulence and en- hance both the mixing and the combustion through buoyancy effects. Our results also indicate that the gravity effects on a transitional / turbulent jet flame is not limited to large-scale flame flickering, and there is a significant impact on small—scale turbulence and mixing as well. Furthermore, the analysis of compositional flame structures suggests that the finite-rate chemistry effects are more significant in finite-gravity conditions than in zero—gravity. In (2) the LES/FMDF method is extended towards multi-step chemistry with realistic thermodynamic properties, such that the predictions could be compared with laboratory flame measurements. In LES / FMDF, the effects of chemical reactions appear in closed forms, allowing for a reliable prediction of complex turbulent reacting flows. The consistency of the Eulerian and the Lagrangian solutions are discussed, and an efl'icient algorithm for parallelization of the hybrid code is presented. Comparisons with the Sandia’s piloted methane jet flames (flame D and F) are performed and good agreements in the case of the near-equilibrium flame D have been achieved with a flamelet-based chemistry model. Finite-rate reduced kinetics in the form of the global l-step, and multi-step kinetics are employed as well. While reasonable predictions have been made in the case of flame D using the 1-step reaction, the extinction predicted by the multi—step methods is more pronounced compared to the experimental measurements. The observations are analyzed and the parameters responsible for this difference are identified. To Negin iv Acknowledgements I would like to take the opportunity to thank the NASA Microgravity Combustion Program and the ONR, who made this work possible. I would like to thank Dr. Farhad A. Jaberi my advisor for his support and scientific insight. I would also like to thank my committee members, teachers, and mentors: Dr. Uday Hegde, Dr. Manooch Koochesfahani, Dr. Tom Shih, Dr. C. Y. Wang, and Dr. Indrek Wichman. I am grateful to Dr. John Foss and Dr. Charles Petty for their valuable insight on turbulent and multi—phase flows. Finally, I would like to thank my wife Negin for her continuous support and patience throughout the years of my PhD research. Table of Contents Abstract ....................................... ii List of Tables .................................... ix List of Figures ................................... 1 Direct Numerical Simulation of Transitional and 'Ihrbulent Buoyant Planar Jet Flames .................................... 1.1 Introduction ................................ 1.2 Governing Equations and Computational Methodology ........ 4 1.3 Results - Non Reacting Jet ........................ 7 1.4 Results - Reacting Jet .......................... 11 1.4.1 T wo-Dimensional Simulations .................. 12 1.4.2 Three-Dimensional Simulations ................. 16 1.4.3 Comparison Between 2—D and 3-D Simulations ......... 26 1.5 Summary and Conclusions ........................ 27 1.6 Figures and Tables ............................ 29 1.7 Nomenclature ............................... 56 Filtered Mass Density Function for the Large-Eddy Simulations of Complex Turbulent Hydrocarbon Flames ........................ 57 2.1 Introduction ................................ 57 2.2 Formulation ................................ 61 vi 2.3 Numerical Solution Procedure ...................... 65 2.3.1 The Lagrangian Monte Carlo Procedure ............ 65 2.3.2 Consistency ............................ 68 2.3.3 Parallelization ........................... 69 2.3.4 Implementation of Chemistry .................. 72 2.3.5 Computational Domain and Parameters ............ 74 2.4 Results for Turbulent Methane Jet Flames ............... 75 2.5 Conclusions ................................ 80 2.6 Figures ................................... 82 2.7 Nomenclature ............................... 101 Large Eddy Simulations of Buoyant and Non-Buoyant Turbulent Jet Flames 104 3.1 Introduction ................................ 104 3.2 Formulation ................................ 106 3.3 Numerical Simulation Procedure ..................... 106 3.3.1 Implementation of chemistry ................... 107 3.4 Results ................................... 108 3.4.1 Physical flow structure ...................... 109 3.4.2 Statistical description ....................... 111 3.5 Conclusions ................................ 116 3.6 Figures and Tables ............................ 118 vii 3.7 Nomenclature ............................... 4 Conclusions and Summary ........................... Bibliography .................................... viii 140 143 List of Tables 1.1 Flow parameters for 2-D simulations. Re = 3000, Ce = 2.5 for all cases except for case 5 which has Re = 6000 ............... 29 1.2 Flow parameters for the 3-D Simulations. Re = 3000, U2 = 0.2 and Ce = 2.5 for all cases ............................ 29 3.1 The chemistry, the form of “energy equation” and the source term used in the direct and flamelet approaches. f is the conserved scalar and G is BRT / 8 f . ................................ 119 3.2 Identification of the simulations ...................... 119 ix List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 1.11 1.12 1.13 Schematic representation of the domain and the coordinate system used in the 3D simulations. The flow is in the upward direction. ISO-surfaces of vorticity in the 3-D nonreacting jet at the final time of the simulation, t/to = 4. The domain shown is 1 < :r/h < 16, —4 < y/h < 4,0 < z/h < 4. .......................... Instantaneous vorticity magnitude contours (0—10) in the 3-D nonre- acting planar jet, at z/h = 2,0 < x/h <16,—4 < y/h < 4. ..... Half-width of the 3-D nonreacting jet along the jet. A straight line, 0.087 (i + 0.23) is fitted to the self-similar portion of the jet growth rate. .................................... Mass entrainment from the lateral boundaries and increase in the axial mass flux. ................................. Comparison of (a) mean axial velocity, (b) longitudinal Reynolds stress, with experimental and DN S data. Our data is from 33/12 = 15. Comparison of (a) lateral Reynolds stress and (b) Reynolds shear stress profiles with experimental and DN S data. Our data is from x/h = 15. Comparison of total “product Integrals” (lower curves) and total ki- netic energy (upper curves) for different Da numbers. ........ Instantaneous temperature contours (1-3.25) for the F r = 20 (left), and F r = 00 (right) 2-D jet flames. The domain shown is 0 < z/h < 24, —4 < y/h < 4. ............................ Comparison of Product Integrals (lower curves) and total kinetic en- ergy (upper curves) for different F 7' numbers. ............. Instantaneous temperature contours (300~2200 K) for the buoyant (left), and non-buoyant (right) 2-D methane jet flames. The domain shown is 0 < :r/h <16,—3.6 < y/h < 3.6. ................... Jet half—width of the 2-D methane jet flames. ............. Total 002 mass fraction in the planes normal to the 3: axis, for the methane jets. ............................... 30 31 32 33 33 34 35 36 1.14 Scatter plot of reaction rate versus mixture fraction for the methane flame at :1: = 15h; (a) buoyant, (b) non-buoyant jets. ......... 1.15 Instantaneous vorticity magnitude (left, 0-3), and temperature (right, 1-3.25) contours in an :1: - y cross section of the 3-D jet flame with Ze = 8, Fr = 00 (case 12). The domain shown is 0 < :c/h < 16, —3.6 < y/h < 3.6. ................................ 1.16 Instantaneous vorticity magnitude (left, 0—4.5), and temperature (right, 1-3) contours in an III-y cross section of the 3-D jet flame with Z 6 = 8, Fr = 40 (case 13). The domain shown is 0 < x/h < 16, —3.6 < y/h < 3.6. .................................... 1.17 Instantaneous vorticity magnitude (left, 045), and temperature (right, 1-3.5) contours in an 23—31 cross section of the 3-D jet flame with Z e = 9, Fr = 00 (case 14). The domain shown is 0 < :c/h < 16, —3.6 < y/h < 3.6. .................................... 1.18 Instantaneous vorticity magnitude (left, 0-7), and temperature (right, 1-3.5) contours in an x—y cross section of the 3-D jet flame with Z 6 = 9, Fr = 40 (case 15). The domain shown is 0 < .r/h < 16, —3.6 < y/h < 3.6. .................................... 1.19 Vorticity magnitude contours in the plane normal to the jet axis at x/h = 15 for case 12 (Ze = 8, Fr = 00) (left, 0-2.5), and case 13 (Ze = 8, Fr = 40) (right, 0-9). The domain shown is —-3.5 < y/h < 3.5,0 < z/h < 4. ............................. 1.20 (a) Time-averaged values of temperature for Z 6 = 8 cases, (b) Time- averaged values of temperature for Z 8 = 9 cases. The overbar repre- sents the time-averaged quantities. ................... 1.21 Jet half-width for reacting cases. .................... 1.22 Total product mass fraction in the planes normal to the 2: axis. 1.23 Profiles of mean velocity for Z 6 = 8. .................. 1.24 Center-line velocity excess for 3-D reacting jets. ............ 1.25 Vertical velocity on the lateral boundaries of reacting jets (positive sign is inward flow). .............................. 1.26 Profiles of longitudinal velocity fluctuations at 23/ h = 15. ....... xi 47 48 48 1.27 1.28 1.29 1.30 1.31 1.32 1.33 1.34 1.35 1.36 2.1 2.2 2.3 2.4 Ratio of turbulence intensity in the finite—gravity jet to the turbulence intensity in the zero-gravity jet on the jet center-line. ......... Scatter plot of temperature versus mixture fraction at :r = 15h, for (a) Ze=8,Fr=oo, (b) Ze=8,Fr=40. ................ Conditional average of temperature versus mixture fraction for Z 6 = 8 cases at x/h=15. ............................. Conditional average of reaction rate versus mixture fraction for the 3-D reacting jets at x/h=15. ......................... W(:c) = f wdydz along the jet axis and area of Z = 0.5 isosurface at t/to = 4, for (a) Ze = 8 cases, (b) Ze = 9 cases. ............ (a) Conditional average of F term in the interval 14 < 113/ h < 15, (b) Conditional average of G term in the interval 14 < :1: / h < 15. Axial variation of the correlation coefficient between the F and G terms for reacting jets. ............................. Comparison of 2—D and 3-D velocity fluctuations at x=15D, Ze=8, F 7' = 00,40. ............................... Comparison of 2-D and 3-D velocity fluctuations at x=15D, Ze=9, F r = 00,40. ............................... Comparison of 2-D and 3-D product integrals for Ze=9 cases. Schematic representation of the Sandia piloted methane jet flame. N oz- zle diameter = 7.2mm, Pilot diameter = 18.2mm. .......... Comparison between Monte Carlo (MC) and finite difference (F.D.) values of the normalized instantaneous temperatures for A3 = Ag, and Q; = 8, 16. .............................. Comparison between Monte Carlo (MC) and finite difference (FUD) values of the normalized instantaneous temperatures for AB = %Ag, and C,» = 8,16. .............................. Speed-up of the parallelized LES/FMDF calculations relative to the 2-processor calculation. ......................... xii 49 50 51 51 52 53 83 83 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 Temperature distribution vs. mixture fraction in a the low-strain (10 / s) laminar opposed jet flame. ....................... Experimental mean and variance of velocity profiles at the jet inlet for flame D. .................................. Isosurface of enstrophy Q = 1 from the flamelet, flame D case. The domain shown is 0 < x/d < 20, -—5 < y/d < 5, —5 < z/d < 5. Isosurface of (T) L = 1500K from the flamelet flame D case. The domain shown is 0 < ar/d < 20, —-5 < y/d < 5, —5 < z/d < 5. Isosurface of (T) L = 1500K from the 1—step flame F case. The domain shown is 0 < :zr/d < 20,—5 < y/d < 5,—5 < z/d < 5. ......... Isosurface of (T) L = 1500K from the 12-step flame D case. The domain shown is 0 < x/d < 20, —5 < y/d < 5, —5 < z/d < 5. ......... Mean axial values of the filtered velocity field for the 1-step and flamelet chemistry in flame D ............................ RMS of the filtered values of axial velocity for the 1-step and flamelet chemistry in flame D ............................ 86 86 88 88 Mean mixture fraction for the 1-step and flamelet chemistry in flame D. 89 RMS of mixture fraction for the 1-step and flamelet chemistry in flame Mean temperature for the 1-step and flamelet chemistry in flame D. . RMS of temperature for the 1-step and flamelet chemistry in flame D. Mean temperature for the 1-step and flamelet chemistry in flame F. . Mean mixture fraction at the jet centerline for the 1-step and flamelet models in flame D. ............................ Mean temperature at the jet centerline for the 1-step and flamelet models in flame D. ............................ RMS of temperature at the jet centerline for the 1-step and flamelet models in flame D. ............................ xiii 90 90 91 2.21 2.22 2.23 2.24 2.25 2.26 2.27 2.28 2.29 2.30 2.31 2.32 2.33 2.34 2.35 Profiles of mean CH4 mass fraction for the 1-step and flamelet models in flame D. ................................ Profiles of mean 02 mass fraction for the 1-step and flamelet models in flame D. ................................ Profiles of mean C02 mass fraction for the 1-step and flamelet models in flame D. ................................ Profiles of mean H20 mass fraction for the 1-step and flamelet models in flame D. ................................ Profiles of mean H2 mass fraction for the flamelet model in flame D. . Profiles of mean CO mass fraction for the flamelet model in flame D. Profiles of mean 0H mass fraction for the flamelet model in flame D. Profiles of mean heat release for the flamelet and 12-step models in flame D. .................................. Scatter plot of temperature vs. mixture fraction for flame D; experi- mental measurements at :13/ D = 7.5, and r/ D = 0.83 (975 samples). . Scatter plot of temperature vs. mixture fraction for flame D; experi- mental measurements at :1:/ D = 15, and T/ D = 0.83 (780 samples). Scatter plot of temperature vs. mixture fraction for flame F; experi- mental measurements at x/ D = 7.5, and T/ D = 0.83 (1550 samples). Scatter plot of temperature vs. mixture fraction for flame F; experi- mental measurements at x/ D = 15, and r/ D = 0.83 (976 samples). Scatter plot of temperature vs. mixture fraction for the 12—step flame D case, from :r/D = 7.5, and r/D = 0.83 (709 samples) ......... Scatter plot of temperature vs. mixture fraction for the 12-step flame D case, from x/D = 15, and r/D = 0.83 (635 samples). ........ Scatter plot of temperature vs. mixture fraction for the 1-step flame D case, from I/D = 15, and r/D = 0.83 (958 samples). ......... xiv 93 93 95 95 96 96 97 97 98 98 99 99 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 Schematic representation of the LES / F MDF simulations performed for buoyant/nonbuoyant methane jet flames. Nozzle diameter = 28.8mm, Pilot diameter = 72.8mm. ....................... 118 Isosurface of vorticity magnitude ((2 = 1) for case 1 (0-g , l-step). The domain shown isO<$/D<20,-5 0'02 + 2H20, in which the chemical reaction rate is obtained as d) = Mg”, p2YpYoexp(-Ta/ T), where B = 5.2 x 1010 (SI units) and Ta = 14906K. Equations (1.1)-(1.4) are integrated using the extended MacCormack method[22]. In addition, compact finite differencing schemes of Le1e[23] is used to approximate the derivatives, and non—reflecting boundary conditions were used both in the inflow[24] and in the outflow[25] boundaries. The longitudinal mean velocity 1'1 (overbar denotes time averaging) at the inflow boundary condition is prescribed by a hyperbolic tangent profile __U1+U2 Ul—Uz lyl"'h/2 u — 2 2 tanh ( 260 , (1.9) where 00 is the initial shear layer momentum thickness, U1 is the maximum jet velocity, and U2 is the coflow velocity. The momentum thickness, 00 is equal to 0.0625 in all simulations. A distribution similar to the one in equation 1.9 is used to describe the species mass fractions at the inlet, with the fuel mass fraction, Ya, decreasing from 1 at the center to 0 at upper and lower boundaries. The distribution is given as 1 1 IyI-h/Z a=——— —— . 1.10 Y. 2 2tanh< 200 ) ( ) In the 2-D simulations, the discrete harmonic perturbation of the 1) component of velocity is utilized to destabilize the jet. The u component of velocity was fixed. The frequency of the harmonic perturbation is the jet’s most unstable frequency[26] based on the initial momentum thickness. A weighting function, f (y), is used to concentrate the perturbations in the shear layers. In the 3-D simulations, homogeneous turbulence with a prescribed distribution, characteristic of those present in an actual turbulent flow is fed through the jet with the mean jet velocity, and all the velocity components are perturbed. The total intensity of the added inlet perturbation is five percent of the jet velocity. A weighting function, g(y), is used to concentrate the perturbation in the shear layers. It should be noted that buoyant instability would still occur if no perturbation is used, as was the case with the simulations of Shu et al.[15] and Davis et al.[13]. However, it is necessary to employ appropriate “turbulence” at the jet inlet to trigger the shear layer instability and to create a more realistic jet growth. 1.3 Results - Non Reacting Jet Three dimensional computations of a nonreacting planar jet is performed to ensure the accuracy of the numerical scheme, the boundary conditions, and the inlet flow con- ditions required to achieve a “fully developed” turbulence. The results are compared with existing experimental data on turbulent nonreacting planar jets of Gutmark and Wygnanski[19] as well as those obtained via DNS by Sarkar et al.[17]. In addi- tion, the overall accuracy of the scheme and the boundary conditions is assessed by examining the conservation of mass and momentum over the whole domain. The cor- responding Reynolds number and Mach number for the nonreacting jet simulations (and additionally for the 3—D reacting simulations) were 3000 and 0.3, respectively. The velocity difference ratio is AU/(Ul + U2) = 0.667, where AU == U1 — U2. The mesh size is uniform and equal to 0.05h. Moreover, the domain size in the 1:, y, and 2 directions were respectively 16h, 12h, and 4h. Figure 1.1 is a schematic representation of the domain and the coordinate system used in the 3D simulations. The resolution and the domain size are selected in such a way that all small and large flow scales are well resolved. ISO-surfaces and contours of vorticity magnitude as obtained by the 3-D nonreacting jet simulation are shown in Figs. 1.2 and 1.3. The vorticity vector is calculated from the instantaneous velocity field at the final time of the simulation (at t/to = 4), to = (2Lz/(U1 + U2) being the time required to traverse the domain length with the mean convective speed. It was seen that the first and second order moments do not change noticeably by increasing the simulation time. All the other reported results are also calculated at t/to = 4, unless otherwise noted. As the mean velocity statistics will indicate, a nearly “self- similar” status is reached at 113/ h z 4, were there is a sudden amplification of the vorticity (Fig. 1.2). A cross section of the 3—D vorticity field, in the :c — y plane at 2/ h = 2 (Fig. 1.3), indicates a break-up of the “continuous vortex sheet” to tube-like large scale vortices at x/h z 5. Further breakup of the vortex tubes to smaller scale turbulence occurs at :r/h a: 10, where the small-scale structures are clearly visible. Analysis of the vorticity field also reveals that the vorticity vector is mainly in the z direction before a: / h z 5. However, the three-dimensional effects gradually increase after 23/ h z 5. Also, after 17/ h z 10 directional dominance in the vorticity vector is no longer significant. Spreading of jets is customarily quantified through the jet half-width 61/2. Half- width is defined as 6 11(x,—1h/3,z) — U2 = % (fl(:1:,0,z) — U2). (1.11) Here, the calculated half-widths are also averaged along the cross-stream (z) direction. The growth of half-width along the axial direction (:5) is shown in Fig. 1.4. The jet thickness is almost constant before a: / h a: 6. However, after this point, the half—width increases at a constant rate. This is the expected behavior for a planar jet when self- similarity is reached[27]. A straight line is fitted to the self-similar section of the half-width plot as 51/2 77 = s. (5; + cl). (1.12) where SI 2 0.087 is slightly less than the average value of 0.1, observed in most planar jet experiments[19]. The value of virtual origin, CI, in our simulation is 0.23. There is a large scatter in the values of virtual origin in the literature. In the self-similar region of a planar jet, the center-line velocity, UC decreases as 23—1/2 , as expected[27]. In order to keep the axial momentum constant, the jet half-width has to grow linearly. The axial momentum flux in our nonreacting jet simulation varies within 2 percent of the initial momentum flux. Furthermore, our results (not shown) indicate that the center-line mean velocity excess along the a:- direction does not significantly change before x/ h z 6.0, but increases linearly after 33/ h z 6.0, which is consistent with the linear growth of the half-width in Fig. 1.4. A straight line is fitted to the linear section of the center-line velocity excess, as AUO=U0-’U2=S2( 1' AU. U. — U2 ’ + 02): (1‘13) h where 52 = 0.24, and C2 = —3.4. The slope 52 is slightly higher than its correspond- ing experimental values. As a result of the jet growth, the jet draws fluid from the surroundings and the mass flux increases accordingly in the axial direction. The change in axial mass flux, Amx is calculated as 'dydz A3,: I" —1=/end, 1.14 m (ffldydzfizo m ,3; ( ) where men, is the vertical mass flux density on the lateral boundaries, which is the source of the increased mass flux in the axial direction. As shown in Fig. 1.5, the integral of the vertical mass flux density matches the axial mass flux increase profile, which is calculated independently. This shows the code’s mass conservation accuracy over long periods of time. The entrainment mam, reaches its peak at :r/h z 8, consistent with the experimental observations of Hussain and Clark[28]. The decline in entrainment is attributed to the breakup of large vortices and the development of the fine-scale turbulence. Self-similarity in the mean velocity profiles is observed at :r/ h > 5, where the initial top—hat velocity profile has been transformed to the profile shown in Fig. 1.6 (a). After :1: / h = 5, the normalized mean axial velocity at all axial locations collapse to those shown in Fig. 1.6 (a). Fhrthermore, there is an almost perfect matching with the experimental data of Gutmark and Wygnanski[19] (taken from the self-similar part) and Ramaprian and Chandrasekhara[29], despite the fact that self-similarity in the experiments is reached at a different axial location. While the first order statistics becomes almost self-similar at 22/ h z 5, self- similarity for the second order statistics is not achieved until farther downstream distances. In our simulation, the trend in the second order statistics suggests an ap- proach to self-similarity, only after :1: / h a: 12. Other DNS and experimental results also indicate that the second order statistics reach a self-similar state much later than the first order statistics. For example, Gutmark and Wygnanski[19] observe that the W profiles reach a self-similar profile at 123/h = 30, while self-similarity for W and W is achieved even farther downstream. In the case of DNS, this distance is short- ened by adding a ‘turbulent” perturbation at the jet inlet. The added turbulence has the spectrum of a developed, solenoidal isotropic field. Figure 1.6 (b) shows the axial velocity fluctuations as obtained by our nonreacting simulation from :1: / h = 15, which is compared to the experimental results of Gutmark 10 and Wygnanski[19], Ramaprian and Chandrasekhara[29], as well as the DNS results of Stanley and Sarkar[17]. Similar comparisons have been made for the v’—2/ and W components in Figs. 1.7 (a) and 1.7 (b), respectively. Evidently, the results obtained from various studies compare very well, despite the differences in the velocity ratio and inlet flow conditions. However, since the flow parameters in the present study are relatively close to the parameters in the simulations of Stanley and Sarkar[17], the results compare best with their results. 1.4 Results - Reacting Jet The results obtained from the simulated 2-D and 3-D jet flames are described in this section. The 3-D-reacting simulations have been performed using the same resolution, perturbations and domain size used for the nonreacting jet. Most of the emphasis would be on the 3-D results, as 2-D simulations do not completely represent the evolution of the turbulent jet flows. However, certain as— pects of a reactive buoyant jet can still be captured by the 2-D simulations, as will be demonstrated later in this paper (in section 1.4.1). Since the 2-D simulations are computationally less expensive, they were used to investigate “near-field” character- istics of the buoyant jet and to obtain an estimate of the flow parameters. Moreover, the accuracy of the 2—D simulations in modeling turbulent reactive jets is investigated. The flow parameters for the 2-D and 3-D simulations are listed in Table 1.6 and Table 1.6, respectively. The axial domain size in the 2-D simulations is 24h, as opposed to 16h in the 3-D simulations. The combustion parameters are changed by varying the Damkohler number, Da, the Zeldovich number, Ze, and the heat release parameter, Ce. The 2-D methane jet flame simulations are not included in tables 1.6 and 1.6, but their jet height and velocity was 3.84cm and 10.98m/s respectively, and fixed chemistry parameters were used. 11 1.4.1 Two-Dimensional Simulations In this section, the main findings from the 2-D simulations are briefly discussed. First, we discuss the variations caused by changing the Damkéhler number, which effectively changes the relative reaction rate, while keeping other parameters constant. In the near field of the jet, the Damk6h1er number effects are somewhat similar to the Zeldovich number effects. The overall response of the flame to the Damkdhler number is shown in Fig. 1.8, where the total kinetic energy and total product for the F r = 40 buoyant cases with different Damkohler numbers (cases 2,6, and 7) are considered. Since the Zeldovich number for these cases is equal, the reaction would start earlier in the higher Damkohler number jet flame. While the product integral for the lowest Damkéihler number flame starts to rise after the other two cases, it will eventually surpass the values for those. Also, at the lowest Damkohler number, highest values of the mean temperature are observed. At a: = 15h, the peak of time-averaged temperature in the Da = 60 case is about 2.6, while in the other two cases, it is about 2.1. This is an unexpected result, since combustion is normally less effective in lower Damkohler numbers. In the Da = 60 case, the temperature rises much more slowly than the two other cases, which results in the jet being less affected by combustion. This slower temperature rise provides better initial mixing between the fuel and the oxidizer, creating a more intense reaction and higher temperature at the downstream locations. The turbulent kinetic energy in the Da = 60 case also exceeds the values for the other cases after 23/ h = 12, and retains a higher growth rate. There seems to be a direct relation between the total turbulent kinetic energy and the total product formation in this case. The effect of gravity on the reaction is studied by varying the Froude number, while keeping the chemistry and flow parameters constant. The non-buoyant jet has a laminar-like structure with the reaction taking place mainly at two distinct and parallel mixing layers. This case has the highest values of the local temperature. 12 However, by increasing the gravity (reducing the Froude number), flame oscillations and flickering increases considerably. The unsteady nature of the 2-D large scale structures and the flame flickering is shown in Fig. 1.9, where the temperature con- tours obtained from the non-buoyant (F r = 00) and buoyant (F r = 20) jets are compared. The overall performance of these cases can also be compared in terms of their product integrals. Figure 1.10 shows that after :1: / h = 17, the product formation in the buoyant jet increases sharply, while the difference between the F r = 40 and F 7' = 20 cases remains small. The total kinetic energy is plotted on the same figure, and shows its direct correlation with gravity. The higher the gravity level, the more significant the turbulent fluctuations become. In the 2-D cases, as seen in table 1.6, two different coflow velocities (0.2 and 0.5 of the jet velocity) have been used to investigate the effect of shear on the buoyant jet flames. In the cases with the coflow velocity being 0.5 of the jet velocity, the combined effects of the shear and the buoyancy are not sufficient to significantly disturb the jet. The jet is very similar to a laminar jet and the fluctuation intensities are an order of magnitude smaller compared to the cases with U2 = 0.2U1. The effect of higher Reynolds number has also been studied by comparing buoyant jet flames with Re = 3000 and Re = 6000. The time averaged values of the velocity, the turbulent stresses and the total product formation indicate that the Reynolds number effects are not significant in the 2-D simulations, for the range of parameters considered. A more realistic jet flame has also been simulated, using diluted methane as the fuel and air as the oxidizer. Realistic thermodynamic properties, such as the specific heat are used. A one-step global mechanism is considered [21], considering the cost of more detailed kinetics. Even with a global mechanism, the spatial and temporal resolution requirements for buoyant hydrocarbon jet flames are very stringent. In our methane jet simulations, a grid of 960 x 600 was used corresponding to a do— 13 main of 18h X 12h. The resulting grid spacing is 0.02h compared to the value of 0.05h for the 3-D simulations. The reaction is started by preheating the jet stream (%5OCH4, %50N2 by volume) to To = 1100K. The coflow has the standard air com- position and temperature (298K). The reaction does not start immediately after exiting from the nozzle, similar to the simulations performed at lower Damk6hler numbers. The Reynolds number based on jet conditions and air properties at the reference temperature (1100K) is 3000. The Froude number is 80, which is twice the value of the other simulations. However, the initial density difference at the nozzle exit is more significant. The effective Zeldovich number is 13.5, but the main differ— ence between this simulation and the previous simulations is the magnitude of the Damkéhler number, which is very large (5.7 x 108). The Damk6hler number con- sidered in the other 2-D and 3-D simulations is much smaller. Nevertheless, these cases exhibit some of the characteristics of hydrocarbon flames at higher Reynolds numbers, such as local extinction. The dimensional jet height and velocity in the methane jet flame simulations are 3.84cm and 10.98m/s, respectively. Due to the limitations of calculations, the grav— itational acceleration was increased to 4 times its standard value, in order to reduce the Houde number and to amplify the buoyant effects at the jet height and velocity considered. Experiments of methane diffusion flames at different gravity levels (1- 59) were performed by Sato, Amagai, and Arai [30, 31]. The global characteristics of the flames (flame length, flickering frequency and shape) were investigated and their dependence on gravity was analyzed. A flickering frequency of about 25 was reported at 49, while the flickering frequency measured in our simulations was about 28. A separate simulation without any perturbations was performed to obtain this frequency. The flame flickering observed without the use of perturbations is symmet- rical in shape, consistent with the low Reynolds number experiments of Davis et al. [13] and Sato et al. [31]. Employing external perturbations would alter the symmetry 14 and frequency of the flickering. The contour plots of temperature for the buoyant and non—buoyant methane jet flames are shown in Fig. 1.11. Evidently, the flame has very distinct boundaries and the inner core of the jet has an almost uniform temperature, because of the reaction speed. Also, the flame flickering is clearly visible in the buoyant case. However, the non-buoyant flame is much wider than the buoyant flame, consistent with experiments of Bahadori et al. [9]. Acceleration of hot products in the buoyant jet results in narrowing of the jet flame. Jet half-width is calculated and shown for the methane jet simulations in Fig. 1.12. The half-width is larger in the buoyant case, consistent with the 3-D results shown in Fig. 1.21. The larger magnitude of half-width for the buoyant jet is due to the flickering motion, despite the wider shape seen in the non-buoyant jet. The growth pattern is consistent with the 3-D results at lower Damkéhler numbers; an initial period of slow growth followed by a more steep growth and a final stage of almost constant but lower growth rate. In the case of the methane flames, the growth is very sudden due to more intense reaction. The half-width drops slightly in the early stages of the buoyant jet, because of the initial acceleration of the jet due its density difference from the coflow. The integrated mass fraction of 002 over the planes normal to the flow direction are shown in Fig. 1.13. In the non-buoyant jet, the initial growth is faster than that in the buoyant jet, since the molecular mixing between the fuel and oxidizer streams is more effective in the slower jet. Analysis of reaction rates at locations after this (:1: / h = 4), show a stronger reaction for the buoyant jet, consistent with the observation of reduced flame lengths at normal gravity, seen in experimental studies [9]. Consistent with the simulations of Katta and Roquemore [32], the scatter plots of temperature and species concentration (not shown), do not reflect a significant change 15 with respect to gravity. However, the magnitude and the amount of scatter in the reaction rate data is significantly larger for the buoyant case (Fig. 1.14), suggesting an increase in mixing and turbulence by gravity. This is consistent with the general behavior seen in our other 2-D and 3—D simulations (see Fig. 1.30 for example). As will be shown below, depending on the level of turbulence, the gravity usually induces higher strain rates and increases the chance for local flame extinction. However, the overall effect of gravity on the flame is very much dependent on the turbulence and the relative flame speed. This issue would be thoroughly investigated in the next section, where we discuss the 3-D simulations, since the 2-D simulations do not fully model the turbulence development in our jet. 1.4.2 Three-Dimensional Simulations The conditions for the 3-D reacting simulations (table 1.6) are somewhat similar to the 3-D nonreacting simulations; since the same perturbations, mesh size, Reynolds number, and velocity ratio have been used for both cases. The location of the peak reaction rate from the nozzle is changed by varying the Zeldovich number. The development of turbulence is affected indirectly by changing the onset of the flame. However, the effect of having different Zeldovich numbers on the reaction rate is compensated by changing the Damk6hler number, such that the same reaction rate is achieved at the adiabatic flame temperature, which is kept constant in turn by keeping the heat release parameter, Ce constant. To investigate the effect of buoyancy (or lack of it), two different levels of gravity are considered. The corresponding Froude numbers are Fr = 00 and Fr = 40 for these cases, which will be referred to as the zero-gravity and the finite-gravity (buoyant) cases, respectively. A higher gravity level (with F 7‘ = 20) is also considered in some of the 2-D simulations. The direction of the gravity force is opposite to the jet flow direction in all simulations. The analysis of the 3-D results are presented in the following order: contour plots 16 and the physical flow structure, the time-averaged statistics, and the compositional structure. Physical flow structure Instantaneous vorticity magnitude contours for the 3-D non-buoyant reacting case 12 are shown in Fig. 1.15. Vorticity contours do not exhibit the vortex sheet roll-ups and strong three-dimensionality seen in Figs. 1.2 and 1.3 for the nonreacting jets, and moving downstream along the jet axis, the vorticity magnitude is actually decreased. A comparison of Figs. 1.15 and 1.16 reveals the overall differences between zero-gravity and finite-gravity jet flames. In the finite-gravity or buoyant flame, unlike the zero- gravity one, the vorticity magnitude increases along the jet and is generally higher as compared to the zero-gravity jet. Peak of the vorticity magnitude is around 4 for the buoyant jet, compared to 3 for the non-buoyant jet. Furthermore, in the finite-gravity jet, there is a strong vortex sheet roll-up at x/ h z 10 and subsequent formation of rather complex 2-D (and to some extent 3-D) vortex structures in contrast to the zero-gravity flame. The larger structures have cross-stream movements that are known as flame flickering. In the 3—D cases 12 and 13, the reaction mainly takes place in the mixing layer zone, and the temperature field coincides with the vorticity field as Fig. 1.15 suggests. In addition, the “temperature layer” (narrow zones of high temperature) is visibly thinner in the buoyant case, due to better mixing and gravity induced acceleration of the high-temperature regions. However, the peak of temperature in the non-buoyant jet, is higher than the peak of temperature in the buoyant jet, which is due to the higher strain rates, more intense convection and stronger finite-rate chemistry effects in the buoyant jet as demonstrated in the next section. The results for cases 14 and 15 are shown in Figs. 1.17 and 1.18, respectively. In these cases, Z e = 9 and the flame / turbulence structures are very different from those 17 shown in Figs. 1.15 and 1.16 for the Ze = 8 cases. In cases with higher Zeldovich number, the maximum reaction rate and heat release occurs farther downstream; allowing the development of vorticity field and partial premixing of the fuel and oxidizer at regions close to the nozzle. As a result, turbulent mixing is more significant in higher Ze cases. Figure 1.17 shows that in the zero-gravity condition, there is a significant growth in the jet and the large/ small scale turbulent structures. However, in cases with Z e = 9, the jet growth rate and instability are not significantly affected by the gravity, (compare Figs. 1.17 and 1.18) unlike the effect seen in the Z 6 = 8 cases. A comparison between cases with Z e = 8 and Ze = 9 at zero-gravity, indicates that in the Z 8 = 8 case, the early increase in temperature stabilizes the shear layer to some extent and clamps the vorticity generation. However, the temperature rise in the Z e = 9 case occurs at locations farther downstream (as will be quantitatively shown in the next section), where the turbulence is more developed compared to the Z e = 8 case. Consequently, in cases with Z e = 9, gravity does not play a leading role in the jet instability as seen in lower Zeldovich number cases, and it acts rather like a “random perturbation function” in the momentum equation. In these cases, the effect of gravity is to accelerate the already “turbulent” jet and to create even higher strain rates, which stretches the reaction regions and makes the high temperature layers even thinner. This argument will be supported by the results presented below in the statistics section. Fig. 1.19 shows the vorticity magnitude contours for buoyant and non-buoyant jet flames at a plane perpendicular to the :1: axis at x/ h = 15 for cases with Z e = 8. Three-dimensional effects exist in both cases, however they are more prominent in the finite-gravity case, and the vorticity magnitude is clearly higher as compared to values in the zero-gravity condition. There is also a substantial difference between turbulent structures in zero- and finite-gravity flames. 18 Statistical description In this section, we discuss the effects of gravity on various turbulent statistics via detailed analysis of our 3-D data. Figure 1.20 shows the time averaged values of the cross-stream temperature at :r/ h = 5 and :r/ h = 15 in the zero-gravity (F r = 00) and finite-gravity (F r = 40) environments. The results for Z e = 8 and Ze = 9 are presented in Figs. 1.20(a) and 1.20(b), respectively. Due to the earlier reaction in the Z e = 8 case, the peak temperature is much higher than that in the case with Z 6 = 9 at a: / h = 5. Consequently, the buoyancy has a more significant effect on velocity and temperature fields in cases with lower Zeldovich number. Figure 1.20 shows that for both Z e = 8, Z e = 9 cases, the peak temperature at x/ h = 15 decreases by gravity. However, this is accompanied by the temperature profiles getting broader by gravity, specially in the Z e = 8 case; suggesting that the extent of combustion and the overall rate of product formation are not necessarily reduced by the gravity, as will be shown later (see Fig. 1.22). Another important difference between the temperature plots in Figs. 1.20(a) and (b) is that while gravity shifts the peak temperature away from the centerline in the case with Z e = 8, it has the opposite effect on the case with Z e = 9. As a result, the maximum mean temperature is closer to the jet centerline in the buoyant jet that has higher Zeldovich number. In the Z e = 9 cases, the reaction develops at a relatively long distance from the nozzle as shown in Figs. 1.18 and 1.20(b). However, the peak temperature in the zero-gravity jet with Z e = 9 eventually reaches the same peak as the zero—gravity jet with Z e = 8. The centerline temperature is higher for the higher Zeldovich number due to stronger mixing in the jet core. In the case with Z e = 8, the buoyancy-induced instability plays a dominant role in the jet growth rate; leading to large scale variations in temperature and more spread in the mean (time-averaged) temperature field and also a significant decrease in the peak of mean temperature. In the case with Z e = 9, the buoyancy affects an already 19 developed turbulent field at all scales; leading to a very complicated effect on the ‘ physical/ compositional structure and the statistical behavior of the jet flame. In this case, the peak of mean temperature shifts towards the centerline despite the fact that the average jet growth rate increases by gravity (see Fig. 1.21). This could be due to the enhanced turbulent mixing and partial premixing of the fuel and oxidizer at the jet core that leads to flow acceleration at the jet centerline. Figure 1.21 shows the growth of half-width in all the simulated 3-D reacting jets. In all cases, the half-width grows almost linearly after x/ h z 7, with the rate of growth being dependent on the magnitudes of Z 6 and F 1‘ numbers. The non-buoyant case with Ze = 8, (case 12) has the lowest jet growth rate, which is consistent with the results in Fig. 1.15; showing a laminar-like vorticity field. The buoyant case with Z e = 8, (case 13) and the non-buoyant one with Ze = 9, (case 14) exhibit almost the same growth rate in half-width. Moreover, the buoyant case with Ze = 9, (case 15), has the highest growth rate. This indicates that the gravity increases the jet growth rate in all cases. However, its effect on turbulent combustion is complex and difficult to predict. This is illustrated in Fig. 1.22, where the axial variation of the product thickness or total product for various cases is shown. While in the cases with Z e = 8 the difference between the zero-gravity and the finite-gravity results are less significant than those in cases with Z e = 9, the gravity seems to decrease the total product at as/h < 15. In the non-buoyant jet with Ze = 8 (case 12), the product formation rate is initially lower, but the predicted average product in this case becomes nearly equal to that of the finite-gravity jet at the end of the computational domain. Moreover, as the results of our 2-D simulations suggest, the predicted values for the finite-gravity jet would exceed the zero-gravity values, if a longer domain is considered. We will explain these results further in the composition section, where the effects of gravity on the compositional flame structure is discussed. 20 The time averaged profiles of the normalized axial velocity at :r/h == 5, 15 for the flames with Ze = 8 are shown in Fig. 1.23. The results for the Ze = 9 cases are similar and are not shown. The mean velocity profiles are not self-similar in either case and as expected the buoyant jets show a slightly wider profile. The variations in the center-line velocity (in terms of the center-line velocity excess) are shown in Fig. 1.24. The non-buoyant jet flames exhibit steady increase in the centerline velocity difference after x/ h > 5. In case 12 (non-buoyant with Z e = 8), the growth rate is linear but the slope is much less than that in the nonreacting jet. Also in case 14 (non-buoyant with Z e = 9), the center-line velocity is almost constant up to a: / h z 7, with a subsequent linear decline between 7 < 33/ h < 10 and a milder decline that continues till x/ h = 15. The intermittent variations in this case are attributed to the nearly isothermal early development phase, followed by a sudden heat release and volumetric expansion phase, and ending with a well developed turbulent phase. Evidently, the evolution of the centerline velocity and the jet growth rates are very much dependent on the combustion parameters. Nevertheless, gravity always tends to diminish the growth in velocity excess due to acceleration of hot combustion products in the jet centerline. The gravity effects appear stronger in case 14 with higher Zeldovich number. The axial variation of the cross-stream velocity, 1'), on the lower boundary (y/ h = —6) of the reacting jets are shown in Fig. 1.25. Negative (outflow) values are caused by the volumetric expansion of the jet due to heat release. In nonreacting jets, the well known entrainment process draws “co-flow” fluid to the jet’s core, and the aver- age cross-stream velocity at the boundary is positive. However, the dominant effect in our non-buoyant reacting jet simulations is the volumetric expansion and 2‘); always takes negative values. In the buoyant jet simulations, stronger turbulence/vorticity production counteracts the effect of volumetric expansion. Consequently, the en- trainment rate is almost zero for the buoyant jets at :r/ h < 4, and slightly positive 21 at 4 < III/h. < 12. After :c/ h = 12, 37, changes sign to negative values due to the diminished vorticity field. The profiles of the normalized longitudinal velocity fluctuations are shown in Fig. 1.26. The shape of the profiles are similar to those for the nonreacting jet with a peak at y / h z 0.8. However, the center-line values are lower in the reacting jets, due to the diminishing effect of combustion on turbulence. While the longitudinal turbulence intensity substantially increases by the gravity in cases with Z e = 8, the difference between finite-gravity and zero—gravity cases is not noticeable when Ze = 9. The results in Fig. 1.26 are consistent with the earlier results shown in the physical flow structure section for Z e = 9 cases. The effect of gravity on the turbulence intensity is illustrated in Fig. 1.27, where the ratio of turbulent intensity of the finite—gravity jet to that of the zero—gravity jet is shown. The turbulent intensity ratio, K E, is defined as I I I I I I KEl/Z = (uu +vv +ww)”:40 7' (quI + Div! + wlwl) ° Fr=oo Figure 1.27 shows the axial variations of K E:/ 2 at the jet centerline for different Zeldovich numbers. For the flame with Z e = 8, the turbulence and E, is enhanced significantly by the gravity. However, the turbulence intensity ratio is much smaller and mostly less than unity in cases with Ze = 9. Despite the fact that this ratio is less than unity at the jet centerline, the average values over planes perpendicular to the jet are larger than unity in all cases. This suggests that the gravity always enhances turbulence. Nevertheless, gravity does have a larger impact on the turbulent energy in under-developed lower Zeldovich number cases considered here. The effect on combustion is less clear, and requires a detailed analysis of the compositional structure of the flames. 22 Compositional flame structure The effects of gravity on the compositional structure of the jet flames are investigated here by looking at the variations of different turbulent quantities in the mixture fraction, Z domain. The scatter plots of temperature for non-buoyant and buoyant jets are shown in Figs. 1.28 (a) and 1.28 (b), respectively. The data are collected from the y-z plane perpendicular to the jet axis at 23/ h = 15. The gravity clearly enhances the finite- rate chemistry effects in the flame as seen by the considerable scatter in Fig. 1.28 (b). The increased scatter is due to the higher strain rates generated by the buoyancy. Our results (not shown) indicate that the gravity indeed increases the strain rates significantly, specially in the Z 6 = 8 case. The calculated temperature between 14 < :r/ h < 15 are conditionally averaged with respect to mixture fraction (< TIZ >) and plotted in Fig. 1.29. The plotted conditional average values (denoted as < TIZ >) are consistent with the scatter plots in Figs. 1.28 (a) and 1.28 (b) and demonstrate that the gravity reduces the peak of the mean conditional temperature profiles without significantly changing the shape of the profiles. The conditional average values of the reaction rate versus mixture fraction ((u'IIZ)) at :r/h = 15 are shown in Fig. 1.30 for all the 3-D reacting jet simulations. At this location from the nozzle, the reaction rate is clearly higher in the buoyant jets. However, to understand the gravity effects on combustion, one must look at the dynamical evolution of the flame and the history of the reaction. To investigate the dynamics of the jet flames and to study the structure of the various reacting jets Considered here, the total reaction rate on planes normal to the :1: axis are computed for all reacting cases. The total reaction rate, W(:I:) 2: frbdydz, (1.15) 23 is defined such that the integral of W(:r) along the :1: axis would indicate the total product formation in the jet. In addition, for each Zeldovich number considered, the area per unit length of the stoichiometric surface (Z = 0.5 isocontours) is computed and is plotted together with W(:r) along the :1: direction in Figs. 1.31(a) and 1.31(b). For the Ze = 8 jets, W(:r) peaks at ar/h z 4.5 (Fig. 1.31(a)), while for the Ze = 9 jets, the peak values are located much farther downstream at :c / h z 8 (Fig. 1.31(b)). Better turbulence and molecular mixing is achieved due to the delayed reaction in the Z e = 9 jets; resulting in a higher peak for W(z) in this case. However, gravity slightly reduces the peak of W(:r) regardless of the magnitude of the Zeldovich number. Due to poor mixing of the reactants in the zero-gravity jet with Z e = 8, W(a:) has an overall decline after it reaches its peak values around 115/ h = 8. In this case, the nearly constant stoichiometric surface area values in Fig. 1.31(a) is consistent with the “laminar-like” structure of the zero-gravity jet in Fig. 1.9. In contrast to the zero-gravity jet, the stoichiometric surface is highly fluctuating in the finite-gravity jet, although it increases on average. The W(2:) term varies similarly and its local maxima coincides with the local minima of the stoichiometric surface area term. This is contradictory to the expectation that the reaction rate is proportional to the reaction surface area. The difference between the W(:I:) profiles in the buoyant and non-buoyant Z e = 9 jets is smaller than the corresponding difference in the Ze = 8 jets (Fig. 1.31(a)). Additionally, the finite-gravity jet exhibits lower values of W(:1:) and the stoichio— metric surface area of the zero-gravity Z e -_.-_- 9 jet fluctuates with higher frequency and smaller amplitude (Fig. 1.31(a)) in comparison to the finite-gravity Z e = 8 jet (Fig. 1.31 (a)). Similar to the Z e = 8 cases, the finite-gravity Ze = 9 jet exhibits larger stoichiometric surface values in comparison with the zero-gravity jet in cases with Ze = 9. However, there is an increase in stoichiometric surface area farther downstream after x/ h a: 12. 24 To better understand the effects of gravity on the reaction, the constituents of the reaction rate expression are separately evaluated and analyzed here by decomposing the reaction rate as proposed by Jaberi et al. [4] in the following form: (12 = DaszAYB exp (—Ze/T) = F x G. (1.16) F G The conditional average of F and G terms are calculated from the data and plotted in Figs. 1.32(a) and 1.32(b) for the interval 14 < az/h < 15. The F term is higher in case 14 (the zero—gravity case with Ze = 9) compared to that in case 12 (zero- gravity case with Z e = 8), which shows a better mixing in the former case. In the Ze = 8 case, the addition of gravity yields a more uniformly distributed F term with considerably higher values. Similarly, the conditional average values of the F term increases by gravity in the Z e = 9 case. However, in this case, the distribution still has a distinct peak at the rich side of the mixture fraction space. The G term follows the trend seen in the conditional averaged temperature distributions. Due to the finite-rate chemistry effects and more significant temperature fluctuations in the higher Zeldovich number jet, the non-buoyant jet shows lower G values as compared to that with lower Zeldovich number. However, the G term values in the buoyant jets are lower due to considerable scatter in the temperature. As Eq.1.16 suggests, in addition to the magnitudes of F and G terms, the corre- lation between these two terms is important. Figure 1.33 shows that the correlation coefficient between the F and G terms is initially high in all cases. However, it quickly drops to very low values in the Z e = 8 case and also soon afterwards in the Z e = 9 case. This suggests the rapid consumption of the reactants and the separation of the high-temperature zones from the “mixing” (F) zones. The finite-gravity Ze = 8 jet has the highest correlation coefficient at far field (:13/ h > 10). This is the very same case that demonstrated an almost uniform distribution of the F term in the mixture 25 fraction space in Fig. 1.32(a). 1.4.3 Comparison Between 2—D and 3-D Simulations Two—dimensional simulations of the jet flames with the same domain size and pertur- bation as the 3-D jets are performed, in order to have a better understanding of the limitations of our 2-D simulations. The 2-D jets seem to have a completely different behavior compared to the 3-D jets at the far-field. Since there is no mechanism for the break-up of the large-scale vortices in the 2-D simulations, large-scale flame flickering is very dominant in the far-field in these simulations. The normalized rms of the velocity fluctuations as obtained by 2-D and 3-D sim- ulations for cases with Z e = 8 are compared in Fig. 1.34. The profiles look similar in the non-buoyant cases, with the velocity fluctuations being stronger in the 3-D case as expected. There is much less difference in the buoyant cases in terms of intensity. However, the profile in the 2-D case does not have a distinct peak, due to the more significant large scale flickering. The turbulent fluctuations from the 2-D simulations with harmonic perturbation (not shown) exhibit similar behavior, while the magni- tudes are larger since the perturbation frequencies match the natural jet instability modes. However, the overall behavior is similar for different perturbations. The rms of the velocity fluctutions for cases with Z e = 9 are shown in Fig. 1.35. There is a rel- atively significant difference between the non—buoyant and buoyant results in 2-D jets when compared to the corresponding results in the 3-D jets. This suggests that the large scale flickering is less pronounced in the jets with more developed and realistic (3-D), small-scale turbulence. This indifference to gravity was seen in the reaction rate profiles for Ze = 9 flames as well (Fig. 1.31(b)). This finding is consistent with the analysis of Bahadori et al. [9], where it was suggested that at sufficiently high Reynolds numbers, the microgravity flame will possibly have characteristics identical to a fully developed turbulent flame in normal gravity. On the other hand, at low 26 Reynolds numbers, where the turbulence is less developed and the flow is almost lam- inar, 2—D simulations are expected to provide results similar to the 3-D simulations or experimental results. For example, the 2—D (axisymmetric) simulations performed by Azzoni et al. [11] and Katta and Roquemore [32] have compared well with their cor- responding experiments, where in both cases the flow Reynolds number was relatively small. The zero-gravity flames at this level of Reynolds numbers are almost laminar as also shown by the experiments of Hegde et al.. While the normal gravity flames become unstable even at lower Reynolds number, they are still far from developed turbulence. Our simulations indicate that the product integral profiles for the 2-D, Ze = 8 cases, where very close to their 3-D counterparts shown in Fig. 1.22. This again suggests that the instabilities in the Z e = 8 flame are of large-scale nature. For the Ze = 9 cases, the product integrals for the 3-D and 2-D flames are shown in Fig. 1.36. The results in this figure reflect the more developed nature of turbulence in the Z e = 9 cases and the fact that the smaller scale turbulence is noticably affected by buoyancy. In this case and for the flames with similar characteristics, a 3-D simulation would be necessary for accurate representation of the flow and combustion. 1.5 Summary and Conclusions Direct numerical simulations (DNS) of various 3-D and 2-D nonreacting and reacting diffusion planar jet flames are performed with a one-step irreversible chemical kinetics model. The Damkohler number considered in these simulations was lower than the typical values in hydrocarbon flames, in order to capture flame-turbulence interac- tions and the finite-rate chemistry effects in zero- and finite-gravity levels, via DNS. Buoyant and non-buoyant 2—D methane flame simulations with realistic thermody- namic properties and a global one-step kinetics model are also performed. Results 27 from the nonreacting 3-D jet simulation are compared to the available experimen- tal and DNS data in the literature. Various statistics such as jet growth rate and Reynolds stresses are shown to be in good agreement with experimental and previous numerical results. The results of our reacting jet simulations are also consistent with previous find- ings and indicate that in the absence of gravity, combustion damps the flow instabil- ity; hence reduces “turbulence production” and jet growth. However, in the “finite— gravity” conditions, combustion generated density variations and buoyancy effects, promote vorticity generation and enhance the turbulent mixing and combustion at nearly all length and time scales. As expected, buoyancy was found to have a desta- bilizing effect on laminar jets, and the formation of large-scale flame flickering was consistent with previous investigations. However, the effect on transitional/ turbulent jet flames is not limited to large-scales, and the turbulence and mixing were affected in the small scales as well. Comparison between the 2-D and 3-D simulations suggests that the effects of com- bustion and gravity on turbulence (and vice versa) cannot be appropriately evaluated by 2—D computational models when the turbulence is more developed. However, the early effects during transition are adequately captured by the 2-D simulations. The gravity effects are also explained via the analysis of the compositional struc- ture of the flame. It is shown that the finite-rate chemistry effects become more significant in the buoyant cases due to the increased turbulence and the higher strain rates. When the turbulence field is damped by the heat release, the addition of grav- ity could help the overall production by destabilizing the flow, resulting in improved mixing and an extended flame surface. On the other hand, when the turbulence field is “strong”, the addition of gravity may induce excessive amounts of strain and cooling of the flame, which in turn adversely affects the combustion. 28 1.6 Figures and Tables Table 1.1: Flow parameters for 2-D simulations. Case Da Fr U2 COCONGO'IIbDONJI—II )—| O 120 120 120 0 120 200 60 120 120 120 20 40 40 40 40 40 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.5 0.5 0.3 except for case 5 which has Re = 6000. Table 1.2: Flow parameters for the 3-D Simulations. Ce = 2.5 for all cases. Re = 3000, Ge = 2.5 for all cases Case Ze Fr Da 11 - oo 0 12 8 00 120 13 8 40 120 14 9 00 160 15 9 40 160 29 Re = 3000, U2 = 0.2 and a v ; Figure 1.1: Schematic representation of the domain and the coordinate system used in the 3D simulations. The flow is in the upward direction. 30 Figure 1.2: ISO—surfaces of vorticity in the 3—D nonreacting jet at the final time of the simulation, t/to = 4. The domain shown is 1 < :r/h < 16, ——4 < y/h < 4,0 < z/h < 4. 31 Figure 1.3: Instantaneous vorticity magnitude contours (0—10) in the 3-D nonreacting planar jet, at z/h = 2,0 < x/h <16,—-4 < y/h < 4. 32 1.0- dug/h. 0.8 - 10 15 Figure 1.4: Half-width of the 3-D nonreacting jet along the jet. A straight line, 0.087 (% + 0.23) is fitted to the self-similar portion of the jet growth rate. . , . , . 0.12 H "‘1 s E 0.08 5 E” <5 E E x 0.04 53 0o 5 1 10 15 z/h Figure 1.5: Mass entrainment from the lateral boundaries and increase in the axial mass flux. 33 I I ' I ' I ' I . I . . (a) . 1 — _. . l 0.8 - _ S d 0.6 — -‘ \ b - _ <1 0.4 - - - U Gutmark - o Ramaprian 0'2 _ --- Stanley _ - — Our DNS ~ 0 a D D O l - J I L l . l r I} r -3 -2 -1 O 1 2 3 9/51/2 f l ' l r l ' l ' l ' 0.3 b (b) o°°°°°go ._ 0° 0 o S 0.2 '- -—I <1 > x; - F :3 0'1 T D Gutmark m f o Ramaprian --- Stanley — Our DNS 0-3 -2 -l O l 2 3 9/51/2 Figure 1.6: Comparison of (a) mean axial velocity, (b) longitudinal Reynolds stress, with experimental and DNS data. Our data is from :r/ h = 15. 34 0.2 - v’v’/AUC 0.1- D Gutmark '3‘ o Ramaprian my --- Stanley 0 — Our DNS °° 1 l 1 l l l l l . l 1 0-3 -2 -1 0 1 2 3 9/51/2 0.02 0.01 <1 0 \\ ‘3 :3 -0.01 -0.02 0 Gutmark - Our DNS _ . l . l . l . J . 1 . 0'03-3 -2 -1 0 1 2 3 9/51/2 Figure 1.7: Comparison of (a) lateral Reynolds stress and (b) Reynolds shear stress profiles with experimental and DNS data. Our data is from 113/ h = 15. 35 — Da=120 /’ 0.3 .9 N deA/ fdA .O y... f mA/ f «M 1 . L . 0 6 12 18 24 :r/h Figure 1.8: Comparison of total “product Integrals” (lower curves) and total kinetic energy (upper curves) for different Da numbers. 36 :.. ‘ [I4 “”8491“ ' ‘ j '. hi» Figure 1.9: Instantaneous temperature contours (1—3.25) for the F 7‘ = 20 (left), and Fr = 00 (right) 2—D jet flames. The domain shown is 0 < :c/h < 24, —4 < y/h < 4. 37 I I I ’ — F1200 /’ 04" — Fr=40 ,’ d -- Fr=20 I l I _ I I 0.3“ /’ x’ ‘1‘. ’ - rs __ _ ’x :* [(u’u’ + v’u’)dA/f dA ’/ I <2 — I / .602 I I Q. I’ ’ ’ I ‘8 ’z , 0.1 _ / ffidA/fdA ‘ 0 1 I 1 l 1 l 1 0 6 12 18 24 x/h Figure 1.10: Comparison of Product Integrals (lower curves) and total kinetic energy (upper curves) for different Fr numbers. 38 Figure 1.11: Instantaneous temperature contours (300—2200 K) for the buoyant (left), and non-buoyant (right) 2—D methane jet flames. The domain shown is 0 < x/ h < 16, —3.6 < y/h < 3.6. 39 1.5 1 I I I fl __ O-g --- 4—g. 1 1 . l . 0 5 10 15 :r/h Figure 1.12: Jet half-width of the 2—D methane jet flames. 0.06 *- 0.04 - f YCO,dA/ f dA 0.02 - - I l l l 00 5 10 15 x/h Figure 1.13: Total C02 mass fraction in the planes normal to the 2: axis, for the methane jets. 40 0.01 . (a) ON ., D -3 00 0.4 0.003 . , . (b) s . , I D .3 . l 1 00 0.2 0.4 2 Figure 1.14: Scatter plot of reaction rate versus mixture fraction for the methane flame at a: = 15h; (a) buoyant, (b) non-buoyant jets. 41 Figure 1.15: Instantaneous vorticity magnitude (left, 0-3), and temperature (right, 1—3.25) contours in an x — y cross section of the 3—D jet flame with Ze = 8, Fr = 00 (case 12). The domain shown is 0 < z/h < 16, -3.6 < y/h < 3.6. Figure 1.16: Instantaneous vorticity magnitude (left, 0-4.5), and temperature (right, 1—3) contours in an 1: — y cross section of the 3—D jet flame with Ze = 8, Fr 2 40 (case 13). The domain shown is 0 < :r/h < 16, —3.6 < y/h < 3.6. 42 Figure 1.17: Instantaneous vorticity magnitude (left, 0—4.5), and temperature (right, 1-3.5) contours in an a: — y cross section of the 3-D jet flame with Ze = 9, Fr = 00 (case 14). The domain shown is 0 < :r/h < 16, —3.6 < y/h < 3.6. Figure 1.18: Instantaneous vorticity magnitude (left, 0—7), and temperature (right, 1-3.5) contours in an a: — y cross section of the 3—D jet flame with Ze = 9, Fr = 40 (case 15). The domain shown is 0 < x/h < 16, —3.6 < y/h < 3.6. 43 Figure 1.19: Vorticity magnitude contours in the plane normal to the jet axis at z/h = 15 for case 12 (Ze = 8, Fr = 00) (left, 0-2.5), and case 13 (Ze = 8, Fr = 40) (right, 0-9). The domain shown is -3.5 <‘y/h < 3.5,0 < z/h < 4. 44 3.0 T l — Fr=oo, x/h=5 — Fr=40, xlh=5 ---— Fr=oo, x/h=15 - — Fr=40, x/h=15 'E~1 2.0 _ \ \ 1~“~- 1.00 2 3.0 r I (b) — FFOO, x/h=5 — Fr=40, x/h=5 - ’ fffff ~ ---- Fr=°°, x/h=15 - - Fr=40, xlh=15 .5. 2.0 1.00 2 9/51/2 Figure 1.20: (a) Time—averaged values of temperature for Ze = 8 cases, (b) Time- averaged values of temperature for Z e = 9 cases. The overbar represents the time- averaged quantities. 45 1.6 . I I I — Fr=oo, Ze=8 ‘ 1,4 — - - Fr=40, Ze=8 — — Fr=oo, Ze=9 I _ ---- Fr=40, Ze=9 I 51/2 Figure 1.21: Jet half-width for reacting cases. 0.25 . 1 I I ' — Fr=oo, Ze=8 i 0 2 _ - - Fr=40, Ze=8 ' — Fr=oo, Ze=9 ---- Fr=40, Ze=9 ’I I I ’I I I ’I I 15 Figure 1.22: Total product mass fraction in the planes normal to the :1: axis. 46 1.0 r , . — Fr=oo, x/h=5 I O 8 _ — Fr=40, X/h=5 _ ' \ ---- Fr=°°, x/h=15 . \\ - - Fr=40, x/h=15 0.6 — 0.4 - 0.2 - 0.00 1 2 9/51/2 Figure 1.23: Profiles of mean velocity for Z e = 8. 1.3 . I I I — Fr=oo, Ze=8 - - Fr=40, Ze=8 12 _ —— FFOO, Ze=9 _ ---- 13:40, Ze=9 0 I 5 I 10 I 15 :r/h Figure 1.24: Center-line velocity excess for 3-D reacting jets. 47 ,-_ 3:. ---------------------- :§.§ 0 -- ’ “a. .\\ \\\ \\\ . \\ \\ g -0.01 \— I; ‘1 ‘ — Fr=oo, Ze=8 -o.02 -— Fr=40, Ze=8 _ — Fr=oo, Ze=9 —--- Fr=40, Ze=9 _ . 1 . l 1 0'030 5 10 15 x/h Figure 1.25: Vertical velocity on the lateral boundaries of reacting jets (positive Sign is inward flow). O 1 l 4 l I l 1 l 1 0 0.5 1 1.5 2 2.5 9/51/2 Figure 1.26: Profiles of longitudinal velocity fluctuations at :r / h = 15. 48 x/h Figure 1.27: Ratio of turbulence intensity in the finite-gravity jet to the turbulence intensity in the zero-gravity jet on the jet center—line. 49 I I I T m I I I I (a) 3.0 — PI 2.0 — ° 1 l 1 1'00 0.2 I I I I I I I I I (b) ,. 2.0 1.0 Figure 1.28: Scatter plot of temperature versus mixture fraction at :1: = 15h, for (a) Ze=8,Fr=oo, (b) Ze=8, Fr=40. 50 3-0 I I I I I I I I I — Fr=oo - - Fr=40 ‘ I ’ - ‘ \ . / \ - I \ I, \ /\ I, \‘ N 2.0 " \\ _ —— \ 5' \ v \ \ I- \ -I \ \ \ 1 l 1 l 1 l 1 I 1 1'00 0.2 0.4 0.6 0.8 1 Z Figure 1.29: Conditional average of temperature versus mixture fraction for Z e = 8 cases at x/h=15. ' l ' l I l ' I — Fr=oo, Ze=8 _ - - Fr=40, Ze=8 _] 0.02 __ Fm, Ze=9 ,’II/ I _ \\ \\ ---_ Fr=40, 26:9 ”1”, \ \\ Figure 1.30: Conditional average of reaction rate versus mixture fraction for the 3-D reacting jets at x/h=15. 51 l r l 1 —" 10XW(X), Fpoo — 10XW(X), Fr=40 2.0 - . _ _ 52.415, Fm” fi . ---- Sz=o.5’ Fr=40 a A .. ll ], v 111‘. (I) [ 1”?” \ . .‘l I ] A u I 1... ', . a \ I, .\ , k 1.0 — r n g — X :I ,I C13 "' l' :I It ’ .I h! .I ‘ _ — —-—..;——5‘_._.a.a‘r'£‘:;v-_' __ I ‘.’_. _____ _. ._ a .1 4 (a) 1 I 1 L + 0.00 5 10 15 x/ h — 1mm, Fr=8o ‘ I — 10XW(X), Fr=40 J] 3.0 . S Fr“, r " " z=o.s’ " : ‘---- 82:05, Fr=40 \. [a 131 1 I." 03:24) — x '1 _ 3 film”: E f i r x l P, U) I "A I‘\ I, S 1.0 '- "l i ‘t’r’ _, / 1 1.05... I 1w *1} 1 If” 53/ 1’ \' \" .- -..- Ll! q (1)) 1 I 1 l 1 0.00 5 10 15 :r/h Figure 1.31: W(:r) = f dIdydz along the jet axis and area of Z = 0.5 isosurface at t/to = 4, for (a) Ze = 8 cases, (b) Ze = 9 cases. 52 r l' I I _ Frzoo, Ze=8 I!” \\ 3 ” - - Fr=40, Ze=8 1| - _ Fr=°°, ZC=9 I,I l. - Fr=40, Ze=9 '1, ~ 2 _ (a) I l 0.08 l ' l ' . — Fr=°°, Ze=8 - - Fr=40, Ze=8 — Fr=°°, Ze=9 0‘0“ Fr=40,Ze=9 Figure 1.32: (a) Conditional average of F term in the interval 14 < :13/ h < 15, (b) Conditional average of G term in the interval 14 < :r / h < 15. 53 1.0 I ' I I l r I 1 l I I w I ff“~ — Fr=oo, Ze=8 1 - - Fr=40, Ze=8 0'8 _ — Fr=oo, Ze=9 '— . ---_ Fr=40, Ze=9 ~ 0.6 - v I, ~ 6 \ L1,“0.4 ~ _ \6 -1 O 2 ’— ‘1\\ / \ I \ - ‘\ ~ 1-. "‘ I J’ / \ “‘ \ x ~- N 0.0 V" W ...... 1 1 1 1 1 . 1 . 1 2 4 6 8 10 12 14 z/ h Figure 1.33: Axial variation of the correlation coefficient between the F and G terms for reacting jets. 0.25 I I I I 9/51/2 Figure 1.34: Comparison of 2-D and 3-D velocity fluctuations at x=15D, Ze=8, F r = 00,40. 54 l r I ' I — F1=oo, 3-D O 2 _ — Fr=40, 3-D __ ' — - Fr=oo, 2-D ----- Fr=40, 2-D S ‘0’ “ ’ s 3 ----- __ s l a’ \ $3 0 l / \ — / \ / \ / \ x 7 \ — x ~~~~~~ \ ........ \ ~ ~ ‘ “‘ . 1 . 1 . 0'00 1 2 3 y/51/2 Figure 1.35: Comparison of 2-D and 3-D velocity fluctuations at x=15D, Ze=9, F r = 00,40. l ' I — Fr=oo, Ze=9, 3-D — Fr=40, Ze=9, 3—D — - Fr=oo, Ze=9, 2-D ---- Fr=40, Ze=9, 2-D 0.2 T as S \ ‘11 .6“. i 0.1 o 5 10 ' 15 :r/h Figure 1.36: Comparison of 2-D and 3-D product integrals for Ze=9 cases. 55 1 .7 Nomenclature The subscript c denotes the centerline properties in this work, while the the subscripts I and I: refer to the internal and kinetic parts of energy, respectively. Also the subscripts a and i are used as species and coordinate indices, respectively. When referring to a specific species, symbols A (main jet composition), B (coflow), and P (product) are used. All variables are expressed in normalized form, unless noted otherwise, and the subscript 0 is used to refer to reference quantities. The symbol 7(— denotes the time averaged of the variable X , while the < X |Z > notation denotes its conditional average with respect to Z. The symbols and nondimensional variables used in this work are summarized as follows: h: jet height p: density 21: velocity U1: jet velocity U2: cofiow velocity P: pressure T: temperature 6: specific total energy Ya: mass fraction of species a 14,: reference dynamic viscosity no: reference thermal conductivity DO: reference mass diffusivity Rea = pouch/110: Reynolds number Mo = uo/V'yRTO: Mach number F r = 11?, / gih: Froude number Pr = poop/n0: Prandtl number Sc = uo/poDo: Schmidt number Z e = E, / RT 0: Zeldovich number (E, is the activation energy) Da = K fpoh/ no: Damkohler number (K j is the reaction rate constant) de: product formation rate -H°: heat of reaction Ce = —H°/(’y — 1)M3cpTo: reaction heat release parameter Q = Cecbp: heat source term R: gas constant 7: specific heat ratio 9,}: stress tensor Sij: strain rate tensor A: dilatation AU: U1 — U2 61/2: jet half width Z: mixture fraction 56 Chapter 2: Filtered Mass Density Function for the Large-Eddy Simulations of Complex Turbulent Hydrocarbon Flames 2.1 Introduction The ability to accurately model turbulent flames is of great practical and theoretical importance. In this regard, the existence of well defined experimental data such as the Sandia’s turbulent jet flames[33], has sparked the investigation of such flames using different models[34, 35, 36, 37]. Both the large-eddy simulation (LES) and Reynolds-averaged simulation (RAS) methods have been used, and more and less successful results were obtained. In the context of RAS, Jones and Kakhi[38] have used PDF/Monte—Carlo modelling to simulate Masri’s flames [39, 40], which are non- premixed methane jet flames. They have observed that at the speed of 41m/s (flame L), there is no significant extinction, but at the speed of 48m/s (flame B), there are some. In this work, turbulence is modelled with standard 11: — e and second moment closures and a four-step reduced kinetics mechanism was implemented. The linear mean square estimation (LMSE) and the coalescence-dispersion models were used for the mixing model. In the extinction-reignition region of flame B, the temperature predicted by LMSE decreases continuously and remains lower than that obtained from the experiments. The coalescence—dispersion model over-predics the experimen- tal data, but yields closer values to experimental data at other flow regions. Xu and Pope[34] have simulated the Sandia’s flames DE, and F with a joint 57 velocity-composition—turbulent frequency PDF methodology. The 16 species aug- mented methane mechanism [41, 42] have been used in the framework of ISAT [43]. Euclidian minimum spanning tree (EMST) [44], which treats the mixing locally in the composition space is used as the mixing model. It is indicated that the IBM or LMSE model incorrectly predicts local and global extinction in the fast chemistry limit. ' Muradoglu et al. [45] have proposed correction algorithms to the Hybrid particle- finite difference method. Flamelet model was used with the IBM model for mixing, which performs well for equilibrium reactions. The results have been validated using experimental data from Masri’s flame L. Landenfeld et al. [35] have coupled a presumed shape-pdf method with a reduced kinetics mechanism, and have used the intrinsic low-dimensional manifolds (ILDM) method [46, 47] to tabulate the kinetic mechanism versus the reaction progress vari- ables. The joint PDF of scalar properties is assumed to be equal to the product of one-variable marginal PDFs. Transport equations for the mean and variance of mix- ture fraction and species mass fractions are solved, together with a 5 PDF assumption for the mixture fraction and a Reynolds stress closure for the turbulence. Generally, the results are found to compare well with the flame D measurements. While most of the investigations have been based on RAS, some researchers have used LES to calculate the Sandia jet flames[36, 48]. For example, Pitsch and Steiner[36] consider the fl-function for the subgrid scale PDF of the mixture fraction and obtain other species from the solution of unsteady flamelet equations. For the chemistry, a 29-step reaction (reduced from OBI-2.11) was used and the results were compared with the flame D data. The primary stumbling block in LES of turbulent combustion has been the lack of closures which can accurately account for the influence of unresolved “subgrid scale” (SGS) fluctuations[49, 50, 51, 52, 53]. These fluctuations are especially important 58 in diffusion flames as they can alter the magnitude of the “filtered” reaction rate by as much as 100%. Thus if neglected, the results can be very erroneous. The importance of subgrid scale (SGS) modeling in large eddy simulation of turbulent flows has been long recognized. This modeling is required, explicitly or implicitly, when simulating flows with a resolution larger than the Kolmogorov scale. The closure problem is associated with the convective transport, where the non-linearity of the N avier-Stokes equation is exhibited. The problem becomes significantly more difficult in turbulent combustion. This is not just due to the additional non-linearities of the chemical source terms, but also due to intricate complexities of turbulence-combustion interactions [54]. The research field of SGS modeling in turbulent combustion is relatively new; it has been only within the past decade or so that such modeling has been the subject of concentrated investigations. However, the extent of research and the rate of progress have been significant. Several reviews are available of the various SGS methodologies currently in use [55, 56, 57, 58]. One such closure in based on the filtered density function (FDF) methodology [59]. This method is the counterpart of the probability density function method which has proven quite effective in Reynolds averaged simulations [55, 60]. The fundamental advantage of the LES/FDF is that it accounts for the effects of chemical reaction and buoyancy exactly, thus allowing a reliable prediction of turbulent flames over different flow regimes ranging from buoyancy-dominated to momentum-controlled. The literature on SGS closures via the FDF method is growing at a relatively fast pace. The scalar FDF is considered in Refs.[61, 62, 63, 64], the scalar filtered mass density function (FMDF) in Refs.[65, 66, 3], the velocity FDF in Ref. [67] and the velocity-scalar FDF in Ref. [68]. In addition, some recent data on experimental validation of FDF are available[69]. In this chapter, the filtered mass density function methodology [2] is employed as a subgrid scale closure in LES of turbulent partially-premixed methane flames. The FMDF represents the joint probability density function of the SGS scalars, and is 59 determined via the solution of its modeled transport equation by a hybrid Eulerian- Lagrangian numerical scheme. In the LES / F MDF method, a better representation of the turbulence is achieved in comparison to RAS methods. In addition, the F MDF method does not assume a particular profile for the PDF of the scalars and again the reaction source terms appear in a closed form in the FMDF equation. The accuracy of the method was previously established for two-dimensional (2-D) and 3—D temporally developing mixing layers and for a 2-D spatially developing planar jet flames by comparison with DNS and experimental data[2, 63]. In this work, we would like to further assess the applicability and the extent of validity of the FMDF method for SGS closure of turbulent round jet flames. The importance of SGS combustion and chemical kinetics models in such flames is well recognized. Here, two different approaches are followed in the implementation of chemistry. In the first approach, all reactive species are obtained by the direct solution of an appropriate kinetics mechanism. Most of our calculations are conducted with a one-step global Arrhenius kinetics model[2l]. In the second approach, only the FMDF of the passive scalar is considered; the temperature and species concentrations are obtained from the solution of steady-state opposed jet diffusion flames at low strain rates via a flamelet table. The flamelet is generated via a detailed (GRI) mechanism. The flame considered here for testing has regions of local extinction. Considering the steady-state flamelet tables used in the second approach, it is expected that extinction would not be predicted correctly. However, it would be useful to evaluate the overall performance of the FMDF/flamelet method for near-equilibrium flames. The flamelet approach considered is based on a single (low) strain rate table. While larger amounts of strain are known to decrease the equilibrium flame temperature, their effect in a turbulent flame is not well understood. Darabhia[70] has indicated that turbulent flames are able to sustain larger strain rates than expected from the steady-state analysis, depending on the time they are subject to those strain rates. 60 The Sandia’s momentum-driven piloted turbulent methane jet flames [33] are con- sidered for partial validation of our computational models. A series of jet flame ex- periments are conducted[33], which are commonly referred to as the flames A to F. Flame A has an almost laminar structure, flame B is transitional and flames C-F are fully turbulent. Flame D has regions of local extinction, while flames E and F even- tually tend towards total extinction. The existence of local extinction in flames D-F provides a good means of assessing the capabilities of the models to predict extinc- tion and reignition. Flames D, E, and F have been studied extensively using different modeling approaches[34, 37]. Two basic flames considered here are the flames D, and F. The geometrical configurations in these two flames is the same, but the jet inlet velocity is varied. In flame D, the fuel jet velocity is the lowest and the flame is close to equilibrium with limited local flame extinction. The jet speed in flame F is twice as fast as flame D, resulting in noticeable non-equilibrium effects. This work is based on Flame D and the available experimental data for this flame are used as the jet inlet conditions. The organization of this chapter is as follows: After discussing the formulation in 2.2, the numerical procedure is discussed in §2.3. In the numerical procedure section, the consistency issues (§2.3.2), along with the parallelization algorithm (§2.3.3), the implementation of chemistry (§2.3.4), and the computational domain and parameters (§2.3.5) are considered. Results from the flames D and F simulations are shown in section 2.4, followed by the conclusions in §2.5. 2.2 Formulation We refer to the nomenclature section for a listing of the parameters. However, it is useful to define the directional quantities here: x = (:r,y, 2) defines the spatial coordinates, with :1: axis parallel to the direction of fuel injection and y, 2 define the 61 plane perpendicular to 1:. The gravity vector is denoted by g, and the normal gravity is denoted by 1G. Large eddy simulation involves the Spatial filtering operation[71] +00 (f (X We = f(X’1t)H(X'1X)dX' (2-1) —00 where H denotes the filter function, (f (x, t))g represents the filtered value of the transport variable f (x, t). In variable density flows it is more convenient to con- sider the Favre filtered quantity, (f (x, t)) L =(p f )g/ (p)g. For spatially & temporally invariant and localized filter functions, H(x’, x) E H (x’ — x) with the properties[71], H (x) = H (—x), and f; H (x)dx = 1. The application of the filtering operation to the fundamental transport equations of reacting flows[54] yields M 3(p>eLL _ _3eR°L :3 (it: (2.5) l me e <11). (52“) + 12% — 36822)) , <11). = Prz. (2.6) 62 The hydrodynamic SGS closure problem is associated with[72] 7;]. = (p)g((u,-uj)L — (u,)L(uj)L) and M? = (p)g((u,—¢Q)L - (u,)L(q§o)L) denoting the SGS stresses and the SGS scalar fluxes, respectively. In reacting flows, an additional model is required for (p30,); = (p) [(30) L denoting the filtered reaction rates. With the assumption of low Mach number, (p53,); = 0. The last term on RHS of Eq. (2.3) represents the effect of gravity on the filtered velocity field and is in a closed form if (p)g is known. It will be shown that both (p) g and (50);, are determined exactly with the knowledge of the F MDF. The FMDF can be implemented in two ways: (i) to consider only the SGS scalar quantities, and (ii) to consider the SGS velocity-scalar quantities. Only the first method is considered in this work. Scalar FMDF: Let ¢(x, t) denote the scalar array. The scalar filtered mass density function (FMDF) is defined as[63, 65, 3, 64] F111», x; ,t) 2 /_ mp1x',t><1w,¢(x',t)1H(x' - x)dx’ (28> <12», ¢(x, t1] = 61¢ — ¢ A at + 6931' — ali< Dal} >l (910081125 pD 8:5,- 623,‘ Id) [FL/p 63 A 8lLlFL _ 3l5a(¢)FLl 3117i 811/101 , (2.10) where (AlB); denotes the filtered value of the variable A, “conditioned” on B, and the hat is used to emphasize the quantities which are dependent only on the scalar field, zle. Q(x,t) :—: Q(¢(x,t)). In this “exact” transport equation (2.10), the last term on the right hand-side represents the effects of chemical reaction and is in a closed form. The second and the third terms on the right-hand side are unclosed and represent the effects of SGS mixing and convection, respectively. The convective flux is modeled by lama/(0)3) Ruth/ill — (Ui>LlFL = —L Mi = _

[Dt 855;“ 3 a=ha ”a (an) The closure for the SGS mixing can be via any of the ones currently in use in RAS/PDF methods[79, 80, 81]. The simplest one is the linear mean square esti- mation (LMSE) model [73], also known as the IBM (interaction by exchange with the mean)[82] 82 a¢a a¢fl A _ a m [,FL/pl — —a¢a lQm(7/}a " <¢a)L)FLl: (2-13) where Qm(x, t) is the frequency of mixing within the subgrid. This frequency is modeled as Qm = Cg((D)L + Dt)/(A§,), where Cg is a model constant and AH is the characteristic size of the filter. This model will be used here but we are well aware of its 64 limitations[79, 80]. Due to space limitations we are unable to discuss this issue here, but per results of extensive previous efforts by us[80, 81, 83, 84, 85, 59, 60, 86, 87, 88], it can be safely indicated that (i) there are currently no “superior” SGS models which are distinctly more accurate, (ii) while the IBM is not quite satisfactory in RAS/PDF, it functions reasonably well in LES/FMDF, and (iii) its numerical implementation is computationally convenient. With these closures, the modeled FMDF transport equation is BFL 5[e((D)L + DOT + (2-14) 8 36,171.] M [girth/Jo: — <¢O>L)FLl _ awn ' where (8);, is the resolved strain rate, E = |(u?)L(u9)L — ((u?)L)y((u?)L)gzl, V; = CRAgx/F, u? = u,- - U.- and U,- is a reference velocity in the 2:,- direction. The subscript 6’ denotes the filter at the secondary level which has a characteristic size (denoted by AH) larger than that of grid level filter. This model is essentially a modified version of that proposed by Bardina et. al[89] wherethe grid and secondary filters are of equal sizes. We have found the performance of this model to be better than the Smagorinsky type closures[90, 91, 92], and we did not even have to implement the “dynamic” methodology[93, 94] for evaluation of the model constants. 2.3 Numerical Solution Procedure 2.3.1 The Lagrangian Monte Carlo Procedure The most convenient means of solving the FMDF transport equation is via the “La- grangian Monte Carlo” procedure [60, 95]. The basis of this procedure is the same as that in recent RAS[96, 97, 45] and LES/FDF[63, 65, 3, 67]. Therefore, here only 65 some of the fundamental properties of the methodology will be described. With the Lagrangian procedure, the FMDF is represented by an ensemble of computational “stochastic elements” (or “particles”) which are transported in the “physical space” by the combined actions of large scale convection and diffusion (molecular and sub- grid). In addition, transport in the “composition space” occurs due to chemical reaction and SGS mixing. In doing so, the notional particles evolve via a “stochastic process,” described by the set of stochastic differential equations (SDEs) [98, 99] dX,(t) 2 mm), t)dt + €(X(t), t)dW,-(t), d¢:(t) = Ra(¢+, t)dt. (2.15) where X,- is the Lagrangian position of the particles, ’D and 8 are known as the “drift” and “diffusion” coefficients, and W,- denotes the Wiener-Levy process [100]. 9b: denotes the scalar value of the particle with the Lagrangian position vector X ,-. Equation (2.15) defines what is known as the general “diffusion” process [99, 100, 101]; thus the PDFs of the stochastic processes (X,(t), ¢:(t)) are governed by the Fokker-Planck equation. A comparison between the standard Fokker-Planck equation corresponding to SDEs (2.15) with the FMDF equation under consideration identifies the parameters of Eq. (2.15). For FMDF 1 allPMlDlL + Dt)l (10): 333i , 8 a \/2((D),, + 0,), 13,-E (111-)1, + (2.16) 7a., 42mm; — (can) + §..(¢+). With this analogy, the FMDF is represented by an ensemble of Monte Carlo particles, each with a set of scalars 455,")(t) = ¢Q(X(")(t), t) and Lagrangian position vector X("). A splitting operation then can be employed in which the transports in physical and compositional domains are treated separately. The simplest means of simulating the spatial transport in Eq. (2.15) is via the Euler—Maruyamma approximation [102]. 66 Higher order numerical schemes are available[103, 104], but one must be cautious in using them for LES since the diffusion term in Eq. (2.15) strongly depends on the stochastic process X(t). The numerical scheme must preserve the It6—Gikhman[105, 106, 107, 108] nature of the process. Based on previous experience[63, 65], the Euler- Maruyamma approximation provides sufficient accuracy for large Damkéhler number flames. For numerical solution of the hydrodynamic field in LES / FMDF, we will use a high-order accurate finite difference procedure which has proven effective for LES[56]. This discretization procedure is based on the “compact parameter” scheme[109] which yields up to 6th order spatial accuracies. All the finite difference operations will be performed on fixed and equally sized grid points. Thus, the filtered values of the hydrodynamic variables are determined on these points. The transfer of information from the fixed finite difference points to the location of the Monte Carlo particles will be conducted via (fourth and second order) interpolation. The SGS empirical “constants” are C1, C3, S ct, Cg. An efficient way of determin- ing these parameters is by consideration of the flow under non-reacting conditions since there would be no additional modelings when going to reacting flows. Based on our experience, in which a variety of different flows (2D and 3D, constant and vari- able density, different chemistry schemes, etc.) are considered, we have determined C; m 0.01, CR z 0.02, Sc, z 0.5 — 0.7, Cg z 3 — 6. Remarkably, the range of some of these values is the same as that typically used in equivalent models in RAS[55]. The magnitudes of the molecular parameters are the same as those typically used for hydrocarbon-air flames[110, 111]. For methane (k/cp)L z 2.58 x 10‘5 ((T)L/298)0'7, Pr z 0.75 and CM is specified through polynomial fits as functions of the temperature. 67 2.3.2 Consistency The equation governing the first subgrid Favre-moment of the scalar (960),.” as ob— tained by integrating Eq. (2.14), would be identical with the filtered Eq. (2.4). Some of the filtered variables (zle. temperature, scalar and density) can also be computed from the Eulerian (finite difference, F.D.) or Lagrangian (Monte Carlo, M.C.) solu- tions. These imply redundancy and require a mathematical consistency between the finite difference and the Monte—Carlo solutions. Moreover, (¢a)g should be indepen- dent from the mixing model and the constant Cd» However, due to the finite size of the Eulerian grid and limited number of Monte Carlo particles, perfect consistency would not be achieved in practice. Muradoglu and Pope [45] have proposed a scheme for eliminating the inconsistency between the Eulerian and Lagrangian fields in their RAS-PDF calculations. The solution of an O.D.E is suggested for enforcing consis- tency in their simulations. However, in this work, the sources of inconsistency are identified and conditions leading to a consistent solution are considered for our final (Sandia’s) jet flame simulations. It was found that based on the flow gradients and the finite difference grid spacing, different values of the ensemble averaging grid size (A E) might be required for achieving very accurate results. Four simulations with values of C4, = Cg/2 = 8,16 and A E = A and AE = 2A have been performed to study the consistency of our LES / FMDF for conditions of Sandia’s piloted jet flames. The reaction is turned off in these simulations but the variable density/ temperature effects are still very important due to strong pilot at inflow. The instantaneous temperature profiles for cases with A E = 2A and A E = A are compared in Figs. 2.2, 2.3 at the same time and location Instantaneous temper- atures, rather than time-averaged temperatures, are compared in order to highlight the differences. The random number generator used in these simulations is seeded with the same value, to minimize the effect of the randomness in the initial/inlet conditions. It is observed that the difference between the Monte-Carlo and the finite 68 difference solutions becomes negligible when Ag = A. Moreover, the particle mixing model coeflicient, which does not appear explicitly in the first moment Favre-filtered equation, does seem to have some minor effects on temperature when Ag = 2A. In this case, the temperature profiles diffuses slightly more when Q, is doubled (compare results for C4, = 16 and Q, = 8 in Fig. 2.2). However, in the case where Ag = A, the effect of Cd, is negligible and virtually undetectable in the time-averaged quantities. For the reacting flame simulations presented in the next section, Ag was chosen to be equal to A in most cases. It was seen that the difference between the finite difference (RD) and Monte-Carlo (MC) solutions was less than 3 percent in the instantaneous temperatures. For the time averaged values, the difference is even lower. 2.3.3 Parallelization As mentioned before, the filtered Eqs. 2.2-2.4 are discretized using the finite difference (F.D) method on an Eulerian 3—D grid system. For FMDF, due to the large number of independent composition variables, a stochastic Monte Carlo approach is adopted. The Monte-Carlo (MC) method used in this work introduces grid-free Lagrangian particles in a finite difference Eulerian grid background. There is no inter-particle interactions, according to Eq. 2.15, suggesting that the method is potentially a good candidate for porting to a a parallel environment. However, the particle transport and mixing require interactions with the FD. field. Moreover, the source term in the the RD. energy equation should be calculated by averaging the particle reaction source terms. Consequently, the MC. calculations can not be performed independent from the FD. calculations. To obtain a good statistical representation, a relatively large number of particles are required for each F .D. grid. Typically, about ten particles are required for each F.D. grid. A typical hybrid (F.D. / M.C.) simulation with ten particles per grid is about 5 times more expensive than its corresponding F.D. simulation with no MC. 69 particles. Hence, the usage of realistic chemistry in these simulations would only be possible by proper application of parallel processing techniques. Two different parallelization schemes are presented below: 1. In the first approach, all processors are solving the same finite difference equa- tions for the whole domain. Consequently no speed-up is gained in the FD. part. However, since the FD. part is not the computationally intensive part, it is theoretically possible to achieve relatively high parallel efficiencies. In the Monte-Carlo part, particles are equally divided among the processors. Since the particles are not transfered between the processors, there is no communica- tion load from the Monte-Carlo calculations and the load distribution is exactly uniform. Inter-processor communications are only required in the averaging process , where a local ensemble average of particles is calculated on each F.D. cell. Since each processor is dealing with particles traversing the whole domain, global communications are required for calculating the summation and broad- casting the result back to the processors. If n f is the number of finite difference grid points and up is the number of processors, the averaging calculation re- quires a summation of np arrays of size n f . Since this addition needs to be done across the processors, 2 x n f x np values have to be transferred between the processors. On distributed memory architecture, this amount of data transfer can be prohibitive, since the network speed is normally much slower than the memory transfer rate. The drawback of this method is the linear increase in the number of operations and communications, with the increase in the num- ber of processors. It is expected that this approach would only be efficient on shared-memory architecture machines and on small number of processors. 2. The second approach is considerably more complex than the first method. How— ever, it is proven to be scalable and much more efficient. In this approach, both 70 the F .D. domain and the MC. particles are distributed among the processors. Minimum number of communications are required when the subgroup of parti- cles are located in the boundaries of the F .D. domain defined for each processor. This concept eliminates the communications required for the averaging process described in the first parallelization method. However, since each processor needs to keep track of the particles located in its own F.D. subdomain, particles that exit the subdomain have to be transferred to its respective processor. For keeping the same load on every processor, it is necessary to have the same num- ber of particles per processor at all times. In the flow that we are simulating (jet flow), and many other flows, there is a dominant flow direction. This direction is used for domain decomposition, since there would be almost the same number of particles in the slabs normal to the main axis. We have found that only about one percent of the particles traverse the boundaries of the CPU domains in each time step, imposing very little communications load. Other than transferring particles to their respective processor, there are a few more operations where communications are required in this approach. Communications are required for calculating the ensemble averages on the boundary points of the slabs, un- less Ag is chosen to be equal to the finite difference grid size and zeroth-order interpolations are used. A larger selection for Ag would increase the communi- cations load, since data from additional planes of the neighboring slabs need to be transferred. Inter-processor communications are also required for interpolat- ing a finite difference quantity on a particle location. Differencing and filtering are the operations in the finite difference domain, where communications will be necessary in this method. A 16 processor cluster with 2 processors per node was used for the benchmarking tests as well as for the final simulations. The processors are Intel Xeon working at the clock rate of 2.8GHz, with a cache size of 512kB. Gigabit networking along with MPI 71 was used for inter-processor communications. A simulation with 1.2 million finite dif- ference grid points and 7 million numerical particles was considered for benchmarking. Single-processor simulations were not possible due to the program structure, thus the speedups reported are relative to a 2 processor simulation. Speedups of 2.4, 8.0, and 15.1 were obtained for the 4, 8, and 16 processor simulations respectively. These results are plotted in Fig. 2.4. In the case of 4 and 8 processor simulations, superlin- ear speedups are obtained, since the problem size for each processor becomes smaller and the higher speed cache is used more efficiently. For the 2, 4, and 8 processor simulations, only one processor per node is being used. However, for the 16 processor simulation, 2 processors per node are used, which have to share the memory and network bandwidth, resulting in a sublinear speedup in this case. It is seen that despite the apparent complexity of this approach, it is much more eflicient than the first method. Turbulent reacting jet simulations with 1.2 million finite difference grid points and 40 million Monte-Carlo elements have been performed successfully on the above hardware. 2.3.4 Implementation of Chemistry In quantitative comparisons with laboratory data a very important issue is associated with the “chemistry of combustion.” The flames are considered in this work under the following two conditions: (i) non- equilibrium, finite rate reduced kinetics, (ii) near equilibrium with realistic kinetics. In (i) the finite rate kinetics effects are modeled via a one-step global mechanism[112], or 12—step reduced mechanisms. In (ii) simulations are conducted with the use of the flamelet library at low strain rates. The flamelet library is generated in a laminar opposed-jet simulation via a detaile chemistry model. In (ii), the transport of the FDF of only the mixture fraction is considered. All other thermo—chemical variables are constructed via the flamelet library. In (i), the transport equation for sensible 72 enthalpy (h, = f]: cp(T’)dT’ is solved 5%thle 3

€LL a __ (9(th (Ji >[ ~ L8L 8331' 3 Mi “"

€Dt axi . (2.18) Moreover, the term (pSahg)g is approximated as (p)g(Sahg)l. In (ii), a transport equation for the filtered value of E = RT is solved (in the flamelet table RT is a function of Z). A transport equation is derived for (E) L, by multiplying the modelled FMDF transport equation by E and integrating over the mixture fraction space. maggot + a

1L : it)? _ fig +

znm11 — LL1 (2.19) where _dE(Z) 0 ~ 1 k 3(E)L G— dZ ’ (JihN LeeD¢ (2.20) The reduced scheme considered in (i) is the 12-step mechanism of Sung et al.[113]. This reduced mechanism is developed from the GRI-MECH 1.2 and involves 16 species (H2, H, 02, OH, H20, H02, H202, CH3, CH4, C0, C02, CHgO, C2H2, C2H4, C2H6, N2). It contains more non-steady-state intermediates than the conventional 4 and 5-steps mechanisms, and it has proven effective in an extensive range of ap- plications including auto-ignition, and laminar flame propagation with strong varia- tions of thermodynamic properties[42]. The 12-step mechanism was also validated in the framework of steady-state laminar opposed jet flames. The 12-step results were found to be almost indistinguishable from the GRI results. For instance the tem- perature profile in the mixture fraction space as calculated by reduced and detailed mechanisms[114] for the conditions of the Sandia flames and the strain rate of 10/3 73 shows an excellent agreement between (Fig. 2.5). 2.3.5 Computational Domain and Parameters Figure 2.1 shows the geometry of the Sandia piloted methane jet flames, and the coordinate system used in our simulations. For these simulations, a finite difference mesh with 160 x 121 x 121 grid points was considered for a domain of 20 x 12 x 12 jet diameters in the x, y, and 2 directions, respectively. The number of Monte-Carlo elements for each FD. grid point was 10 and Ag = AG as mentioned in 2.3.2. The main jet composition is 25% CH4 and 75% air with a Reynolds number of 22400. Flame F has the same parameters, except for the jet speed or Reynolds number, which is doubled. Detailed specifications and measurements of these flames are available on the Sandia’s web—site [115]. The velocity measurements were done seperately by [116], using laser-Doppler anemometry. The main parameters of flame D are as follows: i o Nozzle diameter = 7.2mm a Pilot diameter = 18.2mm 0 Main jet: 25%CH4 and 75%air o Reynolds number: 22400 0 Avg. jet speed: 49.6m/s 0 Peak jet speed: 63.1m/s 0 Main jet temperature: 300K 0 Pilot speed: 11.1m/s 0 Pilot temperature: 1880K 0 Pilot mixture fraction: 0.27 o Co—flow speed: 0.9m/s 74 Non-reflecting boundary conditions [117] are considered for the outlet boundary and zero-gradient conditions are used for the lateral boundaries. For the inflow, non-reflecting boundary conditions are used with prescribed velocity components. It should be noted that temperature composition measurements were not possible at the nozzle exit, and x/d = 1 is the closest distance to the nozzle where data are provided. Moreover, in these near-nozzle locations, less accurate Raman was used rather than Rayleigh scattering. The sensitivity of model predictions to uncertainty in the pilot boundary conditions is an important consideration [118], specially with regards to flame F, which is very close to global extinction. Fig. 2.6 shows the mean and fluctuations of velocity, measured at the nozzle exit [116]. The mean axial velocity was prescribed in the inlet as given in this figure. However, to resolve the sharp velocity gradients, the measured profile was smoothed by an explicit filter. The mean of other velocity components was set to zero at the inlet. The amplitude of the velocity fluctuations was set equal to the measured values. For flame F, velocity measurements were not available and the flame D values are scaled for the inlet conditions. 2.4 Results for Turbulent Methane Jet Flames The isosurfaces of vorticity magnitude, and temperature for flame D and F simulations are shown in Figs. 2.7 - 2.10. These figures represent instantaneous quantities and indicate a transition to turbulence at 23/ D z 10. The temperature contour in the flamelet flame D (Fig. 2.8) case is a continous surface, consistent with the nature of our flamelet approach. The l-step flame F temperature contour (Fig. 2.9) has discontinuities starting from x/ D x 10, indicating local extinction in some regions. In the case of the 12-step flame D case (Fig. 2.10) the isosurface is severely broken from x/D z 16, and it does not exist after as/D z 18. 75 Fig. 2.11 shows the time averaged velocity profiles for the flame D simulation with flamelet chemistry. The values are Favre-averaged and are compared to Favre- averaged data from the experiments. The statistical error associated with the mea- surements of the mean velocity is below 5%, while the accuracy of the RMS measure- ments is about 10% [116]. The agreement between the calculated and the measured filtered values is generally good. For Flame F, there are no measured data available for the velocity profiles. The velocity profiles calculated by the 1-step and flamelet chemistry models are very close, indicating that the chemistry is not substantially affecting the mean flow. The results obtained by the 12-step mechanism simulations are also close to those in Fig. 2.11 and are not shown. The RMS values of axial velocity are shown in Fig. 2.12. At x/D = 7.5, the predicted RMS is lower than the experimental values at r / D values larger than the r/ D of the peak RMS (outer side of the peak). The agreement is better at the inner side of the peak location. Moreover, the flamelet and l-step results are close to each other at the outer side of the peak, while the flamelet simulation predicts higher RMS at the inner side at both cross sections. It is known that heat release generally damps the turbulent perturbations in a jet flame. Since the 1-step chemistry has generally larger heat release compared to the flamelet chemistry, the RMS of velocity is smaller in the 1-step case, specially at the inner part of the peak where the reaction is more important. The higher heat release of the 1-step chemistry is also evident in the mean temperature profiles (Fig. 2.15). The lower intensity of turbulence at 23/ D = 7.5 is due to the smoother inlet velocity profile in the simulations, which effectively delays the growth of the perturbations. Farther downstream at a: / D = 15, the perturbations have grown to the values measured in the experiments. Fig. 2.13 shows the mean mixture fraction profiles from the 1-step and flamelet cases. The trend observed in the mean velocity profiles (Fig. 2.11) can be seen here as 76 well. The mixture fraction profiles shown are obtained from the finite difference values and are consistent with the mean values calculated from the Monte-Carlo values. V The RMS values of mixture fraction are shown in Fig. 2.14. Similar to the RMS of velocity, it is seen that the RMS of mixture fraction is lower in the 1-step case, specially at the inner shear layer zone. However, compared to the velocity RMS, the LES profiles are closer to the experimental profiles :r/ D = 7.5. At :r/ D = 15, the agreement is even better and the peak RMS values match. However, the centerline values do not follow the trend in velocity and are different from the experimental measurements. Favre-averaged mean and RMS of temperature for flame D, are shown in Fig. 2.15 and 2.16. The accuracy of the temperature RMS measurements is not known. The RMS of measured temperature in a flat flame was about 1%, while for mass fractions it ranged betwen 5% and 20% [33]. The predicted mean temperatures are mainly higher than the corresponding experimental values. Moreover, the l—step case has higher temperatures compared to the flamelet case. Considering the fact that the 1—step simulation uses an irreversible reaction, and CO is not included in this mechanism, the agreement between the predicted and experimental data is reasonable, except for the jet core region. This agreement suggests that the reaction is fast around the shear layer. The agreement is better when the flamelet model is used, mainly because the CO is included in the mechanism (Fig. 2.26). Fig. 2.16 shows that the predicted and measured values of the RMS of temperature are generally in good agreement. The RMS of temperature in the l-step case is higher than the flamelet values, which is in contrast to the RMS of velocity (Fig. 2.16) and mixture fraction (Fig. 2.14). The RMS of velocity is lower for the 1-step case due to higher heat release. However, in the case of temperature, the reaction for the l-step case is very fast and as the compositional structure of the flame shows (Fig. 2.35), the l-step flame temperatures are much higher than the the flamelet temperatures in 77 the rich side of the flame. In other words, at any given composition, the temperatures are either below the ignition temperature or they have reached the maximum possible values (Fig. 2.35). Thus the RMS of temperature is affected by the combined factors of velocity and mixture fraction fluctuations on one hand, and the chemical mechanism on the other. Figure 2.17 shows the mean temperature profiles for the flame F simulations. The measured temperatures in flame F are much lower than the flame D; at :1: / D = 15, the peak temperature is about 1100K in flame F compared to 1750K in flame D. There are some differences at :r/ D = 7.5 as well, but not as significant as those at x/ D = 15. Moreover, it is observed that the centerline temperatures have not changed compared to flame D. The flamelet case has higher temperatures compared to the 1—step case, since the flamelet chemistry does not allow for extinction. The 1-step solution is successful in predicting the temperatures of the outer re- gion of flame F, while highly over-predicting the jet core temperatures. Moreover, the predicted temperatures at the centerline are even higher than their flame D counter- parts. The core region of the jet is a fuel rich region, where the l-step reaction does not perform well (Fig. 2.35). This fact combined with the stronger mixing present in flame F, have caused the unrealistic overprediction of temperature at the centerline. Fig. 2.18 shows the axial variations of the mixture fraction along the centerline. Compared to the values predicted by the l-step model, the agreement is much bet- ter with the flamelet model. Both models overpredict the mean mixture fraction. The slower rate of decay in the 1-step simulation can be attributed to higher heat re- lease and volumetric expansion, which counteracts the turbulent mixing and diffusion processes. The predicted rate of velocity decay at the jet centerline (not shown) is lower than the experiment, specially for the 1-step case. The difference could be due to the higher heat release and consequently higher damping effect of the combustion at the 78 centerline (see Fig. 2.19). Fig. 2.19 shows the profile of temperature along the centerline. As expected from the trend in the mixture fraction (Fig. 2.18) and the temperature cross stream profiles (Fig. 2.15), the centerline temperatures overpredict the measurements. The overprediction is higher in the 1-step simulation. The RMS values of temperature at the jet centerline are shown in Fig. 2.20. Again, the RMS is overpredicted by the 1-step model, mainly due to the fast reaction and the nature of this global kinetics model. The agreement is better in the flamelet simulation. The mass fraction of the main reactants, CH4 and 02 are shown in Figs. 2.21 and 2.22, respectively. For CH4, the predictions are good, specially at :r/ D = 7.5. However, for 02, the simulations underpredict the measurements at x/ D = 7.5 at the outer part of the shear layer. Since 02 is originating from the coflow, the lower values suggest that the mixing between the coflow and the main jet / pilot is not as strong as the experiment. This is also confirmed by the lower values of velocity RMS at x/ D = 7.5 (Fig. 2.12). As the mixing becomes stronger farther downstream, the difference between the predicted and measured values reduces (Fig. 2.22, x/ D = 15). The mass fraction of the main products, H20 and 002 are shown in Figs. 2.24 and 2.23, respectively. The agreement between the predicted and measured values for H 20 is acceptable. However, the 002 values are overpredicted in the 1-step simulation, since C0 is not included in the 1-step mechanism. Figs. 2.26 and 2.25 show the profiles of the CO and H2 mass fractions, respectively. These profiles are from the flamelet simulation, as these species are not available in the l-step model. The C0 mass fraction of is considerable relative to the mass fraction of C02, explaining the overpredicted values of C02 in the 1-step simulation. In the case of H2, the computed values are considerably higher than the measured ones. Fig. 2.27 shows the profiles of OH mixture fraction. The OH values match better 79 than the H2 values, slightly underpredicting the peak values in the measurements. The mean heat release from the flamelet and 12-step cases are shown in Fig. 2.28. The a:/ D = 7.5 profiles are relatively close. However, at x/ D = 15, the heat release from the 12-step simulation is substantially lower. This is consistent with the disrupted isosurface of temperature at the same location (Fig. 2.10). Compositional Flame Structure: The experimental scatter plots of tempera- ture vs. mixture fraction for flame D and F are shown in Figs. 2.29 - 2.32. For flame D, as Figs. 2.29 and 2.30 show, the scatter in the a: / D = 15 plot is only slightly higher than in the as/D = 7.5 plot. For flame F, the scatter in the x/D = 7.5 and :c/D = 15 plots are not noticeably different, and at both locations, considerable local extinction exists (Figs. 2.31, 2.32). The scatter plots of temperature as obtained by LES/FMDF with the 12-step mechanism are shown in Figs. 2.33, and 2.34. At 55/ D = 7.5, the values are very close to the flamelet line, while at a: / D = 15, there is a substantial amount of scatter. The scatter predicted by the 12-step chemistry simulation at :r / D = 15 (Fig. 2.34), is substantially larger than the scatter evidenced in the experimental data (Fig. 2.30) Although the l-step global mechanism cannot correctly predict the temperature for the whole range of mixture fraction values, its scatter plot is shown (Fig. 2.35) to clarify the behavior observed for this mechanism in previous figures. 2.5 Conclusions The LES / F MDF method is applied to Sandia’s piloted turbulent methane jet flames. The Sandia flames are chosen for comparison purposes, since extensive experimental data are available for them. Both the flamelet and finite rate kinetics approaches were considered for the simulations. In the kinetics approach, the l-step and 12- step mechanisms were used. The consistency of the Monte-Carlo and finite difference 80 solutions were discussed and a scalable algorithm for parallelization of our hybrid (Eulerian-Lagrangian) methodology is presented. The algorithm was shown to yield superlinear speedups, when single processor nodes are used. The Favre-averaged tem- peratures for two jet Speeds (flames D and F) were compared with the experimental data and reasonably good agreement was observed for the lower jet speed (flame D) with less non-equilibrium effects. Higher degrees of local extinction in the higher jet speed (flame F) causes larger differences between the predicted and measured temper- atures; suggesting that a more accurate finite—rate reduced chemistry model is needed in this case. In the l2-step simulations, it was observed that after a: / D = 10, the flame exhibits substantial extinction. This extinction does not lead to a total blow-out, but rather causes a drop in the mean temperature, and heat release. The main factor that contribute to this problem is considered to be the inadequate resolution of the sharp profiles of the velocity and composition at the inlet. The sharp gradients present at the inlet, compound the problem of an incompatible geometry. Employing larger number of grid points would be prohibitive, in terms of cost. Although the resolution at the domain inlet might not be suflicient to resolve the sharp gradients, the turbulence developed at farther downstream positions is adequately resolved, as shown in the scalar variance profiles. 81 2.6 Figures Figure 2.1: Schematic representation of the Sandia piloted methane jet flame. Nozzle diameter = 7.2mm, Pilot diameter = 18.2mm. 82 S r l I A l I I I ,f/ \\ __ C¢=8.M.C r I”, \ ._ _ C¢=8.F.D 4 ,_ If], \\ _ Car-l6. M.C __ ’1 \ --_- C¢=l6.F.D j \ b3 — — // I 2 — _ 1 i I l L 1 I 1 O 0.5 l 1.5 2 r / D Figure 2.2: Comparison between Monte Carlo (MC) and finite difference (F.D.) values of the normalized instantaneous temperatures for Ag = Ag, and C), = 8, 16. 5 r 1 v //’\\ _ C =8. M.C I ¢ 7 ,1 \ _ _ C¢=8.F.D I’ \ - 4 _ / _ C¢—l6. M.C __, I ---- c,=16.1=.D : . E:13— 4 .— 4” ,r 2 _ _ 1 1 10 l 2 r / D Figure 2.3: Comparison between Monte Carlo (MC) and finite difference (F.D.) values of the normalized instantaneous temperatures for Ag = %Ag, and C4, = 8, 16. 83 I l I I I I I l I I I I I __ —J( 15 [x LES/FMDF . ocfi (‘0 .. \aefi ‘ .565 0,09 '19“ 10 - - Q. :3 g b x 3 f . c“ 53‘ ax’sX\\c\ efiW" \0°°'° 5 - _ l- x -1 1 l 1 I l I l I 4L 1 1 l l 2 4 6 8 10 12 14 16 no. of processors Figure 2.4: Speed-up of the parallelized LES/FMDF calculations relative to the 2- processor calculation. 2000 1500 1000 500 Figure 2.5: Temperature distribution vs. mixture fraction in a the low-strain (10/3) laminar opposed jet flame. 84 Figure 2.6: Experimental mean and variance of velocity profiles at the jet inlet for flame D. 85 Figure 2.7: Isosurface of enstrophy Q = 1 from the flamelet, flame D case. The domain shown is 0 < r/d < 20, —5 < y/d < 5, —5 < z/d < 5. Figure 2.8: Isosurface of (T) L = 1500K from the flamelet flame D case. The domain shown is 0 < :c/d < 20, —5 < y/d < 5, —5 < z/d < 5. 86 Figure 2.9: Isosurface of (T) L = 1500K from the l-step flame F case. The domain shown is 0 < :r/d < 20, —5 < y/d < 5, —5 < z/d < 5. Figure 2.10: Isosurface of (T) L = 1500K from the 12-step flame D case. The domain shown is 0 < x/d < 20, —5 < y/d < 5,—5 < z/d < 5. 87 I I I I C\ O — — x/D=7.5, l-step 60 -‘~\ 0 — xlD=15, l-step ‘ 1 \ ----- fo=7.5, flamelet \ O -——- x/D=15. flamelet ' Cl ‘\ O x/D=7.5 experiment ‘ \ . D xlD=15 experiment \ \ 40 — . _ it 2 \\\ 3 _ \\\\\ O 4 \\ \\ \ \ 20 _ \§\ a \\\€\ .. ‘s\ a o flax O I 0 1 l 10 . z 3 ‘ - 3 O 1 2 3 T/ D Figure 2.11: Mean axial values of the filtered velocity field for the 1-step and flamelet chemistry in flame D. 10 . c1 1 . 1 r [- —— x/D=7.5, l-stcp ‘ D — xID=15.1-step 8 - ----- x/D=7.5, flamelet -‘ — xlD=15.flamelet - O xlD=7.5, experiment ] 0 x/D=15. experiment :5 6 3 53’ 2 4 Q: 2 O 1 L 4 I \.'.——--—---1 0 1 2 3 T/ D Figure 2.12: RMS of the filtered values of axial velocity for the l-step and flamelet chemistry in flame D. 88 I I I — — x/D=7.5, l-step ~ — x/D=15, l-step ----- x/D=7.5. flamelet -] -— x/D=15, flamelet O x/D=7.5, experiment ~ 0 x/D=15. experiment r/D Figure 2.13: Mean mixture fraction for the l-step and flamelet chemistry in flame D. 0.2 r I ‘ ' ' K —- XID=7.S.1’Step _ XID=15. I'Stcp a ----- x/D:7.5,flamelet . 0 — x/D=15.flamelet ‘ o xlD=7.5,experiment A O ‘ n x/D=15. experiment 1 ’16 A ,’ \\ £20.1— ,\ tn ,"0 I \ \‘ 2 a I’ I ‘dx 03 u n I u I O ——-’ O I 0 O I / / / 0 ’ : ‘ 0 l Figure 2.14: RMS of mixture fraction for the l-step and flamelet chemistry in flame D. 89 2000 . 1 . 1 . {:5 —— x/D=7.5,1-step .. {f0 \\ — x/D=15,1-stcp - ’1' \ ----- xlD=7.5.flamclet I : , — x/D=15.flamelet 1500 ~ . “I o xID=7.S. experiment -* I \u D xlD=lS, experiment q ..., I \ -1 A 0 E1, \ 1_ o _. 1000 “P o \\ \ o \\ D 500.. ,0 o \ \ 4 Jo’/0 O \ \ - ' " . 1 1 1 0 ‘1‘"—— 0 2 3 T/ D Figure 2.15: Mean temperature for the l-step and flamelet chemistry in flame D. 5% I I l I I ’ A D — - x/D=7.5, l-step 1 I I o — x/D=15, l-step 400 - I ‘ —-- x/D=7.5,flamelet ~ — xlD=15, flamelet . I 9‘. ‘ 0 o x/D=7.5. experiment ~ I I, x ‘ u xlD=15, experiment ”3300 _ ’1' \I o / \ - A ’I \ ”l. \“‘ z '- I b, \\\\ O , o \\\ \\ D -1 z 200 ,’/ \‘v,’ o \ - B I, \ \ D o: _ . x \ . ll f \\\ \ _ 0 \ _ 100 [,6 ,9 D [- [1, \\\‘ -1 o ‘~\~__ O L I l l o 1\ ;.::--F 0 l 2 3 r / D Figure 2.16: RMS of temperature for the 1-step and flamelet chemistry in flame D. 90 2000 1500 — ' l I — — x/D=7.5. l-step — x/I)=15, l-step ----- xlD=7.5, flamelet —— x/D=15. flamelet 0 x/D=7.5, experiment 0 x/D=15. experiment Figure 2.17: Mean temperature for the 1-step and flamelet chemistry in flame F. l 0.95 - A“ \N/ - 0.9 — P — l-step .. — flameiet 0.85 _ 0 experiment _ 1 L 1 i l 0 5 10 15 a: / D Figure 2.18: Mean models in flame D. mixture fraction at the jet centerline for the l-step and flamelet 1000 V l 1 l I l — l-step _ — flamelet 4 0 experiment 750 - - ’4 .. SOO— 1 l 1 l 1 l 250o 5 10 15 Figure 2.19: Mean temperature at the jet centerline for the l-step and flamelet models in flame D. 300 v I ' I ' l — l-step __ — flamelet . 0 experiment 200 - :3 b a“ E Q: 100 - O ‘ O Figure 2.20: RMS of temperature at the jet centerline for the 1-step and flamelet models in flame D. 92 0.2 _o O U: l Figure 2.21: Profiles of mean CH4 mass fraction for the l-step and flamelet models I' I — — xlD=7.S, l-step — x/D=15. l-step ----- x/D=7.S. flamelet — x/D=15.flamelet 0 fo=7.5. experiment 0 x/D=lS, experiment in flame D. 0.25 0.2( 0.15 ,5 C? E; 0.1 . — — xID=7.5, l-step — xlD=15,l—step ----- x/D=7.5,flamelet 0.05 _ \\\ 0,,1' / — xID=l5,flamelet {v’ I O xlD=7.5. experiment ’ ”' D xlD=15,experiment 0 ' ' ' 0 l 2 Figure 2.22: Profiles of mean 02 mass fraction for the l-step and flamelet models in flame D. 93 r j I I l — — x/D=7.5, l-step _ I' "" \ — x/D=15. l-step "l I \ ----- x/D=7.5, flamelet I \ — x/D=15. flamelet / \ O xlD=7.5. experiment 4 0.1 — / 0 A \ c] xlD=lS. experiment I I” \ ,3 ’ ’ = . g 9 I > I II \\ v0 05 — I ' ° \ - I (I7 \\\ U I \\ l ,1 I O \‘\ b I 9” Q“ I ' , o \ ' U 9” 0 §\ .‘\ . , \ x . ”1'3"“ 1 1 1 01 ‘~~‘..~ \== 0 l 2 3 Figure 2.23: Profiles of mean 002 mass fraction for the l-step and flamelet models in flame D. I I U I I —— x/D=7.5, l-step _ — x/D=15.l-step . ----- x/D=7.5,flamelet , — 1:10:15. flamelet 1" ‘_\a‘ o x/D=7.S, experiment 0.1— ;/o u xlD=15.experiment ‘j ,5 0.. S 0.05 Figure 2.24: Profiles of mean H20 mass fraction for the 1-step and flamelet models in flame D. 94 0.004 I I T I i- — — x/D=7.5, flamelet q xlD=15, flamelet O x/D=7.5. experiment 0.003 "— [\ D x/D=IS, experiment _. Figure 2.25: Profiles of mean H2 mass fraction for the flamelet model in flame D. 0.05 7 I l I l i / 1 — — x/D=7.5, flamelet - f — x/D=IS. flamelet 0.04 " O x/D=7.5. experiment - D x/D=15. experiment 0.03 - _. 2 8 I. d E; 0.02 — _ 0.01 "’ I \ _ I! o \ " ’/ o \\ -4 0% FLO I 1 O (3 fi 3 l 2 3 Figure 2.26: Profiles of mean CO mass fraction for the flamelet model in flame D. 95 0.002 . , . 1 . O _ — — x/D=7.5. flamelet . — x/D=15. flamelet O O xID=7.S. experiment 0.0015 — D x/D=15. experiment — Figure 2.27: Profiles of mean OH mass fraction for the flamelet model in flame D. 150 ' I I I l ----- x/D=7.5, flamelet — x/D=15. flamelet — — x/D=7.5. lZ-step — x/D=15.12-step Figure 2.28: Profiles of mean heat release for the flamelet and 12-step models in flame D. 96 2500 v 1 I I I I r I I . ' I ' o o ' 0. lg... ‘5‘ co 1- o 8 O. O . :- .\ h'fg. ” -1 o o o 0 19 ' $00 0 0° 0 ° ' .~ ’5: o 1500 — ° . v _ ‘10 0 E4 _ O 50 ,1“. . l 0° ‘ ° 0% \z' o o Q‘; .-~\ ‘K. '. . .- __ o ' .' .-. 1 o .0 1000 ° ad“ it ..., 00.0... u ‘00 500 - q 1 L 1 l 1 l 1 l 1 0 0.2 0.4 0.6 0.8 1 Z Figure 2.29: Scatter plot of temperature vs. mixture fraction for flame D; experimental measurements at x/ D = 7.5, and r/ D = 0.83 (975 samples). 25m I I fit I I I 1' I 7 2000 1500 1000 500 Figure 2.30: Scatter plot of temperature vs. mixture fraction for flame D; experimental measurements at x/ D = 15, and r/ D = 0.83 (780 samples). 97 2500 .rr1r1111 2000 1500 1000 500 Figure 2.31: Scatter plot of temperature vs. mixture fraction for flame F; experimental measurements at 33/ D = 7.5, and r/ D = 0.83 (1550 samples). 250011111111- 2000- 00600 0 Och" 3,1191 ‘9 2 o o :ggo o °¢P o 1500— . hi; . v w" - ° 0 q, 8 o I!“ 0° [‘4 0 081°”! 000% Al, _ ° 0° °° ° 0 . . " o 03% 2030?; 1 1000 — . 01:511. ° — °o 00: can 0 2.— " ° o°Jn°0 3° " 500 Figure 2.32: Scatter plot of temperature vs. mixture fraction for flame F; experimental measurements at x/ D = 15, and T/ D = 0.83 (976 samples). 98 2500 1 1 1 1 ' 1 ' I 2000 1500 1000 500 1 1 1 1 1 1 1 1 1 0 0.2 0.4 Z 0.6 0.8 l Figure 2.33: Scatter plot of temperature vs. mixture fraction for the 12-step flame D case, from x/D = 7.5, and r/D = 0.83 (709 samples). 2500 T 1 1 1 1 1 T 1 2000 1500 1000 500 Figure 2.34: Scatter plot of temperature vs. mixture fraction for the 12-step flame D case, from x/D = 15, and r/D = 0.83 (635 samples). 99 250011111111 2000 1500 1000 500 Figure 2.35: Scatter plot of temperature vs. mixture fraction for the l-step flame D case, from :c/D = 15, and r/D = 0.83 (958 samples). 100 2.7 Nomenclature CPO Cp Co : The model in the VSFDF : Specific heat of species a at constant pressure : Specific heat of the mixture at constant pressure : Constant in the VSFDF formulation : Constant in the VSFDF formulation : Constant in the subgrid scalar stress : Constant in the subgrid scalar stress : Constant of the stochastic mixing closure : Molecular diffusion coefficient : Subgrid diffusion coefficient : Drift coefficient in the SDE : Inner jet diameter : Diffusion coefficient in the SDE : Joint scalars filtered mass density function : Joint velocity-scalars filtered mass density function : Froude number : Normal gravity : Gravity vector : ith component of the gravity vector : Filter function : Filter function : Enthalpy : Enthalpy of species a : Enthalpy of formation of species a : Jet momentum : ith component of the flux of scalar a : Thermal conductivity of the mixture : Luminous flame height : Molecular Lewis number : Lateral dimension of the flame : ith component of the subgrid scalar flux of species a : Number of species : Molecular Prandtl number : Universal gas constant : Velocity ratio at the inlet 1‘ = %m t : Reynolds number : Richardson number : SGS Schmidt number : Production rate of species a : Strain rate tensor : Temperature 101 (“$511353 xggs <1§$+E QC!“§$‘3 : Reference temperature : Subgrid scale stresses : Time : Velocity vector : Streamwise velocity at the jet inlet : ith component of the velocity vector : ith component of the stochastic velocity vector : ith component of the reference velocity in MKEV closure : Velocity space : ith component of the velocity vector in composition space : Molecular weight of species a : Wiener-Levy process : Position vector : ith component of the position vector : Lagrangian position of the particles : Streamwise coordinate : Mass fractions of species a : Coordinates defining the plane normal to :1: : Mixture fraction Greek Symbols: ’7,- : Standardized Gaussian random variable A : Grid spacing in LES 6 : Dirac delta function 15,-,- : Kronecker delta AH : Grid level filter width A H’ : Secondary level filter width AU : Velocity difference at the inlet, AU = Um — U1m AP 2 pin _' pf € : SGS dissipation ,u : Molecular viscosity, ,u = pz/ Vt : Subgrid viscosity 9m : SGS mixing frequency 1p : Composition space (I) : Scalar field (pa : The compositional values of scalar Oz 3' : The compositional values of stochastic scalar a p : Density 0 : Number of scalars, a = N3 + 1 Ti, : Molecular stresses Symbols 102 ( )1; : Filtered value ( )L : Favre filtered value ( | )L : Conditional Favre filtered value A : Quantities which depend only on the scalar composition, zle. §(¢,x,t) E S(¢(x,t)) oo : Ambient Subscripts : Flame in : Inner jet k : Discretized (time) level on : Outer jet Superscripts + : Properties of the stochastic Monte Carlo particles (11) : Index of the Monte Carlo particles 103 Chapter 3: Large Eddy Simulations of Buoyant and Non-Buoyant Turbulent Jet Flames 3. 1 Introduction The importance of gravity in diffusion flames is well recognized[119, 120], and it is rather easy to demonstrate the effects: all one needs to do is to visualize the spatial structure of such flames in both “upward” and “downward” fuel injections in nor- mal gravity[121]. Unless the jet velocity is very large, the flame structures will be different in the two configurations. It is also easy to recognize that the reason for this difference is due to mechanisms of “buoyancy-driven” flows. In the case of an upward jet flame, buoyancy results in an upward transport of combustion products. This naturally modifies the entrainment process and the rate of reactant conversion as compared to those in zero or microgravity. While this simple experiment demon- strates the importance of buoyancy in jet flames, the quantification of gravity effects on the physical and “compositional” structures of momentum-buoyancy dominated transitional and turbulent jet flames is not fully established. There have been significant contributions on experimental investigation of buoy- ancy effects in jet flames (ag. Refs.[122, 123, 124, 9, 125, 126, 127]). These efforts have been noticeably enhanced under the Microgravity Combustion Science Program sponsored by NASA[119, 120, 128]. These experiments have been very useful in elu- cidating the spatial flame structure in microgravity, and also in providing valuable data sets for model validations. With these results, it is now understood that modifi- cation of the flame structure by gravity is primarily due to vortex-flame interactions 104 as displayed by the large scale flow structures. In normal gravity, these structures are caused by the inherent instability of the fuel jet coupled with the instability of the buoyancy induced natural convection layer.[123, 125, 127] In reduced gravity levels, the influence of the second instability mechanism is decreased [129], resulting in a flame structure which is different from that in normal gravity. Most previous computational investigations of gravity effects in transitional and turbulent flames were based on “direct numerical simulation” (DN S). Amongst such efforts are the contributions of Kailasanath et al.[130] and Elghobashi et al.[16] The former deals with premixed flames and mostly in laminar (albeit unsteady) flows. The latter deals with DNS of a nonpremixed reacting flow under idealized conditions: homogeneous turbulence, one-step chemical kinetics, low Damkohler numbers, In addition, Bahadori et al.[131, 132] have conducted some test simulations mostly for laminar two-dimensional (2D) diffusion flames. These simulations are for support of their-extensive laboratory experiments which constitute the primary component of their research. In addition to the NASA sponsored research, there are several other efforts devoted to computational investigation of buoyant turbulent diffusion flames. Most of these are via either DNS (e.g. Refs.[124, 126, 133, 134, 135, 136, 137]) or the more conventional “Reynolds-averaged” simulations (RAS)[138, 139, 140, 141]. In this work, we would like to further assess the applicability and the extent of validity of the FDF methods for SGS closure of low-speed methane jet flames. This is conducted by implementation of the joint scalar FMDF for prediction of the gravity effects in partially premixed flames. The importance of gravity in such flames is well recognized[119, 121, 123, 125, 127]. It is now understood that the flame structure is in- directly affected by gravity through the modification of the large scale flow structures. Perhaps not as well recognized is the impact that the gravity has on the behavior of the flame and trubulence at small scales. In normal gravity, the large scale flow struc- tures are caused by the inherent instability of the jet coupled with the instability of 105 the buoyancy induced natural convection layer[125, 127]. In reduced gravity, the in- fluence of the second instability mechanism is decreased, resulting in modifications of the flame structure. Here, we examine some of the influences of gravity on the spatial and the compositional structures of transitional and turbulent hydrocarbon flames via LES / FMDF. We will consider gas-jet diffusion flames which are close representa- tions of those considered experimentally in both reduced gravity[125, 127, 132] and normal gravity conditions[33, 115]. The values of the simulation parameteres such as the Reynolds and Damkohler numbers will be in the same range as those in these experiments. The same goes for the Froude number, and because of computational flexibility a wide range of this number can be considered. It will be shown that the effects of gravity also appear in a closed form in the FMDF transport equation (in both formulations). The FMDF simulations will also facilitate an effective means of assessing the influence of gravity on “the compositional structure” of the flame. 3.2 Formulation The governing equations used for the simulations of this chapter are the filtered con- tinuity, momentum, and scalar equation, discussed in (§2.2). The transport equation for FMDF is solved by the Lagrangian Monte-Carlo procedure (§2.3.1). The Monte- Carlo solution is coupled to the finite difference solution through the reaction source term. 3.3 Numerical Simulation Procedure Details of the numerical simulation procedure can be found in §2.3. Here we briefly discuss the issues related to the chemistry models. 106 3.3.1 Implementation of chemistry As previously mentioned, two different approaches have been taken in the implemen- tation of the chemistry: (i) a direct approach for non-equilibrium flames, and (ii) a flamelet approach for near-equilibrium flames. In (i), fuel oxidation is modeled via finite-rate reduced mechanisms. In (ii), the flamelet library obtained from a laminar counterflow (opposed jet) flame configuration characterizes the reaction. The library is constructed by a detailed (GRI) methane oxidation mechanism. In method (i), Eq. 2.15, is used to determine the composition of the Monte-Carlo elements in FMDF using an appropriate reaction mechanism, while in method (ii), only the conserved scalar equation is calculated via FMDF. Method (i) is more expensive, specially if a complex kinetics mechanism (6.9. a reduced 10—step or 12-step mechanism) is em- ployed. The flamelet tables used in method (ii) are often based on a the most complete kinetics mechanism available. However, finite-rate effects are not represented in this method. In method (i), sensible enthalpy (h,) is computed as one of the scalars in Eq. (2.4), while in (ii), RT equation is calculated. The transport equation of the RT equation is obtained by multiplying Eq. (2.14) by RT, and integrating over the scalar space. This transport equation is different from Eq. (2.4), only in the source term. The source terms used in each method is shown in table 1. The solution of the finite difference scalar equation provides the pressure through the equation of state. The temperature could be obtained by averaging the temperatures of the Monte-Carlo el- ements, instead of using Eq. (2.15). However, as shown in Ref. [2], the pressure field obtained by direct averaging of the Monte-Carlo solution is noisy and would destabi- lize our fully-compressible numerical solver. On the other hand, stable and virtually noise-free solutions are obtained when the temperature (and consequently the pres- sure field) are computed by the finite difference solution of h, or RT equation, while the only interaction with the Monte-Carlo elements is through the chemical source terms. Again, these source terms are presented in table 1. 107 3.4 Results A schematic representation of the simulations performed is shown in Fig. 3.1. Four simulations of a round, piloted, turbulent methane jet were performed using LES / FMDF. Some of the parameters of these simulations are listed in table 3.6. Two different grav- ity conditions, zero-gravity (O-g ) and normal-gravity (1-g ) were considered and both the 1-step and flamelet kinetics were used for each of these gravity conditions. Except for the gravity and the chemistry, all the other flow quantities are the same. For these simulations, a PD mesh with 176 x 141 x 141 grid points was considered for a domain of 22 x 14 x 14 jet diameters in the x, y, and 2 directions, respectively. The number of Monte-Carlo elements for each F.D grid point was 10 and AB = 2A, resulting in about 12 million particles. The main jet composition is 25% CH4 and 75% air with a Reynolds number of 7534. The configuration of the jet is similar to the Sandia’s piloted jet flames [115], and the inlet conditions are similar, scaled relative to the jet velocity. The main parameters are as follows: 0 Nozzle diameter = 28.8mm 0 Pilot diameter = 72.8mm 0 Main jet: 25%CH4 and 75%air 0 Avg. jet speed: 4.13m/s 0 Peak jet speed: 5.26m/s 0 Main jet temperature: 300K 0 Pilot speed: 0.95m/s 0 Pilot temperature: 1880K 0 Pilot mixture fraction: 0.27 o Co—flow speed: 0.075m/s o Reynolds number: 7534 o Froude number: 60 108 The Reynolds number is defined as prUrLr/ur, where Ur, and L,, are the reference velocity and length, respectively, taken equal to the bulk velocity and diameter of the jet. Reference density, p,, and viscosity, u, are the bulk density and viscosity of the jet stream. The Froude number is defined as F 1‘ = Uf/gLr, where g is the gravitational acceleration. Again, it is emphasized here that the mixture fraction is solved for in the Monte-Carlo calculations, and there is no need to calculate it from the mass fractions. However, Bilger’s definition [115] was used to calculate the stoichiometric mixture fraction (Zst = 0.351): (Tu/170,11 '— Y2H)) + (-,,3—C(Yc — Yzcll Z = (mu/1H - Yam + (“Emma — Yzcll (3.1) where YH is the Hydrogen element mass fraction, YC is the Carbon element mass fraction and subscripts 1 and 2 refer to the main jet and the coflow respectively. The atomic weights of H and C are represented by WH and WC, respectively. 3.4.1 Physical flow structure Two important issues are investigated here: (1) the effect of gravity on transitional and turbulent jet flames, (2) the performance of equilibrium (flamelet) and finite-rate 1-step chemistry models in low-speed buoyant and nonbuoyant jet flames. Figures 3.2, 3.3 show the isosurfaces of vorticity magnitude (0 = wad, = 1) for O-g case 1 and 1—g case 3, respectively. In the O—g case (Fig. 3.2) the flow is similar to a laminar jet at the potential core (x/ D < 5) as expected, but after the break up of the vortex ring, it gradually evolves toward a fully developed turbulent jet. The ring like structure seen in the inlet is the vorticity introduced at the interface between the pilot and the coflow. The 1-g jet becomes unstable at distances very close to the jet exit. Unlike the O—g jet, the initial vorticity introduced at the pilot does not disappear quickly and merges with the main jet (better seen in Fig. 3.5). After z/ D > 14 the jet 109 grows rapidly and shows signs of smaller scale turbulence, not seen in Fig. 3.2 for the O—g jet flame. As the results in Figs. 3.2 and 3.3 indicate, when gravity is present, the combustion induced buoyant accelerations assist the production of turbulence, despite the damping resulted from the heat release. The temperature contours (not shown) reflect a similar trend. The DNS results for planar jets in chapter 1 also exhibit similar trends. The animated sequence of vorticity and temperature contours for the 1-g case suggests that the flame flickering originates at distances close to the jet exit. The instabilities propagate downstream in the form of periodic expansions and contractions of the flow, until they exit the domain. The 2—D contours of vorticity magnitude at the y = 0 cross section for cases 1 and 3 are shown in Figs. 3.4, and 3.5. The vorticity field in the O-g case is clearly damped by the heat release, except in some small regions of the flow. In the O-g case 1, the vorticity at the jet inlet is due to the added velocity perturbations. However, the inlet flow perturbations are significantly damped, specially at the pilot region. The 1-g case exhibits a very different trend. As shown in Fig. 3.3, the vorticity from the pilot region merges with that of the jet. At locations far from the inlet, the flow becomes completely “turbulent” and regions of high vorticity magnitude are evident. To better understand the changes in vorticity and its source, the baroclinic term appearing in the vorticity magnitude transport equation [142], B = WlfwlL ' (V(P>t x prlLll- is calculated and compared for the O-g and 1-g cases. The instantaneous values of B exhibit a sporadic pattern in the flow domain. Even the integrated values of this term in planes normal to the flow direction B, = f dedz are not very smooth. Fig. 3.6 shows the cumulative integral of B I plotted versus the axial distance from the nozzle. The quantities used in the definition of B are the resolved quantities and the vorticity vector is calculated from the resolved velocity vector, (w)L = V x (V) L. Naturally, there are subgrid contributions to this term, which are not considered here. It is 110 shown in Fig. 3.6 that the baroclinic vorticity generation term at the early phases of the jet development for the O-g flames and at (117/ D ~ 8) there is some vorticity generation which does not continue afterwards. For the 1—g jets, there is a steady generation of vorticity magnitude starting from very early phases (a: / D ~ 2.5), consistent with Figs. 3.3, 3.5. This steady trend does not continue after (a: / D N 8— 10), and there are periods of generation and destruction at this point. It should be noted that these values are instantaneous values and are not suitable for a quantitative comparison, but they can reflect the trend in vorticity generation. Contours of mixture fraction in case 3 obtained by the finite-difference and the Monte-Carlo particle methods are shown in Figs. 3.7 and 3.8, respectively. Similar results for case 4 are shown in Figs. 3.9 and 3.10. The Monte-Carlo results are ob- tained by ensemble averaging of the particle values. The finite-difference calculations of the mixture fraction are only performed for comparison with the Monte—Carlo re— sults and are not necessary. Also, for the l—step, 1—g flame (case 3), mixture fraction is not required in the particle calculations, and it is included for consistency check. A comparison between the finite-difference and the Monte-Carlo results shows a general agreement and consistency. At farther downstream locations, coinciding with regions of higher vorticity, signs of numerical oscillations are evident in the finite-difference solution (compare Figs. 3.9 and 3.10). The results obtained from the Monte-Carlo particles do not suffer from this problem. 3.4.2 Statistical description A general quantitative analysis is provided in this section, mainly by studying the time-averaged statistics. In order to evaluate the resolution of the finite difference solution, the resolved ( )2 — W2 and the total variance of scalar ((152) — (05)2 are calculated and shown in 111 figures 3.11 and 3.12. The resolved part was obtained by the finite difference solution, and the subgrid contribution was found from the Monte-Carlo particles and added to the resolved part. The profiles from the O-g simulations are close to each other, except at the shear layer where the subgrid contribution seems to be more than 30 percent. This region has the highest gradients of velocity and scalar in the flow. Since the grid used for ensemble averaging is double the size of the finite difference grid in each direction, the differences are most pronounced in regions with large gradients. A ‘good’ LES grid is the one that allows the'subgrid contribution of the turbulent energy to be less than 20% of the total energy. Based on this criterion, the grid size in the turbulent regions of the flow in case 4 is in the acceptable range. As shown in Fig. 3.14, the velocity fluctuations peak around :1: / D = 0.8 in the 1-g case, while the largest subgrid contributions are located around x/ D = 0.4. This also confirms that the higher differences are caused by sharper gradients in the mean profiles, rather than the stronger turbulent fluctuations. For smaller domain sizes, simulations with an ensemble averaging grid size (AB), equal to the finite difference grid were performed and it was observed that the subgrid contribution does not exceed 10%. Moreover, the first order moments such as the heat release were not affected by AB. Other low-order flow quantities are also not expected to change by the choice of the ensemble averaging domain size. Fig. 3.13 shows the mean axial velocity profiles for case 1 and case 3 at x/ D = 10, 18. Case 3 (1-g , l-step) has higher peak velocities compared to case 1 (O-g , 1-step) because of the buoyant acceleration and larger jet growth rate (see Fig. 3.3). The peak velocity does not change between at = 100 and a: = 18D in either of these cases, but the jet width is higher in case 3 at the same interval. Fig. 3.14 shows the axial velocity fluctuations for cases 1 and 3. Evidently, the velocity fluctuations grow substantially between :2: = IOD and x = 18D in the 1-g simulation, while for the same interval, the O-g jet shows negligible growth. This 112 suggests that the heat release has prevented the O-g jet from transition to a fully turbulent jet. Moreover, while the Reynolds numbers are higher in LES (Re = 7532) compared to our planar jet DNS (Re = 3000), the jet does not become a fully turbu- lent jet (even at :r = 18D), unless the buoyancy generated instabilities are added. In addition to the Reynolds number, the chemistry and the inlet boundary conditions have significant effects on the jet instability. For example, the damping of the flow instabilities in our simulations is expected to be affected by the pilot at the jet inlet. Although the pilot has only a small share of the total jet heat release, the presence of a high temperature region in the inlet slows down the development of turbulence in that vicinity, due to an increase in viscosity and changes in the density field. The integrated value of the axial velocity fluctuations over cross sections normal to the flow are shown in Fig. 3.15. The profiles from the O-g simulations (cases 1 and 2) are very close. The profiles from the l-g simulations (cases 3 and 4) are very close as well; indicating that the flow field is not very sensitive to the chemical mechanism. While both the O—g and 1-g simulations show a linear increase in the average RMS values in the flow direction, the rate of increase is substantially higher in the 1-g jets. Fig. 3.16 shows the instantaneous value of total kinetic energy in O-g and 1-g jets as simulated by the 1-step chemistry model (cases 1 and case 3). Both jets follow a similar trend for kinetic energy up to :1: / D as 12, followed by a sudden increase in the 1-g jet. The locations where the kinetic energy rises and falls in the 1-g case coincides with the locations that large scale structures are formed and propagated downstream in the flowm, due to flame flickering. Time averaged profiles of temperature are shown in Figs. 3.17 and 3.18 for the O-g and l-g cases respectively. The temperatures are obtained by time averaging the ensemble averaged particle temperatures. The particle temperatures would provide a more accurate description of the temperature, compared to the one given by the finite difference solution. The difference is partly because of using an ensemble do— 113 main length twice as long as the finite difference grid spacing, and partly due to the numerical errors of the finite difference scheme. Nevertheless, the values obtained from the particle and the finite difference solutions are within a few percent of each other. A comparison between Figs. 3.17 and 3.18 reveals the higher peak of tem- perature at the O-g jets, consistent with the experiments and our DNS simulations (see the results in chapter 1). When the profiles of case 1, and case 2 are compared, it is observed that even at close distances such as x/ D = 5, there are some differ- ences between the temperature profiles. The differences become considerably larger, farther downstream, with the l-step model always predicting higher temperatures. There are substantially smaller differences in the 1-g jets and as figures 3.17 and 3.18 indicate, at x/ D = 5 there is very little difference. At z/ D = 18, the flamelet model predicts slightly higher temperatures at the coflow side. Although different chemical mechanisms have been used for the flamelet and 1-step simulations, it was expected that there would be little difference for “high Damkohler number flames”, such as the ones considered here. The difference between the temperature predictions of the 1-step and flamelet flames would be better explained by comparing the heat release profiles. Figs. 3.19 and 3.20 show the time averaged heat release for the O—g and 1-g simu- lations, respectively. There is a substantial difference between the heat release values in the O-g simulations. This difference is occurring in the fuel side of the flow. In the l-step jet, the centerline heat release is initially low but increases in the downstream direction (x/ D = 18). The flamelet simulation shows comparatively less heat release at the same location. In the 1-g cases 3, and 4, , there are much less differences in the profiles (Fig. 3.20); at 33/ D < 5, there are virtually no differences. There are noticable differences at 23/ D = 10 and :c/ D = 18, but compared to the O-g cases, the agreement is much better. At x/ D = 18, the heat release is larger at the oxidizer side, consistent with 114 higher temperatures seen in Fig. 3.18. Fig. 3.21 shows the integrated values of the heat release along the axial direction. Between the O-g cases, case 1 has higher heat release compared to case 2 before :r/D < 18. At 113/D z 18, cases 1 and 2 have equal heat release (Fig. 3.21), while the surface areas under the curves of Fig. 3.19 indicate a non-equal heat release. The radial profiles of heat release are obtained by integration in the circumferential direction and since the jet is non-symmetric in case 2, it causes a difference in the heat release values. This problem is less pronounced in the l-g jets, since the heat release results shown in Fig. 3.19 as obtained from the flamelet simulation is higher than that of the 1-step simulation. Moreover, the integrated values of heat release in case 3, and 4 are very close to each other for a major part of the domain length, consistent with the results shown in Fig. 3.20. While the differences in the temperature profiles are consistent with the differences in heat release values, there is a fundamental difference between the two approaches in terms of the chemistry. The time averaged mass fraction of the reactants and the products would clarify this difference. For example, Fig. 3.22 shows the time averaged values of the mass fraction of CH4 for cases 1 and 2 at both a: / D = 10 and a: / D = 18 for the O-g jet flames. It is clear that the l-step results are higher than the flamelet ones. Fig. 3.23 shows a similar trend for Lg jet flames, although the difference is smaller than that in the O—g jets. Figs. 3.24 and 3.25 show the 002 mass fraction profiles for the O-g and 1-g flames. Although 002 is a product, it still reflects the trend seen for the fuel CH4. However, for 002 the difference between the 1-step and the flamelet simulations is not small for the 1-g flames. Moreover, the peak 0'02 value is lower compared to that in the 0g jets. Fig. 3.26 shows the mass fractions of CO from the flamelet simulations case 2, and 115 4. The 1—step mechanism does not include CO. The CO values for the O-g simulation are substantially larger compared to the 1-g jet. Figs. 3.27 and 3.28 show the mass flow rate of CH4 for case 1, 2 and case 3, 4 respectively. Very little difference is seen between the flamelet and the 1-step simulations at :c/ D = 18. The integrated values of CH4 mass flow rate are shown in Fig. 3.29. The mass flow rate entering the domain in the 1-g calculations are slightly higher than the O-g calculations, due to the increased density in the inlet, caused by hydrostatic pressure. Despite the higher mass flow rate of CH4 into the domain in the 1-g cases, the exiting mass flow rate is lower, compared to the O—g cases, which is consistent with the higher values of heat release. The difference in heat release between case 1, and 2 is due to considerable production of CO in case 2. For the same reason, since there is much less CO production in case 4, the heat release contours for case 3, and 4 are closer. The increase in heat release in case 4 after a: / h ~ 13 is partly due to faster reaction (Fig. 3.27) and partly because of lower CO production at this stage of the jet. 3.5 Conclusions LES/FMDF of round methane jet flames are conducted to study the physical and compositional structure of transitional and turbulent jet flames and flame-turbulence interactions in zero-gravity (O-g ) and normal-gravity (1-g ) environments. The results of the simulations are consistent with the DNS results and reflect the same trend with regards to changes in gravity conditions. At O-g conditions combustion damps the flow instability and turbulence production, while in the 1-g conditions, combustion generated density variations and buoyancy effects significantly modifies the flow (turbulence) structure. The flamelet and the 1-step global mechanism were used to evaluate the effects 116 of reaction. In the absence of gravity, the two different models predict equal fuel consumption rates, suggesting a near equilibrium flame behavior in the zero-gravity conditions. However, at the presence of gravity, the finite rate chemistry effects become important and the flamelet model is less accurate. 117 3.6 Figures and Tables Figure 3.1: Schematic representation of the LES/FMDF simulations performed for buoyant / nonbuoyant methane jet flames. Nozzle diameter = 28.8mm, Pilot diameter = 72.8mm. 118 "V Method chemistry scalar source term Direct reduced mechanisms hs —d20h2 Flamelet complete mechanism RT ——(p);Qm[( f G) L — (f) L(G) L] Table 3.1: The chemistry, the form of “energy equation” and the source term used in the direct and flamelet approaches. f is the conserved scalar and G is BRT / 8 f . Case gravity chemistry 1 O-g l-step 2 O-g flamelet 3 1-g 1-step 4 lg flamelet Table 3.2: Identification of the simulations. 119 \s Figure 3.2: Isosurface of vorticity magnitude ((2 = 1) for case 1 (O-g , 1-step). The domain shown is 0 < I/D < 20, —5 < y/D < 5, —5 < z/D < 5. \§ Figure 3.3: Isosurface of vorticity magnitude (9 = 1) for case 3 (l-g , 1—step). The domain shown is 0 < x/D < 20, -—5 < y/D < 5, —5 < z/D < 5. 120 -. Figure 3.4: Contours of vorticity magnitude (0-10) for case 1 (O—g , 1-step). The domain shown is 0 < x/D < 20,y = 0, —5 < z/D < 5. 121 Figure 3.5: Contours of vorticity magnitude (0—10) for case 3 (l-g , 1-step). The domain shown is 0 < :c/D < 20,y = 0, -5 < z/D < 5. 122 12 . , . , , I , o. l ,r‘ - _ 8v 'SICP \ . —-—— 0-g,flamelet I.) ‘ I, \ — l-g, l-step . ,- I ‘ / \ 9 - — - l-g. flamelct I \\f I v \I ‘ _ I \l \ ~ I \A‘ - I H I I ] “=2. 6 ~ | ' CQ Lo - / V1 / \ V/ [ 3 ’— ’,/‘\ ________ .1 I f\ _______________________ ' II I . fi‘ fi._\ ’..\\ 1”“ x ’1’ O \11 ' - l ‘T 1 l J l 1 0 5 10 15 20 z/ D Figure 3.6: Cumulative integral of vorticity magnitude generation 8,. 123 Figure 3.7: Contours of mixture fraction from the finite difference solution for case 3 (l—g , 1-step). The domain shown is O < x/D < 20,y = 0, —5 < z/D < 5. 124 Figure 3.8: Contours of mixture fraction from the particles for case 3 (1-g , 1—step). The domain shown is 0 < x/D < 20,y = 0, —5 < z/D < 5. 125 Figure 3.9: Contours of mixture fraction from the finite difference solution for case 4 (1-g , flamelet). The domain shown is 0 < x/D < 20,y = 0, —5 < z/D < 5. 126 Figure 3.10: Contours of mixture fraction from the particles for case 4 (1-g , flamelet). The domain shown is 0 < :r/D < 20,y = 0, —5 < z/D < 5. 127 0.04 0.02 resolved and total variance Figure 3.11: Resolved and total scalar variance for case 2 (O—g , flamelet). resolved and total variance Figure 3.12: Resolved and total scalar variance for case 4 (1-g , flamelet). x=15D, resolved x=18D, resolved x=ISD. total x=l8D, total 128 . x=15D, resolved x=l8D, resolved _ x=15D, total x=18D, total 6 I I I T I ’ ---- 0-g,x=lOD * 5 — 0-g,x=l8D _ —— l-g,x=lOD - — l-g,x=lSD 4 __ As 3 — 3 2 F- 1— O I l 1 FT ‘1- 0 1 2 3 r/D Figure 3.13: Mean axial velocity for case 1 (O-g , 1-step) and case 3 (l-g , l-step). 2 r I I I I ---- Oog. x=10D — O-g, x=18D — — l-g, x=10D — l-g,x=18D 1.5 ’3 A 3 1 v E Q3. 0.5 p "—— \ ~~~~~~ \_ ______ ‘ ffffff t-~-------- ---__-_--- O m l 1 + O 1 2 3 T/D Figure 3.14: RMS of axial velocity for case 1 (O-g , l—step) and case 3 (lg , 1-step). 129 80 I f l l I l I — O-g, l-step —---— 0-g,flamelet / — l-g, l-step / 6O — — — l-g, flamelet Figure 3.15: Integrated values of RMS of axial velocity, for different cases. 130 10 I I I I I l I _ — O-g. l-step - — l-g, l-step 1 O 5 10 15 20 z/D Figure 3.16: Total instantaneous energy (normalized) for case 1 (O—g , 1-step) and case 3 (1-g , 1-step). 131 U I I I l o—o l-step, x=SD — - l-step. x=10D — l-step,x=18D o—o flamelet, x=5D ---- flamelet,x=lOD — flamelet. x=18D I 2000 I r/D Figure 3.17: Time averaged temperature profiles for O-g cases, case 1 (l-step) and case 2 (flamelet). 1 500 I I I I I o—o l-step. x=5D — — l-step, x=10D — l-step, x=18D - o—o flamelet, x=5D --—— flamelet, x=IOD — flamelet. x=18D 1000 - 5 ,‘i la ~ 500 _ 3 Figul‘e 3.18: Time averaged temperature profiles for Lg cases, case 3 (l-step) and case 4 (flamelet). 132 3o . , . , . — - l-step. x=10D _ — l-step, x=18D 25 __ --—- flamelet, x=10D I \ — flamelet,x=l8D l — Figure 3.19: Time averaged heat release for O-g cases, case 1 (1-step) and case 2 (flamelet). 40 I l I l I — - l-step, x=10D P -— l-step,x=18D \ ---- flamelet, x=10D 30 r- /. — flamelet, x=18D ~ / .\ I Figure 3.20: Time averaged heat release for Lg cases, case 3 (1-step) and case 4 flaIIlelet). 133 I j I I I I 8 _. — O-g, l-step I“ -—-- 0—g,flamelet / _ — l-g. l-step / q — — l-g, flamelet // x/D Figure 3.21: Integral of time averaged heat release for all cases. 134 0.2 1 , . , , — — l-step. x=IOD — l-step. x=lBD ---- flamelet, x=10D — flamelet, x=18D 0.15 Figure 3.22: Time averaged mass fraction of CH4 for 0-g cases, case 1 (l-step) and case 2 (flamelet). 0.2 . , . , . — — l-step, x=10D —— l-step. x=18D ---- flamelet, x=IOD 015 g —- flamelel, x=18D I Figure 3.23: Time averaged mass fraction of CH4 for Lg cases, case 3 (l-step) and Case 4 (flamelet). 135 0.2 I I l I I _ —- -- l-step, x=10D — l-step, x=l8D ---- flamelet, x=10D 0.15 _ —- flamelet. x=l8D . / \ . / \ q / - \ ’2 O 1 — / x" \ \ _ 8 l I” \\\ \ >"‘ I I, ‘ V 7 I ’1’ \\\ - I ’I \\\\ _ \ _ 0005 ll ’1’ \\\\\\\ I / \ b- I” \\§ \ d r ’ ,,I’/ “~\\\ 0 TTTT l l 1 1 ~‘_._ ‘ 0 1 2 3 T/ D Figure 3.24: Time averaged mass fraction of C02 for 0-g cases, case 1 (1-step) and case 2 (flamelet). f l I I I ’\ / _ 0.08 "' / \ — — l-step. x-lOD _ \ — l-step,x=l8D _ I \ ----' flamelel, x=10D . I —- flamelet, x=18D 0.06 — I \ _ I I,’ \\ \ Q '- I I, \\ \ " E [I I" \ 00.04 _ \\ _ E" I I,” \\ .. I III \\ \ . I I” \\\\\ 0.02 _ I I, \\\ _, ’ \ / / \\ I ll, \\\ .I r’III \‘\“\ O L l # r“§ ._—__ _L. 0 1 2 3 r/D Figure 3.25: Time averaged mass fraction of C02 for 1-g cases, case 3 (1-step) and C388 4 (flamelet). 136 0.06 . , , I . , ---- O-g. flamelet, x=10D - 005 - I’ \ — O-g, flamelet, x=18D _ — — l-g, flamelet, x=lOD — l-g, flamelet, x=l8D L .— r/D Figure 3.26: Time averaged mass fraction of CO for case 2 (O-g , flamelet) and case 4 (l-g , flamelet). 137 0.2 I I I I I _ — — l-step. x=10D —— l-step, x=18D ---- flamelel,x=10D 0.15 F — flamelet. x=l8D-I i . E? 2;, 0.1 a ,‘i 3 . 3 0.05 — 0 l 2 3 Figure 3.27: Time averaged pal/CH, for 0-g cases, case 1 (l—step) and case 2 (flamelet). 0.2 T I I l I _ —— l-step. x=10D _ — l-step,x=18D L ---- flamelet,x=10D 0.15 — flamelet,x=18Dr \ ,‘i . S D 23, 0.1 — :1 a -I :5 3 0.05 ‘* Figure 3.28: Time averaged [mi/CH, for l—g cases, case 3 (1-step) and case 4 (flamelet). 138 El! .,_._."J~_‘. 0.4'11 4 ‘ I —— O-g, l-step ---- O-g,flamelet —— l-g, l-step F — — l-g. flamelet f

2(U>L(YCH4)Ld3/dz 0.05 ‘ ' 1 10 15 20 x/D Figure 3.29: Integrated values of puYCH‘ for different cases. 139 3.7 Nomenclature Acme-90°C) amass DER”, .Q m '6 EU «23" :0 o : The model in the VSFDF : Specific heat of species a at constant pressure : Specific heat of the mixture at constant pressure : Constant in the VSFDF formulation : Constant in the VSFDF formulation : Constant in the subgrid scalar stress : Constant in the subgrid scalar stress : Constant of the stochastic mixing closure : Molecular diffusion coefficient : Subgrid diffusion coefficient : Drift coefficient in the SDE : Inner jet diameter : Diffusion coefficient in the SDE : Joint scalars filtered mass density function : Joint velocity-scalars filtered mass density function : Froude number : Normal gravity : Gravity vector : ith component of the gravity vector : Filter function : Filter function : Enthalpy : Enthalpy of species a : Enthalpy of formation of species a : Jet momentum : ith component of the flux of scalar a : Thermal conductivity of the mixture : Luminous flame height : Molecular Lewis number ‘ : Lateral dimension of the flame : ith component of the subgrid scalar flux of species a : Number of species : Pressure : Molecular Prandtl number : Universal gas constant : Velocity ratio at the inlet r = %Q“ m : Reynolds number : Richardson number : SGS Schmidt number : Production rate of species a : Strain rate tensor 140 flaw scrs E §§s<§g N§5Haex : Temperature : Reference temperature : Subgrid scale stresses : Time : Velocity vector : Streamwise velocity at the jet inlet : ith component of the velocity vector : ith component of the stochastic velocity vector : ith component of the reference velocity in MKEV closure : Velocity space : ith component of the velocity vector in composition space : Molecular weight of species a : Wiener-Levy process : Position vector : ith component of the position vector : Lagrangian position of the particles : Streamwise coordinate ' : Mass fractions of species a : Coordinates defining the plane normal to x : Mixture fraction Greek Symbols: 7,- : Standardized Gaussian random variable A : Grid spacing in LES 6 : Dirac delta function 5ij : Kronecker delta AH : Grid level filter width AH: : Secondary level filter width AU : Velocity difference at the inlet, AU :2 Um — Um. Ap = pin - pf e : SGS dissipation u : Molecular viscosity, ,1 = pu Vt : Subgrid viscosity Qm : SGS mixing frequency 1p : Composition space ()5 : Scalar field qba : The compositional values of scalar a (I): : The compositional values of stochastic scalar a p : Density 0 : Number of scalars, 0 = N8 + 1 73-]- : Molecular stresses Symbols 141 ( )3 : Filtered value ( )L : Favre filtered value ( I )L : Conditional Favre filtered value A : Quantities which depend only on the scalar composition, zle. §(¢,x, t) E S(¢(x, t)) oo : Ambient Subscripts f : Flame in : Inner jet 1: : Discretized (time) level on : Outer jet Superscripts + : Properties of the stochastic Monte Carlo particles (n) : Index of the Monte Carlo particles 142 Chapter 4: Conclusions and Summary The influence of gravity on turbulent jet flames is investigated via two computational approaches. In the first approach, three-dimensional (3D) direct numerical simula- tions (DNS) of diffusion planar jet flames with an ideal one-step reaction model are performed. DNS of nonreacting planar jets were also performed for validation pur- poses, where good agreement was achieved compared to the available experimental and DNS data. Buoyant and non-buoyant 2-D methane flame simulations with realis- tic thermodynamic properties and a global one-step kinetics model are also performed. The results of our reacting jet simulations are consistent with previous findings and indicate that in the absence of gravity, combustion damps the flow instability; hence reduces “turbulence production” and jet growth. However, in the “finite—gravity” conditions, combustion generated density variations and buoyancy effects, promote vorticity generation and enhance the turbulent mixing and combustion at nearly all length and time scales. As expected, buoyancy was found to have a destabilizing effect on laminar jets, and the formation of large-scale flame flickering was consistent with previous investigations. However, the gravity effects on transitional/ turbulent jet flames are not limited to large-scales, and the turbulence and mixing were affected in the small scales as well. Comparison between the 2—D and 3-D simulations suggests that the effects of combustion and gravity on turbulence (and vice versa) cannot be appropriately evaluated by 2-D computational models when the turbulence is more developed. However, the early effects during transition are adequately captured by the 2—D simulations. The gravity effects are also explained via the analysis of the compositional struc- ture of the flame. It is shown that the finite-rate chemistry effects become more 143 significant in the buoyant cases due to the increased turbulence and the higher strain rates. When the turbulence field is damped by the heat release, the addition of grav- ity could help the overall production by destabilizing the flow, resulting in improved mixing and an extended flame surface. On the other hand, when the turbulence field is “strong”, the addition of gravity may induce excessive amounts of strain and cooling of the flame, which in turn adversely affects the combustion. In the second approach the recently developed large eddy simulation (LES) method- ology, termed (FMDF) is extended towards the simulation of realistic methane flames, such as Sandia’s piloted turbulent methane jet flames. The Sandia flames are cho— 4 sen for validation purposes, since extensive experimental data are available for them. a F Both the flamelet and finite rate kinetics approaches were considered for the 3-D round jet simulations, with the GRI-2.01 for the former and the 1-step and 12-step mechanisms for the latter. The consistency of the Monte-Carlo and finite difference solutions were discussed and a scalable algorithm for parallelization of our hybrid (Eulerian-Lagrangian) methodology is presented. The algorithm was shown to yield superlinear speedups, when single processor nodes are used. The Favre-averaged temperatures for two jet speeds (the so called flames D and F) were compared with the experimental data and good agreement was seen for the lower jet speed (flame D), which has less non-equilibrium effects. Higher degrees of local extinction in the higher jet speed (flame F) causes larger differences between the predicted and mea- sured temperatures; suggesting that a more accurate finite-rate reduced chemistry model is needed in this case. In the multi-step simulations, it was seen that after :r/ D = 10, the flame exhibits substantial extinction, even at low-speed (flame-D) conditions. This extinction does not lead to a total blow-out, but rather causes a drop in the mean temperature, and heat release. The factor that mainly contributes to this problem is identified as the inadequate cartesian grid resolution. Resolving the sharp profiles of the velocity and 144 composition at the inlet is very important. The sharp gradients present at the inlet, compound the problem of an incompatible geometry. Employing larger number of grid points would be prohibitive, in terms of cost. Although the resolution at the domain inlet might not be sufficient to resolve the sharp gradients, the turbulence developed at farther downstream positions is adequateley resolved, as shown in the scalar variance profiles. In the final part of this dissertation, the gravity effects on low-speed round methane jet flames were studied using the LES/FMDF method. The methane jet flames were investigated under O-g and l-g environments. The LES results are consistent with the DNS results and reflect the same trend with regards to changes in gravity conditions. At O-g conditions combustion damps the flow instability and turbulence production, while in the l-g conditions, combustion generated density variations and buoyancy effects significantly modifies the flow (turbulence) structure. The flamelet and the 1-step global mechanism were used to evaluate the reaction source term. It was seen that the fuel consumption rates are equal in the two ap- proaches in the absence of gravity, suggesting a near-equilibrium flame in this case. Due to the lower amounts of extinction in these flames, the chemical mechanisms considered here did not indicate a substantial difference in the results. However, in normal-gravity conditions, the finite rate chemistry effects become important and the flamelet model is not accurate. 145 Bibliography [1] Edelman, R. and Bahadori, M., Effects of Buoyancy on Gas—Jet Diffusion Flames: Experiment and Theory, Acta Astronautica, 13:681—686 (1986). [2] Jaberi, F., Colucci, R, James, S., Givi, P., and Pope, S., Filtered Mass density function for large-eddy simulation of turbulent reacting flows, J. Fluid Mech., 401285—121 (1999). [3] James, S. and Jaberi, F., Large Scale simulations of two dimensional non- premixed methane jet flames, Combustion and flame. [4] Jaberi, F., Livescu, D., and Madnia, C., Characteristics of chemically react- ing compressible homogeneous turbulence, Physics of Fluids, 12(5):1189—1209 (2000). [5] Clemens, N. and Paul, P., Effects of Heat Release on the Near Field Flow Structure of Hydrogen Jet Diffusion Flames, Combustion and Flame, 102:271— 284 (1995). [6] Savas, O. and Golahalli, 8., Flow Structure in Near-Nozzle Region of Gas Jet Flames, AIAA Journal, 24(7):1137—1140 (1986). [7] Rehm, J. and Clemens, N., The Large-Scale Turbulent Structure of Non- premixed Planar Jet Flames, Combustion and Flame, 1162615-626 (1999). [8] Ellzey, J. and Oran, E., Effects of Heat Release and Gravity on an Unsteady Diffusion flame, 23rd Symposium (Intl.) on Combustion/ The Combustion In- stitute, pp. 1635—1640 (1990). [9] Bahadori, M. Y., Stocker, D. P., Vaughan, D. F., Zhou, L., and Edelman, R. B., Effects of Buoyancy on Laminar, Transitional, and Turbulent Gas Jet Diffusion Flames, in Williams, F. A., Oppenheim, A. K., Olfe, D. B., and Lapp, M., editors, Modem Developments in Energy, Combustion and Spectroscopy, pp. 49~66, Pergamon Press, New York, 1993. [10] U.Hegde and Bahadori, M., The transition to turbulance of microgravity gas jet diffusion flames, Combustion Science and Tech., (1994). [11] Azzoni, R., Ratti, S., Puri, I., and Aggarwal, S. K., Gravity effects on triple flames: flame structure and flow instability, Physics of Fluids, 11(11):3449- (1999). [12] Ellzey, J., Laskey, K., and Oran, E., A Study of Confined Diffusion Flames, Combustion and Flame, 84:249—264 (1991). [13] Davis, R., Moore, E., Roquemore, W., Chen, L., Vilimpoc, V., and G053, L., Preliminary results of a numerical-experimental study of the dynamic structure of a buoyant jet diffusion flame, Combustion and Flame, 83:263—270 (1991). 146 [14] Kaplan, C., Baek, S., and Oran, E., Dynamics of an Unsteady Buoyant ethylene jet diffusion flame, AIAA paper no. 93-0109, (1993). [15] Shu, Z., Aggarwal, S., Katta, V., and Puri, I., Flame-Vortex Dynamics in an Inverse partially premixed combustor: the Froude number effects, Combustion and Flame, (1997). [16] Elghobashi, S., Zhong, R., and Boratav, 0., Effects of gravity on turbulent nonpremixed flames, Physics of fluids, 11(10):3123—3135 (1999), 95. [17] Stanley, S. and Sarkar, 8., Direct Numerical Simulation of the Developing region of turbulent planar jets, AIAA, (1999). [18] Stanley, 8., Sarkar, S., and Mellado, J ., A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation, J. Fluid Mech., 450:377—407 (2002). [19] Gutmark, E. and Wygnanski, I., The planar turbulent jet, J. Fluid Mech., 73:465—495 (1976). [20] Givi, P., Model Free Simulations of Turbulent Reactive Flows, Prog. Energy Combust. Sci, 15:1—107 (1989). [21] Bui—Pham, M., Studies in Structures of Laminar Hydrocarbon Flames, PhD Thesis, UCSD, (1992). [22] Kennedy, C. and Carpenter, M., Comparison of Several numerical methods for simulation of compressible shear layers, NASA Technical paper. [23] Lele, S. K., Compact finite difference schemes with spectral-like resolution, J. Comp. Phys, (1992). [24] Poinsot, T. and Lele, S., Boundary Conditions for Direct simulations of Com- pressible viscous flows, Journal of Computational physics, (1992). [25] Rudy, D. and Strikwerda, J ., Boundary conditions for subsonic compressible Navier-Stokes calculations, Computers and Fluids, (1981). [26] Michalke, A., On Spatially Growing Disturbances in an Inviscid Shear Layer, J. Fluid Mech., 23:521—544 (1965). [27] Pope, S. B., Turbulent Flows, Cambridge University Press, Cambridge, UK, 2000. [28] Hussain, A. and Clark, A., Upstream Influence on the near field of a plane turbulent jet, Physics of Fluids, 20(9):1416—1426 (1977). [29] Ramaprian, B. and Chandrasekhara, M., LDA Measurements in plane turbulent jets, ASME J. Fluids Engineering, 107(3):264—71 (1985). 147 [30] Arai, M., Sato, H., and Amagai, K., Gravity effects on stability and flickering motion of diffusion flames, Combustion and Flame, (1999). [31] H.Sato, Amagai, K., and Mari, M., Diffusion flames and their flickering motions related with Froude numbers under various gravity levels, Combustion and Flame, (2000). [32] Katta, V. and Roquemore, W., Simulation of dynamic methane jet diffusion flames using finite rate chemistry models, AIAA J., 36(11) (1998). [33] Barlow, R. and Frank, J ., Effects of Turbulence on Species Mass Hactions in Methane/Air Jet Flames, 27th International symposium on combustion, pp. 1087—1095 (1998). [34] Xu, J. and Pope, S., PDF Calculations of Turbulent Nonpremixed Flames with Local Extinction, Combustion and Flame, 123:281—307 (2000). [35] Landenfeld, T., Sadiki, A., and Janicka, J., A Turbulence-Chemistry Interaction Model Based on a Multivariate Presumed Beta-PDF method for Turbulent Flames, Flow, Turbulence and Combustion, 68:111—135 (2002). [36] Pitsch, H. and Steiner, H., Large-eddy simulation of a turbulent piloted methane/ air diffusion flame (Sandia Flame D), Phys. Fluids, 12(10):2541—2554 (2000). [37] Roomina, M. and Bilger, R., Conditional moment closure (CMC) predictions of a turbulent methane-air jet flame, Combustion and Flame, 125:1176—1195 (2001). [38] Jones, W. and Kakhi, M., PDF modeling of Finite-rate Chemistry Effects in Turbulent Nonpremixed Jet Flames, Combustion and Flame, 115:210—229 (1998). [39] Masri, A., Bilger, R., and Dibble, R., Turbulent Nonpremixed Flames of Methane Near Extinction - Mean Structure from Raman Measurements, Com- bustion and Flame, 71:245—266 (1988). [40] Masri, A., Bilger, R., and Dibble, R., Conditional-Probability Density- Functions Measured in Turbulent Nonpremixed Flames of Methane Near Ex- tinction, Combustion and Flame, 74:267—297 (1988). [41] Sung, C., Law, C., and Chen, J., An Augmented Reduced Mechanism for Methane Oxidation with comprehensive Global Parametric Validation, 27th Symposium on Combustion/ The Combustion Institute, pp. 295—304 (1998). [42] Sung, C., Law, C., and Chen, J ., Further Validation of an Augmented Reduced Mechanism for Methane Oxidation : Comparison of Global Parameters and Detailed Structure, Combust. Sci. and Tech., 156:201—220 (2000). 148 [43] Pope, S., Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation, Combustion theory and modeling, 1:41—63 (1997). [44] Subramaniam, S. and Pope, S., A mixing model for turbulent reactive flows based on Euclidean minimum spanning trees, Combustion and Flame, 1152487— 514 (1998). [45] Muradoglu, M., Pope, S., and Caughey, D., The hybrid method for the PDF equations of Turbulent Reactive Flows: Consistency Conditions and Correction Algorithms, Journal of Computational Physics, 17 2:841—878 (2001). [46] Maas, U. and Pope, 8., Implementation of simplified chemical kinetics based on intrinsic low-dimensional manifolds, in 24th Symposium (Int. ) on Combustion, pp. 103—112, The Combustion Institute, Pittsburgh, PA, 1992. [47] Maas, U. and Pope, S., Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space, Combustion and Flame, 88:239—264 (1992). [48] Vervisch, L., Hauguel, R., Domingo, P., and Rullaud, M., Three facets of turbulent combustion modelling: DNS of premixed V-flame, LES of lifted non- premixed flame and RANS of jet-flame, J. Turbulence, 5(4):2191—2198 (2004). [49] Madnia, C. K. and Givi, P., Direct Numerical Simulation and Large Edddy Simulation of Reacting Homogeneous Turbulence, in Galperin, B. and Orszag, S. A., editors, Large Eddy Simulations of Complex Engineering and Geophysical Flows, chapter 15, pp. 315—346, Cambridge University Press, Cambridge, UK, 1993. [50] Schumann, U., Large Eddy Simulation of Turbulent Diffusion with Chemi- cal Reactions in the Convective Boundary Layer, Atmospheric Environment, 23(8):1713—1726 (1989). [51] Sykes, R. I., Henn, D. S., Parker, S. F., and Lewellen, W. S., Large-Eddy Simu- lation of a 'Ihrbulent Reacting Plume, Atmospheric Environment, 26(14):1713— 1726 (1992). [52] Liou, T. M., Lien, W. Y., and Hwang, P. W., Large-Eddy Simulations of Turbulent Reacting Flows in Chamber with Gaseous Ethylene Injecting through the Porous Wall, Combust. Flame, 99:591—600 (1994). [53] Boris, J. P., Grinstein, F. F., Oran, E. S., and Kolbe, R. L., New Insights into Large Eddy Simulations, NRL Report NRL/MR/4400-92-6979, Naval Research Laboratory, Washington, DC, 1992. [54] Williams, F. A., Combustion Theory, The Benjamin/Cummings Publishing Company, Menlo Park, CA, 2nd edition, 1985. 149 [55] Pope, S. B., Turbulent Flows, Cambridge University Press, Cambridge, UK, 2000. [56] Poinsot, T. and Veynante, D., editors, Theoretical and Numerical Combustion, R.T. Edwards Inc., Philadelphia, PA, 2001. [57] Peters, N., editor, Turbulent Combustion, Cambridge University Press, Cam- bridge, U.K., 2000. [58] Bilger, R. W., Future Progress in Turbulent Combustion Research, Prog. Energy Combust. Sci, 26:367—380 (2000). [59] Pope, S. B., Computations of Turbulent Combustion: Progress and Challenges, in Proceedings of 23rd Symp. (Int. ) on Combustion, pp. 591-612, The Combus- tion Institute, Pittsburgh, PA, 1990. [60] Pope, S. B., PDF Methods for Turbulent Reacting Flows, Prog. Energy Com- bust. Sci., 11:119—192 (1985). [61] Gao, F. and O’Brien, E. E., A Large Eddy Scheme for Turbulent Reacting Flows, Phys. Fluids A, 5(6):1282-1284 (1993). [62] Réveillon, J. and Vervisch, L., Subgrid—Scale Turbulent Micromixing: Dynamic Approach, AIAA J., 36(3):336—341 (1998). [63] Colucci, P. J ., Jaberi, F. A., Givi, P., and Pope, S. B., Filtered Density Func- tion For Large Eddy Simulation of Turbulent Reacting Flows, Phys. Fluids, 10(2):499—515 (1998). [64] Garrick, S. C., Jaberi, F. A., and Givi, P., Large Eddy Simulation of Scalar Transport in a Turbulent Jet Flow, in Knight, D. and Sakell, L., editors, Recent Advances in DNS and LES, Fluid Mechanics and its Applications, Vol. 54, pp. 155—166, Kluwer Academic Publishers, The Netherlands, 1999. [65] Jaberi, F. A., Colucci, P. J ., James, S., Givi, P., and Pope, S. B., Filtered Mass Density Function For Large Eddy Simulation of Turbulent Reacting Flows, J. Fluid Mech., 401:85-121 (1999). [66] Jaberi, F. A., Large Eddy Simulation of Turbulent Premixed Flames via F il- tered Mass Density Function, AIAA Paper 99-0199, AIAA, 1999. [67] Gicquel, L. Y. M., Givi, P., Jaberi, F. A., and Pope, S. B., Velocity Filtered Density Function for Large Eddy Simulation of 'Ihrbulent Flows, Phys. Fluids, 14:1196-1213 (2002). [68] Sheikhi, M. R. H., Drozda, T. G., Givi, P., and Pope, S. B., Velocity-scalar filtered density function for large eddy simulation of turbulent flows, Phys. Fluids, 15(8):2321-2337 (2003). 150 [69] Tong, C., Measurements of Conserved Scalar Filtered Density Function in a Turbulent Jet, Phys. Fluids, 13(10):2923—2937 (2001). [70] Darabiha, N., Transient—behavior of laminar counterflow Hydrogen air diffusion flames with complex chemistry, Combustion science and technology, 86:163—181 (1992). [71] Aldama, A. A., Filtering Techniques for Turbulent Flow Simulations, Lecture Notes in Engineering, Vol. 49, Springer-Verlag, New York, NY, 1990. [72] Rogallo, R. S. and Moin, P., Numerical Simulation of Turbulent Flow, Ann. Rev. Fluid Mech., 16:99—137 (1984). [73] O’Brien, E. E., The Probability Density Function (PDF) Approach to React- ing Turbulent Flows, in Libby, P. A. and Williams, F. A., editors, Turbulent Reacting Flows, chapter 5, pp. 185—218, Springer-Verlag, Heidelberg, 1980. [74] Vreman, B., Geurts, B., and Kuerten, H., Realizability Conditions for the Turbulent Stress Tensor in Large-Eddy Simulation, J. Fluid Mech., 278:351— 362 (1994). [75] Colucci, P. J ., Large Eddy Simulation of Turbulent Reactive Flows: Stochastic Representation of the Subgrid Scale Scalar Fluctuations, Ph.D. Thesis, Depart- ment of Mechanical and Aerospace Engineering, State University of New York at Buffalo, Buffalo, NY, 1998. [76] Eidson, T. M., Numerical Simulation of the Turbulent Rayleigh-Benard Prob— lem Using Subgrid Modelling, J. Fluid Mech., 158:245—268 (1985). [77] Moin, P., Squires, W., Cabot, W. H., and Lee, S., A Dynamic Subgrid-Scale Model for Compressible Turbulence and Scalar Transport, Phys. Fluids A, 3(11):2746—2757 (1991). [78] Lesieur, M. and Metais, 0., New Trends in Large Eddy Simulations of Turbu- lence, Ann. Rev. Fluid Mech., 28:45—82 (1996). [79] Pope, S. B., An Improved Turbulent Mixing Model, Combust. Sci. and Tech., 28:131—145 (1982). [80] Miller, R. S., Frankel, S. H., Madnia, C. K., and Givi, P., Johnson-Edgeworth Translation for Probability Modeling of Binary Mixing in Turbulent Flows, Combust. Sci. and Tech., 91(1-3):21-52 (1993). [81] Jaberi, F. A. and Givi, P., Inter-Layer Diffusion Model of Scalar Mixing in Homogeneous Turbulence, Combust. Sci. and Tech., 104(4-6):249 (1995). [82] Borghi, R., Turbulent Combustion Modeling, Prog. Energy Combust. Sci., 14:245—292 (1988). 151 [83] Madnia, C. K., Frankel, S. H., and Givi, P., Reactant Conversion in Homo- geneous Turbulence: Mathematical Modeling, Computational Validations and Practical Applications, Theoret. Comput. Fluid Dynamics, 4:79—93 (1992). [84] Kosaly, G. and Givi, P., Modeling of Turbulent Molecular Mixing, Combust. Flame, 70:101-118 (1987). [85] McMurtry, P. A. and Givi, P., Direct Numerical Simulations of Mixing and Reaction in a Nonpremixed Homogeneous Turbulent Flow, Combust. Flame, 77:171—185 (1989). [86] Fox, R. 0., Computational Methods For Turbulent Reacting Flows in Chemi- cal Process Industry, Revue De L’Institut Francois Du Petrole, 51(2):215—246 ( 1996). [87] Subramaniam, S. and Pope, S. B., A Mixing Model for 'Ihrbulent Reactive Flows Based on Euclidean Minimum Spanning Trees, Combustion and Flame, 1152487—514 (1998). [88] Subramaniam, S. and Pope, S. B., Comparison of Mixing Model Performance for Nonpremixed Turbulent Reactive Flow, Combust. Flame, 117:732—754 (1999). [89] Bardina, J ., Ferziger, J. H., and Reynolds, W. C., Improved Turbulence Models Based on Large Eddy Simulations of Homogeneous, Incompressible, Turbulent Flows, Department of Mechanical Engineering Report TF-19, Stanford Univer- sity, Stanford, CA, 1983. [90] Smagorinsky, J ., General Circulation Experiments With the Primitive Equa- tions. I. The Basic Experiment, Monthly Weather Review, 91(3):99—164 (1963). [91] Lilly, D. K., On the Computational Stability of Numerical Solutions of Time-Dependent Non-Linear Geophysical Fluid Dynamics Problems, Monthly Weather Review, 93(1):11—26 (1965). [92] Lilly, D. K., The Representation of Small-Scale Turbulence in Numerical Sim- ulation Experiments, in Proceedings of IBM Scientific Computing Symposium Environmental Sciences, pp. 195—210, IBM Form No. 320—1951, 1967. [93] Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., A Dynamic Subgrid- Scale Eddy Viscosity Model, Phys. Fluids A, 3(7):1760—1765 (1991). [94] Germano, M., Turbulence: The Filtering Approach, J. Fluid Mech., 238:325- 336 (1992). [95] Pope, S. B., Lagrangian PDF Methods for Turbulent FLows, Ann. Rev. Fluid Mech., 26:23—63 (1994). [96] Muradoglu, M., Jenny, P., Pope, S. B., , and Caughey, D. A., A Consistent Hybrid-Volume/ Particle Method for the PDF Equations of Turbulent Reactive Flows, J. Comp. Phys, 154(2):342—371 (1999). 152 [97] Xu, J. and Pope, S. 8., Assessment of Numerical Accuracy of PDF/Monte Carlo Methods for Turbulent Reacting Flows, J. Comp. Phys, 152:192—230 (1999). [98] Risken, H., The F okker-Planck Equation, Methods of Solution and Applications, Springer-Verlag, New York, NY, 1989. [99] Gardiner, C. W., Handbook of Stochastic Methods, Springer-Verlag, New York, NY, 1990. [100] Karlin, S. and Taylor, H. M., A Second Course in Stochastic Processes, Aca- demic Press, New York, NY, 1981. [101] Gillespie, D. T., Markov Processes, An Introduction for Physical Scientists, Academic Press, New York, NY, 1992. [102] Kloeden, P. E. and Platen, E., Numerical Solution of Stochastic Differen- tial Equations, Applications of Mathematics, Stochastic Modelling and Applied Probability, Vol. 23, Springer-Verlag, New York, NY, 1995. [103] Milstein, G. N., A theorem on the Order of convergence of mean-square ap— proximations of solutions of systems of Stochastic differential equations, Theor. Prob. Appl, 32:738—741 (1988). [104] Platen, E., An approximation method for a class of Ito processes, Lietuvos Matem. Rink, 21:121—133 (1981). [105] Ito, K., On Stochastic differential equations, Memoirs. Amer. Math. Soc., 4 (1951). [106] Ito, K., Lecture on Stochastic Processes, Tata Inst. of Fundamental Res, Bombay, India, 1961. [107] Gikhman, I. I. and Skorokhod, A. V., Stochastic Differential Equations, Springer-Verlag, New York, NY, 1995. [108] Stratonovich, R. L., Introduction to the Theory of Random Noise, Gordon and Breach, New York, NY, 1963. [109] Carpenter, M. H., A High-Order Compact Numerical Algorithm for Supersonic Flows, in Mortron, K. W., editor, Proc. 12th International Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 371, pp. 254—258, Springer-Verlag, 1990. [110] Solomon, T. H. and Gollub, J. P., Thermal Boundary Layers and Heat Flux in Turbulent Convection: The Role of Recirculating Flows, Phys. Rev. A, 43(12):6683—6693 (1991). 153 [111] Kee, R. J., Rupley, F. M., and Miller, J. A., The Chemkin Thermodynamic Database, Technical Report SAND87—8215B, Sandia National Laboratory, Liv- ermore, CA, 1994. [112] Bui-Pham, M., Studies in Structures of Laminar Hydrocarbon Flames, Ph.D Thesis, UCSD, (1992). [113] Sung, C. J ., Law, C. K., and Chen, J .-Y., An Augmented Reduced Mechanism For Methane Oxidation with Comprehensive Global Parametric Validation, in Proceedings of 27th Symp. (Int. ) on Combustion, pp. 295—303, The Combustion Institute, Pittsburgh, PA, 1998. [114] Mohammad, H., Numerical Simulation of Partiall Premixed Methane Combus- tion in Laminar and Turbulent Flows, Masters Thesis, MS U, (2004). [115] Barlow, R., Sandia/TUD Piloted CH4 / Air Jet Flames, http: / /www. casandia. gov/ T NF /DataA rch/F lameD. html, (2003). [116] Schneider, C., Dreizler, A., and Janicka, J ., Flow field measurements of sta— ble and locally extinguishing hydrocarbon-fuelled jet flames, Combustion and Flame, 135:185—190 (2003). [117] Poinsot, T. and Lele, S., Boundary Conditions for Direct Simulations of Com- pressible Viscous Flows, J. Comp. Phys, 101:104—129 (1992). [118] Xu, J. and Pope, 8., Probability density function calculations of local extinction and no production in piloted-jet turbulent methane/ air flames, in 28th Sympo- sium (Int. ) on Combustion, pp. 133—139, The Combustion Institute, Pittsburgh, PA, 2000. [119] King, M. K. and D., R. H., Overview of the NASA Microgravity Combustion Program, AIAA J., 36(8):1337—1345 (1998). [120] L., U. D. and K., K. M., NASA’s Microgravity Combustion Research Program: Past and Future, Combust. Flame, 116(3):319—320 (1999). [121] Edelman, R. B., Fortune, 0. F., Weilerstein, G., Cochran, T. H., and Haggard, J. E., An Analytical and Experimental Investigation of Gravity Effects upon Laminar Gas Jet-Diffusion Flames, Proc. Combust. Inst., 14:399—412 (1973). [122] Haggard, J. B. and Cochran, T. H., Stable Hydrocarbon Diffusion Flames in a Weightless Environment, Combust. Sci. and Tech., 5:291-298 (1972). [123] Edelman, R. B. and Bahadori, M. Y., Effects of Buoyancy on Gas-Jet Diffusion Flames: Experiment and Theory, Acta Astronautica, 13(11-12):681—688 (1986). [124] Davis, R. W., Moore, E. F., Roquemore, W. M., Chen, L.-D., Vilimpoc, V., and G035, L. P., Preliminary Results of a Numerical-Experimental Study of the Dynamic Structure of a Buoyant Jet Diffusion Flame, Combust. Flame, 83:263—270 (1991). 154 [125] Hegde, U., Zhou, L., and Bahadori, M. Y., The Tfansition to Turbulence of Microgravity Gas Jet Diffusion Flames, Combust. Sci. and Tech., 102:95—113 (1994). [126] Shu, Z., Aggarwal, S. K., Katta, V. R., and Puri, I. K., F lame-Vortex Dynam- ics in an Inverse Partially Premixed Combustor: The Froude Number Effects, Combustion and Flame, 111:276—295 (1997). [127] Bahadori, M. Y., Zhou, L., Stocker, D. P., and Hegde, U., Functional Depen- dence of Flame Flicker on Gravitational Level, AIAA J., 39:1200—1208 (2001). [128] NASA, Microgravity Science and Applications, Program Tasks and Bibliography for FY 1996, NASA Technical Memorandum 4780, 1996. [129] Koochesfahani, M. M. and Frieler, C. E., Instability of Nonuniform Density Free Shear Layers with a Wake Profile, AIAA J., 27(12):1735—1740 (1989). [130] Kailasanath, K., Patnaik, G., and Oran, E. S., Unsteady Multidimensional Simulations of the Structure and Dynamics of Flames, In Proceedings of the Fourth International Microgravity Combustion Workshop [143], pp. 161—166. [131] Bahadori, M. Y., Hegde, U., and Stocker, D. P., Structure of Microgravity Transitional and Pulsed Jet Diffusion Flames, In Proceedings of the Fourth International Microgravity Combustion Workshop [143], pp. 179—184. [132] Hegde, U., Yuan, Z., Stacker, D. P., and Bahadori, M. Y., Characteristics of Non-Premixed ’Ihrbulent Flames in Microgravity, In Proceedings of the Fourth International Microgravity Combustion Workshop [143], pp. 185-190. [133] Ellzey, J. L. and Oran, E. S., Effects of Heat Release and Gravity on an Unsteady Diffusion Flame, in Proceedings of 23rd Symp. (Int. ) on Combustion, pp. 1635—1640, The Combustion Institute, Pittsburgh, PA, 1990. [134] Ellzey, J. L., Laskey, K. J., and Oran, E. S., A Study of Confined Diffusion Flames, Combust. Flame, 84:249—264 (1991). [135] Chen, L. D., Vilimpoc, V., Goss, L. R, Davis, R. W., Moore, E. F., and Roque— more, W. M., The Evolution of a Buoyant Jet Diffusion Flame, in Proceedings of 24th Symp. (Int) on Combustion, pp. 303—310, The Combustion Institute, Pittsburgh, PA, 1992. [136] Katta, V. R., Goss, L. P., and Roquemore, W. M., Effect of Nonunity Lewis Number and Finite-Rate Chemistry on the Dynamics of a Hydrogen-Air Jet Diffusion Flame, Combustion and Flame, 96:60—74 (1994). [137] J iang, X. and Luo, K., Combustion induced buoyancy effects of an axisymmetric reactive plume, Proc. of the Comb. Inst., 28:1989—1995 (2000). 155 [138] Libby, P. A. and Williams, F. A., editors, Turbulent Reacting Flows, Topics in Applied Physics, Vol. 44, Springer-Verlag, Heidelberg, 1980. [139] Pivovarov, M. A., Zhang, H., Ramaker, D. E., Tatem, P. A., and W., W. F., Similarity Solutions in Buoyancy-Controlled Turbulent Diffusion Flame Mod— eling (Turbulent Buoyant Diffusion Flame Modeling), Combustion and Flame, 92:308—319 (1993). - [140] Kollmann, W., Prediction Methods for Turbulent Flows, Hemisphere Publishing Co., New York, NY, 1980. [141] Libby, P. A. and Williams, F. A., editors, Turbulent Reacting Flows, Academic Press, London, UK, 1994. [142] Jaberi, F. A., Livescu, D., and Madnia, C. K., Characteristics of Chemically Reacting Compressible Turbulence, Phys. Fluids, 12(5):1189—1209 (2000). [143] NASA, Proceedings of the Fourth International Microgravity Combustion Work- shop, NASA Conferecne Publication 10194, Cleveland, Ohio, 1997. 156 l[ l lllllll Ill Ml‘ H]: AH K1E.’RSHY RHA‘? 1293 02736 5232