.K... . any. 3:}. 115.3- 3! - B \ .1121), x. — $5, .1: .. mfijflmhfifiau .. .v . {1P , |..' «#391 a; ‘0.)- .. 5:. I, I i‘ITI’ISI ch; . l“. 1| .1533... l... .1155.- I: 0: as...:£..x.. 5:: i. 3...!» ,1 .. a”... .2 , «$7. ‘ I, ‘ I! 1 . .3... 9. 1,33 :31 2.: 3.... 3...... it.“ w! - CIJQWHI to. V¢ L. i 12!. 1%.... l. 7 tin... uri. A 9. 2:9..- .9311I ‘I. l. Alli! ..' :3 3.): 4‘ I t . IHtéog i? This is to certify that the dissertation entitled LARGE-SCALE AND HIGH PERFORMANCE COMPUTATIONS OF COMPLEX TURBULENT REACTING FLOWS presented by ASGHAR AFSHARI has been accepted towards fulfillment . of the requirements for the Ph.D. degree in Mechanical Engineering 7/«fl’fli ’ Major Professor’s Signature 2/02/2005 Date MSU is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State University -.- -o-c-v-|--c-l-n--u--o--t-t- 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 2/05 p:/ClRC/DateDue.indd-p.1 LARGE-SCALE AND HIGH PERFORMANCE COMPUTATIONS OF COMPLEX TURBULENT REACTING FLOWS By Asghar Afshari A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2006 ABSTRACT LARGE—SCALE AND HIGH PERFORMANCE COMPUTATIONS OF COMPLEX TURBULENT REACTING FLOWS By Asghar Afshari A density-based, multi—block, computational model has been developed for large eddy simulation (LES) of reacting and nonreacting, single- and multi-phase, com- pressible turbulent flows in complex geometries and generalized coordinate systems. The spatially-filtered form of the compressible continuity, momentum, energy, and scalar equations are solved together with various subgrid turbulence closures. All spatial derivatives are approximated by a high-order compact differencing scheme and time derivatives are modeled via a low-storage, three-stage, third-order, Runge- Kutta method. Both Smagorinsky and algebraic RNG (renormalization group) sub- grid scale models are employed. The nonreacting, single-phase large eddy and direct numerical simulations (DN S) results for isotropic turbulence, round/ planar jets, and axisymmetric sudden expansion (dump combustor) flows are found to be in good agreement with those obtained via other validated high-order, numerical methods and with the available experimental data. Effects of boundary conditions, subgrid scale model, and various physical and geometrical parameters on the dump com- bustor flows are studied in detail. For this flow, various boundary conditions are examined. The simulated results for dump combustors with and without inlet noz- zle are found to be similar as long as “appropriate” boundary conditions are used. The effects of combustor length and an added end-plate (n—plate) are studied. It is shown that for sufficiently small combustors, the n-plate increases the turbulence, close to inlet. The results also indicate that the high-order compact differencing scheme is an appropriate numerical method for LES and DNS, while the multi-block capability of the scheme enables its application to complex geometries. However, my analysis shows that the first- or second-order upwind schemes are too dissipa- tive for LES and DNS. The results obtained for nonreacting single-phase flows are described. A generalized Lagrangian/Eulerian, theoretical/ numerical methodology is considered for LES of turbulent reacting flows in complex geometrical configurations via the filtered mass density function (FMDF) for subgrid-scale (SGS) combustion closure. The LES/FMDF method has some advantages over conventional methods. A novel numerical algorithm is developed for solving the Lagrangian equations which is much more efficient than similar unstructured grid system algorithms. The La— grangian scheme is coupled with my high-order multi-block flow solver. This allows LES / FMDF to be extended to general coordinate systems. The consistency, conver- gence, and accuracy of the FMDF and the Monte Carlo solution of its equivalent, stochastic differential equations are assessed for different flows. The consistency be- tween Eulerian and Lagrangian fields is established for non-reacting isothermal and non-isothermal flows as well as reacting flows in a dump combustor. The results show good consistency between conventional LES and FMDF method for nonreacting and reacting cases. The results obtained for a reacting flow in an axisymmetric, premixed dump-combustor, as compared with laboratory data, show favorable agreement. The effects of the inlet flow and boundary conditions on the turbulence and combustion within the combustor are investigated by adding a nozzle to the combustor. For the reacting case the results are significantly improved when the non-realistic sudden expansion inlet boundary condition is eliminated. © 2006 Asghar Afshari All Rights Reserved To my parents. ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Farhad A. Jaberi, for enlightening me into the ways of advanced science, turbulence, and combustion, and for his support, commitment, and encouragement that he has given to me. The mentorship he has provided to me will ensure that I achieve my potential and never settle for anything less than the highest quality in my research. I have been fortunate in my studies to also have a co—advisor, Professor Tom I-P. Shih. we started working together from the very beginning of my thesis, and I immediately benefited from his experience and knowledge in numerical analysis, computer simulations, and turbulence. I am very grateful to Professor Shih for his guidance to excellence in research. I gratefully acknowledge Professor Indrek Wichman and Professor C. Y. Wang for the time and effort put into reading this dissertation and for helpful comments and suggestions. I also would like to express my gratitude to Professor John Foss and Professor Charles Petty for very insightful discussions on turbulent and multi-phase flows throughout the course of this study. Thanks to Craig Gunn for editing and proofing the text. I am thankful to Dr. Siamak Hannani at Sharif University of Technology where I gained my first research experience at the graduate level and learned more about the beauty of computational physics. This research would not have been possible without the financial support from the US. Office of Naval Research through Grant number N00014—01-1—0843 under vi the Program Manager Dr. Gabriel Roy. I am greatly in debt to Dr. Roy for his guidance and his support. Computational resources were provided by the Super Computing Institute for Digital Simulation and Advanced Computation at the Uni- versity of Minnesota and by the High Performance Computing Center at Michigan State University. I further owe a great deal of gratitude to my parents for their steadfast support and continuing love. They have shaped me into the person I have become and I would like to dedicate this dissertation to them. vii TABLE OF CONTENTS LIST OF FIGURES x 1 Cold Flow Analysis Using a High Order Multi—Block Method 5 1.1 Introduction ................................ 5 1.2 Governing Equations ........................... 7 1.3 Numerical Procedure ........................... 11 1.3.1 Discretization of the Computational Domain .......... 11 1.3.2 Spatial Discretization ....................... 12 1.3.3 Temporal Discretization ..................... 13 1.3.4 Multi-block Strategy ....................... 14 1.3.5 Boundary Conditions ....................... 14 1.4 Results ................................... 15 1.4.1 Isotropic Turbulence ....................... 16 1.4.2 Planner and Round Jets ..................... 17 1.4.3 Dump Combustor ......................... 17 1.5 Figures ................................... 25 2 Reacting Flow Analysis Using Filtered Mass Density Function Methodology 42 2.1 Introduction ................................ 42 2.2 Governing Equations ........................... 46 2.2.1 Velocity Eulerian Field - Navier stokes equations for a com- pressible fluid flow ........................ 46 2.2.2 Scalar Lagrangian Field - FMDF Methodology for Chemically Reacting Flows .......................... 48 2.3 Numerical Solution ............................ 52 2.3.1 An Eulerian Finite Difference Method for Velocity Field . . . 52 2.3.2 A Lagrangian Monte Carlo Method for Scalar Field ...... 52 2.3.3 Particle Tracking, Interpolation and Averaging ......... 56 2.3.4 Particle Number Density Control ................ 59 2.4 Results ................................... 59 2.4.1 Nonreacting Simulations - Consistency of FMDF ....... 59 2.4.2 Reacting Simulations - Accuracy of LES / FMDF ........ 62 viii 2.5 Figures ................................... 67 3 Conclusions 80 4 Future Work: LES/FMDF of Multi—Phase Turbulent Reacting Flows 83 4.1 Introduction ................................ 83 4.2 Theoretical and Computational Approach ............... 86 4.3 Preliminary Results and Discussions ................... 88 4.4 Proposed Future Research ........................ 90 4.5 Figures ................................... 93 BIBLIOGRAPHY 97 ix 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 1.14 1.15 1.16 1.17 LIST OF FIGURES Energy spectrum vs. wave number, time=10, 15, and 20. ....... 25 Effect of filtering on energy spectrum, time=10. ............ 26 A comparison of energy spectrum obtained via different numerical schemes, time=15. ............................ 26 (a) Contours of instantaneous axial velocity and (b)isosurfaces of vor- ticity magnitude in a turbulent round jet as obtained by LES. . . . . 27 (a) Centerline mean axial velocity vs. downstream location and (b) radial profile of mean axial velocity at x/ D = 7.5 in a turbulent round jet(D is the jet diameter). ........................ 28 (a) Experimental setup and (b) dump combustor test section. . . . . 29 Schematic perspective of (a) the domain discretization and 3-D grid layout and (b) two blocks and interface boundary in a 2-D cross section. 30 Stream lines of mean velocity field at z=0 plane ............. 31 (a) Contours of instantaneous axial velocity at 220 plane and (b) con- tours of vorticity magnitude at different downstream locations. . . . . 31 Distribution of the flow variables for the reference case(o, experi- ment; ------ , LES). .............................. 32 Schematic perspective of (a) the domain discretization, (b) two blocks and interface boundary, for extended domain with nozzle, and (c) Stream lines at z=0 plane for extended domain with nozzle . ..... 33 A comparison of center line temperature in dump combustor with and without nozzle. .............................. 34 (a) Contours of instantaneous axial velocity at z=0 plane and (b)isosurfaces of vorticity magnitude, for the extended domain. . . . . 35 Distribution of the flow variables for the extended case with nozzle(0, experiment;~~, LES). ........................... 36 Distribution of the flow variables for high resolution case(o, experi- ment;~~~, LES). .............................. 37 A comparison of center line velocity in different dump combustor sim- ulations. .................................. 38 Distribution of the flow variables for the case with RNG subgrid scale mode1(o, experiment;-, LES). ..................... 39 1.18 1.19 1.20 1.21 1.22 2.1 2.2 2.3 2.4 2.5 2.6 2.7 A comparison of center line velocity for different dump combustor sim- ulations. In the extended domain case, a nozzle is added to the com- putational domain at inlet. In all other cases there is no inlet nozzle. A schematic picture of the solid walls, including the end-plate, of the combustor .................................. Mean axial velocity contours and streamlines. ............. Distribution of the flow variables for the case with n—plate: dashed-line, n—plate at 3D; dashdot-line, n—plate at 5D; longdashed-line, n-plate at 7D; solid-line, x—lenght=7D without n-plate. .............. Distribution of the flow variables for the case with different inlet tur- bulence: dashed-line, low inlet turbulence; dashdot-line, average inlet turbulence; longdashed—line, high inlet turbulence; solid-line, reference case. .................................... The attributes of LES/FMDF methodology and its LES-FD and F MDF subcomponents ........................... Schematic perspective of (a) experimental setup, (b) dump combustor test section, and (c) the 3-D grid layout. ................ Schematic description of (a) grid ponits, Monte Carlo elements and control volumes around grid points, (b) determining i component of point “q”, (c) determining j component of point “q”, and (d) the two dimensional interpolation scheme. .................... Contours of the instantaneous and time-averaged values of the filtered conserved scalar in an isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values .......................... Contours of the instantaneous and time-averaged values of the filtered temperature in an isothermal non-reacting axisymmetric dump com- bustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time- averaged FD values ............................. Contours of the density in a non-isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF, (a) FD values, and (b) 40 4O 41 41 41 67 68 69 7O 70 MC values (calculated based on weighting of the Monte Carlo elements). 71 Contours of the instantaneous and time-averaged values of the filtered conserved scalar in an non-isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values .......................... 71 2.8 Contours of the instantaneous and time-averaged values of the filtered temperature in an non-isothermal non-reacting axisymmetric dump combustor as obtained by LES /FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values .......................... 2.9 (a) Three—dimensional iso-levels of the instantaneous vorticity magni- tude and Monte Carlo particles, (b) three-dimensional iso-levels of the fuel distribution, (c) two-dimensional planar contours of the fuel, and (d) two-dimensional planar contours of the heat release; in an axisym- metric dump combustor as obtained by LES / FMDF. ......... 2.10 The mean axial velocity contours and stream lines at the z=0 plane. . 2.11 The instantaneous and mean values of the filtered temperature ob- tained from FD and MC data at different downstream locations. . . . 2.12 Contours of the instantaneous filtered temperature in reacting case, (a) FD values (b) MC values. ...................... 2.13 Mean axial filtered velocity and its comparison with experimental data at different downstream locations (0, experiment [reacting flow]; 0, experiment [non-reacting flow]; ~—, LES). ................ 2.14 Mean filtered temperature and its comparison with experimental data at different downstream locations ..................... 2.15 Schematic perspectives of ((a) 3D and (b) 2D) discretized domain for the case with added nozzle. ....................... 2.16 Contours of the instantaneous filtered temperature obtained by three- dimensional LES/FMDF in the reacting case, (a) FD values and (b) MC values, for the case with added nozzle ................ 2.17 Two-dimensional planar contours of the fuel as obtained by three dimensional LES / FMDF, for the case with added nozzle. ....... 2.18 Mean axial filtered velocity and its comparison with experimental data at different downstream locations, for the case with added nozzle (0, experiment [reacting flow]; 0, experiment [non-reacting flow]; .3 LES). 2.19 Mean filtered temperature and its comparison with experimental data at different downstream locations, for the case with added nozzle. 72 73 74 75 76 76 77 77 78 78 78 79 2.20 Distribution of the various turbulent quantities(o, experiment; ..... , LES). 79 4.1 Experimental setup. ........................... 4.2 Liquid Fuel Droplets (big solid points), Monte Carlo particles (small solid points), and Eulerian Grid layout. ................ 4.3 Contours of the instantaneous vorticity magnitude and droplet and Monet Carlo elements distribution (a) 3D and (b) 2D .......... 93 95 4.4 Contours of the instantaneous gas temperature, pressure, and gas species mass fractions as obtained by two-phase LES / FMDF in a spray- controlled premixed reacting dump combustor .............. 96 xiii Nomenclature .> “6| N 2: Ct 81 Filtered pressure Filtered internal energy Filtered total energy Filtered Cartesian velocity components Filtered temperature Constant of the stochastic mixing closure Smagorinsky model coefficient RNG model coefficients Drift coefficient in the SDE Equivalence ratio Inviscid fluxes: chapter 1 and equation (2.1) Joint scalars filtered mass density function Viscous fluxes Filter function Enthalpy Enthalpy of formation of species a Subgrid heat flux Subscripts i, j, k Conventions J acobian ith component of the flux of scalar or 2th component of the subgrid scalar flux of species a Number of species Turbulent Prandtl number Temperature Time Solution vector z'th component of the stochastic velocity vector Coordinate direction in physical domain Lagrangian position of the particles ith component of the position vector Diffusion coefficient in the SDE Position vector Variable numbers Quantities which depend only on the scalar composition, zle. §<¢,x, t) a S<¢> Favre filtered value; equivalent to ( ) L Filtered value Conditional Favre filtered value Symbols af 5 At A6, A17. AC €t,€x,”',Cz #e #8 gs 7.393 t] F avre filtered value; equivalent to ~ Dot product between two vectors Filtering parameter Filtered density Scalar field Composition space time step filter size for LES model Grid spacing in 5, n and C directions Metric coefficients Thermal conductivity Molecular dynamic viscosity Effective viscosity Subgrid scale viscosity for RNG model Turbulent kinematic viscosity SGS mixing frequency The compositional values of stochastic scalar a The compositional values of scalar 0 Density Number of scalars, o = N5 + 1 Subgrid stress tensor The filtered compositional values of scalar a é, 77,C Coordinate directions in transformed domain Superscripts (72) Index of the Monte Carlo particles + Properties of the stochastic Monte Carlo particles S GS Subgrid scale CHAPTER 1 Cold Flow Analysis Using a High Order Multi—Block Method 1. 1 Introduction LARGE EDDY simulation (LES) is a promising computational method for analysis and prediction of engineering problems involving turbulent reacting flows. LES is superior to Reynolds average simulation (RAS) methods in predicting mixing and combustion in unsteady turbulent flows, which are characterized by sustained oscil- lations, and strong coupling between turbulence, the acoustic field, and combustion. Within the past few decades, significant progress has been made in the development and application of LES to complex turbulent flows due to advances in computational power; more accurate numerical algorithms; and high—fidelity, subgrid-scale (SGS) closures. However, an analysis of the effects of spatial discretization errors on LES establishes the need to develop high-order algorithms applicable to general geome- tries such as compact finite difference methods. Compact differencing [1—4] is a well- developed method for structured grids and can be easily extended to higher orders. It has been used to study flows in complex geometries in a curvilinear structured mesh [2, 5, 6]. Some attempts have also been made to develop a high-order, compact method and finite volume method for unstructured grids. Most numerical schemes used in unstructured flow solvers are not appropriate for DNS and LES, although they are appropriate for RAS. For structured grids, the multi—block capability enables the application to complex geometries [2, 7, 8]. A multi-block flow solver can be efficiently coupled with high- order, compact-difference schemes. Since high—order, compact differencing usually has spectral-like resolution, it is an attractive choice for reducing dispersion, anisotrOpy, and dissipation errors associated with low—order spatial discretization. These schemes can be used with explicit or implicit time integration methods. Here, it is noted that in many applications a low-pass spatial filtering is needed to enforce numerical stability on nonuniform grids[2, 4]. An important objective of this work is to develop a high—order model for the large eddy simulation of turbulent reacting flows in complex geometries by using multi- block, compact—differencing schemes and to evaluate the model for different flows. The main objective is to study the turbulent reacting flow in a dump combustor via LES/FMDF. This chapter deals with nonreacting, isothermal flows and describes my efforts on the development and validation of the flow solver for its subsequent application to the dump combustor. For reacting flows, I have developed a new algorithm that is based on the filtered mass density function (FMDF)[9, 10]. This algorithm and the LES / FMDF results obtained for the reacting case will be presented in chapter 2. The flow through an axisymmetric sudden expansion or dump combustor is a well-known and well-studied example of a turbulent, separated flow in practical en- gineering situations. The sudden expansion is one of the most common geometrical configurations for combustors, in which the flames can be stabilized by a stationary recirculation zone. The interest in this flow configuration has led to a significant number of experimental [11—19] and numerical investigations[20—27]. These studies include those which are focused on flame-turbulence interations in a dump combus- tor setup, a subject of chapter 2. These investigations have significantly contributed to our understanding of turbulence, mixing, and combustion in dump combustors; our understanding, however, is not yet complete. Among the various experimental studies, Gould et al’s. [11, 12] is desirable [13] from computational view point, since various turbulence variables in both nonreacting and reacting conditions are mea- sured. Also, unlike most experimental setups , in their experiments a nozzle was used (instead of an incoming pipe) at combustor inlet which makes the specification of the boundary conditions easier. With a nozzle, the flow at the combustor inlet will be relatively flat or uniform. This experimental work was chosen here to validate the computational model. 1.2 Governing Equations The compressible form of the transport equations of mass, momentum, energy, and scalars are cast in strong conservation law form, after introducing a general curvilinear coordinate transformation (x-y-z coordinate system) to (E-n-C coordinate system) [28— 31]. Large eddy simulation involves the spatial filtering operation [32] N +00 f(x,t) = >g = (.00 fHzfi+1 — ¢i—1 4 + a 2 (1.27) I I / 09527—1 + 9% + 0395241: b For an order of accuracy of four, a, a, and b are 211, g, and 0 respectively, as obtained by Taylor series term-matching procedure. At boundary points, higher-order one- sided formulas are utilized which retain the tridiagonal form of the interior scheme. As mentioned above, the numerical noise generated by the growth of numerical error at very high wave number/ frequency modes is removed here by using a filter- ing procedure. This filtering Operation is different than the standard LES filtering operation and it is an important component of the solution algorithm. In LES and DNS, the application of nondissipative spatial schemes leads to a pileup of energy at the smallest scales of the flow and numerical instability. Compact differencing, like other central difference schemes, is a nondissipative scheme. In this work a low-pass, high-order, spatial implicit filtering operator [4, 38] is used for both interior and near- boundary values. For any component of the solution vector (i.e. 45), filtered values (i.e of) are obtained by solving the following tridiagonal system N 074—1 + 45f + af¢i+1 = Z a7n(¢i+n + (”i—n) (1-28) n=0 Here, a f is the filtering parameter, which satisfies the inequality —0.5 < a f S 0.5. 12 Higher values of 01f correspond to a less dissipative filter. The coefficients are derived in terms of a f by performing Taylor- and Fourier-series analyses. Spatial formulas are used at near-boundary points where equation (1.28) cannot be applied. These numerical filtering operations suppress numerical instabilities arising from “very small scale interactions”, mesh nonuniformities, and boundary condition effects. Filtering is performed at the end of every step of R-K and is applied in each direction. For multi-block cases the data are exchanged between blocks after filtering. I have tested both the fourth- and the eighth-order filtering for interior and boundary points and I have found the effects of filtering on the flow statistics to be negligible. 1.3.3 Temporal Discretization All of the time derivatives in the governing equations are approximated by the fol- lowing low storage three-stage Runge—Kutta scheme, which belongs to the strong, stability—preserving (SSP), R—K family [39, 40] (JU)§]1.I)€= (JU)? k +At (aggrfl (1.29) 2 3 n BJU (1) (Jung), = Z(JU)[].,Z+ZI (JU)(_ 31,) +_ jl—At( 6t )ijk (1.30) n 1 n am (2) (JU)i]i1—_—_§(JU)Z(.]2+§ (JU)§fl)C+ §At(— 8t —)ijk (1.31) where am 0F or: 8H (Wt—>ijszijk=—(a—€ +—+79_C_+JS) (1.32) In equations (1.29)- (1.32), U n is known and U n+1 is the unknown, while U (1) and U (2) are the intermediate solutions. 13 1.3.4 Multi—block Strategy In the multi-block strategy, a series of structured grids of varying sizes is constructed in such a way that they fill the entire computational domain and conform to the surface of the geometry of interest. This segmentation of the complete domain into smaller blocks avoids topological problems usually encountered in constructing grids for complex configurations and multiply connected regions. The main idea is to construct a halo of cells that surrounds each block and contains information from cells in the neighboring blocks [7]. This halo of cells, when updated at the appropriate times during the numerical solution procedure, allows the flow calculations inside each block to proceed independently of the others. This approach requires the identification of halo cells adjacent to block boundaries and the construction of a list of halo cells and their internal counterparts in other portions of the global mesh. The requirement for this double halo procedure is due to the need to calculate all of the necessary fluxes for the internal cells of each block without reference to additional cells located outside the block. If spatial derivatives are computed using a second-order, explicit method, the communication between blocks occurs at every stage of the Runge—Kutta scheme in a way that treats the multi-block grid as if there were only one block with no internal boundaries and a time-lag between blocks. For an implicit method, the accuracy will increase with more overlaps. It is to be noted that, for an implicit method, the multi- block is not treated if there is only one block. In the flow solver developed here, data are exchanged between adjacent domains at the end of each stage of the Runge-Kutta scheme, as well as after each application of the high, wave-number, filtering operation (conducted for removal of numerical noise). 1.3.5 Boundary Conditions In this work several types of boundary conditions are considered for various flow simulations. All boundary conditions fit the physical and mathematical/numerical 14 criteria. For isotropic turbulence, periodic boundary conditions are used in every direction. In all jet simulations, characteristic boundary conditions are used for the inflow boundary. Also, zero-derivative boundary conditions are employed for free stream boundaries, and isothermal non-slip condition for the walls. The inlet bound- ary conditions have been extensively examined for the dump combustor simulations. For the case that the computational domain is truncated right at the sudden expan- sion point, the results suggest that by floating the temperature (e.g. using 939% = 0), fixing density, and computing the pressure through the equation of state, the numeri- cal stability is well preserved and improved results are obtained. For the case with an added nozzle all tested inlet boundary conditions lead to the same results. Also for the inlet boundary condition, the experimental data were used to generate turbulent velocity perturbations based on a multi-frequency, random—phase harmonic, forcing function. Characteristic boundary condition [41] coupled with the convective bound- ary condition is used for the outflow boundary in all inhomogeneous (jet and dump combustor) simulations. 1.4 Results The flow solver developed is validated by performing DNS of isotropic turbulence and 2-D jet flows. LES of 3-D jet flows are also conducted to examine the performance of my LES solver in more realistic turbulent flow conditions. The results obtained from these simulations were found to be in good agreement with solutions obtained by using a validated, high-order spectral method and available experimental data. For further validation and application of the computational model, algorithms and multi—block capability, large-eddy simulations of the flow in a dump combustor (axisymmetric sudden expansion) are also performed. In all simulations, the Mach number is less than 0.3, hence the compressibility effects are small. 15 1.4.1 Isotropic Turbulence Direct numerical simulation of isotropic turbulence is the first test case for evaluating the accuracy of the high-order schemes and the effect of low-pass numerical filtering. In these simulations, the initial Reynolds number based on the turbulence intensity and the Taylor microscale is 25, and a 643 mesh is used. Figures ( 1.1- 1.3) show the energy spectrum versus wave number for different nondimensionalized time and numerical schemes. The finite difference (FD) results are compared with results obtained by spectral and finite volume (F V) methods. As observed in figure 1.1, after about six eddy turn-over times (the eddy turn-over time is ~ 3.5), the FD results are still in agreement with the spectral method. This confines the accuracy of my high order compact scheme. Figure 1.2 shows a comparison of the finite difference results obtained by fourth- and eighth-order filtering (a f is the filtering parameter [2, 38]). As mentioned before these high-order filtering are used for removing the numerical noise. It is clear that the application of eighth-order filtering leads to slightly better results. The comparison between the results obtained with the current FD method with those generated via our second-order FV method is illustrated in figure 1.3. In the finite volume method, instead of low pass filtering, a second-order, upwind, numerical dissipation was added to the equations in order to preserve the numerical stability. As shown in figure 1.3, this numerical scheme is too dissipative for direct numerical simulations and cannot predict the energy spectrum correctly. There is a significant discrepancy between the second-order FV method and the spectral method. Similar results were obtained by the FV method with twice the resolution, suggesting that the discrepancy is not solely due to the spatial order of accuracy of the FV method and the upwinding error is also important. 16 1.4.2 Planner and Round Jets As a simple, inhomogeneous-flow model, a two-dimensional jet with Reynolds number of 3000 (based on jet thickness and mean bulk jet inlet velocity) is simulated with my high-order, compact, finite-difference DNS code. The results (not shown here) indicate that the physical features of the flow have been captured by the DNS code. The results [42] also indicate that our DNS data are in good agreement and consistent (both quantitatively and qualitatively) with the available literature data. Turbulent, round, jet flows have been studied by numerous investigators and have been also considered by us for the validation of my flow solver. Simulations of various jets with different inlet flow conditions indicate that the near filed statistics are very sensitive to the inlet flow conditions. Here, I present the three-dimensional LES results for two of the simulated round jets with the same Reynolds number of 5700 (based on jet diameter and mean bulk jet inlet velocity) and different inlet flow conditions. Figure 1.4 shows the contours of instantaneous axial velocity and isosurfaces of vorticity magnitude for the first simulated jet. The comparison with experiment [43] is shown in figures 1.5 for the second simulated jet. The numerical results are in good agreement with the experimental data and correctly capture the physical features of a spatially developing turbulent round jet. 1.4.3 Dump Combustor The focus of this study is on the large-eddy simulations of turbulent flows in an axisymmetric dump-combustor. Gould et a1. [11, 12] have conducted laboratory experiments on a similar configuration for nonreacting and reacting flows. In this chapter, the LES results for the nonreacting flow with Reynolds number of 115,000 (based on the inlet diameter and mean bulk inlet velocity) are compared with the ex- perimental data [11]. Figure 1.6 illustrates the experimental setup and the geometry of the simulated dump-combustor. 17 The dump-combustor simulations are conducted via a fixed Eulerian, multi-block, grid system (figure 1.7). For an axisymmetric configuration centerline, singularity is avoided by using an appropriate, multi-block system. In this simulation I only use two “nodes” that overlap between the two blocks. Additional overlap nodes reduce the interface numerical errors. The numerical results obtained by LES have been validated by comparing them with the experimental measurements of Gould et al. A two-block, O-H grid is used for the computational domain. The center line singularity associated with O—H grids is eliminated by using an H-H grid in the center of the duct. All grids are gener- ated by using transfinite interpolation with controls for clustering, smoothness, and orthogonality [44, 45]. Grid clustering is used in the regions near the walls, in the shear layer, and close to the inlet. Schematic perspectives of the grid are shown in figure 1.7. The resolution for the reference case is 159x47x58 for the interior block, 159x17 x17 for the outer block. The flow streamline as obtained from the LES, mean-velocity field is shown in figure 1.8. The reattachment point predicted by LES is approximately 7.5 step heights, which is slightly lower than the experimental value (8.0 step heights). Figure 1.9 shows the contours of the instantaneous magnitudes of axial velocity and vorticity. Also, in figure 1.10, mean axial velocity, mean radial velocity, the normalized axial turbulent component of the normal stress, and the nor- malized shear stress have all been compared with the experimental data in different down-stream locations. The mean axial velocity is in good agreement with the ex- perimental data. However the velocity is overpredicted in the developing, shear-layer region. The predicted mean, radial velocity agrees quantitatively with experimental data in most of the domain, although it is underpredicted in recirculation zone. As shown in figure 1.10, the mean radial velocity is not predicted as well as the mean axial velocity. This discrepancy could be partly attributed to the experimental data. Mean radial velocity is an order of magnitude smaller than the mean axial velocity 18 over much of the computational domain so there can be some uncertainties in the experimental data, as noted by Gould and Stevenson [11]. Consistent with the exper- iment, the peak values of axial, normalized, turbulent, intensities occur in the shear layer. There is good agreement between numerical and experimental results at down- stream locations after “x/ H = 4” but the experimental values are underpredicted in the near field region. This is due to the inlet boundary condition as it is not a truly turbulent boundary condition. The turbulence becomes more realistic and the turbu— lent intensities are predicted well at various downstream locations. When the nozzle is added to the computational domain, the effect of the ‘nonphysical” inlet-boundary condition is reduced and a better prediction of turbulent intensities are achieved (this issue will be discussed further in the next subsection). The comparison between axial and radial normal stresses indicates that the flow is not isotropic. In summary, the numerical results are in good agreement with the measurements and they capture the trends in the experimental data. INLET AND OUTLET BOUNDARY CONDITIONS In this study the inlet- and outlet-boundary conditions have been extensively ex— amined and tested. This is an important issue in LES and particularly in dump- combustor simulations. The results of my numerical experiments show that the char- acteristic boundary condition coupled with the convective boundary condition for outflow is stable. However, the convective exit velocity must be large enough to force the structures to exit the domain before initiating any numerical instability. The con- vective velocity is set to be the maximum velocity at the boundary in all simulations. This might make the results unreliable at the outlet, but it is necessary to preserve the stability. At the outlet boundary, the pressure is computed via the characteristic boundary condition [1]. For the inlet boundary condition, experimental data were used to generate tur- 19 bulent perturbations for all velocity components, and the mean velocities are set identical to the experimental values. The numerical methodology needs a boundary condition for one of the thermodynamic variables (temperature, density or pressure) while the other variables will be fixed or computed through the equation of state. Various boundary conditions that fit the physical and mathematical/ numerical crite- ria are used. Basic features of various tested boundary conditions at inlet(INBC) can be summarized as follows: 0 INBC (1): T=cte, characteristic boundary condition is used to compute the density [1] and pressure is computed through the equation of state. 0 INBC (2): T=cte, the density is allowed to float via pseudo-characteristic . . a _ a? _ . boundary condition ,55 — O or a—zlg — 0, and pressure 1S computed through the equation of state. 0 INBC (3): T=cte, the pressure is allowed to float (via pseudo—characteristic boundary condition 32 = 0 or 92% = 0 and the density is computed through ‘I x 81: 3 the equation of state. 0 INBC (4): The temperature is allowed to float at the inlet (ET: 2 0), the pressure was kept constant and the density is computed through the equation of state. 0 INBC (5): The temperature is allowed to float at the inlet ("3% = 0), the density is kept constant, and the pressure is computed through the equation of state. The simulations indicate that IN BC 1, 2, and 3 are causing some numerical insta- bility and error. In fact for these boundary conditions and all other tested boundary conditions that are based on fixed temperature, the instability problem persists re- gardless of the way that the density and the pressure are computed. When the computational domain is cut off right at the the sudden expansion point and the 20 temperature is fixed, the imposed boundary condition cannot adjust to the changes in the interior points and INBC (1)-INBC (3) become numerically unstable. However, INBC (4) and INBC (5) (in which the temperature is allowed to float at the inlet ) do not exhibit any instability. In this study INBC (5) is used, which yields good results consistent with previous observations [46]. To clarify the source of the instability, a nozzle was added to the computational domain before the sudden expansion point, and the data were compared with the case without nozzle. A three-block, O-H grid is used for the computational domain. The resolutions for this case are 207x 17x58 , 207x 17 x 17 and 159x32x58 for three different blocks (the resolution is the same in radial and azimuthal direction as com- pared to the reference case). Schematic perspectives of the gird layout and streamlines are shown in figure 1.11. Reattachment is predicted as for the reference case. In figure 1.12, the mean center-line temperatures obtained by simulations with or with- out the nozzle are compared. There is good agreement between the data obtained with the nozzle and that without the nozzle and INBC (5). This suggests that the boundary condition used for temperature in the reference case is accurate and adjusts correctly to the flow inside the dump combustor. The results also indicate that the temperature continuously decreases through the nozzle, where the density change is negligible. Fixing the temperature at the inlet for the case without the nozzle (refer- ence case) causes significant errors, as expected. The same results were obtained for the case with the nozzle using the characteristic boundary at the inlet and using the pseudo-characteristic boundary condition for the pressure or the density [1]. Figure 1.13 shows the contours of instantaneous axial velocity and vorticity mag- nitude isosurfaces in the dump combustor with the inlet nozzle. Also, in figure 1.14 the mean axial velocity, the mean radial velocity, normalized axial turbulent normal stress, and normalized shear stress obtained by LES have all been compared with the experimental data in different downstream locations. The predicted velocity profile 21 at “:r/ H = 0.2” is in good agreement with the experimental data. As expected, the nozzle generates a flat profile with a very thin shear layer at the combustor inlet. The mean axial velocity predictions are almost the same as the previous (reference) case, where better agreement with the experiment is observed for the mean radial velocity and the various turbulence intensities. As mentioned earlier, the negative effects of a “non-physical inlet boundary condition” on the downstream flow in the dump combustor is minimal when a nozzle is added at the inlet. GRID INDEPENDENCY A grid-sensitivity test is performed to asses the dependency of the results to grid resolution. Various grid layouts with different resolutions are tested. The results for high resolution (159x61 x66 and 159x19 x 19) case are compared with experimental data in figure 1.15. A low-resolution simulation is also performed and is referenced to the low-resolution case in the subsequent discussion. A comparison between the results shown in 1.15 with those in 1.10 indicates that although the “grid effects” are not very significant; there are some small discrepancies between the results obtained with various resolutions / grid layouts. However, the discrepancies are not just due to resolution, rather they are mainly due to the SGS stress model. The model coeflicients in the Smagorinsky closure are proportional to the cell size and can easily lead to the overestimation of the turbulent viscosity in the regions with larger cell sizes. To show that a part of the discrepancy between low resolution and high resolution results is due to SGS model, a new simulation with a different Smagorinsky model coefficient was conducted. In this simulation, Cd in equation (1.21) was adjusted so that CdA2 becomes nearly the same for a grid point in shear layer for both the reference and low resolution cases. As observed in figure 1.16, the center line velocity as obtained with the adjusted coefficient is very close to that for the reference case. Without the coefficient adjustment, the center line velocity in low resolution case does not decay 22 as fast as the reference case, simply because the local filter size in the Smagorinsky model is larger. In other words, the turbulent viscosity is over estimated in the shear layer by the Smagorinsky model, resulting in slower growth in the turbulence and a slower decay rate of the center-line velocity. To further investigate the effects of the SGS model, an algebraic renormalization group (RNG) SGS closure [37] is implemented. The results are compared with those obtained from the Smagorinsky closure. In the model, the coefficients, Crng = 0.157 and C = 100 are obtained from the RN G theory. In the highly turbulent regions, the filtering operation results in a very high subgrid viscosity compared with the molecular viscosity #593 >> It, and the effective viscosity becomes the same as the SGS where, lie g #sgs- In these regions, the models based on the RNG theory yield results not very different from the Smagorinsky model if a different model constant is used. In weakly turbulent regions, the argument of the Heaviside function is negative and the effective viscosity is equal to the molecular viscosity. The RN G SGS model in this limit correctly yields zero SGS viscosity in low Reynolds number flows without any adjustments. A comparison with experimental data is shown in figure 1.17. The results indicate improvement with respect to the reference case. Finally, in figure 2.20 the mean, axial, center-line velocity for all the cases is compared. The discrepancies occur mainly because of the Smagorinsky model, which is more effective for the low resolution case (155 x 35 x 42 and 155 x 13 x 13). OUTFLOW CONDITIONS AND INLET TURBULENCE EFFECTS The effect of the outflow condition and inlet-flow turbulence is very important in practical situations. These influences are studied in this subsection. An n-plate is added to the outlet of the dump combustor, which confines the outgoing flow to an area identical to the cross—sectional flow area at the sudden expansion point at the combustor inlet. Figure 1.19 Shows a schematic picture of the dump combustor with 23 the added n-plate. Numerical data are obtained for three different combustor sizes (x—lengths of 3D, 5D and 7D) and the results are compared with the reference case, which has no n-plate and an x-length of 7D. Figure 1.20 Shows mean axial-velocity contours and streamlines for the case with n-plate and x—length=3D. This picture presents the effect of the n—plate on the mean-flow structure and the recirculation zones. Recirculation zones are formed at the corner of the n-plate and the combus- tor wall. The added n-plate changes the flow structure and increases the turbulent intensities. To examine the effect of the n—plate and combustor length on the tur- bulence statistics, the results are compared in 1.21. The comparison of normalized axial turbulent normal stress and normalized shear stress indicates that by adding the n-plate and by reducing the combustor length, the turbulent intensities close to the combustor inlet increases. This forces the flow to develop faster. The effect of incoming flow turbulence on the flow field down stream of the sudden expansion point in the dump combustor is examined by increasing the turbulent intensities at the nozzle inlet. Different quantities for different cases are compared in figure 1.22. By increasing the inflow turbulence, the shear layer grows faster downstream of the expansion point. It is interesting to note that for sufficiently small combustors, the n-plate has the same effect as increasing the inlet turbulence. 24 1.5 Figures - ‘ ~-= .3 '--::=..:. ~ \""=-~-..- 1o li= 8\'.\ e .\..\l’"e \0 “‘e. a: <2” 10-5 — . '3 x I." N ‘9 x- . a 7 ‘ ' a 3 0 [1110.9 _ o SpectraL time=10 ‘ ho ”a FD, 8Lh order filten'ng, time=10 9‘ 5.x o SpectraL time=15 o ‘ H — — - FD, 8zh order f11tering,time=15 ’5 $ 10' - 0 Spectral. tiIIIe=20 '3’.“ -— ----- — FD, 8:11 order filtering, time=20 0"; O C 10-13 _ .. . . I I 1 I I I I I IIIIIIlIIIIlIIIIlIII 5 10 15 20 25 30 k Figure 1.1. Energy spectrum vs. wave number, time=10, 15, and 20. 25 103 105 :- 0 Spectral ‘w A E — ----- — 4'n order filterng (0.50.495) ‘9 S110, ;- 8m order filtering (0.50.49) Q. E I. 108 b ‘g b g r . . ' 10 .h 6' JD 10 .. O. o l 1 1 l 1 1 L4 1 llllllllllllllllll‘ 10 15 20 253‘0 Figure 1.2. Effect of filtering on energy spectrum, time=10. 10“ - \o‘ 10(5 - \\ C. x ‘5 ’s, x a ‘0. A ‘9 S109 __ 0 Spectral ‘9 _ ....... FD (4.“h 0A), 4m order filtering ‘ _1 ”EL-” 1:1)(4‘h OA), 8"" order filtering b...“ ..— 1010 L _ _ — Fv ( 2nd 0A) Q‘o, 0 ."'“\ O O 10"2 - ° . O l l 1 l 1 I l I I lLJLJlllllllllllll£ 5 10 15 20 25 30 Figure 1.3. A comparison of energy Spectrum obtained via different numerical schemes, time=15. 26 (b) Figure 1.4. (a) Contours of instantaneous axial velocity and (b)isosurfaces of vorticity magnitude in a turbulent round jet as obtained by LES. 27 0.8.....,.. .EXP 0.6 - ' - ’1 ‘3 Mean Axial Velocity p .h :5 0.2 r a a - I n l 0 1 n l n I . -l.25 -0.75 —025 0.25 0.75 1.25 Radial Position (a) 1 L ‘ a cam “ . 49 LBS , W 0.9 - J is - o g 0.8 - It - E * ‘ . J a 07 2 0.6 - s - 05 L . ° C 0.4 1 l 1 L 1 0 5 10 15 Axial Position 0)) Figure 1.5. (a) Centerline mean axial velocity vs. downstream location and (b) radial profile of mean axial velocity at x/ D = 7 .5 in a turbulent round jet(D is the jet diameter). 28 0.914m it 0.676 J 0.762 m Honeycomb Flow R Strighteners /_ Quarts Test Section DOP Seeder 5:] Contraction [Extension Duct (a) L=596.9 mm _I II ”III/II/I/III/IlllIIIIII/IIIIIIII/III]!III/III/III/III/IIIWM (b) Figure 1.6. (a) Experimental setup and (b) dump combustor test section. 29 '9 "ash \\\\\\* .» . ‘6‘- \'\\\\\\‘3‘ (“A \l’. fiswxfifmfi “til“ “‘1‘ r . L‘R‘wfis‘fitv‘t‘hfi’ . as , \\ II‘“ \ \ \\ t gwgs‘s‘; \nfihfl‘f‘a ‘39 “V. \\ “ “‘25 " \\$§\\:“‘\\ ., '1 . “949284139? 2* '8‘ \ \ f" ‘9 \\\\ IA“ Iii ‘9‘“ fits“ (It‘lfig I '- \ s \‘ \ . y \\ \‘ \\ - . I::,-'Z~.~..:.:.~\.§N§§\“‘f\“‘\\\\\\ u, o 0 “\“w‘ \\\\“ “““ ‘el‘t‘xyi [t\\\‘|‘l‘|‘“ I :1 “I. Ili‘uuuumm ff. (“lullllllllllll .‘ a“ ‘| 3.. 5’ I as " ~ ‘ ’ ‘ .~ 6“ - 699 ==- ~ ’ av -- ’6 ’I" z, / I 0.: Q J i . 1,, ‘a 710’ I" l'. I ll ,1, ’I ’I/ ’< ’I/ I \‘\\‘ -h Q“ [’6’i6 e3? -52: *s‘s‘ ’s’v” ’ is 1; @z’ / ”I ’0’, (b) Figure 1.7. Schematic perspective of (a) the domain discretization and 3-D grid layout and (b) two blocks and interface boundary in a 2—D cross section. 30 Figure 1.8. Stream lines of mean velocity field at z=0 plane. n ‘rrw'; (b) Figure 1.9. (a) Contours of instantaneous axial velocity at z=0 plane and (b) contours of vorticity magnitude at different downstream locations. 31 u . M41 I 111 MU... .0000...- u - ”II-II... p L w w to... m a . one... 8 8 ... 8 IIIIIIIIII. .IOII Ico- .... ... . (a) 22/ Mean axial velocity Mean radial velocity (d) 32 Normalized shear stress 0)) Normalized axial turbulent normal stress (C) I. an CC. 3 K... 3 .. ...-... 3 -.fi..:. till!!! 4111 till. fl» 2 A 2 no 2 ...- ~’ ... .. I”... I'f' ,1 Larry - 11-: 1%.)! Distribution of the flow variables for the reference case(o, experiment;~~~-~, 0 0.62 ._ i .85 .w 2. nu AU 0 . | u . 0 .36. . 0 «I 4 ll] 0 5 l 0 5 l “1.0 5 ... .m: o. .mu 0. m “me o. Figure 1.10. LES). Figure 1.11. Schematic perspective of (a) the domain discretization, (b) two blocks and interface boundary, for extended domain with nozzle, and (c) Stream lines at z=0 plane for extended domain with nozzle . 33 300 299 N (O (D 297 296 295 Mean temperature 294 293 : Z \ h— \ . _'_- ‘\ Reduced domain (reference case) - - " ''''' Extended domain (with nozzle) . , . ; P .— l l l l I l l l l I l l l l l l 5 1 0 1 5 Figure 1.12. A comparison of center line temperature in dump combustor with and without nozzle. 34 (b) Figure 1.13. (a) Contours of instantaneous axial velocity at z=0 plane and (b)isosurfaces of vorticity magnitude, for the extended domain. 35 12 10 Mean axial velocity (6!) 12 10 0.2 .‘1 Mean radial velocity I no.0 . no.0- (b) 12 10 xii-l: 0. 2 Normalized axial turbulent normal stress :oooncccooIIIOIo- I. CC... 0.. I .00... on... . .OFOII coon-oo- ICIICI C... II... I. I... III 11 I o c coo-co- £7“...- .ié L)... - u (C) 10 o 0.62 Normalized shear stress ((0 Figure 1.14. Distribution of the flow variables for the extended case with nozzle(o iment; , LES). exper- 9 36 12 10 ... .0... Mean axial velocity 12 10 7. 0.2 1%! ...-II. 0 ll... . .. I- ‘l ‘ all. OIIOCOI III [11 .3. D. I." .- ill!) 0.5 ’ Mean radial velocity 12 10 .- on .0. m In... I .m .... on... m .m X ...... m m 0 on N I no. .500? 5r - a. O 0 P 0 (J. l 0 (C) 12 10 ‘1‘i‘1“‘ o 0.62 Normalized shear stress (60 Figure 1.15. Distribution of the flow variables for high resolution case(o, experiment; , LES). 37 Reduced domain (referaice case) ’_ . Exp. data (Gould et. a1) : - -- - Low resolution I . - ----- Adjusted Smagorinsky model coefficient (low resolution) 1.1 9.0.0.0 O)\lCD(O-| .0 on oIIIIllIII'IIII'IIII'IFIIIIITIIIllIIIIII Mean axial velocity .0 rs .0 w .0 N Figure 1.16. A comparison of center line velocity in different dump combustor Simulations. 38 12 IO P0000000. Mean axial velocity (a) Mean radial velocity I no.0 . 86. (b) 12 10 I. Eco-Ice. :00... llfilUllooe-nle ..... L»?! . Normalized axial turbulent normal stress 0 0.65 (C) 12 10 0.2 xi 1 0“ 0 0.62 Normalized shear stress (01) Figure 1.17. Distribution of the flow variables for the case with RNG model(e, experiment; ——, LES). subgrid scale 39 Reference case 1.1 '- - ----- Extended domain with nozzle : - - - High resolution : . “' """"" RNGmodel 1 .. . Exp. data (Gould et. al) I . \ \ — — — Low resolution I- , ' 0.9 :- 1' 3* - - .5 : .1 o 0.8 _- l '_‘ . Q) I l > _ . _‘ 0.7 :- i .2 : ’ 25 0.6 i- f I I a : i g) 0.5 _- I. 2 E I 0.4 _- .t : 'l 0.3 :1 I/ 0.2 _ L 1 l I l I I I I I I J l l l M l I 0 5 10 15 Figure 1.18. A comparison of center line velocity for different dump combustor simula- tions. In the extended domain case, a nozzle is added to the computational domain at inlet. In all other cases there is no inlet nozzle. \ I 21) “D n-plate D l/- l X-length Figure 1.19. A schematic picture of the solid walls, including the end-plate, of the com- bustor. 40 Figure 1.20. ~ H 1 2 4 5 2 4 s V i ' I - l 1 z i ‘. x R, {._ i , \t \ \t \t; g t» K . ‘ \\\\\ \:\ .> a“: V\ L— \\.N ‘ ~\’ \\\, \ . ‘1: \\\ \ “; 05} ,A) ‘ .. ‘\ ‘\ . \nt - \Ct ,} >) ’\ ‘ _ . \ ’ ,1," xi, ‘1’“. r / I; I , I" .I/ /’_' ‘ .’ , // 9 ’ ll: I I, ' ,1” ‘5‘ {/4 I :1 1g .1 I t I] ’ i l _- , r .' f. l 0.02 0.04 0.06 0.05 0.005 0.02 0.03 0.03 Normalized axml turbulent normal stress Normalized shear stress (a) (b) Figure 1.21. Distribution of the flow variables for the case with n-plate: dashed-line, n-plate at 3D; dashdot-line, n—plate at 5D; longdashed—line, n-plate at 7D; solid-line, x- 1enght=7D Without n-plate. m- 1 3 6 10 K H 1 3 6 10 C 1 i y r u ' 0 I . ‘ \ \ ‘> ‘. \ I ': { \ ‘5 x . \ . i \ I ;. , x i ‘\ ‘K ‘\ \‘\ \\ i I". De». I ' “:\ ix _ \ \ \ \ l I!" t i '0': \i. g” 1. ‘4 .4" >- ‘ \~. 1: ‘~-. ‘_\) _ $\ \ ”t x. ‘~\ \ l! ‘ 0 S- “5.),” - \ \ x - v - I! 0.5 ..',-—/ - \ ; - u . Mr I ,I t I' e I! 4 I ,- L” ,’ ‘v' | ,’ ,’ 3 t ‘(i/ I/ I] l I, I." . l ,i /' ,~’ I’,-'.' /’ ill "I I L’ .’ ' / V l 0 03 0.00 0 00 0.03 1 0.000 0.03 0.02 0.000 Normalized axral turbulent normal stress Normalized shear stress (a) (b) Figure 1.22. Distribution of the flow variables for the case with different inlet turbulence: dashed-line, low inlet turbulence; dashdot-line, average inlet turbulence; longdashed-line, high inlet turbulence; solid-line, reference case. 41 CHAPTER 2 Reacting Flow Analysis Using Filtered Mass Density Function Methodology 2.1 Introduction THE OVERALL performance of a typical combustion system (its efficiency, stabil- ity, and compactness) is dependent on the coupled effects of the fuel/ chemistry, the geometry, the input / output flow conditions, etc. Normally, it is difficult to predict the combustor behavior under various operating conditions. Also, it is costly and time consuming to develop / test a new system via direct (trial and error) experimentation. On the other hand, the LES method provides valuable, time-dependent spatial data that can be used for development of new laboratory-scale combustors and the design of reliable control strategies for them. LES is considered as a new and promising com- putational methodology for analysis and prediction of engineering problems involving turbulent reacting flows. LES is superior to RAS (Reynolds average simulations) predicting mixing and unsteady reaction rates, flame-acoustic wave interactions, and other turbulent combustion phenomena. In the past few decades, LES of turbulent reacting flows has been the subject of widespread investigations [9, 47-59]. More recently, Colucci et a1. [9] and Jaberi et a1. [10] have developed a methodology termed the “filtered density function” (FDF), 42 or alternatively FMDF, based on an idea originally proposed by Pope [60]. The funda- mental property of the FDF or F MDF is that it explicitly accounts for probabilistic subgrid scale (SGS) scalar fluctuations. Colucci et a1. [9] developed a transport equation for the FDF in constant density flows in which the effects of unresolved convection and subgrid mixing are modeled consistent with those in “conventional” LES procedures. FMDF is an extension of FDF to variable density flows. One of the advantages of FMDF for turbulent reacting flows is that chemical reaction and other one-point processes are in closed form. The encouraging results generated by FMDF warrant its extension and application to more complex flows. In previous RAS/PDF modeling investigations, two different directions have been pursued; (i) the application to complex configurations, and practical combustion- device geometries [61, 62]. The numerical methodologies developed for this purpose are so far not appropriate for LES / FMDF due to accuracy issues and several other aspects discussed in chapter 1. Also, in most of the previous PDF modeling efforts unstructured grids have been used, in which searching and locating the Monte-Carlo elements, interpolation, and averaging procedures, are very expensive and not suit- able for LES; (ii) the extension and application of FMDF for the treatment of various variable density reacting flows with complex kinetics [63]. FMDF has been applied to some non-premixed, premixed and partially-premixed gaseous turbulent flames [63—67]. These studies have been limited to relatively simple geometries (e.g., ax- isymmetric and planner jets, often with simple one-step kinetics). This work deals with the development and application of a reliable and affordable LES/ FMDF methodology suitable for calculation of turbulent combustion in com- plex geometries. The LES/FMDF model is implemented via an efficient Lagrangian- Eulerian numerical scheme. High-order Eulerian finite-difference schemes are used and an efficient algorithm has been developed to search and locate particles in a multi—block hexahedral—structured grid system. In this work, an H-H and O-H grid 43 is used for simulations of a dump combustor, although a similar method can be used for other structured grid type systems regardless of the complexity of the system. The sudden expansion or dump combustor is one of the most common geometrical configurations for combustors, in which the flames can be stabilized by a stationary recirculation zone. A significant number of experimental [1 1—19] and numerical in- vestigations, [20—27] have been carried out on this flow configuration for a a better understanding of turbulence and combustion in dump combustors. The experimental work of Stieglmeier et al. [14] on axisymmetric sudden expansions with different dif- fusers indicates that the diffuser geometry strongly influences the flow field inside the dump combustor. Ahmed and Nej ad measured velocity and low-frequency combustor pressure oscillations in a ramjet dump combustor [15]. The data showed substan- tial differences between reacting and non-reacting cases. They also indicated that the intensity and frequency of the oscillations are dependent on the inlet velocity, combustor length, and equivalence ratio. Cole and Glauser [16] did hot-wire veloc- ity measurements in an axisymmetric sudden expansion recirculation zone. They performed a kinetic energy balance between the convection, production, turbulent diffusion, and dissipation rates. Lieuwen and Zinn [17] investigated the accuracy of the common experimental practice that estimates the unsteady pressure in the flame region of an unstable combustor from pressure measurements along the combustor wall. They also conducted computational studies to quantify the difference between the acoustic pressure at the wall and the flame. Their results showed that the duct pressure at the flame and the wall typically differ in magnitude by a factor of 5 and in phase by about 10 degrees. Chuang et al. [20] did RAS calculations for a hot flow in a sudden expansion dump combustor with swirl. They showed that RAS with the standard It — 5 model yields a poor prediction of the mean velocity, particularly the tangential velocity. They suggested it is necessary to use a Reynolds stress model. Guo et a1. [21] using a RAS model studied inlet swirling effects on sudden expansion 44 flow, with the focus on the unsteady flow behavior. They concluded that imparting swirl causes the mean flow to become unstable and oscillatory. They predicted pre- cessing vortex core phenomenon in their numerical simulations. Zhang et al. [22] applied an improved stochastic separated flow model to study particle-laden sudden expansion flows. The improved model requires fewer particles and a lower computing time as compared with traditional stochastic separated flow models. Chuang et al. [23] performed RAS predictions of the swirling flow of a dump combustor with a cen- tral V—gutter flameholder and six side-inlets. Their results indicated that for a fixed swirl strength, the length of central recirculation zone is decreased when the V-gutter angle is increased. The combustor outlet velocity in the reacting flow case is higher than for cold flow because of heat release effects. Stone and Menon [24] used LES to calculate combustion dynamics in a realistic, swirling-flow, dump combustor. They investigated the heat release effect on combustion dynamics, vortex-flame interaction, and vortex breakdown mechanisms. For a certain choice of system parameters (com- bustor length and flow-rate), they found that heat release has an attenuating effect on the combustion dynamics through a reduction in the turbulence intensity. Schliiter et al. [25] applied their integrated LES / RAS methodology to simulate swirling dump combustor flows. They proposed a method based on a virtual body force to impose Reynolds-averaged velocity fields near the outlet on an LES flow domain in order to take downstream flow effects (computed by a RAS) into account. Huang et al. [26] performed LES of turbulent combustion in a lean-premixed swirl-stabilized combus- tor. They identified and quantified several physical processes responsible for driving combustion instabilities, including the mutual coupling between acoustic wave mo- tion, vortex shedding, and flame oscillations. They also studied the mechanisms of energy transfer from chemical reactions in the flame zone and acoustic motions in the combustion chamber. Moin and Apte [27] used a second-order, unstructured flow solver to simulate a swirling, particle-laden cold flow in the coaxial geometry corre- 45 spending to the experiments of Sommerfeld and Qiu [18]. They also performed LES of evaporating droplets in a coaxial, non-swirling jet corresponding to the experiments of Sommerfeld and Qiu [19]. These investigations have significantly contributed to our understanding of turbulence, mixing, and combustion in dump combustors; our un- derstanding, however, is not yet complete. Among the various experimental studies, Gould et al’s. [11, 12] is desirable [13] from computational view point, since various turbulence variables in both nonreacting and reacting conditions are measured. In this work, the LES / FMDF methodology with a new high-order finite-difference (FD) flow solver and Lagrangian Monte Carlo (MC) particle method is used to nu- merically simulate turbulent combustion in the dump combustor experiment of Gould et al. [11, 12]. Local mass and volume consistency requirements between Lagrangian and Eulerian parts were assessed. Also, for isothermal, non-isothermal, and reacting cases, the filtered conserved scalar and temperature fields as obtained via MC and FD flow solvers are compared to establish consistency. Validation is done by comparing the reacting LES / FMDF results with experimental data [11, 12]. 2.2 Governing Equations As mentioned above, in the LES/ F MDF methodology a combinations of Eulerian and Lagrangian equations are solved together for velocity and scalar (temperature and mass fractions) fields. These are described below in separate sections. 2.2.1 Velocity Eulerian Field - Navier stokes equations for a compressible fluid flow The compressible and filtered form of the conservation equations of mass, and mo- mentum in generalized coordinate can be written as a as ac“; 3H - aJU—i-E‘Jr'a—n-‘FEZ—JS (2.1) 46 where J is the determinant of the transformation Jacobian, F, G & H represent the inviscid fluxes and F u, G; & H yare the viscous fluxes. The variables and quantities have been introduced and discussed in chapter 1. Thus, only the differences with the equations presented in that paper, will be addressed in this section. Theoretically, it is not required to solve the standard Eulerian energy and scalar equations in the LES / FMDF methodology since they can be obtained from FMDF. However in order to avoid some of numerical problems arising from coupling of Eulerian and Lagrangian schemes, the Eulerian energy equation is also solved by standard methods. This equation is similar to and can be described by equation (2.1). The energy equation is solved for A ~ NS 0,, A a2+62+w2 Azp E—Zhaqsa =p' i+ (2.2) The first four components of the source term vector S, are zero and the fifth component, which is the reaction source term in energy equation, is defined here as Se 2 —pk fAB exp(-—Ea / RT), where k f is the pre—exponential factor, Ea is the activation energy, and A, B denote the mass fractions of fuel and oxidizer. The source term is computed from the FMDF data and is closed. A one step chemistry model [68] is considered for propane combustion in air, viz., C3H8 + 5(02 + 3.76N2) => 3002 + 4H20 +18.8N2 (2.3) —E where the reaction rate is u') 2: BCCBHSCO26RT with B = 9.9 x 107[m3.mol_1.s-1] as the reaction coefficient, E = 30[kcol.m0l_1] is the effective activation energy, and R is the universal gas constant. CCBH8 and 002 are concentrations of the fuel propane and oxygen. A conserved scalar equation is also solved by the traditional Eulerian (finite difference) method in order to establish consistency between the La- grangian and Eulerian flow solvers. Species mass fraction equations in Eulerian field 47 are expressed by the following equation a a . we (we

€6L+ (PMngfWalL =_ (a;i>g— ax: +g, a=1,2,---,0 (2.4) This equation can also be written in the form of equation (2.1). 2.2.2 Scalar Lagrangian Field - FMDF Methodology for Chemically Re- acting Flows With ¢(x, t) as the scalar array, we can define the FMDF, denoted by F L, as +00 men) s / p< [¢,¢1 a II we. — ¢a(x, a} (26) 0:1 where 6 denotes the delta function and 11: denotes the composition domain of the scalar array. The term C [qb, ¢(x, t)] is the “fine-grained” density [69, 70], and equation (2.5) implies that the FMDF is the mass weighted spatially filtered value of the fine-grained density. The integral property of the FMDF is such that +00 +00 I I I / FL(t/J;x,t)d1/2 = / pg 5 (2-8) 48 Equation (2.8) implies (2') For 00.0 = c, we mm = c (29) (ii) For Q(X, 1t) E Q(¢(X,t)), gi= QM) (2-10) +00 (ii/2'.) Integral property : / (Q(x, t)]¢>eFL(’l/J;X, t)dz/J —00 = (p(X.t)>g(Q(X,t)>L (2-11) where c is a constant, and Q(¢(x, t)) E Q(x, t) denotes the case where the variable Q can be completely described by the compositional variable ¢(x, t) E [(251, 052, . . . ,qbg]. From these properties, it follows that the filtered value of any function of the scalar variables (such as p E fi[¢(x,t)] and Sa 5 Sa[¢(x,t)] ) is obtained by integration over the composition space. It is noted that the mass weighted conditional filtered mean reduces to the conditional filtered mean [9] when the density can be completely expressed in terms of the compositional variables. By applying the method developed by Lundgren [71], Pope [72], and O’Brien [69] to the 4) equation, a transport equation is obtained for the fine-grained density [9]. The transport equation for F L(”¢/J;X, t) is obtained by multiplying the equation for the fine grained density by the filter function G (x, — x) and integrating over x’ space. The final result after some algebraic manipulation is drum»)+atgFL(w;x,t>i _a_ __1_6J.-“ at (9332' (91/20 “(1’) 82:,- |¢>€FL(1/);x,t)] 8[Sa(¢>FL(¢;x,t)1 — 8% (2.12) This is an exact transport equation for the F MDF which can be solved by standard (Eulrian) methods. However, Eulerian methods are very expensive due to added dimensions. In practice only Lagrangian methods can be used. The last term on the right-hand side of this equation is a closed-form chemical reaction term. The 49 unclosed nature of SGS convection and mixing is indicated by the conditional, filtered values. These terms are modeled in a manner consistent with Reynolds averaging and conventional LES in non—reacting flows. The convection term is decomposed as follows: (UiIWgFL = (UDLFL + [(UZ'IWe - (ui>LlFL (2-13) where the second term on the right—hand side denotes the influence of SGS convective flux. This term is modeled as dude/w — L1FL = —’Yta(Fg::p>€) (2.14) The advantage of the decomposition [equation (2.13)] and the subsequent model [equation (2.14)] is that they yield results similar to that in conventional LES [73, 74]. For example, the first Favre moments corresponding to equations (2.13) and (2.14) (“Wan = L<¢0>L + lL ’ (ui)L(¢a)Ll (2-15) (p>gla>L- L1= — v. 5(3ij (2.16) The term within brackets in equation (2.15) is the generalized scalar flux. The closure adopted for the SGS mixing is based on the linear mean square estimation (LMSE) model [69, 75], also known as the IBM (interaction by exchange with the mean) [76] area em J- rel—22:”) 8 +37%) [QmU/Ja—Wall, )F L] (2.17) where Qm(x, t) is the “frequency of mixing within the subgrid”, which is not known a prion. This frequency is modeled as Om = 09(7 + 7t)/((p)gA%;) here but other models can be used. For the first term on the right-hand side of equation (2.17), an 50 additional minor assumption is made 3 (9 F A 8 (9 F p 3:,- 62:2- 02:,- 8:12,- This approximation is not necessary for the treatment of F MDF . It is only adopted to establish consistency between the FMDF and conventional LES. With these ap- proximations, the modeled FMDF transport equation becomes BFL aliuilLFLl _ 8 ’37 + T — 5;; (7 + 7t) 6 F 3 _(_Ié_::.fl6_) + 5% [omwa — <¢a>L>FLl Elgar—fill 01/10 This equation may be integrated to obtain transport equations for the SGS moments. (2.19) The equations for the first subgrid Favre moment, ($0) L, and the generalized subgrid variance, 0% = will) L — <¢(a)>% are, respectively, a a “2' 01 a ((9)3; >L) +8(

gL) :ai%[(7+7t)0[4] +

geUc21+ 2(Ple (<¢(a)s(a)>L — <(..)>LL) (2.21) where the subscripts in parenthesis are excluded from the summation convention. These equations are identical to those which can be derived by filtering equation (2.4) directly, and employing consistent closures for the subgrid flux and the dissipation. In direct moment closure formulation, however, the terms involving (5(1)); remain unclosed. This indicates that for non-reacting flows the FMDF method is consistent with standard LES methods at the first moment level (and possibly second, if the second moment equation is solved by standard methods). 51 2.3 Numerical Solution As mentioned, the LES / FMDF methodology is based on two different mathematical models, namely (i) the conventional LES equations for the velocity field, and (ii) the FMDF for the scalar field. The numerical methods used for each part are solved very different. These methods are discussed in detail below. 2.3.1 An Eulerian Finite Difference Method for Velocity Field The numerical solution procedure, discretization of the computational domain, and the temporal & spatial derivatives for finite difference part of LES / FMDF have been addressed in chapter 1. For any scalar quantity, such as a metric, flux component or flow variable, the spatial derivative of the variable in the transformed domain is obtained using fourth-order compact difference formula [3]. A tridiagonal system is solved using the Thomas algorithm. The derivatives of both inviscid and viscous fluxes are obtained by first forming the fluxes at the nodes and subsequently using a fourth- order compact differencing formula. The numerical noise generated at very high wave number/ frequency modes is removed by using an implicit filtering procedure. The filtering Operation suppresses numerical instabilities arising from “very small scale interactions”, mesh nonuniformities and boundary condition effects. All the time derivatives in governing equations are approximated by using a low-storage three— stage Runge-Kutta scheme. 2.3.2 A Lagrangian Monte Carlo Method for Scalar Field The scalar filtered mass density function methodology is employed as a subgrid scale (SGS) closure in large eddy simulation. The scalar field, 45 E 450, a = 1, 2, . . . ,N3+1 (for N3 species mass fractions and enthalpy) is obtained from the joint scalar FMDF [60], The FMDF transport equation is solved using a “Lagrangian Monte Carlo” procedure [77]. In this procedure, the FMDF is represented by an ensemble of 52 particles which are transported in the physical space by the combined actions of large scale convection and diffusion (molecular and subgrid). 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) [63, 77]. While it is potentially, or eventually, possible to simulate F MDF exclusively via the effective Monte Carlo methods described below, the most practical procedure is “hybrid” methods in conjunction with “deterministic” schemes, such as finite dif- ferences, finite volumes, spectral, and other methods. There has been significant progress in the develOpment of these high-accurate methods. A hybrid method would make use of this accuracy, which is not yet achievable in MO. Moreover, the influence of the MOS dispersion and statistical errors is less significant than in hybrid methods. These issues are investigated in detail in the context of RAS/FDF [78, 79] and con- stitute a major element of the computational procedure in LES/FMDF [63, 80, 81]. The operational procedure is described in the next subsection. The MC particles are distributed randomly and are free to move anywhere within the domain. This transport is Lagrangian, thus the solution is free of the constraints associated with typical FD simulations on fixed grid points. The attributes of conventional LES and F MDF are shown in figure 2.1. The flowchart also shows the redundant variables and the necessary consistency between the two computational parts of the LES / FMDF flow solver. In the Lagrangian Monte Carlo procedure [70], equation (2.19) is solved indirectly via equivalent stochastic (diffusion) equations. Here, the spatial transport of the FMDF is represented by the general, diffusion process governed by the following stochastic differential equation [82, 83]. dX - ”l (t) = Di(X(t), t)dt + E(X(t), t)sz-(t) (2.22) where X,- is the Lagrangian position of a stochastic particle, D,- and E are the “drift” 53 and “diffusion” coefficients respectively, and Wz- denotes the Wiener process [84]. The drift and diffusion coefficients are obtained by comparing the Fokker-Plank equation corresponding to equation (2.22) with the spatial derivative terms in the FMDF transport equation [ (2.19)], 1 304-70. M 3152' (2°23) E —: \/2(7 + we», D.- s L + < The subgrid mixing and reaction terms are implemented by altering the compositional makeup of the particles (we? dt = «met — ((150) L) + sad/2+) (2.24) where ¢lf = ¢Q(X(t), t) denotes the scalar value of the particle with the Lagrangian position vector Xi- The solutions of equations (2.22) and (2.24) yield the same statistics as those obtained directly from the solution of FMDF transport equation according to the principle of equivalent systems [70, 85]. Lagrangian methods have been used for simulation of a wide variety of stochastic problems [86] and have benefited significantly from modern developments in SDE solver technology [87]. They have been the primary means of solving the PDF in RAS [70, 77, 88—90] and thus far the primary method of choice for solving FMDF in LES. Typically, the method is implemented by representing the FMDF by an ensemble of Np particles. These particles carry information pertaining to their positions, x(n) (t), velocities, u(n)(t), and scalar values, ¢(n)(t), n = 1,...,Np. This information is updated via temporal integration of the modelled SDEs (2.13). The simplest means of performing this integration is by the Euler-Maruyamma approximation [91]. For the SDE d2: = Ddt + BdW, this approximation advances the position of the nth particle :rn(t) from time level t k to t k +1 according to $?(tk+1) = 159%) + D?(tk)At + Bn(tk)(At)1/2(zn(tk), i = 1, 2, 3 (2.25) 54 where (in (t k) are independent standardized Gaussian random variables. This formu- lation preserves the Markovian character of the diffusion processes [92, 93] and fa- cilitates affordable computations. Higher-order numerical schemes are available [91], but one must be cautious in using them for LES [9]. Since the diffusion term in the SDEs depends on the stochastic processes, the numerical scheme must be consistent with It6—Gikhman [94, 95] calculus. Equation (2.25) exhibits this property. Statistical information, e.g. filtered values, at any point are obtained by consider- ing an ensemble of NE computational particles residing within an ensemble domain of side length A E centered around the points. This is necessary because, with probabil- ity one, no particle will coincide with the point considered [85]. For reliable statistics with minimal numerical dispersion, it is desired to minimize the size of the ensemble domain and maximize the number of MC particles [70]. In this way, the ensemble statistics would tend to the desired filtered values. A larger ensemble domain (AE) decreases the statistical errors but may increase the dispersion errors, which mani- fests itself in “artificially diffused” statistical results. This compromise between the statistical accuracy and dispersive accuracy as pertaining to Lagrangian Monte Carlo schemes implies that the optimum magnitude of A E cannot, in general, be specified a priori [70]. This limitation does not diminish the capability of the procedure but exemplifies the importance of the parameters that govern the statistics. Transfer of information from the grid points to the MC particles is accomplished by interpola- tion. The transfer of information from the particles to the grid points is accomplished using ensemble averaging as described above. With a hybrid method as such, some of the quantities are obtained by MC, some by FD, and some by both. That is, there is a “redundancy” in the determination of some of the quantities. In general, all of the equations for the filtered quantities can be solved by FD, in which all of the unclosed terms are evaluated by MC. This process can be done at any filtered moment level [81]. 55 Presently, the most serious issue associated with the F MDF is its computational cost. In most cases, LES / FMDF is more expensive than the conventional LES meth- ods (e. g. finite difference). This overhead is expected considering all of the SGS statis- tical information that LES / FMDF provides in comparison to other schemes. Here, a new approach is described in which the computational cost of LES/FMDF for general geometries is minimized. This approach is described below. With the development of highly parallelized computational schemes, it is predicted that LES/FMDF will distinguish itself as a major practical combustion prediction tool. 2.3.3 Particle Tracking, Interpolation and Averaging An effective method has been developed to search and locate Monte Carlo elements or particles in a structured, multi-block Eulerian grid system. This method makes FMDF feasible for LES of flows in complex geometries. The search and locate 0p- eration for Monte Carlo elements is an expensive Operation when an unstructured Eulerian grid system is used. However by dividing a complex geometry into different blocks and using structured grids for each block, the cost of this operation can de- crease dramatically. As shown in figure 2.2 for the special geometry considered in this study (dump combustor) searching and locating can be done for outer and interior blocks, which are both structured. This algorithm can be used for any complex geom- etry by considering an appropriate formulation with respect to the non-uniform but still structured grids. It is important to mention that using a structured, multi-block, high-order flow solver is not only advantageous in terms of accuracy and simplicity, it has extra benefit of being computationally much less expensive when FMDF is used. Unlike an unstructured grid system, where particles are tracked element by e1- ement, for structured grids particles are tracked in three different grid line directions separately. Figure 2.2 illustrates the experimental setup and the grid system. In order to 56 locate the particle, the first step is to determine the block in which the particle is found, so each particle has an index showing the corresponding block. The particles on the common boundary are considered on the outer block. Figure 2.3(a) shows the grid lines, Monte Carlo particles, and imaginary control volumes formed around the grid points. For simplicity, a two dimensional picture is considered, even though the procedure is three-dimensional. Solid black and green lines show the finite difference mesh and the control volumes, respectively, while solid circles are the Monte Carlo elements. The control volumes around the grid points are formed by connecting the mid-points of grid-line segments. In figure 2.3(a) the hexagon “abode f” is formed around point “0”, which is the finite difference grid point. The particles located inside hexagon “abode f” are considered in ensemble averaging, and point “0” is a good representative for the center of the hexagon cell. For interpolation, finite difference grid points are used. As shown in figure 2.3(a), if particle “p” is located inside the quadrilateral “qrst” the values at the corners will be used. To perform ensemble averaging and interpolation for any particle, the hexagon “abode f” and the quadrilateral “qrst” should be determined. As the procedure is the same to determine the hexagon “abode f” and the quadrilateral “qrst”, here I only describe the method to determine quadrilateral “qrst”. When the block in which particle “p” is found has been specified, if 3 indexes of C‘ I, point q are determined, the quadrilateral “qrst” is specified. In other words, we “ ” only need to determine i,j, and k (components of the coordinates) of point q . w” 1 To determine the component of “q”, the auxiliary vectors 57), st], and fl are defined (figure 2.3(b)). A is the normal vector to st and could easily be formed [the direction of .4 could be either way]. With this, G(q) is defined as follows ((9 E dot product between two vectors): C(Q) = (870 G ANS” G ) (226) If C(q) is positive, the particle “p” and point “q” lie on the same side of st, meaning “p” is located at the east-side of st. We can assume the sweeping start point to be i = 1 and that sweeping is performed in the direction of increasing “i”. The first grid point with G’ > O is “q”. A similar procedure is used in other 3 grid-line w” 1 directions as shown in figure 2.3. As the component of “q” is determined, the “k” component could be specified by searching in a “grid-plane”. By knowing the “i” and “k” components, “j” is determined only by sweeping in a segment formed by grid lines. To reduce the cost of the searching procedure, the sweeping start points are determined by previous values (which are stored for each particle) and particle possible movement direction, which is based on the velocity of the particle. The above numerical scheme is applicable to any curvilinear grid. If in a block, the grid lines are straight lines, the searching procedure becomes even more less time consuming, even though this does not limit the algorithm. A similar algorithm has been successfully implemented for a different section of the combustor, including its inflow nozzle. To interpolate the properties needed in particle locations, a (coupled) linear inter- polation is applied in three grid-line directions. A two-dimensional schematic picture is shown in figure 2.3(d) to describe the method. As shown in this figure, par- ticle “p” is located in the arbitrary quadrangle “qrst”. Lines “ab” and “mn” are placed so that if there were grid lines crossing the particle position “p” in the grid system, they would be “ab” and “mn”. This is accomplished by considering trans- formation equations between physical and transformed domains. In the transformed domain,quadrangle “qrst” will be a square and the grids are uniform. The velocity values are interpolated at “a”, “b” , “m” and “n”, based on the values at the corner of quadrangle “qrst”. The velocity for point p is computed by an interpolation on “ab” and “mn”. 58 2.3.4 Particle Number Density Control Ensemble averaging of particles is employed in a way that there is an acceptable number of particles in ensemble domain. As particles move in physical space, the ensemble averaging domain is chosen so that Np,min S NP 5 Np,ma$. This means that if the number of particle is lower than the minimum number specified, a bigger ensemble domain is considered, or if it exceeds a maximum value, a smaller ensemble domain is considered. This approach is an alternative for cloning and annihilation in which property distributions are not exactly preserved. Different values of Np,min were examined and it was insured that Np,min was sufficiently high. 2.4 Results The focus of this work is on the LES of turbulent flow in a dump-combustor or axisymmetric sudden expansion configuration. Gould et al. [11, 12] have conducted various laboratory experiments on similar configuration for both nonreacting and reacting flows. In this section, the performed consistency checks between Eulerian and Lagrangian models are discussed for reacting and non-reacting cases. The numerical data ob- tained by LES/FMDF for the reacting case are compared to experimental data afterward[12]. The simulations for the reacting case are performed for the refer- ence case. The resolutions for the reference case is 159x47><58 for the interior block, and 159x17 x 17 for the outer block. 2.4.1 Nonreacting Simulations - Consistency of FMDF The objective of the results presented in this section is to demonstrate the consistency of the FMDF and its Lagrangian Monte Carlo scheme with the finite—difference and conventional LES models. Both isothermal and non-isothermal, non-reacting cases 59 are studied. The filtered temperature and conserved scalar variables as calculated by both LES-FD and FMDF-MC models. The consistency for the non-isothermal, reacting case will be discussed in ‘reacting simulations’ subsection. The non-reacting, isothermal flow in an axisymmetric, dump combustor is similar to cold-flow simulations presented in chapter 1. However, here I solve an additional equation for the conserved and passive scalar. This case examines mixing without high temperature and density variations. At the initial time, the scalar value is set to 1.0 for core region, to zero in the region surrounding the core. The scalar value is 1.0 for the incoming flow and as shown, as time increases, the incoming flow (with mass fraction of 1.0) later dominates the computational domain. For this case, it is possible to compute filtered conserved scalar and temperature fields from the Monte Carlo particles and finite difference grid points, even though the temperature variations are small. Contours of the instantaneous and time-averaged conserved scalar variable as obtained by LES-F D and FMDF-MC are compared in figure 2.5. The results in this figure show how the mixing process is captured by the two models similarly showing an excellent consistency. The time-averaged values show that consistency between the two methods exists at all times despite the turbulence. As expected, there is no significant temperature/ density variations in the “isother- mal” non-reacting flow. The consistency between LES-FD and FMDF-MC tempera- ture fields can still be well established. Figure 2.5 shows the instantaneous and mean temperature fields as computed by Eulerian finite difference and Lagrangian Monte Carlo methods. Even though, the maximum temperature difference in the flow is only 3 percent of initial (or incoming) flow temperature, the trends are similarly captured by both methods, again showing the consistency of the F MDF methodology and its numerical solver with conventional LES methods. To further establish the consistency and the accuracy of the LES-F D and FMDF- MC methods, a more complex non-isothermal non-reacting case is also considered. In 60 this case, the flow is initiated with non-uniform density and temperature distributions. The temperature is set to be 2T0 (T 0 = 300K) for the core region as well as incoming flow, while it is To in the region surrounding the core. Initially, the pressure has a uniform distribution, so the density distribution is the opposite of the temperature distribution. Figure 2.6, shows a comparison between filtered densities computed via two dif- ferent models. The density from FMDF is calculated based on the weighting of the particles located in the ensemble domain (control-volume cell)by averaging over the particles. The LES-FD density is directly computed and obtained from the grid points. To check the consistency between the two fields, the instantaneous density fields are compared in figure 2.6. has been calculated from the discrete points (par- ticles). Acceptable agreement is observed; however, the values obtained from the MC particles seem to be noisy. It is to be noted that the values of Np (number of particles) perviously used in non-reacting simulations and the consistency check, are less than those used for reacting simulations to magnify the error. The results clearly show that the mean velocity field used to advect the particles satisfies the mean continuity equation. Figure 2.7 shows a comparison between the MC and FD scalar values. In the non- isothermal case, the scalar initialization and inlet conditions are the same as those in the isothermal case, but the influence of variable density on mixing is important. When the incoming flow has a high temperature and mixing is taking place at the shear layer, the Reynolds number drops locally in the layer and thus the flow evolution is different than for the isothermal case. Here the small vortical motions are damped by the high viscosity of the hot fluid and the modification of the instability modes due to density variations. Nevertheless, the data indicates a high level of consistency between the results obtained by the FD and MC methods. The instantaneous and mean temperature contours in a 2D plane as computed by Eulerian LES-FD grids and 61 Lagrangian FMDF-MC particles are compared in Figure 2.8. The results presented in this figure are similar to the others in this section and clearly indicate that my Lagrangian FMDF flow solver and models are fully consistent with our Eulerian models as long as there is no reaction. This suggests that the finite difference and Monte Carlo solvers are both reliable and accurate. 2.4.2 Reacting Simulations - Accuracy of LES / FMDF Large-eddy simulations of reacting turbulent flow in a dump-combustor configura- tion (axisymmetric sudden expansion) are presented in this section. In the reacting simulations, a lean premixed propane-air mixture with F=O.5 is ignited by a high- temperature initial condition. Initially, a region with the thickness of H / 2 (where H is the step hight), adjacent to the side-wall of the combustor is set to be high- temperature products. The flow Reynolds number, based on inlet diameter and mean bulk inlet velocity, is 115,000. A one—step global reaction model[68] is used for propane-air combustion. Figure 2.9 shows the 3D iso—levels of vorticity magnitude and the Monte Carlo particle distribution, and the 2D plane contours of fuel and heat release. As figures 2.9(b) and 2.9(c) indicate, a large amount of unburned fuel and oxidant exit the domain which is consistent with the experiment mentioned by Gould et al. [12]. By increasing the inflow turbulence or swirling flow at inlet, more fuel will burn and an improved combustion efficiency is achieved. Further investigation was done on turbulent combustion in other configurations with the same LES / FMDF computational model and code. The results of these simulations are not presented here. Figure 2.10 shows contours of mean axial velocity and stream lines. The recir- culation zone (defined by the streamline with zero mean axial velocity) is found to be thinner and shorter in the reacting-flow case, due to the reaction and volumetric expansion in the radial direction. To examine the consistency between Eulerian and Lagrangian models for the re- 62 acting case, instantaneous and mean (time-averaged) filtered temperatures, obtained by LES-F D and FMDF-MC, are compared at different downstream locations. As indicated before, in any single simulation (reacting and nonreacting) the filtered tem- perature is a redundant variable that can be obtained from both the Monte Carlo particles and finite difference grids. The values obtained by these two methods should be the same in theory. As observed in figures 2.11, there is indeed good agreement between the two methods. Also figure 2.12 shows contours of instantaneous temper- ature at the 220 plane as obtained by FD and MC data. These results are calculated after the flow has reached the statistically stationary condition, which again indicates good agreement between MC and FD. Figures 2.13 and 2.14 show the comparison between LES / FMDF and experimen- tal results for mean axial velocity and mean temperature at different axial locations. As the results suggest, the central core region upstream of a: / H = 5 in the reacting- flow case is similar to that in the cold-flow case because the temperature in this region is close to the inlet temperature. In the reacting-flow case, the core flow region stays cool and the center-line velocity doesn’t decay as rapidly as in the cold flow due to the low pressure gradient at the center-line. Higher mean axial velocities in the share layer are found in the reacting case due in part to volumetric expansion. It is interesting to note that the reacting-flow center-line velocity does not decay nearly as soon as does the cold-flow center-line velocity. The relatively low temperature at this location suggests that this result is not caused by volumetric, sudden expansion effects alone and might arise from some problems at the inlet boundary. Figure 2.14 shows that the core region stays relatively cool and that the maximum temperatures are achieved in the recirculation zone and developing boundary layer downstream of the reattachment point. As Figure 2.13 indicates, the predicted mean axial velocity is higher than the measured shear layer values. Again this could be caused by boundary conditions 63 at the inlet. A nozzle was added to the computational domain, right before the combustor, consistent with experiment, to examine the accuracy of the boundary conditions. Schematic perspectives of a fixed Eulerian multi-block gird system are shown in figure 2.15. As for the case with the added nozzle, the Monte-Carlo elements are only considered after the sudden expansion point, and the FMDF equation is solved for the combustor region; the consistency between Eulerian and Lagrangian data (figure 2.16) indicates that the F MDF results are reliable. The mass fraction of the species is set as the incoming flow condition in the nozzle, since the Monte Carlo elements are injected at the sudden expansion point. The number of MC elements and their weighings are adjusted based on the velocity and density at the sudden expansion point. The contours of fuel in figure 2.17 show that the flow has a similar behavior to the case without the nozzle, and that the reaction is taking place in the same way. Although, the time-averaged and statistical values determine the effects of the added nozzle, precisely. Figure 2.18 shows a comparison with the experiential data for mean axial velocity. As seen, the mean axial velocity distribution is improved in shear layer. This figure indicates that the discrepancy between predicted and measured values in the shear layer is due to non-realistic inlet boundary condition for the case without the nozzle. Figure 2.19 shows the mean temperature comparison with experimental data, for the case with the added nozzle. The comparison indicates that adding the nozzle to the computational domain improves the predicted values for mean temperature as well. The maximum, normalized temperature achievable in this simulation is bounded by the adiabatic flame temperature and is approximately 5. The apparent asymmetry for temperature on the center-line at the downstream locations x/H=12 may be due to large scale structures, or it may be caused by the inherently non-symmetric nature of the heated horizontal cylinder. The results presented for the case with the nozzle are obtained for a larger computational domain in the x—length=9D ( where D is the diameter of the nozzle at the sudden expansion point). This was done to reduce 64 the effect of the outlet boundary condition on the numerical results. However, the effect of the outlet boundary does not seem to be very significant as long as the flow and pressure oscillations are allowed to exit the computational domain. My analysis indicates that the high mean axial velocity and high turbulent intensities predicted for the case without the nozzle are caused by the inlet boundary condition, not outlet boundary condition. Overall, my numerical results are consistent with the experiment and exhibit some interesting physical features that are summarized in the following. The data indicates that the axial center line pressure gradient in the reacting flow case is approximately four times less than in the cold flow, consistent with an axial momentum balance. This is due to increased velocity (and temperature) in the rapidly growing shear layer, which preserves the central core region. Figure 2.20 shows normalized axial turbu- lent normal stresses and normalized fluctuating temperatures. The turbulence level is found to be lower in the reacting flow case over most of the flow field, possibly indicating either reduced production and/or increased dissipation of the turbulent kinetic energy. The locations of the peak turbulence levels are shifted radially out- ward, consistent with the location of the dividing streamlines defining the thinner recirculation zone in the reacting flow. The reacting-flow field was more isotropic than the cold flow. The results indicate that the maximum Reynolds stress in the reacting flow is an order of magnitude lower than in the cold-flow. The production of turbulence generated by shear is much lower than in the cold flow. Low shear in re- acting flow occurs by an increased turbulent dissipation rate (due to the temperature dependent viscosity) and by the volumetric expansion in reacting flow. The dynamic viscosity increases by a factor of approximately 3 at the elevated temperature existing in the combustion chamber. Likewise the kinematic viscosity increases by a factor of approximately 15 (since density decreases by a factor of 5) in the combustion cham- ber. Increased turbulent dissipation may reduce eddy coalescence and structures in 65 reacting flows and shear layers. In order to see the influence of n-plate in reacting simulations, similar to the cold flow simulations, an n—plate is added to the outlet of the dump combustor at different locations from inlet. The added n-plate decreases the turbulent intensities in core and surrounding regions, on the contrary to cold flow simulations. This is because the added n—plate confines the hot products and increases the temperature at the surrounding region. Thus, as a effect high viscosity in that region the turbulent intensities decrease. 66 2.5 Figures (p) 3 Conventional LES g Solver (LES-FD) (“ilL A :: 5(5) — 2:};111250) E (TlL /€/ “ “”3 W "L- £21] (a...) L (conserved scalar) o.a.--- IW—I‘O 1 r » 1 (We [ FMDF Solver ] (m " < (080:) 8 ($0) L J Figure 2.1. The attributes of LES/FMDF methodology and its LES-FD and FMDF subcomponents. 67 0.914 m 0.676 J 0.762 m J Honeycomb Flow Suighteners : Quarts Test Section 7" l l 1 [IL I] l U U U L? DOP Seeder 5:1 Contraction Extension Duct (a) Q. \ s \ 0‘ \\ \i ”3' fi sa‘o \ \Q‘gfi‘i‘l‘“ tit] \ . Figure 2.2. Schematic perspective of (a) experimental setup, (b) dump combustor test section, and (c) the 3-D grid layout. 68 (b) |i(q) w - . a, 999'" direction (C) component of point 5‘ ” q Figure 2.3. Schematic description of (a) grid ponits, Monte Carlo elements and control volumes around grid points, (b) determining i component of point “q”, (c) determining j , and (d) the two dimensional interpolation scheme. 69 (C) (d) Figure 2.4. Contours of the instantaneous and time-averaged values of the filtered con- served scalar in an isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values. (C) (d) Figure 2.5. Contours of the instantaneous and time-averaged values of the filtered tem- perature in an isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values 70 ( (b) Figure 2.6. Contours of the density in a non-isothermal non-reacting axisymmetric dump combustor as obtained by LES / FMDF, (a) FD values, and (b) MC values (calculated based on weighting of the Monte Carlo elements). (C) Figure 2.7. Contours of the instantaneous and time-averaged values of the filtered con— served scalar in an non-isothermal non-reacting axisymmetric dump combustor as ob- tained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time- averaged FD values, and (d) time-averaged FD values. 71 (C) (d) Figure 2.8. Contours of the instantaneous and time-averaged values of the filtered tem- perature in an non-isothermal non-reacting axisymmetric dump combustor as obtained by LES/FMDF (a) instantaneous FD values, (b) instantaneous MC values, (c) time-averaged FD values, and (d) time-averaged FD values. 72 ((1) Figure 2.9. (a) Three-dimensional iso—levels of the instantaneous vorticity magnitude and Monte Carlo particles, (b) three-dimensional iso-levels of the fuel distribution, (c) two- dimensional planar contours of the fuel, and (d) two-dimensional planar contours of the heat release; in an axisymmetric dump combustor as obtained by LES / FMDF. 73 Figure 2.10. The mean axial velocity contours and stream lines at the z=0 plane. 74 1 ~ g 1 r ’ j : ] l/r 4 h ‘,....:_:’,‘ I— L I!) 0‘5 * ......;.:_::___,__,.S' 015 . “I ”-....- My,” / I [// ------ FD ' _ - - _ - - m ~"‘ ~ MC ] MC ll '5 l l l! 1 1 4 o 2 4 o 2 4 m° T /TO (a) (b) 1 — 1 P I l , , r "I \’\/ (’1’ \‘__ iwmh“ " b- -... »;‘—-4"."-‘ ' "-... ‘J‘J'N “H h ‘ I] 0 5 [I I I 0 5’ ‘x/fij/ ..1 . _ '*-- v‘ ‘ 1”,” .l ------ 1'1) ’,,/ I] ________ MC 7 ------ H) V " ————- MC LI] H; I] I! .] .: L 1 1 12 k 1 k l o 2 4 0 2 4 m0 T no (C) (d) 1 — I . 1 r 'I s / J r” s. #7,.) ‘ ‘ \— [l’ _ : t- ,I /’ __.,. ”-5.15; I 0.5 » l’,,.;-;—;»~.--~.-~..r.r.r-—r o 5 _ I / r ------ a [i ] / . \]\ ------ 11) 7 y f:‘ - - _ ~_ ML / _. “mung- ‘ - I 1 PP ~ ‘ T ‘ 1 l 1 A A 1 o 2 4 0 2 4 Tfl‘0 Tin/To Figure 2.11. The instantaneous and mean values of the filtered temperature obtained from FD and MC data at different downstream locations. 75 (b) Figure 2.12. Contours of the instantaneous filtered temperature in reacting case, (a) FD values (b) MC values. X/ll’l’:| 1 o - 8 . O 8 ' 8 g?- 8 H . c- C? 0.5 - dim 8 . ' g 0 g; s 1 o 1 ' . - Mean ax1al veloc1ty Figure 2.13. Mean axial filtered velocity and its comparison with experimental data at different downstream locations (0, experiment [reacting flow]; 0, experiment [non-reacting flow]; ——, LES). 76 O '6. - y-L. l: .... r/ R1 .° 'sII on Mean temperature different downstream locations. Figure 2.14. Mean filtered temperature and its comparison with experimental data at case with added nozzle. (b) Figure 2.15. Schematic perspectives of ((a) 3D and (b) 2D) discretized domain for the 77 (b) Figure 2.16. Contours of the instantaneous filtered temperature obtained by three- dimensional LES/FMDF in the reacting case, (a) FD values and (b) MC values, for the case with added nozzle. 9516'?! ’S-éfi-f'f'é‘ifi'f'ifi'fig 12'; : > Figure 2.17. Two-dimensional planar contours of the fuel as obtained by three—dimensional LES/FMDF, for the case with added nozzle. [x/l-l‘TI l_ n 12 15 O O O ' O O 02‘- O O t . 8 o s - 0° " _ 0 § .' a :' O O E a - 1 ) 1 ' . . Mean ax1al velocuy Figure 2.18. Mean axial filtered velocity and its comparison with experimental data at different downstream locations, for the case with added nozzle (0, experiment [reacting flow]; 0, experiment [non-reacting flOW]; "'4 LES). 78 x5} 1* l 2 3 4 0 Mean temperature Figure 2.19. Mean filtered temperature and its comparison with experimental data at different downstream locations, for the case with added nozzle. 1 - 3 _ 5 12 E I I I I I I I I I I I I I I '. 'u '. I I I I I I I I .0 ' e 3 1 6 0.05 Normalized axial turbulent normal stress (a) 1 Normalized fluctuating temperature ('0) Figure 2.20. Distribution of the various turbulent quantities(e, experiment;~—-—, LES). 79 CHAPTER 3 Conclusions Chapter 1 deals with the development and validation of a high-order, compact, finite- difference, numerical methodology for large eddy simulation (LES) of compressible turbulent flows in complex geometries and its applications to the numerical investi- gation of cold flow in a dump combustor. The numerical schemes often used for flows in complex configurations are based on low (second) order, finite-volume schemes. These schemes are not usually appropriate for LES. The numerical scheme developed here includes high-order, compact differencing formulas and multi-block strategy cou- pled with high-order, implicit, spatial filtering. The flow solver is validated via direct and large—eddy simulations of isotropic turbulence and round/ planar jet flows. The computational model and algorithms are further validated by carrying out an exten- sive investigation of the flow in a dump combustor. The analysis includes a detailed investigation of the effects of boundary conditions, the subgrid scale model, inlet tur- bulence, an added inlet / outlet nozzle and various physical / geometrical parameters. The results were found to be in agreement with those obtained via other validated, high-order, numerical methods and with the available experimental data. This indi- cates that my high-order multi-block flow solver is a quantitatively accurate and reliable numerical scheme for turbulent-flow simulations in complex geometries. It is also shown that adding a nozzle to the computational domain reduces the “non— physical” effect of the inlet boundary condition upstream of the dump combustor. 8O This is consistent with experiments. Investigations of the effects of various inlet boundary conditions indicate that for sufficiently small combustors an n-plate cause more turbulence in the combustor. A similar effect on the flow is observed when the inlet turbulence is increased. The computational model has been extended for LES of turbulent reacting flows by using the filtered mass density function (FMDF) methodology as a subgrid-scale (SGS) scalar/ combustion model [10]. The LES / FMDF methodology and the results obtained for a reacting dump com- bustor are discussed in chapter 2. A Lagrangian/Eulerian numerical methodology is developed and implemented for large eddy simulation of turbulent reacting flows in general geometries. The subgrid-scale combustion model is based on the filtered mass density function approach that is extended here for LES of turbulent combus- tion in complex geometries and applied to the flow in a dump combustor. An efficient algorithm is developed for Lagrangian Monte Carlo F MDF numerical scheme equa- tions, which makes the application of FMDF to complex geometries feasible. The applied scheme offers much higher computational efficiency than the algorithms used in unstructured-mesh flow solvers. The new approach to the Lagrangian Monte Carlo solution of the FMDF in structured grid systems, coupled with a previously devel- oped high-order, finite-difference, multi-block flow solver, is a promising tool as the results presented in this paper suggest. The consistency, convergence, and accuracy of the FMDF and the Monte Carlo solution stochastic system are assessed. The con- sistency between Eulerian conventional finite difference LES and Lagrangian Monte Carlo FMDF methods for non-reacting isothermal and non-isothermal flows as well as reacting flows in a dump combustor are established. The results show that LES—FD and FMDF-MC results are consistent in all test flow conditions. Additionally, the LES / F MDF results for a reacting flow in a dump combustor show favorable agree- ment with the laboratory data of axisymmetric dump combustor. However, it is possible to improve the numerical predictions by adding an inlet nozzle component 81 to the computational domain. With this, the effects of a non-realistic inlet bound- ary condition at the sudden expansion point which are important in reacting sim- ulations, are eliminated. The results clearly indicate that the Eulerian/ Lagrangian LES / FMDF methodology is an affordable, consistent and reliable methodology for re— acting flows in complex geometries with any structured, multi-block grid system. The results also indicate that LES / FMDF is quantitatively accurate for three-dimensional, time-dependent flows with large-density variations and strong turbulence-combustion interactions. The third (last) part of this work deals with the SGS modeling and simulation of multi-phase transport and spray combustion. In the next chapter, as a future work, preliminary results via a new two-phase Lagrangian-Eulerian-Lagrangian large eddy simulation methodology are presented and the challenging issues will be addressed. 82 CHAPTER 4 Future Work: LES/FMDF of Multi-Phase Turbulent Reacting Flows 4. 1 Introduction THE PERFORMANCE (efficiency, stability, and compactness) of a combustor in ad- vanced air-breathing propulsion systems is dependent on the coupled and complicated effects of various parameters such as the fuel spray, the fuel / chemistry, the geometry, and the input/output flow conditions. Normally, it is very difficult to predict the combustor behavior under various operating conditions, and it is extremely costly and time consuming to develop/ test a new system via direct (trial and error) experi— mentation. The developed robust and affordable large-eddy simulation methodology can be potentially used for the prediction of multi-phase combustion and unsteady turbu- lent flow evolution in realistic configurations. As mentioned earlier, LES provides detailed time-dependent, spatial data that are very valuable for the development and evaluation of new combustors and the design of reliable control strategies for them. It has been shown that the LES/FMDF is capable of capturing the interactions among turbulence, combustion, and spray and has shown to be applicable to a variety of turbulent flames (slow, fast, non—premixed, premixed, and partially pre- 83 mixed) [10, 64, 65, 96-104]. In this research, the complex interactions among turbu- lence, combustion and spray in realistic combustion systems (e.g., a spray-controlled lean premixed dump combustor) is modeled and simulated via a new two-phase Lagrangian-Eulerian-Lagrangian large-eddy simulation methodology. The spray is modeled with a Lagrangian model [105—108] , which allows two-way mass, momen- tum and energy coupling between phases [109, 110]. The subgrid gas-liquid combus— tion is based on the two-phase filtered mass density function method that has several advantages over conventional methods. The LES/FMDF is employed in conjunc- tion with non-equilibrium and equilibrium reaction models and reduced and detailed chemical kinetics mechanisms. The Eulerian gas-phase part of the solver is based on a generalized high-order multi-block finite-difference method that was used for the computation of compressible turbulent flows in complex geometries in pervious chapters. The proposed LES / FMDF methodology is built on a robust mathemati- cal/ computational framework and can be continuously improved and applied to in- creasingly more complex systems without having to change its basic framework. Ex- pectedly, the LES / FMDF calculations are more intensive than similar RANS compu- tations but are completely feasible, especially for the relatively (geometrically) simple laboratory combustors to be considered in this work. For a systematic assessment of LES / FMDF and its submodels, various single— and multi-phase turbulent reacting and nonreacting flows have been simulated [109—121]. However, the model is not complete and fully accurate in its present form. There are several issues that need to be investigated before we would be able to use the model as a predictive tool for realistic combustion/ propulsion systems; the most important issues are: 0 Chemical kinetics, flamelet approach vs. direct ODE solver, appropriate re- duced kinetics model 84 0 Improved subgrid mixing model, differential diffusion e Reliable and efficient thermal radiation models for droplet and gas phases 0 Group combustion model 0 Improved droplet model that accounts for the internal droplet motion & tem— perature gradient and droplet shape variations 0 Improved wall models for the droplet and gas phases 0 Correct implementation of inflow and outflow conditions for realistic (unsteady) operating conditions 0 Improved spray models for LES, droplet breakup and collision models 0 Improved spray models to capture the interface evolution in secondary breakup and dispersed regions 0 Add a Lagrangian stochastic droplet breakup model to the LES/FMDF flow solver for simulation of dense spray region Presently, the inlet / initial droplet size and velocity distribution needs to be spec- ified a priori from experiment, because the primary breakup spray regime is not simulated. However, a droplet-droplet collision model will be added, similar to that in Refs. [122—124] , to the existing spray model. The collision model has been previously tested for the LES of moderately dense droplet-laden or particle-laden isothermal turbulent flows [124, 125]. I also propose to add a Lagrangian stochastic droplet breakup model to my LES / FMDF flow solver for the simulation of the dense spray region. The model is consistent and easily adaptable to the current Lagrangian droplet flow solver and has been used for LES of cold sprays in diesel engines at var- ious pressures [126]. There are some recent Eulerian-Eulerian LES models in which the effects of the primary jet breakup formation are taken into account [127, 128]. 85 These models might be able to capture the interface evolution in the primary breakup region but are not practical for the secondary breakup and dispersed regions of the spray in a realistic situation. 4.2 Theoretical and Computational Approach My mathematical/ computational approach here is based on a recently—developed, Lagrangian-Eulerian-Lagrangian LES/FMDF methodology. In this methodology, the “resolved” carrier gas velocity, and pressure (Eulerian) fields are obtained by solv- ing the filtered form of the compressible Navier-Stokes equations with ”standard” numerical methods. However, the liquid-droplet phase and scalar (temperature and species mass fractions) fields are both obtained by Lagrangian models. The droplet evaporation and combustion effects are incorporated via a FDF-based methodology, developed for a liquid-gas, reacting system. There is a two-way coupling between the phases and all the Lagrangian and Eulerian fields. Governing equations and numerical procedures for LES / FMDF have been presented for single—phase flows in chapters 2 and 3, and the droplet field is solved via a Lagrangian method for a two-phase system. The LES/FMDF can be implemented via two combustion models: (1) a finite rate, reduced chemistry model for non-equilibrium flames, or (2) a near equilibrium model employing detailed kinetics. In (1), a system of nonlinear ordinary differential equations (ODEs) is solved together with the FMDF equation for all of the scalars (mass fractions and enthalpy). Finite-rate chemistry effects are explicitly and “ex- actly” included in this procedure since the chemistry is closed in the formulation. In (2), the LES/FMDF is employed in conjunction with an equilibrium fuel-oxidation model. This model is enacted via “flamelet” simulations , which consider a laminar counterflow (opposed jet) flame configuration. At low strain rates, the flame is usually close to equilibrium. Thus, the thermo—chemical variables are determined completely 86 by the “mixture fraction.” A flamelet library is coupled with my LES / F MDF solver in which the transport of the mixture fraction is considered. It is useful to emphasize here that the PDF of the mixture fraction is not “assumed” a priori. Rather, it is calculated explicitly via the FMDF. The LES / FMDF / flamelet solver is computation- ally less expensive than that described in (1); thus it can be used for more complex flows. For two-phase flow calculations only non-equilibrium models based on (1) are used. The filtered Eulerian carrier-gas equations are solved together with diffusivity- type closures for the SGS stress and the SGS scalar flux terms via “standard” finite difference and finite volume schemes. For the droplet phase, a stochastic velocity model is considered by which the residual or subgrid velocity of the carrier fluid at the droplet location is constructed. The combined large- and small—scale fluid velocity is then used for calculations of droplet location and velocity. The most convenient means of solving the F MDF transport equation is via the “Lagrangian Monte Carlo” procedure. 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 subgrid). In addition, changes in the “composition space” occur due to chemical reaction, SGS mixing, and droplet evaporation. All of these are implemented via a “stochastic process” described by the set of stochastic differential equations (SDEs)[1]. These SDEs are fully consistent with the original FMDF transport equation. The Lagrangian FMDF represents the gas scalar and energy fields and is used to evaluate the local values of the temperature, density, and species mass fractions at the droplet location. Therefore, no SGS model is needed for the droplet temperature and mass equations. The discretization procedure of the carrier fluid is based on the “compact para- meter” finite difference scheme, which yields up to sixth order spatial accuracy, and 87 a finite-volume collocated grid method. In the latter, the time differencing is based on a third-order low storage Runge—Kutta method. Once the fluid velocity, density and temperature fields are known, the dr0plet transport equations are integrated via the second-order Adams-Bashforth scheme. The evaluation of the fluid velocity at the droplet locations is based upon a second order accurate Lagrangian interpola- tion scheme. The locations and sizes of the droplets vary in different simulations. Currently, the droplet size and velocity at injection location are specified a priori. However, work is in progress to incorporate a stochastic droplet breakup model in my LES solver. 4.3 Preliminary Results and Discussions As mentioned above, the LES / FMDF flow solver is made of several different coupled modules and various Eulerian / Lagrangian numerical schemes. The gas-phase solver is based on a generalized Cartesian (hexahedral) cell applicable to complex configura- tions. It has a multiblock capability and can be used for the computation of industrial combustors with complex geometries. The LES models and the numerical algorithms are tested by comparing the results obtained via this solver for a relatively simple uniform grid, with those generated with my high-order finite-difference and spectral DNS codes. Also, for overall validation of the computational model, the numerical results have been compared with the experimental data. Some recent studies on the numerical simulation of dispersed two-phase turbulent flows in complex geometries also indicate that LES predictions are more reliable than those by RANS [129, 130]. The two-phase reacting LES/FMDF methodology has been applied to a laboratory-scale spray-controlled dump combustor [131]. In this combustor, a lean premixed preheated air-fuel flame is controlled by injection of a relatively small amount of liquid fuel. The experimental setup is shown in figure 4.1. This ex- periment with and without n—decane fuel droplet evaporation and combustion has 88 been simulated, with realistic flow and dr0plet parameters. Figure 4.2 shows the three interacting fields: (1) the Eulerian finite difference field describing the gas dynamic variables; (2) the grid-free Lagrangian Monte Carlo field describing gaseous species and temperature through FMDF; and (3) another La- grangian field describing liquid-fuel droplets and spray. The two-phase LES / F MDF has been used for studying the effects of the fuel injector and other controlling para- meters on the droplet evaporation and combustion process in this combustor. Sample results are presented here in figure 4.3, where the contours of the carrier gas vorticity and droplet & Monet Carlo elements distribution are shown. The results are consistent with the experiment and indicate that the sprays tend to decrease the pressure oscillations in the combustor. Figure 4.4 shows the instantaneous (isolevel) values of the gas temperature, pressure, evaporated fuel, carbon dioxide, water vapor, and oxygen mass fractions in the combustor. I have conducted many simulations and the preliminary results are promising; but this part of work has been considered as “future work”, since it needs much/ careful consideration. The main observation are summarized in the following: The results indicate that close to the combustor inlet the droplet concentration is localized and significant, evidently and expectedly. However, depending on the spray angle and the locations of the injectors, the peak droplet concentration moves closer to the wall and is less localized at the downstream locations. The maximum concentration of the evaporated fuel droplet may or may not coincide with the high temperature and chemically active premixed flame zones. In general, the spray can have negative or positive effects on the reaction in terms of flame temperature and speed. However, close to the inlet the spray usually decreases the reaction rate. This is caused by the cooling effect of the cold droplets and droplet evaporation on the carrier gas field. The spray also has a dampening effect on the pressure field and leads to less pressure oscillations, which is consistent with the experiment. At 89 the downstream locations, the evaporated fuel eventually increases the equivalence ratio of the lean premixed mixture and leads to higher reaction rates and increased temperatures. The effect on the pressure field is complicated and is dependent on the added fuel distribution. In the case studied here the added fuel concentration is not very localized and is not causing significant dilatational fluid motions and pressure oscillations. There is a good correlation between the evaporated fuel and temperature increase caused by the spray, as expected. It was observed that the active droplets have two-way interactions with the carrier gas. But passive droplets have one-way interactions and are being influenced only by gas temperature/velocity field while they have no influence on the carrier gas. Evidently there are noticeable differences between the passive and active droplet cases. The active droplets can attain higher temperatures because of more effective combustion and extra heat release induced by the evaporated fuel. The effect on the droplet velocity seems to be smaller. 4.4 Proposed Future Research Future research can be directed toward the following three foci: (1) Development and validation of a high-fidelity LES-based computational model for numerical simulation of the combustor and related inlet/ outlet components. This may be carried out in three different phases. First, one can consider the unsteady turbulent reacting flows in the combustor. Second, the input / output components (nozzles and diffusers) may be modeled. Third, the integrated simulation of all components might be consid- ered. (2) Detailed simulation of a laboratory combustor can be conducted and the results can be compared with experimental data. The LES is too expensive for the design and optimization of the combustor or prOpulsion unit but complements the experiment and may be used to test ideas in a time / cost effective manner and to sim- ulate conditions that are not possible to generate in the laboratory such as unsteady, 90 high/ low pressure, and highly turbulent inputs. (3) Detailed integrated simulation of the entire propulsion (ramjet) unit is the long-time goal of the proposed research. The vision is to develop a virtual prototype that replicates a simple ramjet engine on the computer. This work deals with the development of a reliable and affordable LES model for turbulent reacting flows with and without multi-phase flow effects such as dr0plet evaporation. Within the past few years, I have extended, implemented, and validated a novel subgrid-scale combustion model (termed the filtered mass density function) for numerical simulations of turbulent combustion, which was discussed in the pre- vious chapters. The FMDF has been successfully applied to various non-premixed, premixed, and partially-premixed turbulent flames and was extended to two—phase reacting systems. The two-phase formulation of LES/FMDF involves the two-way mass / momentum / energy coupling between phases and droplet transport, evaporation and combustion. The LES / FMDF model is implemented via an efficient Lagrangian- Eulerian-Lagrangian numerical scheme. The efforts need to be continued in the following directions: 0 Development and application of LES to gaseous combustion in various configu- rations with reduced and detailed kinetics models 0 Development and application of LES to flows involving sprays with: — droplet dispersion — droplet breakup and coalescence — droplet evaporation and combustion 0 Comparison of the LES results with experimental and DNS data whenever pos- sible 91 Furthermore, direct numerical simulations of “simple” turbulent flames may also be considered with the main objectives being the understanding of the underlying processes and assessing the subgrid LES models. Further validation can be done by comparing the LES results with the experimental data, whenever possible. In support of this research, the intend is to conduct fundamental investigations for better understanding of flow physics and the development of improved closures for some of the subcomponents of the LES / FMDF. This will be partially based on direct numerical simulation of reacting turbulent flows with a particular emphasis on flow complexities (such as realistic chemical kinetics, flame-wall interactions, and complex geometries etc.). The execution of such a plan is well-justified if the final goal of practical use is kept in mind. Currently, there are very limited data available for turbulent combustion. Development of an extensive DNS data base is one of the objectives of my proposed plan. Validation of the simulated results via comparison with experimental data will also be an essential part this future research. 92 4.5 ‘Figures l liquid magi-ID H2 C2H4 tuel (b) Figure 4.1. Experimental setup. 93 Figure 4.2. Liquid Riel Droplets (big solid points), Monte Carlo particles (small solid points), and Eulerian Grid layout. 94 Figure 4.3. Contours of the instantaneous vorticity magnitude and droplet and Monet Carlo elements distribution (a) 3D and (b) 2D. 95 (0) 0101722 and 002 (d) CioH22. 02, H20, and 002 Figure 4.4. Contours of the instantaneous gas temperature, pressure, and gas species mass fractions as obtained by two-phase LES/FMDF in a spray-controlled premixed reacting dump combustor. 96 [ll [2] [4] [5] [7] [8] l9] BIBLIOGRAPHY Lele, S. K., “Compact Finite Difference Schemes With Spectral-Like Resolu- tion,” Jornal of Computational Physics, Vol. 103, 1992, pp. 16—42. Visbal, M. R. and Gaitonde, D. V., “High-Order Accurate Methods for Complex Unsteady Subsonic Flows,” AIAA Journal, Vol. 37, No. 10, 1999, pp. 1231— 1239. Visbal, M. R. and Gaitonde, D. V., “Very High-Order Spatially Implicit Schemes For Computational Acoustics on Curvilinear Meshes,” Journal of Com- putational Acoustics, Vol. 9, No. 4, 2002, pp. 1259—1286. Gaitonde, D. V. and Visbal, M. R., “A further Development of a Navier-Stokes Solution Procedure Based on High-Order Formulas Simulation,” AIAA paper 99-0557, Jan. 1999. Visbal, M. R. and Rizzetta, D. P., “Large-Eddy Simulations on Curvilinear Grids Using Compact differencing and Filtering Schemes,” Journal of Fluid Engineering, Vol. 124, 2002. Visbal, M. R. and Gaitonde, D. V., “On the Use of High-Order F mite-Difference Schemes on Curvilinear and Deforming Meshes,” Journal of Computational Physics, Vol. 181, 2002. Jameson, A. and Baker, T., “Solutions of the Euler Equations for Complex Configurations,” AIAA paper 83—1929, 1983. Jixian, Y., Jameson, A., and Juan, J. A., “Development and Validation of a Massively Parallel Flow Solver for Turbomachinery Flows,” Journal of Propul- sion and Power, Vol. 17, No. 3, 2001, pp. 659—667. Colucci, P. J ., Jaberi, F. A., Givi, P., and Pope, S. B., “Filtered Density Emc- tion For Large Eddy Simulation of Turbulent Reacting Flows,” Phys. Fluids, Vol. 10, No. 2, 1998, pp. 499—515. 97 [10] Jaberi, F., Colucci, R, James, S., Givi, P., and Pope, S., “Filtered Mass Density Function for Large Eddy Simulation of Turbulent Reacting Flows,” Journal of Fluid Mechanics, Vol. 401, 1999. [11] Gould, R. D. and Stevenson, W., “Investigation of Turbulent in an Axisymmet- ric Sudden expansion,” AIAA Journal, Vol. 28, No. 2, 1990, pp. 276—283. [12] Gould, R. D., Stevenson, W., and Doyle, H., “Simultaneous Velocity and tem- perature Measurements in a Premixed Dump Combustor,” AIAA Journal, Vol. 10, No. 5, 1994, pp. 639—645. [13] Dellenback, P. A., Metzger, D. E., and Neitzel, G. P., “Measurment in Turbulent Swiriling Flow Through an Abrupt Axisymmetric Expansion,” AIAA Jornal, Vol. 26, No. 6, 1988, pp. 669—681. [14] Stieglmeier, M., Tropea, G, Weiser, N., and Nitsche, W., “Experimental In- vestigation of the Flow Through Axisymmetric Expansions,” Journal of Fluid Engineering, Vol. 111, 1989, pp. 464—471. [15] Ahmed, S. A. and Nejad, A. S., “Premixed, Turbulent Combustion of Axisym- metric Sudden Expansion Flows,” International Journal for Heat and Fluid Flow, Vol. 13, No. 1, 1992, pp. 15—21. [16] Cole, D. R. and Glauser, M. N., “Flying Hot Wire Measurements in an Ax- isymmetric Sudden Expansion,” Journal of Experimental Thermal and Fluid Science, Vol. 18, 1998, pp. 150—167. [17] Lieuwen, T. and Zinn, B. T., “On the Experimental Determination of Com- bustion Process Driving in an Unstable Combustor,” Combustion and Science Technology, Vol. 157, 2000, pp. 111—127. [18] Sommerfeld, M. and Qiu, H.-H., “Experimental Studeies of Spray Evaporation in Turbulent Flow,” International Journal for Heat and Fluid Flow, Vol. 19, 1998, pp. 10—22. [19] Sommerfeld, M. and Qiu, H.-H., “Detailed Measurments in a Swirling Particu- late Two—Phase Flow by a Phase-Doppler Anemometer,” International Journal for Heat and Fluid Flow, Vol. 11, No. 1, 1991, pp. 20—28. [20] Chuang, S.-H., Lin, H.-C., Tai, F.-M., and Sung, H.-M., “Hot Flow Analysis of Sweirling Sudden Expansion Dump Combustor,” International Journal for Numerical Methods in Fluids, Vol. 14, 1992, pp. 217—239. [21] Guo, B., Langrish, T. A. G., and Fletcher, D. F., “Simulation of Turbulent Swirl Flow in an Axisymmetric Sudden Expansion,” AIAA Journal, Vol. 39, No. 1, 2001, pp. 96—102. 98 [22] Zhang, H. Q., Chan, C. K., and Lau, K. S., “Numerical Simulation of Sudden Expansion Particle Laden Flows Using an Improved Stochastic Separated Flow Model,” AIAA Journal, Vol. 40, 2001, pp. 89—102. [23] Chuang, S.-H., Yang, C.-H., and Wu, N .-J ., “Predictions of Swirling Flow in Sudden Expansion Dump Combustor with Flameholder Side Inlet using Two Step Combustion Model,” International Journal of Numerical Methods for Heat and Fluid Flow, Vol. 9, No. 7, 1999, pp. 764—787. [24] Stone, C. and Menon, 8., “Large Eddy Simulations of Isotropic Decaying With Structured and Unstructured Grid Finite Volume Method,” High performance computing, grand challenges in computer simulations, seattle, WA, 2001. [25] Schliiter, J. U., Pitsch, H., and Moin, P., “Consistent Boundary Conditions for Integrated LES/RANS Simulations: LES Outflow Conditions,” AIAA paper 02—3121, 2002. [26] Huang, Y., Sung, H.—G., Hsieh, H.-Y., and Yang, V., “Large Eddy Simulation of Combustion Dynamics of Lean Premixed Swirl Stabilized Combustor,” AIAA Journal, Vol. 19, No. 5, 2003, pp. 782—794. E [27] Moin, P. and Apte, V., “Large-Eddy Simulation of Realistic Gas Turbine Com- bustors,” AIAA Journal, Vol. 44, No. 4, 2006, pp. 698—708. [28] Pulliam, T. H. and Steger, J. L., “Implicit finite-difference simulation of three- dimensional compressible flow,” AIAA Journal, Vol. 18, No. 2, 1980, pp. 159. [29] Vinokur, M., “Conservation equations of gasdynamics in curvilinear coordinate systems,” Jornal of Computational Physics, Vol. 14, 1974, pp. 105. [30] Steger, J. L., “Implicit finite-difference simulation of flow about arbitrary two— dimensional geometries,,” AIAA Journal, Vol. 16, No. 7, 1978, pp. 679. [31] Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer, McGrawHill, New York, 1984. [32] Aldama, A. A., Filtering Techniques for Turbulent Flow Simulations, Vol. 49 of Lecture Notes in Engineering, Springer-Verlag, New York, NY, 1990. [33] Dailey, L. D., Simons, T. A., and Pletcher, R. H., “Large Eddy Simulations of Isotropic Decaying With Structured and Unstructured Grid Finite Volume Method,” Being reviewed for presentation at the ASME winter annual meeting, Nov. 1996. 99 [34] VanderVan, H. S. and Stone, C., “A Family of Large Eddy Simulation (LES) filters with Nonuniform Filter widths,” Physics of Fluids Journal, Vol. 7, No. 5, 1995, pp. 1171—1172. [35] Won-Wook, K. and Suresh, M., “LES of Turbulent Fuel/Air Mixing in a Swirling Combustor,” AIAA paper 98—0200, Jan. 1999. [36] Suresh, M. and Stone, 0, “Simulation of a Fuel-Air Mixing Combustion in a Trapped-vortex Comubstor,” AIAA paper 2000—0478, Jan. 2000. [37] Yakhot, A., Orszag, S. A., Yakhot, V., and Israeli, M., “Renormalization Group Formulation of Large-Eddy Simulation.” Journal of Scientific Comput- ing, Vol. 1, 1986, pp. 1—51. [38] Visbal, M. R. and Gaitonde, D. V., “Pade Type High-Order Boundary Fil- ters for the Navier-Sotkes Equations,” AIAA Journal, Vol. 38, No. 11, 2000, pp. 2103—2112. [39] Shu, C. W. and Osher, S., “Efficient Implementation of Essentially Non- Oscillatory Shock-Capturing schemes,” Journal of Computational Physics, Vol. 77, 1988. [40] Gottlieb, S., Shu, C. W., and Tadmor, E., “Strong Stability-Preserving High- Order time discretization Methods,” SIAM Review, Vol. 43, No. 1, 2001, pp. 89— 112. [41] Poinsot, T. and Lele, S., “Boundary Conditions for direct Numerical Simula- tions of Compressible Viscous Flows,” J. Comput. Phys, Vol. 101, 1992. [42] Afshari, A., Shotorban, B., Mashayek, F., Shih, T.-P., and Jaberi, F. A., “Devel- opment and Validation of a Multi-Block Flow Solver for Large Eddy Simulation of Turbulent Flows in Complex Geometries,” AIAA paper 2004-657, January 2004. [43] Gillandt, I., Fritsching, U., and Bauckhage, K., “Measurement of Phase Inter- action in Dispersed Gas/ Particle Two-Phase Flow,” Int. Journal of Multiphase flow, Vol. 27, 2001. [44] Shih, T.-P., Bailey, R., Nguyen, H., and R.J.Roelke, “Algebraic Grid Genera- tion for Complex Geometries,” International Journal for Numerical Methods in Fluids, Vol. 13, 1991, pp. 1—31. [45] Steinthorsson, E., Shih, T.-P., and R.J.Roelke, “Enhancing Control of Grid Dis- tribution in Algebraic Grid Generation,” International Journal for Numerical Methods in Fluids, Vol. 15, 1992. 100 [46] Rudy, D. H. and Strikwerda, J. C., “Boundary Conditions For Subsonic Com- pressible Navier-Stokes Calculations,” Computers and Fluids, Vol. 9, 1981. [47] McMurtry, P. A., Menon, S., and Kerstein, A. R., “A Linear Eddy Sub-Grid Model for Turbulent Reacting Flows: Application to Hydrogen-Air Combus- tion,” Proceedings of 24th Symp. (Int. ) on Combustion, The Combustion Insti- tute, Pittsburgh, PA, 1992, pp. 271—278. [48] Menon, S., McMurtry, P. A., and Kerstein, A. K., “A Linear Eddy Sugbrid Model of Turbulent Combustion,” Large Eddy Simulations of Complex En- gineering and Geophysical Flows, edited by B. Galperin and S. A. Orszag, chap. 14, Cambridge University Press, Cambridge, UK, 1993, pp. 287—314. [49] Gao, F. and O’Brien, E. E., “A Large Eddy Scheme for Turbulent Reacting Flows,” Phys. Fluids A, Vol. 5, No. 6, 1993, pp. 1282—1284. [50] Madnia, C. K. and Givi, P., “Direct Numerical Simulation and Large Edddy Simulation of Reacting Homogeneous Turbulence,” Large Eddy Simulations of Complex Engineering and Geophysical Flows, edited by B. Galperin and S. A. Orszag, chap. 15, Cambridge University Press, Cambridge, UK, 1993, pp. 315—346. [51] Frankel, S. H., Adumitroaie, V., Madnia, C. K., and Givi, P., “Large Eddy Sim- ulations of Turbulent Reacting Flows by Assumed PDF Methods,” Engineering Applications of Large Eddy Simulations, edited by S. A. Ragab and U. Piomelli, ASME, FED-Vol. 162, New York, NY, 1993, pp. 81—101. [52] Cook, A. W. and Riley, J. J ., “A Subgrid Model for Equilibrium Chemistry in Turbulent Flows,” Phys. Fluids, Vol. 6, No. 8, 1994, pp. 2868—2870. [53] Fureby, C. and Lofstrom, C., “Large-Eddy Simulations of Bluff Body Stabilized Flames,” Proceedings of 25th Symp. (Int. ) on Combustion, The Combustion Institute, Pittsburgh, PA, 1994, pp. 1257-1264. [54] Branley, N. and Jones, W. P., “Large Eddy Simulation of a Turbulent Non- Premixed Flame,” Proceedings of the Eleventh Symposium on Turbulent Shear Flows, Grenoble, France, 1997, pp. 21.1—21.6. [55] Cook, A. W., Riley, J. J., and Kosaly, G., “A Laminar Flamelet Approach to Subgrid-Scale Chemistry in Turbulent Flows,” Combust. Flame, Vol. 109, No. 3, 1997, pp. 332—341. [56] Mathey, F. and Chollet, J. P., “Large Eddy Simulation of Turbulent Reac- tive Flows,” Proceedings of the Eleventh Symposium on Turbulent Shear Flows, Grenoble, France, 1997, pp. 16.19—16.24. 101 [57] DesJardin, P. E. and Frankel, S. H., “Large Eddy Simulation of a Turbulent Nonpremixed Reacting Jet: Application and Assessment of Subgrid-Scale Com- bustion Models,” Phys. Fluids, Vol. 10, No. 9, 1998, pp. 2298—2314. [58] Jaberi, F. A. and James, S., “A Dynamic Similarity Model for Large Eddy Sim- ulation of Turbulent Combustion,” Phys. Fluids, Vol. 10, No. 7, 1998, pp. 1775— 1777. [59] Réveillon, J. and Vervisch, L., “Subgrid-Scale Turbulent Micromixing: Dynamic Approach,” AIAA J., Vol. 36, No. 3, 1998, pp. 336—341. [60] Pope, S. B., “Computations of Turbulent Combustion: Progress and Chal- lenges,” Proceedings of 23rd Symp. (Int. ) on Combustion, The Combustion Institute, Pittsburgh, PA, 1990, pp. 591—612. [61] Suramaniam, S. and Haworth, D. C., “A Pdf Method for Turbulens Mixing and Combustion on Three-Dimensional Unstructured Deforming Meshes,” Submit- ted to Inernational Journal of Engine Research, 24 March 2004. [62] Zhang, Y. Z. and Haworth, D. C., “A General Mass Consistency Algorithm for Hybrid Particle/Finitevolume PDF Methods,” Jornal of Computational Physics, Vol. 194, 2004, pp. 15—193. [63] 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., Vol. 401, 1999, pp. 85—121. [64] Mehravaran, K., Mathematical Modeling and Numerical Simulations of Chem- ically Reacting Turbulent Jets, Ph.D. thesis, Michigan State University, June 2005. [65] Afshari, A., Almeida, T., Mehravaran, K., and Jaberi, F. A., “Large Scale Simulations of Turbulent Combustion and Propulsion Systems,” Proceedings of the seventeen ONR propulsion meeting, June 2004. [66] Givi, P., “Sugrid Scale Modeling in Turbulent Combustion: A Review,” AIAA paper 2003—5081, July 2003. [67] Givi, P., “Filtered Density Function for Subgrid Scale Modeling of Turbulent Combustion,” AIAA Journal, Vol. 44, No. 1, 2006, pp. 16—23. [68] YU, H., LUO, L.-S., and GIRIMAJI, S. S., “Scalar Mixing and Chemical Reac- tion Simulations Using Lattice Boltzmann Method,” International Journal of Computational Engineering Science, Vol. 3, No. 1, 2002, pp. 73—78. 102 [69] O’Brien, E. E., “The Probability Density Function (PDF) Approach to Reacting Turbulent Flows,” Turbulent Reacting Flows, edited by P. A. Libby and F. A. Williams, chap. 5, Springer-Verlag, Heidelberg, 1980, pp. 185—218. [70] Pope, S. B., “PDF Methods for Turbulent Reacting Flows,” Prog. Energy Com- bust. Sci, Vol. 11, 1985, pp. 119—192. [71] Lundgren, T. S., “Model Equation for Nonhomogeneous Turbulence,” Phys. Fluids, Vol. 12, No. 3, 1969, pp. 485—497. [72] Pope, S. B., “The Probability Approach to Modeling of Turbulent Reacting Flows,” Combust. Flame, Vol. 27, 1976, pp. 299—312. [73] Germano, M., “Turbulence: The Filtering Approach,” J. Fluid Mech., Vol. 238, 1992, pp. 325—336. [74] Salvetti, M. V. and Banerjee, S., “A priori Tests of a New Dynamic Subgrid- Scale Model for Finite-Difference Large-Eddy Simulations,” Phys. Fluids, Vol. 7, No. 11, 1995, pp. 2831—2847. [75] Dopazo, C. and O’Brien, E. E., “Statistical Treatment of Non-Isothermal Chem- ical Reactions in Turbulence,” Combust. Sci. and Tech, Vol. 13, 1976, pp. 99— 112. [76] Borghi, R., “Turbulent Combustion Modeling,” Prog. Energy Combust. Sci, Vol. 14, 1988, pp. 245—292. [77] Pope, S. 3, “Particle Method for Turbulent Flows: Integration of Stochastic Model Equations,” J. Comp. Phys, Vol. 117, 1995, pp. 332—349. [78] 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, Vol. 154, No. 2, Sep 1999, pp. 342—371. [79] Muradoglu, M., Pope, S. B., and Caughey, D. A., “The Hybrid Method for the PDF Equations of Turbulent Reactive Flows: Consistency Conditions and Correction Algorithms,” J. Comp. Phys, Vol. 172, 2001, pp. 841—878. [80] Gicquel, L. Y. M., Givi, P., Jaberi, F. A., and Pope, S. B., “Velocity Filtered Density Function for Large Eddy Simulation of 'Ilirbulent Flows,” Phys. Fluids, Vol. 14, No. 3, 2002, pp. 1196—1213. [81] 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, 2003, in press. 103 [82] Risken, H., The F okker-Planck Equation, Methods of Solution and Applications, Springer-Verlag, New York, NY, 1989. [83] Gardiner, C. W., Handbook of Stochastic Methods, Springer-Verlag, New York, NY, 1990. [84] Karlin, S. and Taylor, H. M., A Second Course in Stochastic Processes, Acad- [85] [86] [87] [88] [91] [92] [93] [94] [95] emic Press, New York, NY, 1981. Pope, S. B., “Lagrangian PDF Methods for Turbulent FLows,” Ann. Rev. Fluid Mech., Vol. 26, 1994, pp. 23—63. Grigoriu, M., Applied Non-Gaussian Processes, Prentice-Hall, Englewood Cliffs, NJ, 1995. Kloeden, P. E., Platen, E., and Schurz, H., Numerical Solution of Stochastic Difi‘erential Equations through Computer Experiments, Springer-Verlag, New York, NY, corrected second printing ed., 1997. Haworth, D. C. and POpe, S. B., “A Second-Order Monte Carlo Method for the Solution of the Ito Stochastic Differential Equation,” Stochastic Analysis and Applications, Vol. 4, No. 2, 1986, pp. 151—186. Haworth, D. C. and Pope, S. 8., “Monte Carlo Solutions of a Joint PDF Equa- tion for Turbulent Flows in General Orthogonal Coordinates,” J. Comp. Phys, Vol. 72, No. 2, 1987, pp. 311—346. Xu, J. and Pope, S. 8., “Assessment of Numerical Accuracy of PDF/Monte Carlo Methods for Turbulent Reacting Flows,” J. Comp. Phys, Vol. 152, 1999, pp. 192—230. Kloeden, P. E. and Platen, E., Numerical Solution of Stochastic Differential Equations, Vol. 23 of Applications of Mathematics, Stochastic Modelling and Applied Probability, Springer-Verlag, New York, NY, 1995. Billingsly, P., Probability and Measure, John Wiley and Sons, Inc., New York, 1979. Gillespie, D. T., Markov Processes, An Introduction for Physical Scientists, Academic Press, New York, NY, 1992. Ito, K., “On Stochastic differential equations,” Memoirs. Amer. Math. Soc., Vol. 4, 1951. Gikhman, I. I. and Skorokhod, A. V., Stochastic Difierential Equations, Springer-Verlag, New York, NY, 1995. 104 [82] Risken, H., The Fokker-Planck Equation, Methods of Solution and Applications, [83] [84] [85] [86] [87] [88] [89] [90] [91] [92] [93] [94] [95] Springer-Verlag, New York, NY, 1989. Gardiner, C. W., Handbook of Stochastic Methods, Springer-Verlag, New York, NY, 1990. Karlin, S. and Taylor, H. M., A Second Course in Stochastic Processes, Acad- emic Press, New York, NY, 1981. Pope, S. B., “Lagrangian PDF Methods for Turbulent FLows,” Ann. Rev. Fluid Mech., Vol. 26, 1994, pp. 23—63. Grigoriu, M., Applied Non- Gaussian Processes, Prentice—Hall, Englewood Cliffs, NJ, 1995. Kloeden, P. E., Platen, E., and Schurz, H., Numerical Solution of Stochastic Difierential Equations through Computer Ezperiments, Springer-Verlag, New York, NY, corrected second printing ed., 1997. Haworth, D. C. and POpe, S. B., “A Second-Order Monte Carlo Method for the Solution of the Ito Stochastic Differential Equation,” Stochastic Analysis and Applications, Vol. 4, No. 2, 1986, pp. 151—186. Haworth, D. C. and Pope, S. B., “Monte Carlo Solutions of a Joint PDF Equa- tion for Turbulent Flows in General Orthogonal Coordinates,” J. Comp. Phys, Vol. 72, No. 2, 1987, pp. 311—346. Xu, J. and Pope, S. B., “Assessment of Numerical Accuracy of PDF/Monte Carlo Methods for Turbulent Reacting Flows,” J. Comp. Phys, Vol. 152, 1999, pp. 192—230. Kloeden, P. E. and Platen, E., Numerical Solution of Stochastic Differential Equations, Vol. 23 of Applications of Mathematics, Stochastic Modelling and Applied Probability, Springer-Verlag, New York, NY, 1995. Billingsly, P., Probability and Measure, John Wiley and Sons, Inc., New York, 1979. Gillespie, D. T., Markov Processes, An Introduction for Physical Scientists, Academic Press, New York, NY, 1992. Ito, K., “On Stochastic differential equations,” Memoirs. Amer. Math. Soc., Vol. 4, 1951. Gikhman, I. I. and Skorokhod, A. V., Stochastic Difierential Equations, Springer-Verlag, New York, NY, 1995. 104 [96] [97] [100] [101] [102] [103] [104] [105] [106] James, S. and Jaberi, F. A., “Large Scale Simulations of Two—Dimensional Nonpremixed Methane Jet Flames,” Combust. Flame, Vol. 123, 2000, pp. 465— 487. Sheikhi, M. R. H., Drozda, T. G., Givi, P., Jaberi, F. A., and Pope, S. H, “Large Eddy Simulation of a Turbulent Nonpremixed Piloted Methane Jet Flame (Sandia Flame D),” Proceedings of the Combustion Institute, Vol. 30, 2005, pp. 549—556. J aberi, F. A., “Large Eddy Simulation of Turbulent Premixed Flames via Fil- tered Mass Density Function,” AIAA paper 99-0199, January 1999. Garrick, S. C., Jaberi, F. A., and Givi, P., “Large Eddy Simulation of Scalar Transport in a Turbulent Jet Flow,” Recent Advances in DNS and LES, edited by D. Knight and L. Sakell, Vol. 54 of Fluid Mechanics and its Applications, Kluwer Academic Publishers, The Netherlands, 1999, pp. 155—166. Zhou, X. Y. and Pereira, J. C. F., “Large Eddy Simulation (2D) of a Reacting Plan Mixing Layer Using Filtered Density Function,” Flow, Turbulence and Combustion, Vol. 64, 2000, pp. 279—300. Venkatramanan, R., Cook, D., and Pitsch, H., “Hybrid LES/FDF Simulation of a Non-Premixed Bluff-Body Stabilized Flame,” Bulletin of the American Physical Society, Vol. 49, No. 9, 2004, pp. 57—58. Glaze, D. J ., Frankel, S. H., and Hewson, J. C., “Non-Premixed Turbulent Jet Mixing Using LES with the FMDF Model,” 2004, pp. 79, in 30th International Symposium on Combustion, Abstracts of Work-In-Progress Posters, Pittsburgh, PA. Raman, V., Pitsch, H., and Fox, R. O., “Consistent Hybrid LES-FDF For- mulation for the Simulation of Turbulent Combustion,” 2004, pp. 231—241, in Annual Research Briefs, Center for Turbulence Research, NASA Ames Research Center, Stanford University, CA. Lu, L., Ren, Z., Raman, V., Pope, S. B., and Pitsch, H., “LES/FDF/ISAT Computations of Turbulent Flames,” 2004, pp. 283—294, in Proceedings of the 2004 Summer Program, Center for Turbulence Research, NASA Ames Research Center, Stanford University, CA. Faeth, G. M., “Mixing, Transport and Combustion in Sprays,” Prog. Energy Combust. Sci., Vol. 13, 1987, pp. 293-345. Sirignano, W. A., Fluid Dynamics and Transport of Droplets and Sprays, Cam- bridge University Press, 1999. 105 [107] Crowe, C., Sommerfeld, M., and Tsuji, Y., Multiphase Flows Droplets and Par- ticles, CRC Press, 1998. [108] Loth, E., “Numerical approaches for motion of dispersed particles, droplets and bubbles,” Prog. Energy Combustion Sci., Vol. 26, 2000, pp. 161—223. [109] Mashayek, F. and Pandya, R. V. R., “Analytical Description Particle/Droplet- laden Turbulent Flows,” Prog. Energy Combustion Sci, Vol. 29, 2003, pp. 329— 378. [110] Bellan, J., “Perspectives on large eddy simulations for sprays: Issues and solu- tions,” Atomization and Sprays, Vol. 10, 2000, pp. 409—425. [111] Multiphase Flows Droplets and Particles, R. T. Edwards, Inc., 2001, J. B. Geurts, editor. [112] Pope, S. B., Turbulent Flows, Cambridge University Press, Cambridge, UK, 2000. [113] Ciofalo, M., Large Eddy Simulation: A Critical Survey of Models and Applica- tions, Vol. 25, Academic Press, New York, NY, 1994. [114] Oran, E. S. and Boris, J. P., Numerical Simulation of Reactive Flows, Cam- bridge University Press, New York, NY, 2001, second edition. [115] Givi, P., “Subgrid Scale Modeling in Turbulent Combustion: A Review,” AIAA Paper 2003-5081, 2003. [116] Jaberi, F. A., Mashayek, F., Madnia, C. K., Taulbee, D. B., and Givi, P., Advances in Chemical Propulsion: Science to Technology—chapter 9: Analytical Description of Turbulent Reacting Flows, pp. 149-164, CRC Press LLC, Boca Raton, FL, 2002, Roy, G. D., editor. [117] Fox, R. 0., Computational Models for Turbulent Reacting Flows, Cambridge University Press, Cambridge, UK, 2003. [118] Peters, N., Turbulent Combustion, Cambridge University Press, Cambridge, UK, 2003. [119] Veynante, T. P. D., Theoretical and Numerical Combustion, R. T. Edwards, Inc., Philadelphia, PA, 2001. [120] Menon, S., “Subgrid Combustion Modeling for Large-Eddy Simulations,” Iner- national Journal of Engine Research, Vol. 1, No. 2, 2000, pp. 209—227. 106 [121] Jaberi, F. A., Madnia, F. C. K., and Givi, P., Handbook of Numerical Heat Transfer—chapter 5: Large Eddy Simulation of Heat and Mass Ti‘ansport in Turbulent Flows, John Wiley and Sons, 2005, in press, Edited by Minkowycz, W. J., Sparrow, E. M. and Murthy, J. Y. [122] Schmidt, D. P. and Rutland, C. J ., “A new Droplet Collision Algorithm,” Jornal of Computational Physics, Vol. 164, 2000, pp. 62—80. [123] Post, S. L. and J .Abraham, “Modeling the Outcome of Drop-Drop Collisions in Diesel Sprays,” Int. Journal of Multiphase flow, Vol. 28, 2002. [124] Gorokhovski, M. and Chtab, A., “Hypothetical Heavy Particles Dynamics in LES of Turbulent Dispersed Two-Phase Channel Flow,” 2003, Annual research briefs, Center for Turbulence Research, Stanford University, CA. [125] Apte, S. V., Mahesh, K., and Lundgren, T., “A Eulerian-Lagrangian Model to Simulate Two-Phase/ Particulate Flows,” 2003, Annual research briefs, Center for 'Ihrbulence Research, Stanford University, CA. [126] Apte, S. V., Gorokhovski, M., and Moin, P., “LES of atomizing spray with stochastic modeling of secondary breakup,” Int. Journal of Multiphase flow, Vol. 29, 2003. [127] Shen, L. and Yue, D., “Large—Eddy Simulation of Free-Surface Turbulence,” Journal of Fluid Mechanics, Vol. 440, 2001. [128] Herrmann, M., “Modeling Primary Breakup: Three-Dimensional Eulerian Level Set / Vortex Sheet Method for Two-Phase Interface Dynamics,” 2003, Technical report, Center for Turbulence Research nnual Briefs, Stanford University, CA. [129] Apte, S. V., Mahesh, K., Moin, P., and Oefelein, J. C., “Large-eddy simulation of Swirling Particle-Laden Flows in a Coaxial-Jet combustor,” Int. Journal of Multiphase flow, Vol. 29, 2003. [130] Lakehal, D., “On the modeling of multiphase turbulent flows for environmental and hydrodynamic applications,” Int. Journal of Multiphase flow, Vol. 28, 2002. [131] Hsu O., Pang, B. and Yu, K., “Active Control Studies for Advanced Propulsion: Temperature Dependence of Critical Fuel Flux,” Proceedings of the Fourteenth ONR propulsion Meeting, edited by G. Roy and F. Mashayek, 2001, pp. 19—24, Chicago, IL. 107 111111111111