v :1) . :1 r? n. .2... 6.5% I ,rhvaw . k} u i. 31.3%, r i?.m$$..¥¥3 5.71.3 flillj V ,. at. hr: n .. awn.» I II iv‘ ’1‘? .v ..,._\<....- . C... «3.5. . 1., i an .: )1... .IO! v $ .v .. . . . . (51:21.... ‘ ti 7.23: 15...: .:Li\.i . .Itxaia .a. ‘ '- {.3511 .9. .qrzta A. .. : 53?.33‘11 .13.... 5.4 b L :2 . h.w....b..ztt4 .1}... . 1.5;: s: a??? .3; \‘)‘I . $3 iflhJflXS 70 .H. ‘ . , A . . a. fi...§$m§mwaw ”g: Jam, .. l masts LIBRARY ," Michigan State ‘ T University This is to certify that the dissertation entitled MODELING AND SIMULATION OF TURBULENT MULTIPHASE FLOWS presented by Zhaorui Li has been accepted towards fulfillment of the requirements for the Doctoral degree in Mechanical Engineering «OK/Z MaWrofe? or 5 Signature ogéé'fl lg, :24“: Date MSU is an Affirmative Action/Equal Opportunity Employer 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 5/08 K:/Prolecc&Pres/ClRC/DateDueiindd MODELING AND SIMULATION OF TURBULENT MULTIPHASE FLOWS By Zhaorui Li A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Mechanical Engineering 2008 ABSTRACT MODELING AND SIMULATION OF TURBULENT MULTIPHASE FLOWS By Zhaorui Li The atomization of liquid spray in turbulent reacting and non-reacting flows usually occurs in two successive steps, i.e., primary breakup and secondary breakup. In the primary breakup region, the evolution of the interface between the phases is usually complex and very difficult to model. In the secondary breakup region, the average droplet size and volume occupied by the droplets are relatively small but the number of droplets is usually very significant. In this study, we use different mathematical and numerical models for different regions of the spray. For dense spray simulations, a coupled Lagrangian interface- tracking and Eulerian level set method is developed and implemented. In this method, the interface is identified based on the locations of notional particles and the geometrical information concerning the interface and fluid properties are obtained from the level set function. The level set function maintains a signed distance function via the particle- based Lagrangian re-initialization technique. Numerical simulations of several ‘standard interface-moving’ problems and two-fluid laminar and turbulent flows are conducted to assess this new hybrid method. The results of our analysis indicate that the hybrid particle-level set method can handle interfaces with complex shape change, and can accurately predict the interface values without any significant mass loss or gain. The results obtained for isotropic two-fluid turbulence via the new particle-level set method are validated by comparison with those obtained by the ‘zero Mach number’, variable- density method. The two-way interactions between the turbulent velocity field and the interface are studied by the particle-level set method. Extensive analysis of vorticity and energy equations indicates that the destabilization effect of turbulence and stability effect of surface tension on the interface motion and interface’s effect on turbulence are strongly dependent on the density ratio and Weber number. For dilute spray simulations, a robust and efficient Eulerian-Lagrangian- Lagrangian mathematical/numerical LES model is employed. This is based on the filtered mass density function (FMDF) methodology and is applicable to two-phase turbulent reacting flows with two-way mass, momentum and energy coupling between phases included. In the LES/FMDF methodology, the “resolved” carrier gas velocity field is obtained by solving the filtered form of the compressible Navier-Stokes equations with a high-order finite difference scheme. The sugrid species, energy and combustion are modeled with the two-phase scalar FMDF transport equation, which is solved by a Lagrangian Monte Carlo method. The liquid droplet/spray is simulated with a non- equilibrium Lagrangian model and stochastic SGS closures. The two-way coupling is implemented through series of source/sink terms. The two-phase LES/FMDF is employed for systematic analysis of turbulent combustion in the double swirl spray burner and spray-controlled dump combustor for various flows and spray parameters. The effects of fuel type, spray/injection angle, mass loading ratio, droplet size and its distribution, fuel/air composition, wall, and other parameters on the combustion and turbulence are investigated. ACKNOWLEDGMENTS Completing this work would not have been possible without the support of and inspiration from a large number of people. First of all I would like to thank my advisor, Professor Farhad Jaberi for enabling me to carry out this research. I gratefully acknowledge his advice and guidance, continuous encouragement, and patient support throughout this work. It has been a very enlightening and fruitful experience to work with Professor Jaberi. I am also grateful to members of my PhD dissertation committee, Professor Tom Shih, Professor Guowei Wei, Professor Mei Zhuang and Professor Andre Benard. I deeply appreciate their time and efforts and their valuable comments and suggestions. I also wish to thank Professor Ashwani Gupta at the University of Maryland for valuable discussions and suggestions. I would like to acknowledge the financial support from the US National Science Foundation under grant NSF-0325760 & CTS-OO92665 and the US Office of Naval Research under grant NOOOl4-01-l-084 (Program Manager: Dr. Gabriel Roy). Computational resources are provided by the High Performance Computing Center (HPCC) at the Michigan State University, and the National Center for Supercomputing Applications (NCSA) at the University of Illinois at Urbana. Many thanks should also go to my parents. Their love and encouragement are the backbone of my persistence in life and my achievements. Most of all I would like to thank my wife, Lijuan. Her endless love and support are the sources of my dedication and inspiration throughout my PhD study. Table of Contents LIST OF TABLES - -- vii LIST OF FIGURES - viii Chapter 1. A Hybrid Lagrangian-Eulerian Particle-Level Set Method for Numerical Simulations of Two-Fluid Turbulent Flows _ - - 1 1.1 Introduction ....................................................................................................... 1 1.2 Mathematical Formulation ................................................................................ 4 1.2.1 Interfacial Particle Level Set Equations ................................................... 5 1.2.2 Zero-Mach Number Variable Density Equations .................................... 7 1.3 Numerical Solution of IPLS Equations ............................................................. 9 1.3.1 Temporal Discretization ........................................................................... 9 1.3.2 Spatial Discretization ............................................................................. 10 1.3.3 The Interface Algorithm ........................................................................ 12 1.4 Numerical Solution of ZMA Equations .......................................................... 16 1.5 Results and Discussions .................................................................................. 18 1.5.1 Single-Fluid Isotropic Turbulence ......................................................... 18 1.5.2 Interface Tracking via IPLS Method ..................................................... 19 1.5.3 Rising Bubble ......................................................................................... 24 1.5.4 Two-Fluid Isotropic Turbulence ............................................................ 25 1.6 Summary ......................................................................................................... 29 1.7 Figures ............................................................................................................. 31 1.8 References ....................................................................................................... 47 Chapter 2. Turbulence-Interface Interactions in a Two-Fluid Homogeneous Isotropic Flow - 52 2.1 Introduction ..................................................................................................... 52 2.2 Results and Discussions .................................................................................. 54 2.2.1 Numerical Accuracy .............................................................................. 55 2.2.2 Physical Structures of Turbulence and Interface ................................... 58 2.2.3 Statistical Analysis ................................................................................. 59 2.3 Summary and Discussions .............................................................................. 70 2.4 Table ............................................................................................................... 73 2.5 Figures ............................................................................................................. 74 2.6 References ....................................................................................................... 93 Chapter 3. LES/FMDF Methodology for Two-Phase Turbulent Reacting Flows.... 98 3.1 Introduction ..................................................................................................... 98 3.2 LES/FMDF for Two-Phase Turbulent Reacting Flows ................................ 106 3.2.1 Governing Equations for Carrier Gas Field ......................................... 108 3.2.2 Two-Phase FMDF Equations for Scalar Field ..................................... 112 3.2.3 Lagrangian Equations for Droplet Field .............................................. 116 3.3 Numerical Solution ...................................................................................... 119 3.3.1 An Eulerian Finite Difference Method for Velocity Field ................... 120 3.3.2 Lagrangian Model for Dilute Spray ..................................................... 122 3.3.3 Lagrangian Monte Carlo Solution of FMDF ....................................... 122 3.4 Consistency of Eulerian LES and Lagrangian FMDF Solutions .................. 126 3.5 Summary ....................................................................................................... 129 3.6 Figures ........................................................................................................... 132 3.7 References ..................................................................................................... 140 Chapter 4. Applications of Two-phase LES/FMDF Model to Spray Combustion . 149 4. 1 Introduction ................................................................................................... 149 4.2 Results and Discussions ................................................................................ 150 4.2.1 Double Swirl Spray Burner .................................................................. 151 4.2.2 Spray-Controlled Dump Combustor .................................................... 155 4.3 Summary ....................................................................................................... 157 4.4 Figures ........................................................................................................... 159 4.5 References ..................................................................................................... 173 vi LIST OF TABLES Table 2-1. The Specifications of DNS Cases .................................................................... 73 vii LIST OF FIGURES Images in this dissertation are presented in color Figure 1-1. A schematic view of the interface, Eulerian grid and Lagrangian particles, illustrating the re-initialization technique in the particle level set method. (i,j) are the coordinates of the Eulerian grid points and [P represents the nearest particle index number. ..................................................................................................................... 31 Figure 1-2. Turbulent statistics calculated by different methods: (a) enstrophy; (b) turbulent energy spectral density at t=4.0; (c) enstrophy dissipation rate; and ((11) Reynolds number based on integral scale and turbulent energy ............................... 32 Figure 1-3. Turbulent statistics obtained by incompressible Pseudo—Spectral and modified IPLS method: (a) enstrophy dissipation rate and (b) Reynolds number based on integral scale and turbulent energy. In the modified IPLS method all the derivatives are calculated by second-order numerical methods .................................................. 33 Figure 1-4. The interface movement as predicted by various methods for the Zalesak problem: (a) level set method without reinitialization; (b) level set method with reinitialization; and (c) IPLS method ........................................................................ 34 Figure 1-5. The interface and particles’ distribution in a vortex flow (vortex stretching problem) as obtained by different methods at t=0.0, 3.0; (a) level set method without reinitialization; (b) level set method with reinitialization; and (c) IPLS method. In Figures (a) and (b), NX=NY=200 and in Figure (c) NX=NY=128. A total of np=2000 particles is used in the IPLS method .......................................................... 35 Figure 1-6. A collapsing circle interface with NX=NY=100, np=1000 as predicted by the IPLS method: (a) the interface at different times and (b) the contours of level set function at t=0.175, the contour levels from outer to interior are (I) = 0.1, 0.05, 0.0, - 0.05, — 0.1 ................................................................................ 36 Figure 1-7. Temporal variation of the interface for: (a) collapsing C-shape and (b) collapsing star. In Figure (a), t=0.0, 1.0, 2.0, 3.0, 4.0 from outer to interior, respectively. In Figure (b), t=0.0, 0.05, 10.0, 0.2, 0.4 from outer to interior, respectively. In both figures, NX=NU=200, np=2200 ............................................. 36 Figure 1—8. Signed distance function for the collapsing C—shape interface at t=2.0 as obtained by: (a) level set method without reinitialization and (b) IPLS method ...... 37 Figure 1-9. Temporal variation of the rising bubble as predicted by the IPLS method for thecase with Re=100,We=200,/i=1:5and77:1:200 .................................. 38 Figure 1-10. Temporal variation of the rising bubble as predicted by the IPLS method for the case with Re =100,We 2 200,1 =1:100and 77 =1: 200 ............................. 39 viii Figure 1-11. Rising bubble as predicted by the IPLS method with different grid resolutions at different times for the case with Re =100,We = 200 ,1 = 1 : 100 and 77 = 1 : 200 ................................................................................................................ 40 Figure 1-12. The predicted vorticity contours and interface in a two-fluid isotropic turbulent flow by the IPLS method: (a) without particle addition algorithm and (b) with particle addition algorithm. ............................................................................... 40 Figure 1-13. Initial contours of (a) density and (b) vorticity for an oil-water two-fluid isotropic turbulent system. ........................................................................................ 41 Figure 1-14. Density contours at different times as obtained by IPLS and ZMA methods in an oil-water two-fluid isotropic turbulent system. Figures (a)—(d) are obtained by the ZMA method and Figures (e)-(f) by the IPLS method. In Figures (a) nd (b), Pr=0.7, and in Figures (c) and (d) P1260. Figure (a), (c) and (e) show the results at t=1.0 and Figures (b), (d) and (f) show the results at t=3.0. In all figures LO/Lyzl/IO,Re=1000,77=200:1,and A=0.8:1. ...................................... 42 Figure 1-15. The vorticity contours as obtained by ZMA and IPLS methods in an oil- water two-fluid isotropic turbulent system at different times. Figures (3) and (b) are obtained by the ZMA method and Figure (0) and (d) by the IPLS method. Figure (a) and (c) show the results at t=1.0 and Figures (b) and ((1) show the results at t=3.0. In all figures, LO/Ly zl/lO, Re=1000, 77=20021,and Z=O.8:1 ..................... 43 Figure 1-16. Turbulent statistics calculated by different methods: (a) enstrophy; (b) turbulent energy spectral density at t=3.0; (c) enstrophy dissipation rate; and (d) Reynolds number based on integral scale and turbulent energy. .............................. 44 Figure 1-17. Interface contours in a two-fluid (oil-water) turbulent flow as predicted by the IPLS method. In Figures (a) and (b) We = 100 and in Figures (c) and (d) We 2 20. Figures (a) and (c) show the results at t=2.0 and Figures (b) and ((1) show the results at t=4.0. In all figure, 77 = 200:1, L0 /Ly 2: 1/10 and Re = 1000. ....... 45 Figure 1-18. Interface contours in a two-fluid (oil-water) turbulent flow as predicted by the IPLS method. In Figures (a) and (b) We = 100 and in Figures (c) and (d) We = 20. Figures (a) and (c) show the results at t=2.0 and Figures (b) and ((1) show the results at t=4.0. In all figure, 77 = 200:1, L0 IL), z1/6 and Re =1000. ........ 46 Figure 2-1. Instantaneous vorticity contours for case 3 at t=2.0 for various grid resolutions and interface thicknesses.(a) e = 2Ax, 256x 256 grid resolution; (b) e = 3Ax, 256x 256 grid resolution; (c) e = 4Ax, 256x 256 grid resolution; and (d) e = 4Ax, 512x512 grid resolution. ......................................................................... 74 Figure 2-2. Temporal variations of the total interface surface areas for various interface thicknesses and various resolutions. ......................................................................... 75 Figure 2-3. Comparison between different forms of baroclinic terms and surface tension terms in vorticity equation for various interface thicknesses: (a)-(b) e = 2Ax; (c)—(d) e = 3Ax; and (e)-(f) e = 4Ax. The results are plotted for point C or particle C on the interface. Grid resolution is 256 X 256 ..................................................................... 76 Figure 2-4. Temporal variations of baroclinic terms and surface tension terms for different resolutions .................................................................................................. 77 Figure 2-5. Vorticity contours and interface t0pologies in case 3 at different times: (a) t=0.0; (b) t=0.5; (c) t=1.5; and (d) t=2.5 ................................................................... 78 Figure 2-6. Instantaneous vorticity contours and interface topologies at t=2.0 for different cases: (a) case 0; (b) case 1; (c) case 3; (d) case 4; (e) case 5; and (f) case 7 ........... 79 Figure 2-7. Temporal variations of mean enstrophy for various cases and in different regions of flows: (a) case 3, fluids 1 and 2 and interface regions; (b) various cases, fluid 1 region; (c) and (d) various cases, interface region ........................................ 80 Figure 2-8. Temporal variations of different terms in the mean enstrophy equation in different regions for various cases: (a) viscous dissipation term III in different regions in case 1; (b) baroclinic term 11 in interface region; (c) surface tension term IV in interface region; and (d) viscous dissipation term 111 in interface region ........ 81 Figure 2-9. Two-dimensional spectral density function of (a) turbulent kinetic energy and (b) enstrophy at t=2.0 ................................................................................................ 82 Figure 2-10. Temporal variations of the total interface areas for different cases. (a) various density ratios and (b) various Weber numbers ............................................. 82 Figure 2-11. Temporal variations of mean turbulent kinetic energy in various regions for (a) case 1 and (b) case 3 ............................................................................................ 83 Figure 2-12. Temporal variations of mean turbulent kinetic energy in interface region for different cases. (a) various density ratios and (b) various Weber numbers .............. 83 Figure 2-13. Temporal variations of different terms in the mean turbulent kinetic energy equation in various regions for case 3: (a) surface tension work; (b) pressure work; (c) viscous dissipation; and (d) viscous stress transport terms ................................. 84 Figure 2-14. Temporal variations of different terms in the mean turbulent kinetic energy equation in interface region for cases with different density ratios: (a) surface tension work; (b) pressure work; (c) viscous dissipation; and (d) viscous stress transport terms .......................................................................................................... 85 Figure 2-15. Temporal variations of different terms in the mean turbulent kinetic energy equation in interface region for cases with different Weber numbers: (a) surface tension work; (b) pressure work; (c) viscous dissipation; and (d) viscous stress transport terms .......................................................................................................... 86 Figure 2-16. Temporal variations of interfacial kinetic energy and its two-components in different cases: (a) case 1; (b) case 3; (c) case 6; and (d) case 5 .............................. 87 Figure 2-17. Temporal variations of normalized interface area change rate and mean surface tension work in interface regions with various density ratios: (a) case 1; (b) case 2; (c) case 3; and (d) case 4 ............................................................................... 88 Figure 2-18. Temporal variations of normalized interface area change rate and mean surface tension work in interface regions with various Weber numbers. (a) case 5; (b) case 6; (c) case 3; and (d) case 7 ............................................................................... 89 Figure 2-19. Temporal variations of normalized viscous dissipation term of mean turbulent kinetic energy and normalized mean enstrophy in interface region: (a) case 1; (b) case 4; (c) case 5; and ((1) case 6 ..................................................................... 90 Figure 2-20. Temporal variations of different terms in the vorticity equation on particle A for different cases: (a) case 1; (b) case 3; (c) case 5; and (d) case 7 ......................... 91 Figure 2-21. Temporal variations of different terms in the vorticity equation on particle B for different cases: (a) case 1; (b) case 3; (c) case 5; and ((1) case 7 ......................... 92 Figure 3-1. A block diagram showing different components of the LES/FMDF and its Lagrangian-Eulerian-Lagrangian flow solver ......................................................... 132 Figure 3-2. Contours of the instantaneous vorticity magnitude, droplets, Monte Carlo particles and grid layout in a spray-controlled dump combustor as obtained by LES/FMDF ............................................................................................................. 133 Figure 3-3. The attributes of LES/FMDF methodology and its LES-FD and FMDF subcomponents ........................................................................................................ 134 Figure 3-4. Contours of instantaneous filtered temperature in the premixed dump combustor without fuel spray, (a) FD values and (b) MC values ........................... 134 Figure 3-5. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations in the premixed dump combustor without fuel spray as obtained by LES/FMDF .......................................................................................... 135 Figure 3-6. Contours of instantaneous filtered fuel mass fraction in the premixed dump combustor with fuel spray, (a) FD values and (b) MC values ................................ 136 xi Figure 3-7. Contours of instantaneous filtered temperature in the premixed dump combustor with fuel spray as obtained by LES/FMDF, (a) PD values and (b) MC values ...................................................................................................................... 136 Figure 3-8. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations as predicted by LES/FMDF in the premixed dump combustor with fuel spray ....................................................................................... 137 Figure 3-9. Contours of instantaneous filtered temperature in the double swirl spray burner, (a) FD values and (b) MC values ............................................................... 138 Figure 3-10. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations in the double swirl spray burner as predicted by LES/FMDF ............................................................................................................. 139 Figure 4-1. A schematic of double swirl-stabilized spray burner and its fuel injector....159 Figure 4-2. Contours of the instantaneous filtered temperature and vorticity in the double swirl spray burner, (a) with enclosure and (b) without enclosure ........................... 160 Figure 4-3. Radial variations of instantaneous filtered temperature at different axial locations for both enclosed and unenclosed double swirl spray burner .................... ................................................................................................................................. 161 Figure 4-4. The 3D instantaneous vorticity ISO-surface and fuel droplets’ positions in the double swirl spray burner, (a) with enclosure and (b) without enclosure ................ 161 Figure 4-5. Contours of the instantaneous filtered temperature and vorticity in an unenclosed double swirl spray burner with different mass loading ratios, MLR and inlet air temperature Tin = 500K , (a) MLR=0.04 and (b) MLR=0.1 ..................... 162 Figure 4-6. Radial variations of instantaneous filtered temperature at different axial locations for an unenclosed double swirl spray burner with different mass loading ratios, MLR and inlet air temperature Tm = 500K ................................................. 163 Figure 4-7. Contours of the instantaneous filtered temperature and vorticity in an unenclosed double swirl spray burner with different mass loading ratios, MLR and inlet air temperature Tm = 300K , (a) MLR=0.04, (b) MLR=0.1, and (Cb) MLR=0.3 ................................................................................................................. 164 Figure 4-8. Radial variations of the droplet temperature in an unenclosed double swirl spray burner for (a) uniform droplet size distribution and (b) nearly Gaussian droplet size distribution ....................................................................................................... 165 Figure 4-9. Radial variations of the droplet mass in an unenclosed double swirl spray burner for (a) uniform droplet size distribution and (b) nearly Gaussian droplet size distribution .............................................................................................................. 166 xii Figure 4-10. Contours of the instantaneous filtered values of temperature and vorticity in an unenclosed double swirl spray burner with (a) spray angle of a = 200 and (b) spray angle of a = 600. .......................................................................................... 167 Figure 4-11. Radial variations of the instantaneous filtered values of temperature in an unenclosed double swirl spray burner with different spray angles a at (a) axial location ofx/ R = 2 and (b) axial location ofx/R = 4 ........................................... 168 Figure 4-12. Contours of the instantaneous filtered values of filtered temperature and fuel mass fraction in a spray-controlled dump combustor for different droplet sizes, (a) SMD=30um; (b) SMD=45um; (c) SMD=60um; and (d) SMD=90um .................. 169 Figure 4-13. Radial variations of filtered temperature and fuel mass fraction of the gas mixture in a spray-controlled dump combustor for different droplet sizes (SMD=20, 35, 60, 90) at two axial locations of x/D=2, 4 ......................................................... 170 Figure 4-14. Radial variations of filtered temperature, the filtered fuel mass fraction and the RMS pressure fluctuations of the gas mixture in a spray-controlled dump combustor for different injection angles ,6 at two axial locations of x/D=2, 6 ....... 171 Figure 4-15. Instantaneous Iso-surface of the filtered pressure of the gas mixture in a spray-controlled dump combustor for different pulsed spray injection frequencies ([1 and duty cycles D ......................................................... _ .......................................... 172 xiii CHAPTER 1 A Hybrid Lagrangian-Eulerian Particle-Level Set Method for Numerical Simulations of Two-Fluid Turbulent Flows 1.1 Introduction Two-fluid flows have been widely studied due to their importance in variety of engineering and scientific problems such as the spray atomization, oil-water separation, ink jet printing and coating [1-3]. One of the important issues in the study of two-fluid flows is the stability and evolution of the interface between the fluids [4-6]. Experimentally, it is often difficult to visualize and quantify the interface motion in transitional and turbulent flows. This is due to intrinsic features of the interface, such as the near-zero thickness and complex topological changes. Numerical methods offer three advantages over the experimental methods: (1) numerical simulations can provide wealth of data on almost everything, (2) various physical processes can be turned on and off at will in numerical simulations [7], (3) it is relatively convenient to numerically study the effects of various fluid/flow parameters. Despite these advantages and extensive efforts that have been devoted to the development of reliable methods, the numerical study of two-fluid turbulent flows has been somewhat limited. A challenge is that in most two- fluid systems the physical properties (e.g. density and viscosity) change abruptly across the interface, creating a hard-to-capture “discontinuity” in the system. Another challenge is associated with the breakup and coalescence of interfaces; processes which are very much dependent on the fluid properties as well as the kinematics of the flow. Generally, the interface motion is computed by (i) “tracking” and (ii) “capturing” methods. Examples of the interface tracking methods are the boundary integral [8-9] and front-tracking or surface-marker [4, 10] methods. In the boundary integral method, only the points on interface are discretized, making it relatively efficient in locating the interface. However, once the interface experiences merging, breaking and other significant topological changes, the boundary integral method often becomes hard to implement [1.1]. The advantage of front-tracking or surface-marker method is that the interface can be tracked explicitly and accurately via marker particles, whose number can be increased to improve the accuracy. In comparison with the interface-capturing method, the computational time needed to resolve the interface accurately is less in the front-tracking method since the interface is represented by Lagrangian particles. However, the front-tracking methods need special efforts like particle reseeding/eliminating to deal with stretching, disconnection and reconnection of the interface. Also, it is difficult to calculate the surface tension and curvature from the Lagrangian particles [12]. In contrast to interface-tracking methods, interface-capturing methods implicitly locate the interface based on an indicator function. There have been two popular interface- capturing methods, one is the Volume-of—Fluid (VOF) method, and the other is the level set (LS) method [13]. In the VOF method, the indicator function is the volume fraction which is governed by a convection equation. This method can handle relatively complex topological changes and is implemented in many commercially available codes [14-16]. The VOF method is usually implemented in two steps, advection and interface reconstruction. The former updates the volume fraction and the latter constructs an approximated interface. The piecewise linear interface reconstruction algorithms that are used in most VOF schemes have lirrrited accuracy and suffer from discontinuity problem [16-18]. Another problem with the VOF method is that it is usually difficult to calculate the interface curvature from the volume fraction. To address some of these problems, higher order interface reconstruction algorithms are suggested [7,19]. In contrast to VOF methods, LS methods do not need interface reconstruction. In the LS method, the interface is modeled as a smooth interface function. Also in this method, the interface always has a fixed thickness to alleviate the numerical oscillations caused by discontinuities in fluid properties across the interface. Therefore, the interface is described as a signed distance function at any time to keep a uniform and constant front thickness. Like VOF method, no additional procedure is needed in the LS method to model the complex topological changes of the interface. Also, the interface related curvature and surface tension can be directly obtained from the level set function. Due to these advantages, the LS method has been widely adopted in numerical simulation of multi-fluid flows [20-24]. However, there are some problems with the LS method, the most important one being the unphysical mass loss/gain. To alleviate this problem, a “re- initialization” procedure is proposed in which the signed distance function is adjusted during the simulation [20,24-26]. The re-initialization procedure cannot completely remove mass loss/gain problem and as shown below generates erroneous results in two- fluid turbulent flows. Hybrid schemes have been suggested to improve the accuracy of VOF and LS methods by combining them [15,27]. Most numerical methods developed for VOF and LS methods so far are based on the low order finite volume and/or finite difference schemes. Examples are the variable density project method of Bell et al. [7,24,28,29] for structured grids and the artificial compressibility matrix-free implicit dual-time scheme of Zhao et al. [30] for unstructured grids. These methods can usually handle large density and viscosity ratios but are not fully applicable to two-fluid turbulent flows due to low order numerical accuracy. Motivated by the fact that the level set function can provide smooth geometrical information about the interface and the Lagrangian particle method can accurately capture the interface location without any mass loss/gain, we have developed and tested a coupled (Lagrangian) interfacial particle and (Eulerian) level set method that takes the advantages of both methods. The new method, termed the interfacial particle level set (IPLS), is different than the other hybrid particle-LS and particle-VOF methods [31-34] and is successfully applied here to various laminar and turbulent flows. The effect of turbulence on the interface and the effect of interface on the turbulence are studied via IPLS method for various flow/interface parameters. This chapter is organized as follows: in Section 1.2, the two-fluid “incompressible” particle-level set equations and the zero Mach number form of the Navier-Stokes (ZMA) equations are presented and compared. The numerical solution procedure for IPLS and ZMA along with the re-initialization technique for IPLS is described in sections 1.3 and 1.4. Section 1.5 discusses the IPLS and ZMA results for some representative two-fluid problems. The main summaries or conclusions are presented in section 1.6. 1.2 Mathematical Formulation Two different sets of equations are presented in this section. The first one corresponds to the interfacial particle level set (IPLS) method and the second one describes the variable-density zero Mach number (ZMA) method. The ZMA method is considered here for the assessment of the IPLS method. The IPLS and ZMA equations and the numerical solution procedures employed for solving them are different. However, these two methods should yield similar results under special limiting conditions, at least in theory. It is shown below that the results generated by them for a two-fluid isotropic turbulent flow are indeed similar. This indicates the consistency and the accuracy of both methods. 1.2.1 lnterfacial Particle Level Set Equations In the IPLS method, instead of treating the two fluid regions separately, a coupled solution is sought via one set of equations and “appropriate” representation of the interface. The non-dimensional equations describing such a two-fluid incompressible and immiscible system include the continuity and momentum equations [35]: -o V-U=O (1.1) LU:3140“.)+_LV'<2#<¢>D>_LW__L;, (1.2) a: pm) Re M) We no Frly where U and p represent the velocity and pressure fields, respectively. Equations (1.1) and (1.2) are non-dime‘nsionalized with the dimensionless parameters being the Reynolds number Re = plLU 0 / #1 , the Weber number We = plLUg /0' , and the Froude number Fr 2 U8 /gL , where ,0,- and ,ui are the density and dynamic viscosity for fluid 1' and L and U0 are the reference length and velocity, respectively. In equation (1.2), the 5'2 vector denotes the gravity direction and D=(VU+(VU)T)/2 is the rate of deformation tensor. Gravity and surface tension parameters are represented by g and 0, respectively. Here, the fluid density p, the viscosity ,u , and the curvature coefficient k are all functions of the level set function 4!) , mm = 4 + (IT/0115“”) #(¢)=77+(1-77)H8(¢) (1-3) K(¢) =V-(V¢/|V¢ll where/i = p2 / ,0] is the density ratio and 77 = #2 / ,ul is the viscosity ratio. The level set function itself is defined as: > 0, if E fluidl $06,!) = = 0,556 interface (1.4) < 0,55 6 fluid2 In the IPLS method, the evolution of the level set function is described by the standard level set equation 5’9 = —(U -V>¢ (1.5) at and the interface movement is obtained via the Lagrangian particle equation 3” U*<‘ ) <16) : x . dt p where Ftp is the location of the particle p and (7*(5c'p) is the fluid velocity at the corresponding particle position. In theory, the Heaviside function is defined as: O ¢O However, to alleviate the numerical difficulty due to discontinuity of the fluid properties across the interface, the interface is given a nominal thickness 8comparable to the grid size and equation. (1.7) is modified as: 0 ¢<-e Hg(¢)= (1+¢/e+sin(7r¢/e)/7r)/2 Mgr: (1.8) 1 ¢>e 1.2.2 Zero-Mach Number Variable Density Equations For numerical simulations of low speed variable density flows, the compressible ’ equations may be solved in the zero Mach number limit. For this, the acoustic field which imposes stringent Courant—Friedrichs-Lewy (CFL) restriction on the time step is removed from the compressible equations [36]. In this limiting condition, the non- dimensional ZMA equations take the following form: 23+V-(pl7)=0 (1'9) 3: 3(10U)+V_(pg0)=_vp+_l_v.r (1.10) at RC = (IT — 1 l at+( ) RePrp (.11 ) ( ) _ 0 pT_p (1.12) where the viscous stress tensor is defined as: Z=fl[20‘%(v'0)l] and L is the unit second-order tensor. In equations (1.9)-(1.12), p0, T are the thermodynamic pressure and temperature, respectively. Note that, p0 is uniform in space and its non-dimensional value is unity. The Prandtl number in temperature equation is defined as Pr: C py/ K , where K is the thermal conductivity, C p is the specific heat at constant pressure, and ,u = ,u(T) = T" is the dynamic viscosity. It should be emphasized here that ZMA and IPLS equations describe the flows in different fluid systems; while ZMA is for an ideal gas, IPLS represents the flow in two incompressible (liquid or gas) fluids, separated by the interface. However, ZMA and IPLS equations become similar for some limiting conditions. When Pr—>oo, the energy (temperature) equation in the ZMA formulation (equation (1.11)) becomes the same as the interface equation in the IPLS formulation (equation (1.5)). Also in this limit, equations (1.1) and (1.2) can be derived from equations (1.9) and (1.10) by ignoring the surface tension and gravity terms. Thus, under the above limiting conditions, ZMA and IPLS equations are the same and should yield similar results for velocity and density, even though the temperature in the ZMA equations might not have any physical meaning. This suggests that in theory the ZMA method is applicable to any two-fluid flow problem. However, in practice this method is limited to “low” density ratios and is not suitable for liquid-liquid or liquid-gas systems. Nevertheless, we are using the ZMA method here for partial validation of the IPLS method. In the numerical solution of ZMA equations presented below, the density and dynarrric viscosity are not constant in each fluid region as they are in the IPLS method. But the variations are found to be very small and negligible in the tested conditions. 1.3 Numerical Solution of IPLS Equations The numerical scheme for the velocity and pressure calculations in the IPLS method is based on a variable density projection algorithm [35,37]. In this algorithm, the momentum equation (1.2) is written as: ELLEELW (1.13) a: pm) where Wz—(U-Vfl7+iv.(2’u(¢)D)——l_k(¢)VH(¢)—_l_i . Re pm» We p<¢> Frlrl By taking the divergence of equation (1.13) and using the continuity equation, one can readily derive the Poisson equation for pressure p , v.(&_):V.W (1.14) PM) Once the pressure field p is determined via equation (1.14), U is updated from the following equation: 294-3% (1.15) 8! mm 1.3.1 Temporal Discretization The time differencing in the IPLS equations is based on the second-order Adams-Bashforth scheme [38] which yields the following difference equations: ¢Il+i :¢N _dtx(%(ljll 'V)¢II __;_(Un—l _V)¢Il-‘l) (1.16) -: [1+1 _ '3" +th 20* '1 N)_ll—]*(-° Il‘l)) (117) AP —x p (2 (11,, 2 xp . t7" =U" +dtx(%VI/" —%W"—l) (1.18) -°* 011+1 _U -VPn+l dt p11 (1.19) With the known interrnedial velocity, [7 * and the time step,dt, equation (1.19) is solved together with the following Poisson equation forU "+1 and P"+1 at time step n+1with an iterative method proposed by Yu [39]: n+1 v.(VP )=(;ll_)v.(j* (1.20) n l ,0 1.3.2 Spatial Discretization To minimize the numerical oscillation near the interface, which is caused by the discontinuity of the density in the Poisson equation, a second-order finite volume method is used to discretize the left hand side (LHS) of equation (1.20). The final form of the discretization scheme can be summarized in the following (two dimensional) compact form: VP . u I : n u u I u u I O ._ I l 1.21 [v (7)1“ Ax(P),+1’J +Bx(P),_LJ +Cx(P),,J+1+Ex(P),,J_1 FX(P),,J ( ) ,0 where, A— l _ 2 _ 2 n. .- 2 n. . n.. Ax )0 l+l/2,j Ax (,0 I+1,} +10 1,1) 8_ 1 _ 2 szpni—1/2,j sz(p"i, j + p"i—1.j) C_ 1 _ 2 (1.22) Ay2pni,j+l/2 A)’2(p"i.j+1 + pniJ) E: 1 2 Ayzp”i, j—1/2 Ay2(p"i, j + p"‘i.j—1) F = (A + B + C + E) 10 For low to moderate density ratios, both standard multigrid method (MG) [40] and multigrid preconditioned conjugate gradient method (MGCG) [41] yield quick convergence when they are applied to the above 5-stencil linear system (equation (1.21) and (1.22)). However. for larger density ratios, the convergence rate of the MG method is relatively slow [42]. Thus in this work we have used the MGCG method to solve the Poisson equation. For the convection term in the interface equation, the following second-order upwind scheme is applied, [(0 vwrid- =(u,-,j +|u,-,j|X3¢,-,j —4¢,-_1,j+¢,-_2,j)/(4Ax)+(u,-,j —|u,-,jl) Xi— 3¢i,j +4¢i+1,j -¢I+2,j)/(4Ax) + (Vi, j + In, jlx3¢i, j -4¢i,j—1 + ¢i, j—2)/(4Ay) + (Vi, j -Ivz', 1|) >< (— 3a, j + 4¢.-, j+1 — a, j+2)/(4A)’) (1.23) The convection term in the momentum equation is solved by a sixth-order central compact scheme proposed by Garnet [43]. Here, the continuity equation (1.1) is employed to rewrite the convection term, (U -V)U as: ((7 -V)17 = f, + g, (1.24) 2 uv where f = [u j , g =[ 7] and the subscripts x and y imply partial differentiation in x LIV V and y directions. Taking f x as an example to show the details of compact scheme, fl- is the sixth-order approximation of ( f x ) ,- and satisfies: 1 ., , 1 ., l4 fi+1-fi-1 1fi+2-fi—2 3f’l f’ 3f”1 9 2m 9 4m ( ) The viscous term is rewritten as 11 _l_V-(2#D)_ 1 {ZwaquDywwyu+va>)} (126) Re p<¢> ' Re M) szMDva Dxtmbyu + va>) where Dx and Dv derivatives are again approximated by a sixth-order central compact scheme. The divergence term in equation (1.20) and the pressure gradient term in equation (1.19) are also calculated by similar compact schemes. The curvature in the surface tension term is rewritten into a new form suggested by Osher [44], 3 kw» = V ' (V¢/|V¢l) : (¢x2¢yy _ 2¢x¢y¢xy + ¢y2¢xx)/IV¢| (1-27) and then approximated by a fourth—order central difference scheme. To stabilize the high order scheme, the unphysical high wavenumber numerical oscillations in velocity and level set function are removed by a sixth-order accurate filter function [45], A. . . d C b aft-1 + fi + 0fi+1= aft +5(fi+3 + fI—3)+§(fi+2 + fi—2) +-2-(fi+1 + fi—l) (1-28) with a -—l—(ll+10a) b —i(15+34a) c —i(-3+6a) d —i(1—2a) 16 ’ 32 ’ 16 ’ 32 where fr.’ is an unfiltered variable at grid point x,- and fi denotes the filtered value. In this work, the filtering effect is very small and is confined to the smallest flow scales with a = 0.475. 1.3.3 The Interface Algorithm In order to keep a uniform and constant interface thickness in the level set method [24], the level set function, e is usually set as a signed distance function close to the . interface. However, it is often difficult to properly maintain this distance function when 12 the interface undergoes significant changes. To resolve this problem, Sussman er al. [24] proposed a “re-initialization” procedure for maintaining e as a proper signed distance function. In this procedure, an auxiliary equation, 3% = Sign(¢0)(1—]V¢|) , (1.29) is solved for the “corrected steady-state” value of the interface function. In equation (1.29),Sign(¢)0) = $0 “M63 +62 , where 19 is a small number comparable to grid size and (2)0 , (I) are the level set functions before and after re-initialization. As mentioned in reference [46], equation (1.29) may introduce significant error into the calculations as the interface is not normally located on the grid points. Even with the re-initialization procedure, the level set method still generates considerable unphysical mass loss/gain [24]. Despite various re-initialization procedures and constraining techniques proposed to alleviate this problem [20,25,26], there still remains considerable numerical error in the level set method. Also, it is often difficult to obtain a pseudo steady-state solution for equation (1.29) in turbulent flows, where the interface geometries are complex and undergoes significant topological changes. In contrast to the particle-level set method described in reference [32], in our IPLS method, the position of interface is explicitly and accurately represented by a set of Lagrangian particles. The level set function provides smooth geometrical information concerning the interface, while the particles are used for keeping the level set function as a signed distance function in the interface region without using any constraining techniques. The interface algorithm and the basic features of the IPLS method can be described in three main steps: 13 (1) Calculate the terms in the right hand side (RHS) of equations (1.5), (1.6) and (1.13) by .9 using the variables at present and previous time steps (i.e.U",¢",xz,U"_1,¢"_1, —. — . . "* ... . . . 1;; 1). Note that the partrcle velocrty U (x2) rs obtarned from the surroundrng Eulerian velocity U" by a fourth-order interpolation scheme. -:H+l Iii-I, ¢n+l and "p (ii) The primary variables at next time step (i.e. U ) are obtained from equations (1.16)- (1.20). (iii) After certain number of iterations, the values of level set function ¢"+l are corrected by using a new re-initialization technique based on the Lagrangian particles’ - ~ + locanons r" l . p to maintain the level set function as a signed distance function. The re- initialization procedure is described in details below with the help of Figure 1-1. For correcting the level set function¢ at an Eulerian grid point, like grid (i,j) in the Fluid 1 side in Figure 1-1, the particle that has the minimum distance to this point is found first . This particle is represented with index IP in Figure 1-1. The indexes of Lagrangian particles resided on the interface are monotonically increased in an initially predetermined direction (anti-clockwise in Figure 1-1) along the interface. So, for the particle with index IF, the tangential direction vector "1' of interface at the particle location can be approximated as? = (3,1)“ — 361,34 )/|.Y',p+l - EULII . The vector starting from particle IP and heading to the node (i,j) is given by E = 5&(i,j) — 52”)“. In the two- dimensional space, the corrected (9,7,1: is the cross product of f and 11¢); j =fo, where the magnitude of (15,7,1- (i.e. h in Figure 1-1) is the approximate distance of node 14 (i,j) to the interface. The positive sign of Q; j indicates that the node (i,j) is in Fluid 1 side and the negative sign of ¢i+l,j—l means node (i+l,j-1) is in the Fluid 2 as shown in Figure 1-1. The advantage of this technique is that the corrected ¢i, j automatically maintains the right signed distance function and less number of particles is needed. Note that, the computational cost of particles’ transportation and re-initialization are linearly dependent on the number of particles. The IPLS method described here has some similarities with the hybrid Lagrangian-Eulerian method proposed by Ceniceros and Roma [47]. In the later, the interface is represented by piecewise linear functions and polygons containing a series of nearby points are used to find the minimum distance. Hybrid particle level set methods are also proposed by Enright et al.[32] and Hieber and Koumoutsakos [33]. However, in these works the particles are uniformly distributed in a narrow band around the interface and are only used for adjusting the zero level set function. For the IPLS method, the particle addition and removal procedures are necessary as the particle number density may become too low or two high in the “stretching” and “compressing” regions of the interface especially in turbulent flows. Inadequate number of particles makes the re-initialization erroneous and cause numerical oscillations as discussed, in Section 1.5.4. Too many particles cause unnecessary computational cost. The particle addition and removal algorithm implemented in the IPLS method are straightforward and are based on the distance between neighboring particles. Taking particles in Figure 1-1 as example, once the distance between IP-l and IP particles becomes larger than an upper threshold, the indexes larger than IP-l is increased by l. and a new particle with index IP is added in the middle location between particles IP-l 15 and IP+1. On the contrary, if the distance between particles 1P-1 and [PH becomes smaller than a lower threshold, the particle IP is removed (or is given a zero weighting) and those indexes larger than IP is reduced by 1. In practice, it is not necessary to implement the particle addition/removal procedure at every time step. 1.4 Numerical Solution of ZMA Equations In this chapter, the predictor-corrector integration algorithm of Najm [48] is used for solving the ZMA equations (1.9)-(1.12). The algorithm is briefly described here. First, the equation of state (equation (1.12)) is used to rewrite the time rate of change of density as: a_p_ = —E—qr— (1.30) a: T 3. Subsequently, equation (1.30) is used together with equation (1.11) to calculate the local time derivatives of (GT/8t)" and (dp/dt)" in the predictor step. The second- order Adams-Bashforth scheme is implemented to evaluate the predicted density ,0* and the Intermedrate velocrtyui, whrle predrcted temperature T rs obtarned from equation (1.12), * n n n—l .___'0 _'0 23(2) __l_{§'2) (1.31) dt 2 a: 2 at *A ' —1 p u,- —p"u{' _§{_ 9(Wi”1)+ifl]n _1[_m.23’2]n (132) dz — 2 dxj Re dxj 2 dxj Re dxj T = (1.33) 16 The intermediate (dynamic) pressure is determined by inverting the Poisson equation, 2 * 1 a(paicuhi) 1 * n n— V P =— + 3 —4 + (II ax, 2dt( p '0 p 1 ) (1.34) Finally. in the projection step, the predicted velocity u: is determined as: >1: =1< *. >1< ,0 u,- —p u; :_BP (135) dt axl' . In the corrector step, the new corrected values for density are calculated via a second-order quasi Crank-Nicolson integration scheme to stabilize the oscillation caused by the variable density effects, n+1 n '1 ** p -p __1. 92. in.) . .11.) +1. ) In equation (1.36), (GT/8t)” and (dp/dt)W are calculated from equations (1.11) and n+1 (1.30) based on the “predicted” values. By replacing p*, T*,P* , and It: with ,0 , T”+1 , P"+1 and ui’IH, equations (1.32)-(1.35) are employed again for the corrected values of T"+1 , P"+1 and uin+l. The Poisson equation is solved by the spectral method and Fast Fourier Transform (FFT). Note that, there is no discontinuity in the Poisson equation which is solved in the predictor and corrector steps. The above predictor-corrector integration scheme improves the maximum computable density ratio to increase [36]. In this work, all spatial derivatives are calculated by hi gh-order spectral/FFT schemes. 17 1.5 Results and Discussions Before applying the IPLS method to a two-fluid turbulent problem, its reliability and accuracy are first established for simpler “standard” test problems. In section 1.5.1, a single-fluid isotropic turbulent flow is simulated via IPLS and ZMA methods by setting the fluid properties the same and by ignoring the surface tension term. The ability of IPLS method in capturing the interface evolution in various moving—interface problems is demonstrated in section 1.5.2. In section 1.5.3, the rising of a low density bubble in a high density fluid is simulated by the IPLS. The results obtained by the IPLS and ZMA methods for a two-fluid isotropic turbulent flow are presented in section 1.5.4. All turbulent simulations are conducted in a square domain with uniform grids and grid resolutions of 256x256 and 512x512. The turbulent kinetic energy spectrums in Figures 1-2(b) and 1-16(b) and all the other results not shown in this chapter indicate that for both single- and two-fluid flows thetturbulence is well resolved with 256x256 grid points. The results shown below are obtained with lower (256x256) grid points to establish the performance of the IPLS method and its interface tracking algorithm in the limiting conditions. 1.5.1 Single-Fluid Isotropic Turbulence In this section, simulations of an isotropic decaying turbulent flow with periodic boundary conditions are conducted to establish the accuracy of the numerical models for single-fluid flows before applying them to two-fluid flows. The initial turbulence in all of our isotropic turbulent flow simulations is the long time solution of an initially random solenoidal velocity field with Gaussian spectrum [49]. In addition to ZMA method, spectral-based incompressible and compressible methods are also employed for further 18 assessment of the IPLS method. The compressible simulation is conducted at low Mach number (Ma=0.05). In the following, several important turbulent statistics as obtained by different methods are compared. This includes the mean turbulent energy k, = (lliul‘)/ 2, the mean enstrophy, £2 = (wiwi>/ 2 and the enstrophy dissipation rate, 771 = (,u/p)<|foul>, whereL = k,“2 /7711/3 and ReL are turbulent macro length scale and the corresponding Reynolds number. Figure 1-2 shows that the predicted long time statistics by the IPLS, ZMA and spectral methods (compressible and incompressible) are virtually identical. This indicates the reliability and the accuracy of both ZMA and IPLS methods in single-fluid turbulent flows. Note that in the IPLS method, the pressure (Poisson) equation is solved with a second-order finite volume scheme to avoid the numerical instability due to discontinuity in density in two-fluid calculations. Nevertheless, the IPLS results remain accurate as long as the other spatial derivatives are approximated with high-order difference schemes. This is confirmed in Figure 1-3, where it is shown that the IPLS method loses its accuracy when the spatial derivatives are calculated by second-order schemes instead of sixth or forth-order schemes. The IPLS results discussed below are all based on high order space differencing. 1.5.2 Interface Tracking via IPLS Method In this section, the IPLS method is used for prediction of the interface movement in two types of problems. In the first one, the velocity field, I7 = 17(55) is independent of the interface and is externally imposed. The so-called Zalesak and vortex stretching problems belong to this type. The other type of problem considered here is the one in 19 which the velocity field is self-generated and is directly dependent on the interface. An example is a closed interface that moves inward with a local velocity which is proportional to the local curvature of the interface, i.e.l7 = —C1KJ\7 (C1 is a constant), where A7 =V¢/|V¢i is the unit normal direction vector, and K=V - (IV) is the curvature. ,, ‘6 The so called “collapsing circle, collapsing C-shape” and “collapsing star-shape” problems are among the second type of problems which are considered here. These problems are more difficult to model or simulate than the first type of problems in which the velocity field is not affected by the interface. 1.5.2.1 Zalesak Problem In the so-called Zalesak problem [50], the interface is a two-dimensional rotating slotted disk that is moved with a specified velocity field. This problem has been considered by many investigators for testing of the interface-tracking numerical models and is considered here for preliminarily assessment of the Lagrangian interface tracking algorithm in the IPLS method. Figure 1-4 shows the IPLS results for the Zalesak problem. The computational domain size is 100x100 units and the slotted circle which has the radius of 15 units is initially centered at (x = 50,y = 75) point. The width and length of the slot are 5 and 25, respectively. The interface is revolved by a pre—specified velocity field defined as u = (7r/314)x(50—y),v = (7r/314)x(x—50), where u and v are the x and y components of the velocity vector, respectively. The level set function, is set as (D < 0 inside the disk, and a > 0 outside the disk. The interface moves from time t = 0 to I = 628 to its initial position after a complete revolution of the disk. Figure 1-4 shows the interface (¢=O) contours at initial (t=0) and final (t=628) times as obtained by the IPLS and standard level set methods. Evidently, the IPLS predictions in 20 Figure 1—4(c) are very accurate as the final and initial interfaces are virtually identical and there is no measurable numerical error even though the grid resolution is moderate in these calculations. Unlike IPLS method which perfectly maintains the interface shape at all locations, the standard level set method is not so accurate and causes significant interface shape change (Figure 1-4). This indicates that the IPLS method is superior to standard level set methods as it can handle complex interface topology with high level of accuracy. It is to be noted that the level set simulations are conducted with no “reinitialization” technique. The level set results obtained with sophisticated reinitialization techniques [26,46] (not shown here) are better than those shown in Figure 1-4(a) but are still considerably less accurate than those predicted by the IPLS method. 1.5.2.2 Vortex Stretching Problem Turbulent flows are usually characterized by complex vortical motions with various length and time scales. Consequently, the interface is expected to experience significant changes in a turbulent environment. In some cases, the distance between two adjacent interfaces or the characteristic length scale of a fluid filament moving in another fluid could become comparable or smaller than the grid size. This could create some difficulties for the standard grid-based methods. The vortex stretching problem involves the stretching and folding of a fluid filament in an isolated vortex and was first introduced by Bell et al. [28] to assess the performance of numerical methods. Here we consider this problem to access the accuracy and reliability of the IPLS method. For this, an initially circular interface is moved according to a pre-specified velocity field: u=—sin2(7zr)xsin(272y),v=sin2(7zy)xsin(27zx). The computational domain size is 21 chosen to be 1x1 unit and the initial circle is centered at (050,075) with radius of 0.15. Theoretically, the circle should be stretched into a long, thin fluid filament which progressively wraps itself around the center [32]. Figure 1-5 shows the predicted interface contours ((7) = 0) and the Lagrangian particles by the IPLS method. Evidently, the interface evolution is more accurately predicted by the IPLS (Figure l-5(c)) than by the standard level set method with and without reinitialization (Figures 1~5(a) and (b)). In fact, the IPLS is able to capture the smallest scales in the interface which are comparable to or even smaller than the mesh size without the necessity of refining the grids near the interface [35]. With the Lagrangian scheme, the filament size is not lirrrited by the grid size in the IPLS method and can be much smaller than the grid size. Accurate results are expected as long as the velocity field is well resolved and the interpolation scheme is sufficiently accurate. 1.5.2.3 Collapsing Circle The collapsing circle is an interesting interface moving problem that has an analytical solution. Here, the tangential velocity is zero on the interface and the normal velocity is a simple function of the interface curvature, K, V" =dr/dt=-K=-l/r (1.37) By integrating the above equation, we get the interface evolution equation, r2 = r02 — 2t , where r is the radius of the circular interface that has the initial radius of r0 = 0.6. Figure 1-6 (a) shows that the predicted interface via IPLS method match with the analytical solution very well at all times. Theoretically, the level set function¢ cannot stay as a perfect signed distance function because the magnitude of interface 22 advection velocity increases with 1/r and the distance between any two «1 contours increases in time in a collapsing circle. However, Figure l-6(b) shows that the level set function, ,9 remains a proper signed distance function in the IPLS method. 1.5.2.4 Collapsing C-shape and Star-shape To further assess the IPLS method, the collapsing of a C-shape and a star-shape interface under the influence of its curvature are also studied. The initial C-shape is made out of a 6x6 square and is shown in Figure 1-7(a). The initial star—shape interface is shown in Figure l-7(b) and is described by the following five-pointed star function: y(s) = 2.0 + 0.6 sin(5 >< 2m)(cos(2m), sin(2m)) + O.9(cos(27[s), sin(27zs)), s 6 [0,1] where s is a number between 0 and 1. The surface advection velocities for both problems are l7 = —IdV, where KIS the curvature and 1V is the unit normal vector as defined before. Figure 1-7 shows that for both problems the interface collapses into a circle due to curvature dependent velocities. Also, in both problems, the interface experiences a relatively complex evolution that is not easy to capture with the standard interface tracking algorithms [51]. Nevertheless, the IPLS method is able to accurately predict the interface movement and its evolution from a C-shape or a star-shape initial form to a circle. The results in Figure 1-8(b) for the level set function are consistent with those in Figure 1-7 and indicate that the IPLS method maintains the level set function, ¢ as a signed distance function very well. In contrast, without re-initialization the standard level set method is not accurate and yields erroneous results for the level set function (Figure l-8(a)). 23 1.5.3 Rising Bubble In this section, the temporal evolutions of a two-dimensional bubble [30] in a fully filled container for various density ratios are studied via IPLS method. Initially, the bubble is a circle with radius R and the computational domain size is 4Rx8R. The density ratio,/i and the viscosity ratio, 77are defined in section 1.2.1, where the fluid 2 is air inside the bubble and the fluid 1 is the liquid in the container. All the other parameters are set to be same as those considered in the numerical simulations of Zhao et al. [30], e.g. Fr : 1.0, Re = 100 and We = 200. Figures .1-9 and 1-10 show the temporal evolution of the rising bubbles with the density ratios of /I =1 :5 and}. = 1 : 100, respectively. In both figures, the viscosity ratio is 77:1:200 and the grid resolution is 256x512. Expectedly, the bubble is driven upward by the buoyancy force and the larger the density ratio, the faster the rising speed would be. The larger pressure gradient on the lower surface of the bubble induces a jet of liquid that pushes the lower surface upward and causes the formation of vortex sheet around the lower surface. Owing to the presence of surface tension, the lower surface cannot fully penetrate and break the upper surface and the bubble takes a mushroom shape. At later times, the vortex sheet around the lower surface generates two symmetrical bubble skirts which folds inward. Eventually, the two bubble skirts detach from the bubble under the intensive stretching caused by the vortex sheet. Note that, although our results for density ratio of /I = 1:5 are similar to those in Zhao et al. [30], there is a noticeable difference between the results obtained with smaller and larger density ratios. Figure 1-10 shows that the shedding bubbles can continue to break up into two separate smaller bubbles at later times, which is a reasonable and physically possible 24 phenomenon that may not be well captured by low-order numerical methods. Also, unlike some of the other numerical methods which are not able to predict shedding bubbles and sometimes generate (unphysical) asymmetric results, the IPLS predictions are accurate and symmetric. These results further confirm that the H’LS method is able to accurately predict the two-fluid flows. In Figure 1-11, the rising bubbles as predicted by the IPLS method with different grid resolutions are compared at different times. At earlier time (t=3.0), there is virtually no difference between the interfaces predicted with different resolutions. At later time (t=6.0), only the smallest scales of the broken bubble are not adequately resolved in the lower resolution case and most of the interface including the two shedding bubbles are similarly predicted in both cases. This indicates that the employed high resolution grid system used in Figures 1-9 and 1-10 is adequate. 1.5.4 Two-Fluid Isotropic Turbulence The two-fluid turbulent flow simulations considered in this chapter are two- dimensional. We are well aware of the differences (and similarities) of material- element/interface transport in two-dimensional and three-dimensional turbulent flows. However, even a two-dimensional two-fluid isotropic turbulent flow is a very challenging flow to model and simulate. We consider this flow as a relatively simple two-fluid turbulent flow that has some important features of more realistic three- dimensional inhomogeneous flows. Here, not only the interface experiences complicated changes (e. g. stretching, folding and breaking) due to turbulence, but the turbulence itself is affected by the interface in a two—way coupling manner. To correctly capture the interface movement in turbulent flows and the subsequent effects of the interface on turbulence, a particle addition/removal algorithm is 25 implemented in the IPLS method. This is an inexpensive algorithm which is based on the distance between neighboring particles. When the distance becomes larger than an upper threshold, new particles are added and particles’ indexes are adjusted accordingly. Alternatively, when the distance between particles becomes smaller than a lower threshold value, particles are removed and again their indexes are adjusted to maintain a proper order. The importance of particle addition algorithm is illustrated in Figure 1-12, where the contours of the vorticity magnitude and the interface, obtained by the IPLS method, in two-fluid isotropic turbulent flow are shown. Evidently, without the particle addition algorithm, the number of particles is inadequate in the “stretched” regions of interface which may cause some numerical oscillations as described in section 1.3.3. However, the IPLS method with particle addition algorithm is able to accurately predict the interface evolution in a two-fluid turbulent system. As shown in Figure 1-13(a), the two-fluid system considered in this chapter is an oil layer that is surrounded by water. The oil thickness is changed by changing L0 / Ly , where L0 and L), are the initial thicknesses of oil layer and domain size in the normal (y) direction, respectively. A fully developed isotropic turbulence, illustrated by the vorticity contours in Figure 1-13(b), is implemented as the initial velocity field and periodic boundary conditions are applied in all directions. This two-fluid flow is simulated by the IPLS and ZMA methods. In the IPLS method, the density and viscosity ratios are set to fixed values of A = 0.8 :1 and 77 = 200:], while the other parameters are varied. In the ZMA method, the desired viscosity ratio is obtained from ,u=,u(T)=T" and equation (1.12). 26 As indicated before, in the IPLS method, the interface can be identified by either particles’ locations or the zero level set (¢=0) contours. In contrast, in the ZMA method. the interface is identified by monitoring the density itself. This is illustrated in Figures l-l4(a)-(d), where the results obtained by the ZMA method for two different Prandtl numbers are shown. Figures 1-14(a) and (b) show that for Pr = 0.7 the initial sharp density (or temperature) gradient is smoothed and the minimum density increases from 0.8 to 0.87 due to relatively high values of the thermal diffusivity coefficient in the energy equation. In contrast, the initial density jump across the interface is well maintained when Pr =60 and the thermal diffusion is low (Figures 1-l4(c) and (d)). Also for Pr = 60, the results obtained by the ZMA method are reasonably close to those predicted by the IPLS method when Weber number is large or surface tension is small (Figures 1-14(e) and (f)). Furthermore, Figure 1-15 shows that the computed turbulent vorticity contours via IPLS and ZMA methods are very similar; indicating the consistency of these two methods and the accuracy of their numerical solution procedures. Further evidence is provided in Figure 1-16, where several statistical quantities, obtained from the ZMA and IPLS data are compared. It is to be emphasized here that the governing equations and the numerical solution procedures in the ZMA and IPLS methods are very different. Nevertheless, the interface evolution, the vorticity field and all the flow/turbulence statistics, predicted by these two methods, are similar when the Prandtl number in the ZMA equations and the Weber number in the IPLS equations are sufficiently large. This supports the argument we made in section 1.2.2 that when the surface tension and the thermal diffusion are negligible, both IPLS and ZMA methods should yield similar results. Note that, since the 27 numerical approximation in the ZMA method is based on a high-order Pseudo-spectral scheme, the similarity of the IPLS and ZMA results indicates that the IPLS method is indeed a hi gh-order numerical method for two-fluid turbulent calculations. The results in Figures 1-17 and 1-18 confirm that the IPLS method is able to accurately predict the two-fluid turbulence and interface evolution for a variety of flow parameters. Figure 1-17 shows the interfaces for various values of the Weber number. It is clear that for a higher Weber numbers, the interface becomes more “active” and less “stable” as the surface tension decreases and the interface shape becomes very complex. In fact, in some cases, the interface is broken up into fluid ligaments or droplets. Nevertheless, the IPLS method is able to capture the complicated interactions between the turbulence and the interface. The effects of the oil thickness or L0 / Lyon the interface are shown in Figure 1-18. A comparison between the results in Figures 1-17 and 1-18 indicates that for thinner oil layer ( L0 /Ly = 1/10 ), the interface is broken up quickly into several portions and the oil ligaments entraps in the surrounding water. In this case, the overall effect of the turbulence on the interface is indeed substantial. In contrast, for the thicker oil layer (L0 /Ly = l/ 6), the large distance between the two neighboring vorticity fields in water reduces the global effect of turbulence and the pinching effect of the turbulence on the oil layer becomes more important. This is shown in Figure 1-18, where it is observed that some small ligaments or droplets are pinched off from the thick oil layer by small-scale turbulent motions, while the large scale motions are not able to completely break up the layer. The formation of ligaments and droplets in a liquid-gas system is studied by Boeck er al. [52], who discuss the effects of Weber number and Reynolds number in details. Again the results in Figures 1—17 and 1-18 demonstrate the ability of the IPLS method to 28 capture the complicated features of the interface evolution in two-fluid isotropic turbulent flows. 1.6 Summary A hybrid Lagrangian-Eulerian interfacial particle level set (IPLS) method is developed for numerical simulation of two-fluid turbulent flows. In this method, the interface is identified based on the locations of notional (Lagrangian) particles and the geometrical information concerning the interface and the fluid properties are obtained from the (Eulerian) level set function. The Lagrangian particles can accurately represent the interface evolution without any numerical mass loss or gain and the level set function provides smooth geometrical information concerning the interface and its effects on the flow. The new methodology has been assessed by applying it to several standard interface-moving and two-fluid laminar problems. It has also been used for simulation of two-fluid isotropic turbulent flows under various flow/fluid conditions. The numerical results are evaluated by monitoring the mass conservation law, the turbulence energy spectral density function and the consistency between Eulerian and Lagrangian components of the method. It is shown that the IPLS method can handle interfaces with substantial topological complexity, and accurately predict the interface evolution even at the “sub—grid range” as long as the velocity field is well resolved. The accuracy of the Navier-Stokes flow solver in the IPLS method is established first by comparing the results obtained by the IPLS method for a single-fluid isotropic turbulent flow with those obtained via high-order incompressible and compressible Pseudo-Spectral numerical schemes. The IPLS results for a two-fluid isotropic turbulent flow are also compared with those predicted by a high-order zero-Mach number variable density flow solver. 29 Despite the differences in the equations and the numerical schemes, we have found that the IPLS method generates virtually identical results with other methods. The similarity of the results confirms the accuracy of the IPLS flow solver in single—fluid and two-fluid turbulent systems. Analysis of the velocity statistics and vorticity/interface contours in turbulent flows indicate that the destabilization effect of turbulence and the stability effect of surface tension on the interface motion are strongly dependent on the fluid density and viscosity ratios. As expected, the interface becomes more stable and less “active” as the surface tension increases or the Weber number decreases. The turbulence can completely break up a thin layer of oil in water. However, for an initially thick layer, turbulence can only pinch off small ligaments or droplets from the layer. Despite the complex and sometimes very significant effects that the turbulence and the interface have on each other, it is shown that the IPLS method can correctly capture the interface evolution and the effect of turbulence on the interface in all the cases considered in this study. 30 1.7 Figures 1 \2 If” \0\9\ ‘1’ 7 In IVA 1? . . (1.1-1) (1+1.1-1) e96 (iii-2) (i+17j-2) Figure 1-1. A schematic view of the interface, Eulerian grid and Lagrangian particles, illustrating the re-initialization technique in the particle level set method. (i, j) are the coordinates of the Eulerian grid point and IP represents the nearest particle index number. 31 400 300 012t3 4 5 Figure 1-2. Turbulent statistics calculated by different methods: (a) enstrophy; (b) turbulent energy spectral density at t=4.0; (c) enstrophy dissipation rate; and (d) 350 _ ReL_ 250 I Q (d) Spectral_Incomp ...... Spectral_Comp —o— ZMA IPLS O 1 1 1 1 l 1 L A MEL L #4 41 I 1 1 A l l A A 1 1 200 A A ‘ L Reynolds number based on integral scale and turbulent energy. 32 400 l ( b) : Spectral_Incomp 350 - “ —-O— IPLS_Low oder Re L — 300 I 250 h 200 1 A A 1 LLL 1 A I l l 1 A L l L L x l 1 A 1 1 A Figure 1-3. Turbulent statistics obtained by Incompressible Pseudo-Spectral and modified IPLS method: (a) enstrophy dissipation rate and (b) Reynolds number based on integral scale and turbulent energy. In the modified IPLS method all the derivatives are calculated by second-order numerical methods. 33 100 100 f (a) I (b) 80- , 80; _ , f \ _ \ t 60- ‘ 60~ : y: Y. t=0.0 - 40— 40- — ------- t=628.01 . , 20; 20- OF 1 11.11 0-444411 ILL.rrl._r. 0 20 40 x 60 80 100 0 20 40 x 60 80 100 100 ' (C) 80 ' I 60* 40 20 I I 0 J l A J J A J 4 #9. A A l A J l 0 2o 40x 60 80 100 Figure 1-4. The interface movement as predicted by various methods for the Zalesak problem: (a) level set method without reinitialization; (b) level set method with reinitialization; and (c) IPLS method. 34 0.8; 0.6:- v . 0.4e 0.2; 00 * . 1 1 0.8 l , t ,’ 0.6 "' l v » ' 0.4 - Particels at t=0 . Particels at t=3 0-2 _' Interface at t=0 — Interface at t=3 OLALL111_L 14 441111 0 0.2 0.4 X 0.6 0.8 1 Figure 1-5. The interface and particles’ distribution in a vortex flow (vortex stretching problem) as obtained by different methods at t=0.0, 3.0; (a) level set method without reinitialization; (b) level set method with reinitialization; and (c) IPLS method. In Figures (a) and (b), NX=NY=200 and in Figure (0) NX=NY=128. A total of np=2000 particles is used in the IPLS method. 35 2 2 I (a) . (b) 1.5; 1.5;~ C r Y: Y: . I 0.5- 0.5r O~.1.1L...11...l.11. OIIAIJrrrLIAMIIIIJII 0 0.5 x 1 1.5 2 0 0.5 x1 1.5 2 Figure 1-6. A collapsing circle interface with NX=NY=100, np=1000 as predicted by the IPLS method: (a) the interface at different times and (b) the contours of level set function at t=0.175, the contour levels from outer to interior are ¢ = 0.1, 0.05, 0.0, — 0.05, — 0.1 10 OPA.IIIA_II.L11..I o 2 4x6 Figure 1-7. Temporal variation of the interface for: (a) collapsing C-shape and (b) collapsing star. In Figure (a), t=0.0, 1.0, 2.0, 3.0, 4.0 from outer to interior, respectively. In Figure (b), t=0.0, 0.05, 0.1, 0.2, 0.4 from outer to interior, respectively. In both figures, NX=NY=200, np=2200. 36 10 10 i (a) i (b) 8L 8L 6- 6;- YZ v. 4; 4- 2- 2L 0 l l l l 044. l J 44L O 2 4x6 8 10 0 2 4XG 8 10 Figure 1-8. Signed distance function for the collapsing C-shape interface at t=2.0 as obtained by: (a) level set method without reinitialization and (b) IPLS method. 37 ” t=0.0 * t=1 .0 ' t=2,0 ’ t=3.0 6~ 6f 6- 6— I r f : .. s. E . S - E 4- 4- 4 4 C 1 gel 1 C 1 124444“ ....1....1L...1.... 14.4....1.... 0 1 2 3 4 O 1 2 3 4 CO 1 2 3 4 C0 1 2 3 4 XIR XIR x/n xm 8 8 8 8 ~ t=4.0 . - t=5.0 ~ t=6.0 . t=7.0 l- . . . 6e 6- 6. 5- e- e- E E- 4- 4r 4F 4" I I C l 2" 2r Z-U U 2- C Z I i CPUUIHUIHHLLH CF....IL...LL.LL1. ....I....1....I. P...-J.L..r....1.... O 1 2 3 4 0 1 2 3 4 GO 1 2 3 4 C0 1 2 3 4 X/R XIR X/R XIR Figure 1-9. Temporal variation of the rising bubble as predicted by the IPLS method for the case with Re = 100, We = 200, xi = 1:5 and 77 =1: 200. 38 * t=0.0 si- mi- \ >1. 4)— 2- 0....1....l Llrrr O 1 2 3 X/R 8 f t=4.0 6!- m. \ >.. 42 1 2h 0 ILLLIL l Lrlraaa 0 1 2 3 XIR 8 8 ~ t=1.0 ~ t=2.0 6' 61- C” m» 4- 4.. : t 2- 2- t _ 0.11.1...111...l... G-lull+4qlgul,l,_,, 0 2 4 O 1 2 3 XIR x13 8 8 1 t=5.0 - t=6.0 6- 6- ' l rc~ me §r ;~ 4" 4)— ..b 0) .(l i) C-.. rel“ .1 cl .1. 1......1 0 2 4 0 1 2 3 XIR XIR 8 * t=3.0 f 6.— E». >1- 41- t 21— 01.141...l..-.1..1. 0 1 2 3 XIR 8 ‘ t=7.0 6i E). >r 4r- 2- ./ \ c....r....1..-.1.... 0 1 2 3 XIR Figure 1-10. Temporal variation of the rising bubble as predicted by the IPLS method for the case with Re = 100, We = 200, 1 = 1:100 and 77 = 1: 200. 39 r t=3.0 t=3.0 U“ t=6.0 0’ 1:6.0 . 256X512 _ 128X256 25BX512 128x256 s~ 6— s- - e» _ m: I: >. E- ;- :— 4— 4- 4- 4— l _ . - 2- 2— 2- Q P) 2— 6 P) C l l l — I I l ”L l l I h I I....I._._._. 0 1 2 3 4 G0 1 2 3 C0 1 2 3 c0 1 2 3 4 XIR XIR X/R XIR Figure 1-11. Rising bubble as predicted by the IPLS method with different grid resolutions at different times for the case with Re=100, We =200, /i. =1:100 and 7721:200. ‘ 1 . \ 7r \\ 3‘\\ \\" \~ W J34 /T\\‘ ,5 0 ‘1V 2? / \// ll 1 I .2\ g; ’l IA “€le .. , .,. gs‘. (A) r \\ \ \» \\ J Figure 1-12. The predicted vorticity contours and interface in a two-fluid isotropic turbulent flow by the IPLS method: (a) without particle addition algorithm and (b) with particle addition algorithm. 40 Fluid 2 12x3 4 5 6 Figure 1-13. Initial contours of (a) density and (b) vorticity for an oil-water two- fluid isotropic turbulent system. 41 (a) me we W @F" A\ Figure 1-14. Density contours at different times as obtained by IPLS and ZMA methods in an oil-water two-fluid isotropic turbulent system. Figures (a)-(d) are obtained by the ZMA method and Figures (e)-(f) by the IPLS method. In Figures (a) and (b) Pr =0.7, and in Figures (c) and (d) Pr = 60. Figures (3), (c) and (e) show the results at t=1.0 and Figures (b) , (d) and (f) show the results at t=3.0. In all figuresLO/Lv zl/IO, Re=1000, 77=200:1,and 2:082]. 42 ear; ] I :2 f;— x". J. \_ /" — 245‘s??? Figure 1-15. The vorticity contours as obtained by ZMA and IPLS methods in an oil-water two-fluid isotropic turbulent system at different times. Figures (a) and (b) are obtained by the ZMA method and Figures (c) and (d) by the IPLS method. Figures (a) and (0) show the results at t=1.0 and Figures (b) and (d) show the results at t=3.0. In all figures L0 IL, =1/10,Re =1000,77 = 200:1 and xi = 0.8 : 1. 43 1 0" EIK) 1 0'6 1 0'8 2 1.1.11.1...l..11 10.10 0 1 t 2 3 4 100 400 ( d ) ------ ZMA with Pr=0.7 1 —— ZMA with Pr=60 ‘, A IPLS with We=100 300 REL. 200 0....LAILILLII.III.1 1OO,J,,1,,.,|.,.LLLJ.. 0 1 12 3 4 o 1 t2 3 4 Figure 1-16. Turbulent statistics calculated by different methods: (a) enstrophy; (b) turbulent energy spectral density at t=3.0; (c) enstrophy dissipation rate; and (d) Reynolds number based on integral scale and turbulent energy. 44 (a) (b) @125“ KIN A < ( (e) (d) 8‘)... Figure 1-17. Interface contours in a two-fluid (oil-water) turbulent flow as predicted by the IPLS method. In Figures (a) and (b) We = 100 and in Figures (0) and (d) We = 20. Figures (a) and (c) show the results at t=2.0 and Figures (b) and ((1) show the results at t=4.0. In all figures,77=200:1, LO/Ly =1/10 and Re =1000. 45 (a) (b) (C) (d) W Figure 1-18. Interface contours in a two-fluid (oil-water) turbulent flow as predicted by the IPLS method. In Figures (a) and (b) We = 100 and in Figures (0) and (d)We = 20. Figures (a) and (c) show the results at t=2.0 and Figures (b) and (d) show the results at t=4.0. In all figures, 77:20021, L0/Ly z1/6 and Re=1000. 46 1.8 [1] [2] I3] [4] [5] l6] [7] [8] I9] [10] [11] [12] References Li, X., “Mechanism of Atomization of a liquid jet,” Atomization and Sprays Vol. 5, 1995; PP. 89-105. Mamora, D. D, Sutadiwiria G., “Analytical mode of distillation of oil at the steam-oil interface,” The seventh UNITAF? international conference on heavy crude and tar sands, Beijing, China, Oct. 1998, pp.27-30. Torpey, P.A., “Prevention of air ingestion in a thermal ink—jet device,” Proceeding of the 4th international congress on advances in non-impact print technologies, Springfield, VA, Mar.1998. Esmaeeli, A. and Tryggvason, G., “An inverse energy cascade in two- dimensional low Reynolds number bubble flows,” Journal of Fluid Mechanics, Vol. 314, 1996, pp. 315-330. Herrmann, M., “Modeling primary breakup: A three-dimensional Eulerian level set/vortex sheet method for two-phase interface dynamics,” Annual Research Briefs, Center for turbulence research: Stanford, 2003. 1313.185- 195. Reitz, R. D. and Bracco, F. V., “Mechanisms of breakup of round liquid jets,” Encyclopedia of fluid mechanics, Vol. 3, 1986, pp.233-249. Puckett, E. G., Almgren, A. S. and Bell, J. 8., etc, “A high-order projection method for tracking fluid interfaces in variable density incompressible flows,” Journal of Computational Physics, Vol. 1 30, 1997, pp.269-292. Boulton-Stone, J. M. and Blake, J. R., “Gas bubbles bursting at a free surface,” Journal of Fluid Mechanics, Vol.254, 1993, pp.437-466. Strain, J., “A boundary integral approach to unsteady solidification," Journal of Computational Physics, Vol.85, 1989, pp.342-389. Unverdi, S. O. and Tryggvason, G., “A front-tracking method for viscous, incompressible, multi-fluid flows,” Journal of Computational Physics, Vol.100, 1992, pp.25-37. Sussman, M. and Smereka, P., “Axisymmetri free boundary problem,” Journal of Fluid Mechanics, Vol. 341, 1997, pp.269-294. Helmsen, J. J., “A comparison of three-dimensional photolithography simulators,” Ph.D. thesis, U. C. Berkeley, 1994. 47 [13] [14] [15] [15] [17] [18] [19] [20] [21] [22] [23] [24] Osher, S. and Sethian, J. A., “Fronts propagating with curvature- dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, Vol. 79, 1988, pp. 1 2-49. Mashayek, F. and Ashgriz, N., “Advection of axisymmetric interfaces by volume of fluid method,” International Journal for Numerical Methods in Fluids, Vol. 20, 1995. 991337-1361 . Mashayek, F. and Ashgriz, N., “A hybrid finite element - volume of fluid method for simulating free surface flows and interfaces,” International Journal for Numerical Methods in Fluids, Vol. 20, 1995, pp.1363-1380. Scardovelli, R. and Zaleski, 8., “Direct numerical simulation of free-surface and interfacial flow,” Annual Review of Fluid Mechanics, Vol. 31, 1999, pp.567-603. Ashgriz, N. and P00, J. Y., “FLAIR: flux line-segment advection and interface reconstruction,” Journal of Computational Physics, Vol. 93, 1991, pp.449-468. Youngs, D. L., “Time-dependent multi-materials flow with large fluid distortion,” in Numerical Simulation for Fluid Dynamics, Morton, K. W. and Baines, M. J. (eds). Academic Press: New York, 1982, pp.274-285. Pilliod, J. E. and Puckett, E. G., “Second-order accurate volume-of-fluid algorithms for tracking materials interface,” Technical Report No. LBNL 40744, Lawrence Berkeley National Laboratory, 1997. Chang, Y. C., Hou, T. Y., Merriman, B., and Osher, S., “A level set formulation of Eulerian interface capturing methods for incompressible fluid flows,” Journal of Computational Physics, Vol. 124, 1996, pp.449- 464. Chen, S., Merriman, B., Osher, S. and Smereka, P., “A simple level set method for solving Stefan problem,” Journal of Computational Physics, Vol. 135, 1997, pp.8-29. Mulder, W., Osher, S. and Sethian, J. A., “Computing interface motion in compressible gas dynamics,” Journal of Computational Physics, Vol. 100, 1992, 139209-228. Russo, G. and Smereka, P., “A level-set method for the evolution of faceted crystals,” SIAM Journal on Scientific Computing, Vol. 21, 2000, pp.2073-2095. Sussman, M., Smereka, P. and Osher, S., “A level set approach for computing solution to incompressible two-phase flow,” Journal of Computational Physics, Vol. 1 14, 1994, pp. 1 46-1 59. 48 [25] [25] [27] [28] [29] [30] [31] [32] I33] [34] [35] [36] Hou, T. Y., Li, Z., Osher, S. and Zhao, H., “A hybrid method for moving interface problems with application to the Hele-Shaw flow,” Journal of Computational Physics, Vol. 134, 1997, pp.236-252. Sussman, M., Fatemi, E., Smereka, P. and Osher, S., “An improved level set method for incompressible two-phase flows,” Computers and Fluids, Vol. 27, 1998, pp.663-680. Sussman, M., “A second order coupled level set and volume of fluid method for computing the growth and collapse of vapor bubbles,” Journal of Computational Physics, Vol. 187, 2003, pp.110-136. Bell, J. B., Colella, P. and Glaz, H. M., “A second-order projection method for the incompressible Navier-Stokes equations,” Journal of Computational Physics, Vol. 85, 1989, pp.257-283. Bell, J. B. and Marcus, D. L., “A second-order projection method for variable-density flows,” Journal of Computational Physics. Vol. 101, 1992, pp.334-348. Zhao, Y., Tan, H. H. and Zhang, B., “A high-resolution characteristics- based implicit dual time-stepping VOF method fro free surface flow simulation on unstructured grids,” Journal of Computational Physics, Vol. 183, 2002, pp.233-273. Aulisa, E, Manservisi, S. and Scardovelli, R., “A mixed markers and volume-of—fluid method for reconstruction and advection of interfaces in two-phase and free-boundary flows,” Journal of Computational Physics, Vol. 188, 2003, pp.611-639. Enright, D., Fedkiw, R. Ferziger, J. and Mitchell, 1., “A hybrid particle level set method for improved interface capturing,” Journal of Computational Physics, Vol. 183, 2002, pp.83-116. Hieber, S. E. and Koumoutsakou, P., “A Lagrangian particle level set method,” Journal of Computational Physics, Vol. 210, 2005, pp.342-367. Raad, P. E. and Bidoae, R., “The three-dimensional Eulerian-Lagrangian marker and micro cell method for the simulation of free surface flows,” Journal of Computational Physics, Vol. 203, 2005, pp.668-699. Sussman, M. and Almgren, A, 8., “An adaptive level set approach for incompressible two-phase flows,” Journal of Computational Physics, Vol. 148, 1999, pp.81-124. Cook, A. W., and Riley, J. J., “Direct numerical simulation of a turbulent reactive plume on a parallel computer,” Joumal of Computational Physics, Vol. 129, 1996, P9263283. 49 [37] [38] [39] [40] [41 l [42] [43] [44] [45] [45] [47] [48] Almgren, A, 8., Bell, J. B. and Szymczak, W. G., “A numerical method for the incompressible Navier-Stokes equations based on an approximate projection,” SIAM Journal on Scientific Computing, Vol. 17, 1996, pp358- 369. Kim, J. and Moin, P., “Application of a fractional-step method to incompressible Navier-Stokes equations,” Journal of Computational Physics, Vol.59, 1985, pp.308-323. Yu, X., “Iterative pressure passion equation method for solving the unsteady incompressible N-S equations,” Chinese Journal of Numerical Mathematics and application, Vol. 23, 2001, pp.447-456. Briggs, W. L., “A multigrid tutorial,” SIAM, Philadelphia, PA, 1987. Tatebe, 0., “The multigrid preconditioned conjugate gradient method,” In 6th Copper Mountain Conference on Multigrid Methods, Copper Mountain, CO, April 1993. Alcouffe, R. E., Brandt, A., Dendy, J. E. and Painter, J. W., “The multigrid method for the diffusion equation with strongly discontinuous coefficient,” SIAM Journal on Scientific Computing, Vol. 2, 1981, pp.430-454. Garnet, L., Ducros, F., Nicoud, F. and Poinsot, T., “Compact finite difference scheme on non-uniform meshes, application to direct numerical simulations of compressible flows,” International Journal for Numerical Methods in Fluids, Vol. 29, 1999, pp.159-191. Osher, S. and Fedkiw, R, Level set methods and dynamic implicit surfaces, Applied Mathematical Sciences, Springer: New York, 2002. Lele, S. K., “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, Vol. 103, 1992, pp.16-42. Dupont, T. F. and Liu, Y., “Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function,” Journal of Computational Physics, Vol. 190, 2003, PP.311-324. Ceniceros, H. D. and Roma, A. M., “A multi-phase flow method with a fast, geometry based fluid indicator,” Journal of Computational Physics, Vol. 205, 2005, [39391-400 Najm, H. N., Wyckoff, P. S. and Knio, O. M., “A semi-implicit numerical scheme for reacting flow. I Stiff chemistry,” Journal of Computational Physics, Vol. 143, 1998, pp.381-402. 50 [49] [50] [51] [52] Jaberi, F. A. and James S., “Effects of chemical reaction on two- dimensional turbulemce,” Journal of Scientific Computing, Vol. 14, 1999, pp.31-72. Zalesak, S. T., “Fully multidimensional flux-corrected transport, algorithms for fluid,” Journal of Computational Physics, Vol. 31 , 1979, pp.335-362. Ong, B., “Level set Framework for fronts propagating with curvature dependent speeds,” Level Set Research-Curvature Motion, 2002. http://mathqeek.us/research/papers/curvature/curvature.pdf. Boeck, T., Li, J., Lopez-Pages, E., Yecko, P. and Zaleski, S., “Ligament formation in sheared liquid-gas layers,” Theoretical and Computational. Fluid Dynamics, Vol. 21, 2007, pp.59-76. 51 CHAPTER 2 Turbulence-Interface Interactions in a Two-Fluid Homogeneous Isotropic Flow 2.1 Introduction Liquid atomization, the process of liquid jet break up into droplets in sprays, is a typical two-fluid flow problem which has been extensively studied due to its importance in numerous applications, which range from combustion in aircraft engines to spray forming and coating in material processing. Fundamental understanding of liquid atomization is indispensable to achieve the best performance in all of these applications which require atomization to produce desired droplet size and velocity distributions. All atomizers force the liquid to break up by disrupting the balance between the consolidating effects of surface tension, liquid viscosity and the aerodynamic forces created by relative motion between the gas and liquid phases at their interface [1]. When the flow is laminar, a vibration in the liquid jet or an external disturbance could cause atomization as predicted by the classical Rayleigh [2] and Kelvin-Helmholtz [3-4] break up theories. However, when the flow is turbulent, the break up of the liquid jet or sheet can take place from the rapid deformation in the liquid-gas interface created by the turbulent structures; the average relative motion between the phases may not be that critical. The atomization in turbulent flows usually takes place in two successive steps as defined by Faeth [5]. The initial part of liquid jet/sheet break up into large ligaments and 52 droplets takes place in the so-called primary breakup zone, where the large-scale liquid jet/sheet has strong interactions with the gas at their interface. Secondary breakup refers to the break up of the initially formed ligaments and droplets into smaller droplets until the stable droplet sizes are achieved. In the region followed by the secondary break up zone, the average droplet size and the volume occupied by the droplets are small. The droplets in this region, referred to as dilute spray region, are often treated as point particles [6], which are dispersed by turbulent fluid motions. In the primary break up region, the interactions between the turbulence in both fluids occur through a very sharp interface. 80, besides turbulence, the interface’s location, topology and curvature have to be precisely simulated to assure accuracy. However, the accurate prediction of interface in turbulent flows is extremely difficult since it involves highly convoluted interface and the breakup of liquid ligaments or droplets. Furthermore, the turbulence near the interface is affected by the interface. For example, the normal component of the turbulent kinetic energy is found to be redistributed into the tangential component through pressure fluctuations [7] or dissipated by the surface tension. Also, the deformed interface could generate high- intensity, small-scale vort'ical motions. The importance of liquid atomization has motivated extensive number of theoretical, numerical and experimental investigations [1, 8-12]. Though much has been learned, the understanding is mainly qualitative in the form of experimental correlations. Theory has been limited to idealized situations and progress has been hampered by the difficulties involved in making detailed experimental measurements. 53 Based on current computing capacity, it is not possible to simulate the turbulence and the liquid interface during the entire atomization process via three-dimensional (3D) direct numerical simulation (DNS) methods. It is also very challenging to conduct large eddy simulation (LES) of the entire spray. Therefore, it is not surprising that most of previous efforts on modeling and numerical simulations of atomization have considered the flow to be laminar [13-23] and there are limited studies on two-fluid turbulent flows [24-31]. There are, however, several numerical (LES and DNS) studies which consider droplets’ dispersion, evaporation and combustion in the dilute spray region with point particle methods [32-45]. The main objective of this work is to investigate the two-way interactions between the turbulence and highly convoluted interface in a two-dimensional incompressible homogeneous isotropic flow by a recently developed two-fluid DNS method [46]. In this method, the interface is tracked by the hybrid Lagrangian-Eulerian interfacial particle level-set (IPLS) algorithm. The effects of density ratio, viscosity ratio and Weber number on the turbulence and interface are studied in details. 2.2 Results and Discussions In this section, the numerical data obtained by DNS for several different two- fluid isotropic turbulent flows are analyzed for detailed investigation of turbulence- interface interactions. The effects of density ratio and Weber number on the evolution of interface and turbulence are also studied. The parameters of simulations considered in this chapter are listed in Table 2-1. The reference Reynolds number, as defined in chapter 1, is 100 in these simulations and the gravity is omitted. 54 The initial vorticity field in both fluids is a fully developed 2D isotropic turbulent flow and the interfaces are two straight lines separating the fluid 2 from fluid 1 (Figure 2- 5(a)). Unless otherwise mentioned, the initial conditions for velocity and interface are the same in all cases listed in Table 2-1. Most simulations are carried out on uniform grids with a resolution of (256x256). In order to keep the turbulence decay rate in both fluids the same for different density ratio it, the viscosity ratio r] is adjusted so that the Reynolds numbers remain the same. In case 0, the flow is essentially the same as that in single-fluid isotropic decaying turbulent flows. Casesl-4 are considered to investigate the effects of density ratios, while cases 3, 5-7 are for studying the effects of Weber number or the surface tension. 2.2.1 Numerical Accuracy Before using the DNS data for investigation of turbulence-interface interactions, the numerical accuracy of the computational model has to be established. This includes finding the minimal interface thickness 6‘ for the given resolution. Figures (from Figure 2-1 to 2-4) show the results for case 3 for various interface thicknesses and grid resolutions. The results for other cases are similar and not shown. Despite the differences in interface thickness and resolution, Figure 2-1 shows that the instantaneous vorticity contours as strongly affected by the interface are virtually identical for all grid resolutions and interface thicknesse. The temporal variations of total interface area are also similar for different interface thicknesses and resolutions as shown in Figure 2-2. All the other turbulence and interface statistics such as the decay rate of mean turbulent kinetic energy in the interface zone and the turbulent energy spectrum are found to be the same for lower and higher resolutions and different 8. 55 To better assess the numerical accuracy of the model, the vorticity equation along the interface is considered: 0ch g , 1 Dr =l(“"V) 1,. leX(Vp/p<¢>)1,.+figth(V-(2p<¢>D)/p<¢>)1,, (7) (f1) ‘ v . 1 (III) (2.1) - W [V x (K(¢)VH(¢)/p(¢))lp (iii) The subscript p in equation (2.1) denotes the interface as marked by the interfacial particles and c?) = V x0 is the vorticity vector. The terms I, II, IH, and IV in the RHS of equation (2.1) are the vortex stretching, baroclinic vorticity generation, viscous and surface tension terms, respectively. Obviously, the term I is zero in two-dimensional flow. Terms II and IV are only important in the interface region and automatically vanish as one move away from the interface. By employing equations (1.3) and (1.8), the terms 11, IV can be written into different formats (labeled 11a and Na ) as shown below, (Ila) = [5(¢)(1—1/,1)V¢pr/p2(¢)Jp 1 (2.2) (W a > = —l6(¢)V¢> 6‘ The expressions for terms Ila and Na in equations (2.2) could be further simplified by using equations (1.3), (1.8) and (2.3), knowing that ¢=O at interface for any given interfacial particle, 56 14aa—0 (Hb)s- (vaxvp) 2 P iii+lét (24) IV E————V V ‘ b) £We(/i+l)( “X “P where terms IIb and IVb are inversely proportional to the interface thickness 8 . Mathematically, IIb and IVb are the same as H and IV. A noticeable advantage of IIb and IVb over II and IV is that the spatial derivative of density, which sharply changes in space, is avoided. This makes the numerical values of IIb and IVb to be much closer to the exact values. Here, we compare these terms to evaluate the numerical error caused by the derivative of density across interface. Figure 2-3 shows that the numerical values of terms 11 and IV are noisy and very different than those of IIb and IV b when interface thickness 8 is too small. These results are obtained for Lagrangian particle C on the interface (point C in Figure 2-1). The results for other points on the interface exhibit similar trends and are not shown here. As the interface thickness increases, the magnitudes of baroclinic and surface tension terms are decreased with almost the same rate and the temporal variations of terms 11 and IV become smooth and nearly the same as those of terms IIb and IVb, respectively. We found the numerical error associated with sharp density change across interface to be negligible when the interface thickness 8 is about 4Ax or bigger. Figure 2-4 shows that at a given interfacial (particle) location, the values of baroclinic and surface tension terms in equation (2.1) are roughly doubled when the resolution is increased from (256x256) to (512x512) while the ratio of interface thickness 8 to grid resolution Ax is kept constant. This indicates that the variations of 57 different variables along the interface’s normal direction, such as the curvature and pressure gradient are similar when properly scaled with interface thickness. Basedon the results discussed in this section, the resolution (256x256) and interface thickness of 4Ax are considered to be sufficient and are employed in all simulations. 2.2.2 Physical Structures of Turbulence and Interface Examination of flow structure, primarily in the region close to the interface, reveals several interesting features of turbulence-interface interactions and is considered in this section. The evolutions of vorticity and interface for case 3 are shown in Figure 2- 5. Evidently, the two interface planes (or lines) between fluids 1 and 2 experience significant topological changes due to turbulent stretching and compression by large scale vortical fluid motions in both fluids. As expected, the interactions between turbulence and interface are not important in the regions away from the interface, where the vorticity magnitudes in both fluids are shown to decrease monotonically in time, consistent with the decay rate of two-dimensional isotropic turbulence [47-48]. However, in the regions close to the interface, particularly in the low density (fluid 2) side, larnellar like vortical structures are generated with dimensions much larger in the tangential direction of the interface than the normal direction. The vorticity magnitude in this region is much higher than those in other flow regions. The instantaneous vorticity contours and the geometrical shapes of interfaces for different density ratio and Weber number are shown in Figure 2-6. Obviously, at locations away from interfaces, the vorticity contours are identical in all cases. However, in the vicinity of interfaces, the intensity of interface-generated vorticity, as measured by 58 vorticity magnitude is enhanced with the increase in density ratio (Figure 2-6(b)-(d). This is expected since the vorticity generation terms in the vorticity equation are proportional to density ratio. On the contrary, for the cases with the same density ratios, as shown in Figure 2-6(a), (b), (e) and (f), the magnitudes of the generated vorticity are decreased as the Weber number increases. As Figure 2-6 suggests, the influence of density ratio and Weber number on the interface is such that the total interface area/length is decreased as the density ratio increases and or the Weber number is decreased. Again, this is expected as the vorticity generation terms are directly dependent on the inverse of Weber number. In Figures 2-1, 2-5, and 2-6, several sample particles, labeled as particles A-E, are considered for studying the dynamics of interface and turbulence, which would be discussed in the next section. Also in next section, we will examine the turbulence structure and the geometrical variations of interface by considering several statistical quantities, such as the mean vorticity, the mean kinetic energy, and the interface area. 2.2.3 Statistical Analysis To make a clear distinction between various flow regions, the computational domain is divided here into five regions based on the density and distance to the interface. The mean value of a given variable in these regions is calculated by the followings relations. (a)l = (al ¢ > 26‘) (at/)2 = (al d < —28> .= <2-s> (Ct/>11 = (a| 0 < p S 28) (0')” =11 and( >12, respectively. As the results in Figure 2-5 and 2-6 suggest, one of the important flow variables that is significantly affected by the interface is the vorticity. Figure 2-7 shows the temporal variations of the average vorticity magnitude or enstrophy,£2 = (Ti-(DI 2. The mean enstrophy plots for case 3 in Figure 2-7(a) indicate that the vorticity generation is very significant in the interface region ( |¢>| s 25 ) when density ratio is significant. However, away from the interface (M > 25) in the fluid 1 and fluid 2 sides, the temporal variations of mean enstrophy as shown in Figure 2-7(a) and 2-7(b), are not affected by the interface and are the same for different density ratios and Weber numbers. Close examination of the results in Figure 2-7(a) indicates that the vorticity is generated mainly in the low-density side of the interface; this is due to the fact that the density ,0 appears as the denominator in the enstrophy generation terms (terms II and IV) in equation (2.6). The vorticity generation is more significant at higher density ratio and lower Weber number as shown in Figure 2-7(c)-(d). The results in Figure 2-7 are in agreement with those shown in section 2.2.2. To further explain the above results and to identify the mechanisms responsible for generation and destruction of mean enstrophy, the transport equation for mean enstrophy is considered below: 60 92 __ ~ ~ ~ 2 07 Dr —w-D-w+g)-(Vp(¢)pr)/p (¢3+Eg-v><(vtummy/2(a) (1) (f1) * (III) 4 (D (2.6) _ W: - (V K(¢) x V¢)5(¢) / PM) L (111) Terms 1, II, III, and IV in the RHS of equation (2.6) represent the “vortex stretching”, the “baroclinic”, the “viscous dissipation” and the “surface tension” effects on mean enstrophy, respectively. The vortex stretching term I is zero in two-dimensional flows. The contributions of the last three terms are discussed below. Away from the interface, terms 11 and IV are negligible in both fluid 1 and fluid 2 regions, and the mean enstrophy are monotonically decreasing in these regions (Figure 2— 7(b)) due to negative term III (Figure 2-8(a)). Figure 2-8(a) also shows that the magnitude of mean viscous dissipation term in the interface region (111) is larger than I those in the regions away from the interface ((III)land ([102). Figure 2-8(b) indicates that the contribution of mean baroclinic term (“)1 to the mean enstrophy in the interface region proportionally increases with the density ratio. The average contribution of surface tension term in the interface region (IV) is small and negative at the beginning I but becomes much larger and positive later on as shown in Figure 2-8(c). In contrast to baroclinic and surface tension terms, Figure 2—8(d) shows that the sign of mean viscous dissipation term is always negative while its magnitude always increases with the density ratio. Another observation is that the evolutions of baroclinic and surface tension terms are contrary to each other, i.e. the baroclinic term increases early on and decreases later on but the surface tension term decreases at the beginning then increases. Also the 61 positive peak of baroclinic term at early times is accompanied with the negative peak of surface tension term at nearly the same time. This indicates that in the interface region, the surface tension opposes the baroclinic effects. Since the vorticity is mainly generated in the interface region by the baroclinic and surface tension terms, the physical scale of the generated vorticity field, at least in the direction normal to the interface, is smaller than those already exist in fluids 1 and 2 sides. This is not only observed in the vorticity contour plots in Figure 2-1, 2-5 and 2-6, but is also evident in Figure 2—9, where it is shown that the effects of interface on the turbulence kinetic energy and enstrophy spectrums primarily appear at large wavenumbers or small scales. The total surface area (or length) of the interface is an important variable in quantification of turbulent effects on the interface. The results in Figure 2-10 are consistent with the geometrical shapes of interface shown in Figure 2-6, indicating that the total interface area, A has a tendency to increase in time as the density ratio decreases or the Weber number increases. However, for smaller Weber number, A either increases slowly or does not increase at all. At long time, A tends to not increase so much or even decrease. Another important turbulent quantity that is affected by the interface is the turbulent kinetic energy, ek = LII-u,- / 2. The temporal variations of mean turbulent kinetic energy in different regions of the flow are shown in Figure 2-11. The results in this figure indicate that the decay rates of mean turbulent kinetic energy in both fluids become almost identical with those in single-fluid isotropic homogeneous turbulence 62 [47-48]. However, the mean turbulent kinetic energy in the interface region is decaying faster than those in fluid 1 and 2 sides. Figure 2-12 also demonstrates that the decay rate of mean turbulent kinetic energy in the interface region is varied with the density ratio and Weber number. For all cases considered in this figure, the initial decay rate of mean turbulent kinetic energy is significant. At later time the decay rate decreases for the case 4 with higher density ratio and the kinetic energy may even increase (Figure 2-l3(a)). Figure 2-12(b) shows that for case 5 with smaller Weber number, the mean turbulent kinetic energy decays faster at early time in comparison with other cases, but it exhibits irregular and oscillatory behavior at later times. Comparing the evolution of kinetic energy in Figure 2-12 with that of the interface area in Figure 2-10, one can conclude that with the decrease in mean turbulent kinetic energy, the interface area increases or vice versa. This implies an inverse relationship between the interface area and mean turbulent kinetic energy in the interface region. To understand the mechanisms responsible for this relationship, the transport equation of the instantaneous turbulent kinetic energy is considered here: Ifl _V-(pl7) V-(2#(¢)17-D) _2#(¢)D:D _6(¢)x(¢) = + —— V . (git?) (2.7) Dr p<¢> \ Re pm) 1 , Re pm) J . Weprm g (1) (11) (1h) (1W Terms 1, II, III, and IV in the RHS of equation (2.7) represent the pressure work, the energy transport by viscous stress, the viscous dissipation and the surface tension work, respectively. The average values of the above four terms in different regions of the flow for case 3 are shown in Figure 2-13. Obviously, the viscous dissipation is always negative and the work done by the surface tension force (term IV) is zero away from the interface. 63 The results in Figure 2-13 indicate that away from the interface, in fluid 1 and 2 sides, the viscous dissipation term is the most significant term and the contributions of the other terms are negligible. However, the mean turbulent kinetic energy in the interface region is affected by all terms. The viscous transport term seems to be the least important term in comparison with other terms. Also at early times, the magnitudes of pressure and surface tension terms are larger than that of viscous dissipation term. The positive contribution of term I implies that the turbulent kinetic energy is increased by pressure work near to interface. In turn, the surface tension work is negative at early times and tends to cancel the effects of pressure work. Later on, the surface tension work turns positive, while the pressure work is nearly zero. At these times, the positive effect of surface tension on turbulent kinetic energy is cancelled by the viscous dissipation term. The early (rapid) decay of the mean turbulent kinetic energy in the interface region is shown in Figure 2-11 to be mainly due to negative contribution of surface tension and to a lesser extend the viscous dissipation terms. The later one is somewhat reduced by the positive pressure work. Figures 2-11 and 2-13 also show that the decay rate of the turbulent kinetic energy and the contributing terms to its evolution are nearly the same in fluids 1 and 2 away from the interface. The effects of density ratio and Weber number on various terms in equation (2.7) in the interface region are shown in Figures 2-14 and 2- 15, respectively. Figure 2-14 indicates that the magnitudes of surface tension, pressure and viscous dissipation terms all increase with the density ratio and the contribution of viscous transport term remains smaller than those of the other three terms as expected. Figure 2- 14 also shows that, for the case with large density ratio (i.e. case 4), the surface tension 64 work will turn from a large negative value to a large positive value, while the pressure work remains to be positive or close to zero. At early times, when the contributions of the surface tension and viscous dissipation terms are both negative, the magnitude of viscous dissipation term increase slowly because both terms are balanced by the pressure work. However, once the surface tension work turns positive, the magnitude of viscous dissipation term increases sharply to take over the role of pressure to counteract the effect of surface tension in the interface region. , Figure 2-15 shows that the magnitudes of surface tension and pressure works are both increased by decreasing the Weber number, or by increasing the relative effect of surface tension force, as expected. While at low Weber number, there is a significant variation and oscillation in the mean turbulent kinetic energy and the terms contributing to its evolution (Figure 2-15) in the interface region, they seem to vary much less at higher Weber numbers. The oscillations in surface tension term cause the oscillations in the mean turbulent kinetic energy and other terms in its transport equation. Nevertheless, the magnitude of surface tension work is higher than that of pressure work, while they seem to be inversely proportional. For cases with smaller Weber number or bigger surface tension, the magnitude of viscous dissipation term increases rapidly at early times then becomes very small later on. The magnitude of viscous transport term is somewhat correlated but is far smaller than that of pressure work. The faster decay of mean turbulent kinetic energy for the cases with smaller Weber number (Figure 2-12(b)) is caused by the large surface tension work and viscous dissipation which are both negative at early times. 65 To further better understand the effects of surface tension, the temporal variations of turbulent kinetic energy along the interface are considered in Figure 2-16. These results are obtained by averaging the turbulent kinetic energy and its components normal and tangential to the interface over all particles as described in the following equations: NP _. _. ep = 2m, -vp /2)/NP n=l NP em = Z((V,,)P -(V,,)p /2)/NP (2.8) n=l NP ep, = Z((V,)P -(V,)p/2)INP, n=l where 17,), (V,,) and (V,)p are the fluid velocity, the normal component of the fluid [2 velocity, and the tangential component of the fluid velocity at particle p, respectively. NP is the total number of interfacial particles. Figure 2-16(a) shows that in case 1 the normal component of turbulent kinetic energy decays faster than its tangential component, where the surface tension work is negative at all times. However, for case 3, the results in Figure 2-l6(b) indicate that the normal component of energy increases at later time when the net surface tension work at interface becomes positive (Figure 2-l3(a)). It is shown in Figures 2—l6(c) and (d) that at early time the normal component of turbulent kinetic energy decays faster than the tangential component and then oscillates at later times consistent with the surface tension oscillation shown in Figure 2-15(a). Nevertheless, in all cases the tangential component of kinetic energy seems to be less sensitive to the surface tension work in comparison to the normal component. To explain why the normal component of the turbulent kinetic energy is more sensitive to surface tension than the tangential component, the transport equation for 6,, is 66 considered below: 2c; : u,,0 .213 _U” (Viral), pm» + \7(2710201?!~D)-2u(¢)D:V/II Dr Dr M) M) Re p<¢> _ 6<¢>Kr¢>v - (a?) (29> Wep(¢) ’ (v) where (2,, = U ,, -U ,, / 2 is the instantaneous normal component of turbulent kinetic energy at Eulerian grid point and U” = U - N is the normal velocity component. IV = V¢/IV¢| is the unit normal vector at interface. Note that, term V in equation (2.9) is the same as term IV in equation (2.7), which indicates that the surface tension work only changes the normal component of the instantaneous and average turbulent kinetic energy in the interface region. This confirms that only the normal component of kinetic energy can be directly affected by the surface tension in the interface region. It has been stated before that there is a direct relation between the total interface area and the mean turbulent kinetic energy in the interface region. This is confirmed in Figure 2-17 and 2-18, where the temporal variations of the rated interface area change 2 and the mean surface tension work change (3) in the interface region are plotted for different cases. The variables 2, O are normalized and are defined as: t=z/ , MM” (2.10) 9:949le , where 12% and®=<—i:/¢—I)0K(g:lv-(¢l7)> . e I 67 Figures 2-17 and 2-18 clearly show that the rated interface area change is inversely proportional to the mean surface tension work in the interface region, regardless of density ratio and Weber number. In other words, when the value of mean surface tension work is positive, the rate of interface area change is negative and vice versa. Also, when ,7? increases, 9 decreases with the same rate and vice versa. Figure 2- 17 also shows that at higher density ratio, the rate of interface area change rapidly decline shortly after its early rapid increase. This explains why the total interface area in the case with higher density ratio is generally smaller than that in the case with lower density ratio (Figure 2-10(a)). Figure 2-18 shows that with an increase in Weber number, the mean surface tension work oscillates less and mostly attains negative values, while the rate of interface area change decreases slowly and stays mainly positive. This explains the larger total interface area values in cases with higher Weber number (Figure 2-10(b)). The destabilization effect of turbulence and the stability effect of surface tension on the interface as observed in Figures 2-6 and 2-10 can be explained through the turbulent kinetic energy, primarily its normal component. For this explanation, we consider the interface evolution at various stages. At the first stage, the flat interface is stretched rapidly by the nearby turbulent strain field, resulting in an early rapid and significant increase in total interface area (Figure 2-10, 2-17 and 2-18). This is accompanied with a loss in turbulent kinetic energy, especiallyiits normal component, primarily due to the surface tension which opposes the increase in interface area. The lower the Weber number and the higher the density ratio, the more the turbulent kinetic energy of the surrounding flow is damped (Figures 2-12, 2-14(a), 2—15(a)). At the second 68 stage of interface evolution, turbulent stretching effect is not as strong as that in the first early stage since the surrounding turbulence loses most of its kinetic energy, particularly its normal component. Although the total interface area continues to increase, the rate of interface area change is considerably lower in the second stage. In this stage, the turbulent kinetic energy’s loss due to surface tension is decreasing as shown in Figures 2- l7 and 2-18. For the cases with higher density ratio or lower Weber number, the total interface area decreases at later time. The main reason for the decrease in interface area is that the strain field or turbulent kinetic (especially its normal component) is becoming too weak to overcome the strong surface tension force generated by large curvature. During the time that the interface area decreases, the surface tension is so significant that may even increase the turbulent kinetic energy in the interface region as Figures 2—12, 2- 17 and 2-18 indicate. By increasing the kinetic energy of the surrounding turbulent flow, the surface tension counteracts the further decrease in interface area, such that it again makes the interface more stable. Figure 2-19 shows that in the interface region, the normalized value of viscous dissipation )7 in the mean turbulent kinetic energy equation is proportional to that of mean enstrophy—£3. These normalized quantities are indicators of the turbulence or strain field at small scales and are defined as: i7 = 747ij a 2 I911 /|Q|I,Max > . In the interface region, the increase in enstrophy is I (2.11) _ 2#(¢)D : D where y = < Re i0(¢) primarily due to small-scale turbulent motions generated by the baroclinic and surface tension terms (Figures 2-7 to 2-9). 69 For better understanding of the turbulence/interface dynamics, the vorticity equation for several interfacial particles (described in equation (2.1)) is analyzed. Two representative interfacial particles, i.e. particle A and Particle B in Figure 2-6 are selected. The behaviors of terms 11, III, and IV in equation (2.1) for these two particles are shown in Figures 2-20 and 2-21. Note that the initial locations of the particles are the same in different cases. Figures 2-20 and 2-21 indicate that in the absence of baroclinic term in case 1, the viscous term counteracts the surface tension effect with a comparable strength. As density ratio increases, the magnitudes of surface tension and baroclinic terms are both increased while the magnitude of viscous term is decreased. The contributions of surface tension term and baroclinic term seem to be comparable but opposite to each other. With the decrease in Weber number (or increase in surface tension force), the magnitudes of both surface tension and baroclinic terms increase, as expected. The contribution of viscous term is negligible in the vorticity equation in comparison with those of baroclinic and surface tension terms (Figures 2-20(b) and 2- 21(b)). 2.3 Summary and Conclusions The interactions between turbulence and interface are investigated via data generated by direct numerical simulations of two-dimensional, two-fluid incompressible and immiscible isotropic turbulent flow. The interface motion is simulated with a recently developed hybrid Lagrangian-Eulerian interfacial particle level-set (IPLS) method, where the interface is explicitly represented by the positions of massless particles and the geometrical information concerning the interface is provided by the 70 level-set function (D which is kept as a signed distance function with the particle-based re-initialization process. The destabilization effect of turbulence on the interface is shown to be primarily due to stretching by large-scale vortical motions which increase the total interface area depending on the density ratio and Weber number. The deformed interface not only creates a complicated surface tension force distribution, but also misaligns the pressure and density gradient vectors, causing the baroclinic effect to enhance. The baroclinic term counteracts the influence of surface tension in both momentum and vorticity equations. The generated small-scale high-intensity turbulence by the baroclinic and/or surface tension terms causes the strain rate and viscous dissipation of turbulent kinetic energy in the interface region, especially in the low density side, to become much larger than those away from the interface region. Away from the interface, the turbulence behaves similar to single—fluid isotropic turbulent flows. The surface tension force opposes the increase in interface area by the neighboring turbulent strain field. The stability effect of surface tension on the interface is mostly associated with intense damping of neighboring turbulent flow’s kinetic energy, primarily its normal component. However, in some cases, the surface tension force acts oppositely and increases the kinetic energy of neighboring turbulent flow to oppose the decrease in interface area. The rate of interface area change is found to be inversely proportional to the mean work done by the surface tension force in the interface region. Examinations of the density ratio and Weber number effects on the turbulence and interface indicate that with an increase in density ratio or a decrease in Weber 71 number, the total interface area and the mean turbulent kinetic energy in the interface region are both decreased, while the mean enstrophy in the interface region is increased. In the kinetic energy equation, the magnitudes of pressure work, surface tension work and viscous dissipation of turbulent kinetic energy increase as the density ratio increases or the Weber number decreases. The magnitudes of baroclinic, surface tension and viscous dissipation terms in the mean enstrophy equation are increased by an increase in density ratio. However, the magnitude of viscous term in the vorticity equation decreases as density ratio increases. At early stage of interface evolution, as the interface is stretched rapidly by the nearby turbulent strain field, the turbulent kinetic energy, especially its normal component in the vicinity of interface decreases. This is primarily due to surface tension which opposes the increase in interface area. Later on when turbulent stretching become weaker, the rate of increase in interface area decreases and even in some cases the interface area starts to decrease. At this time, the surface tension starts to increase the turbulent energy in the nearby region of interface. By increasing the kinetic energy of the surrounding turbulent flow, the surface tension effectively opposes the decrease in surface area of the interface. Close to interface, when the contribution of surface tension and viscous dissipation terms in the turbulent kinetic energy are both negative, the viscous dissipation remains moderate because both terms are balanced by the pressure term. However, once the surface tension work turns positive, the magnitude of viscous dissipation term increases sharply to take over the role of pressure to counteract the effect of surface tension in the interface region. 72 2.4 Table Table 2-1. The Specifications of DNS Cases Case No. A T] We 0 1 1 00 1 1 1 10 2 5 5 10 3 10 10 10 4 20 20 10 5 10 10 2 6 10 10 5 7 10 10 2O 73 2.5 Figures é . / , ‘Fkfi is C (QmeH ii I) ‘\If(\\\\::‘ ’ :t Figure 2-1. Instantaneous vorticity contours for case 3 at t=2.0 for various grid resolutions and interface thicknesses. (a) a = 2Ax, 256x256 grid resolution; (b)£ = 3Ax, 256x256 grid resolution; (0) e = 4Ax, 256x256 grid resolution; and (d) 6‘ = 4Ax, 512x512 grid resolution. 74 18 Ifj] 17:- 16} A; 15: e=2Ax, 256 x256 : _._._._ e=3Ax, 256 x256 14 _ _.._..._ g=4Ax, 256 x256 _ —--.-°- €=4Axs 512 X512 13. l 12’s...r....r....r....r.... O 0.5 1 t 1.5 2 2.5 Figure 2-2. Temporal variations of the total interface surface areas for various interface thicknesses and various resolutions. 75 0 0.5 1 t 1.5 2 2.5 0 0.5 1 t 1.5 2 2.5 0"‘\ IIb \ — II Figure 2-3. Comparison between different forms of baroclinic and surface tension terms in vorticity equation for various interface thicknesses: (a)-(b)e = 2Ax; (c)-(d) a = 3Ax; and (e)-(f) e = 4Ax. The results are plotted for point C or particle C on the interface. Grid resolution is 256x256. 76 16 (a) . b) o p 11 e=4AX,256x256 i IV 8=4AX.256><256 ( ---——--— 11/2 e=4AX, 512x512 12» 1W2 8=4AX.512><512 - l r ‘ ./ _\ . q 11" [I .‘ -41. 8: A“ i“ t i ," _ “ i ‘ I 4r '2 r’ i. 8 44 1 H I r r l g l c , A 1 l l . l I x x 0 05 1 t 1.5 2 25 0 05 1 t 15 2 25 Figure 2-4. Temporal variations of baroclinic and surface tension terms for different grid resolutions. 77 Figure 2-5. Vorticity contours and Interface topologies in case 3 at different times: (a) t=0.0; (b) t=0.5; (c) t=1.5; and (d) t=2.5. 78 Figure 2-6. Instantaneous vorticity contours and interface topologies at t=2.0 for different cases: (a) case 0; (b) case 1; (0) case 3; (d) case 4; (e) case 5; and (f) case 7. 79 1O _ — (9)1 —-O— Caseo 8— ----2 2h Case1 : -...e...- (g2)i ; ----- oasez _ _6_ <9)“ % -......-.._.... Cases 6: —O-— (9)12 1.5:- :1 . (9)1 I 1 _. 0.5 6. (e) 1°» (d) —o—Caseo : —O——-Cese5 Figure 2-7. Temporal variations of mean enstrophy for various cases and in different regions of the flow: (a) case 3, fluids 1 and 2 and interface regions; (b) various cases, fluid 1 region; (0) and (d) various cases, interface region. 80 ,. (bl ,"' l. ’ \ " I ‘\ —-o—Case1 5_ ,I \ ——Case2 P, \\ ----- 08303 —-- (111)] - - - - (III): Figure 2-8. Temporal variations of different terms in the mean enstrophy equation in different regions for various cases: (a) viscous dissipation term III in different regions in case 1; (b) baroclinic term II in interface region; (c) surface tension term IV in interface region; and (d) viscous dissipation term III in interface region. 81 113-01 113-00 1E-03 lE-OZ l E-OS lE-04 lE-O7 15-09- ----- Gas“ * £2(K) ..... Case1 ’ ——0— Case2 lE-08 F —0— Case2 IE-ll_ --o-— Case3 —e—Case4 --o--Case3 15-13- IE-IO— —9— ““4 lE-lS . .....2. - ... . 1E-12 . .2. . ... . 1E+00 lE+01 K 1E+02 1E+00 1E+01 K 1E+02 Figure 2-9. Two-dimensional spectral density function of (a) turbulent kinetic energy and (b) enstrophy at t=2.0. 25 25 (a) j (b) . —e— Case1 - —O-— Case2 A I A A A A O 0.5 1 t 1.5 2 2.5 0 0-5 1 11-5 2 2-5 Figure 2-10. Temporal variations of the total interface areas for different cases: (a) various density ratios and (b) various Weber numbers. 82 0.18 (€101 - - - - 2 0.15, 0.12’ ' I 0.09} 0.06 '— 0.03 :- I I I A I O 0.5 1 t 1.5 2 2.5 )— O I L A l l l I A l I Figure 2-11. Temporal variations of mean turbulent kinetic energy in various regions for (a) case 1 and (b) case 3. - (a) (b) 0-1 _ Case 1,),=1 0'1 Case 5, We=2 _ ----- Case 3.}.=10 ; ----- Case 6, We=5 . —0— Case 4.A=20 - Case 7, We=20 0.12; 0-12. (€101: (3101: 0.08} 0.08: F _ 0.04- 0.04_ 0* 2 2 *4 .2 2 0’ ‘ t ' O 0.5 1 t 1.5 2 2.5 0 0.5 1 t 1.5 2 2.5 Figure 2-12. Temporal variations of mean turbulent kinetic energy in interface region for different cases: (a) various density ratios and (b) various Weber numbers. 83 0.8_ 0.8 (a) (b) 0.6 0.6 YY'II 0.4: Ifl'TjF'rlj 0.2: -0.2__ : —— 1 '04 ~ - - - ‘ <1V>2 I —0— (IV)[ -O.6 _- E —— (1101 .016; ---- <111>2 0.02 ____ (“>2 ' ; + (1101 : —o— <11>1 _O.2P. .11....L....1.2.._l....lIO.O4 1.11.1.1....l....l.2..l 0 0.5 1 t 1.5 2 2.5 o 05 1 t 1,5 2 25 Figure 2-13. Temporal variations of different terms in the mean turbulent kinetic energy equation in various regions for case 3: (a) surface tension work; (b) pressure work; (0) viscous dissipation; and (d) viscous stress transport terms. 84 j (a) .0 ts. ' r (1)1) 02 1 ~ Case 1 ----- Case 2 -0 4 .. —0— Case 4 12 2 J I l l 2 I l l l -0.08 (IIDI: 0.12: 016' )- r- 02 LALL ‘AlLlAALA AAAA AAlJAlJl l4l‘l - 0 O 0.5 1 t 1.5 2 2.5 O 0.5 1 t 1.5 2 2.5 Figure 2-14. Temporal variations of different terms in the mean turbulent kinetic energy equation in interface region for cases with different density ratios: (a) surface tension work; (b) pressure work; (c) viscous dissipation; and (d) viscous stress transport terms. 85 ca“ 5 0.04 L -O.16 - ..... Cases Cases —0— Case7 i """ 03995 : ~ —0— Case7 -02 4‘ ‘ ‘ ' 2 14 - L. A . .I_22_. . I . . . '0-06 . . . . I . .2 . I 2. . . I. . . . I . . . . 0 0.5 1 t 1.5 2 2.5 0 0.5 1 t 1.5 2 2.5 Figure 2-15. Temporal variations of different terms in the mean turbulent kinetic energy equation in interface region for cases with different Weber numbers: (a) surface tension work; (b) pressure work; (0) viscous dissipation; and (d) viscous stress transport terms. 86 Figure 2-16. Temporal variations of interfacial kinetic energy and its two- components in different cases: (a) case 1; (b) case 3; (c) case 6; and (d) case 5. 87 Figure 2-17. Temporal variations of normalized interface area change rate and mean surface tension work in interface regions with various density ratios: (a) case 1; (b) case 2; (c) case 3; and (d) case 4. 88 Figure 2-18. Temporal variations of normalized interface area change rate and mean surface tension work in interface regions with various Weber numbers: (a) case 5; (b) case 6; (c) case 3; and (d) case 7. 89 12_ L2 . (a) ~ (b) - _\ : [I \ 0.8— \ O.8~ I ‘ . \ _ I \ : ‘~ t , ‘ O4,- \“‘ 1: III 1 ' F. ......... 0-4- __,—-" 03 0- Figure 2-19. Temporal variations of normalized viscous dissipation term of mean turbulent kinetic energy and normalized mean enstrophy in interface region: (a) case 1; (b) case 4; (0) case 5; and (d) case 6. 9O lllALLALJJ—AQ¥JIAALJL1‘1 h-AAA.llAAAr.lLLJAIA‘AA'rAAA O 0.5 1 t1.5 2 2.5 O 0.5 1 11.5 2 2.5 Figure 2-20. Temporal variations of different terms in the vorticity equation on particle A for different cases: (a) case 1; (b) case 3; (0) case 5; and (d) case 7. 91 2.5 0.5 1 t 1.5 2 2.5 Figure 2-21. Temporal variations of different terms in the vorticity equation on particle B for different cases. (a) case 1; (b) case 3; (0) case 5; and (d) case 7. 92 2.6 [5] [6] [7] [8] [9] [10] [111 [12] [13] References Lefebrvre, A. H., Atomization and Sprays Hemisphere Publishing, New York, 1989. Rayleigh, L., “On the instability of jets,” Proc. London Math. 800., Vol. 10, 1878, pp.4-12. Meiron, D. l., Baker, G. R. and Orszag, S. A., “Analytic structure of vortex sheet dynamics. Part 1. Kelvin-Helmholtz instability,” Journal of Fluid Mechanics, Vol. 114, 1982, pp.283-298. Atsavapranee, P. and Gharib, M., “Structures in stratified plane mixing layers and the effects of cross-shear,” Journal of Fluid Mechanics, Vol. 342, 1997, pp.53-86. Faeth, G. M., “Structure and atomization properties of dense turbulent sprays,” Twenty- Third Symposium (International) on Combustion, Orleans, France, July 1990, pp.1345-1352. Bellan, J., “Perspectives on Large Eddy Simulation for sprays: Issues and solutions,” Atomization and Sprays, Vol. 10, 2000, pp.409-425. Komori, S., Nagaosa, R. and Murakami, Y. et al., “Direct numerical simulation of three-dimensional open-channel flow with zero-shear gas- quuid interface,” Physics of Fluids A, Vol. 5, 1993, pp.115-125. Reitz, R. D. and Bracco, F. V., “Mechanisms of atomization of a liquid jet,” Physics of Fluids, Vol. 25, 1982, pp. 1 730-1 742. Sirignano, W. A., Fluid dynamics and transport of droplets and spray, Cambridge University Press, 1999. Lin, S. P. and Reitz, R. D., “Droplet and spray formation from a liquid jet,” Annual Review of Fluid Mechanics, Vol. 30, 1998, pp.85-105. Lasheras, J. C. and Hopfinger, E. J., “Liquid jet instability and atomization in a coaxial gas stream,” Annual Review of Fluid Mechanics, Vol. 32, 2000, pp.275-308. Crowe, C. T., Troutt, T. R. and Chung, J. N., “Numerical models for two- phase turbulent flows,” Annual Review of Fluid Mechanics, Vol. 28, 1996, pp.11-43. Ashgriz, N. and P00, J. Y., “FLAIR: flux line-segment advection and interface reconstruction,” Journal of Computational Physics, Vol. 93, 1991, pp.449-468. 93 [14] [15] [151 [17] [18] [19] [20] [21] [22] [23] [24] Mashayek, F. and Ashgriz, N., “Advection of axisymmetric interfaces by volume of fluid method,” International Journal for Numerical Methods in Fluids, Vol. 20, 1995, 901337-1361. Pilliod, J. E. and Puckett, E. G., “Second-order accurate volume-of—fluid algorithms for tracking materials interface,” Journal of Computational Physics, Vol.199, 2004, pp.465-502. Unverdi, S. O. and Tryggvason, G., “A front-tracking method for viscous, incompressible, multi-fluid flows,” Journal of Computational Physics, Vol.100, 1992, pp.25-37. Shyy, W., Udaykumar, H. S., Rao, M. M. and Smth, R. W., Computational Fluid Dynamics with Moving Boundaries, Taylor 8. Francis, Washington, 1996. Anderson, D. M., McFadden, G. B. and Wheeler, A. A., “Diffuse-interface methods in fluid mechanics,” Annual Review of Fluid Mechanics, Vol. 30, 1998, Pp.139-165. Sussman, M., Smereka, P. and Osher, S., “A level set approach for computing solution to incompressible two-phase flow,” Journal of Computational Physics, Vol. 1 14, 1994, pp.146—159. Chang, Y. C., Hou, T. Y., Merriman, B., and Osher, S., “A level set formulation of Eulerian interface capturing methods for incompressible fluid flows,” Journal of Computational Physics, Vol. 124, 1996, pp.449- 464. Enright, D., Fedkiw, R. Ferziger, J. and Mitchell, l., “A hybrid particle level set method for improved interface capturing,” Journal of Computational Physics, Vol. 183, 2002, pp.83-116. Scardovelli, R. and Zaleski, 8., “Direct numerical simulation of free-surface and interfacial flow,” Annual Review of Fluid Mechanics, Vol. 31, 1999, pp.567-603. Jafari, A., Shirani, E. and Ashgriz, N., “An improved three-dimensional model for interface pressure calculations in free-surface flows,” International Journal of Computational Fluid Dynamics, Vol. 21, 2007, pp.87-97. Banerjee, S., “Modeling considerations for turbulent multi-phase flows,” in Engineering Turbulence Modeling and Experiments, Rodi, W. and Ganic, E. N. (eds). Elsevier: New York, 1990, pp.831-866. 94 [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [35] Lopez de Bertodano, M., Lahey, R. T. and Jones, O. C., “Development of a k-c model for bubbly two-phase flow,” Journal of Fluid Engineering, Vol. 116, 1994, pp.128-134. Hong, W. L. and Walker, D. T., “Reynolds averaged equation for free- surface flows with application to High-Froude-Number jet spreading,” Journal of Fluid Mechanics, Vol. 417, 2000, pp.183-209. Alaibegovic, A., “Large eddy simulation formalism applied to multiphase flows,” Proceedings ofASME FEDSM’01, New Orleans, LA, 2001. Herrmann, M., “Modeling primary breakup: A three-dimensional Eulerian level set/vortex sheet method for two-phase interface dynamics,” Annual Research Briefs, Center for turbulence research: Stanford, 2003, pp.185- 195. Shirani, E., Jafari, A. And Ashgriz, N., “Turbulence models for flows with free surfaces and interfaces,” AIAA paper 2004-1283, Jan. 2004. Jafari, A. and Ashgriz, N., “Turbulence generated primary breakup: spreading of a liquid jet,” ILASS America, 17‘” Annual Conference on Liquid Atomization and Spray Systems, Arlington, VA, May 2004. Carrara, M. D. and DesJardin, P. E., “A filtered mass density function approach for modeling separated two-phase flows for LES I: mathematics formulation,” International Journal of Multiphase Flow, Vol. 32, 2006, pp.365-383. Crowe, C. T., “Review - numerical models for dilute gas-particle flows,” Journal of Fluid Engineering, Vol. 104, 1982, pp.297-303. Crowe, C. T., Sommerfeld, M. and Tsuji, Y., Multiphase Flows with Droplets and Particles, CRC Press, 1998. Jaberi, F. A., “Temperature fluctuations in particle-laden homogeneous turbulent flows,” International Journal of Heat and Mass Transfer, Vol. 41, 1998, P114081 -4093. Mashayek, F. and Jaberi, F. A., “Particle dispersion in forced isotropic low Mach number turbulence,” International Journal of Heat and Mass Transfer, Vol. 42, 1999, pp.2823-2836. Miller, R. S. and Bellan, J., “Direct numerical simulation and subgrid analysis of a transitional droplet laden mixing layers,” Physics of Fluids, Vol. 12, 2000. 90650-671. 95 [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] Mashayek, F., “Numerical investigation of reacting droplets in homogeneous shear turbulence,” Journal of Fluid Mechanics, Vol. 405, 2000, DIM-36. Loth, E., “Numerical approaches for motion of dispersed particles, droplets and bubbles,” Progress Energy Combustion Science, Vol. 26, 2000, pp.161-223. Ling, W., Chung, J. N. and Crowe, C. T., “Direction numerical Simulation of a two-way thermally coupled droplet-laden mixing layer,” Journal of Fluid Mechanics, Vol. 437, 2001, pp.45-68. Apte, S. V., Gorokhovski, M. and Moin, P., “LES of atomizing spray with stochastic modeling of secondary breakup,” International Journal of Multiphase Flow, Vol. 29, 2003, pp.1503-1522. Salman, H. and Soteriou, M., “Lagrangian simulation of evaporating droplet sprays,” Physics of Fluids, Vol. 16, 2004, pp.4601-4622. Menon, S. and Patel, N., “Subgrid modeling for simulation of spray combustion in large-scale combustors,” AIAA Journal, Vol. 44, 2006, pp.709-723. Mahesh, K., Constantinescu, G., Apte, S., laccarino, G., Ham, F., and Moin, P., “Large-eddy simulation of reacting turbulent flows in complex geometries,” Journal of Applied Mechanics, Vol. 73, 2006, pp.374-381. Almeida, T. and Jaberi, F. A., “Direct numerical simulations of a planar jet laden with evaporating droplets,” International Journal of Heat and Mass Transfer, Vol. 49, 2006, pp.2113-2123. Almeida, T. and Jaberi, F. A., “Large eddy simulation of a particle-laden turbulent jet flow with a stochastic two-phase subgrid closure,” International Journal of Heat and Mass Transfer, in press, 2007. Li, Z., Jaberi, F. A. and Shih, T. l.-P., “A hybrid Lagrangian-Eulerian particle-level set method for numerical simulation of two-fluid turbulent flows,” International Journal for Numerical Methods in Fluids, Vol. 56, 2008, pp.2271-2300. Herring, J. R., Orszag, S. A., Kraichnan, R. H. and Fox, D. G., “Decay of two-dimensional homogeneous turbulence,” Journal of Fluid Mechanics, Vol. 66, 1974, pp.417-444. Brachet, M. E., Meneguzzi, M., Politano, H. and Sulem, P. L., “T he dynamics of freely decaying two dimensional turbulence,” Journal of Fluid Mechanics, Vol. 194, 1988, pp.333-349. 96 [49] [50] [511 [52] [53] [541 Almgren, A. 8., Bell, J. B. and Szymczak, W. G., “A numerical method for the incompressible Navier-Stokes equations based on an approximate projection,” SIAM Journal of Scientific computing, Vol. 17, 1996, pp.358- 369. Sussman, M., Almgren, A. S. and Bell, J. B, etc, “An adaptive level set approach for incompressible two-phase flows,” Journal of Computational Physics, Vol. 148, 1999, pp.81-124. Kim, J. and Moin, P., “Application of a fractional-step method to incompressible Navier-Stokes equations,” Journal of Computational Physics, Vol. 59, 1985, pp.308-323. Lele, S. K., “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, Vol. 103, 1992, pp. 1 6-42. Osher, S. and Fedkiw, R., Level set methods and dynamic implicit surfaces, Applied Mathematical Sciences, Springer, 2002. Tatebe, O., ”The Multigrid preconditioned conjugate gradient method,” In 6’“ Copper Mountain Conference on Multi-grid Methods, Copper Mountain, CO. April 4-9, 1993. 97 CHAPTER 3 LES/FMDF Methodology for Two-Phase Turbulent Reacting Flows 3.1 Introduction As stated before, turbulent two-phase reacting flows are important to a large variety of applications including spray combustion, pollutant dispersion in atmosphere, heat and mass transfer in chemical reactors, material processing, among many others. This chapter is focused on the LES of dilute two—phase turbulent reacting flows. The extent of research and the rate of progress in the field of modeling and large-scale simulations of two-phase turbulent flows have been significant. In fact, the span of the developments has been much more than can be described, with justice, in a short introduction. Therefore, this section is focused on selected issues concerning the droplet/particle transport in turbulent flows and dilute spray combustion. Some of the challenges associated with LES of multiphase flow are also discussed. This discussion is not exhaustive as the physical and the mathematical complexity of multiphase transport do not allow a “complete” coverage. All three simulation strategies, namely direct numerical simulation (DNS), large- eddy simulation (LES), and Reynolds-averaged simulation (RAS), have been widely employed in computational analysis of turbulent two-phase flows. Despite all of the progress in supercomputer technology, and even considering the (most optimistic) projected rate of progress of this technology, DNS will remain primarily as a basic research tool. This should not imply a “shortcoming” of the method. On the contrary, 98 DN S has proven very effective in capturing many physical phenomena which could have not been studied by other means. It has also been very useful in assessing the performances of many of the closures in RAS and LES. RAS has been the method of choice for practical applications and perhaps this trend will continue for many years to come. This is the case even with proven non-universality of the currently employed RAS closures for two-phase turbulent (non) reacting flows. LES has experienced significant growth in two-phase flow simulations within the past decade or so. This methodology makes use of many of the pleasing features of both RAS and DNS. It is more reliable than the former, and is less computationally intensive than the latter. The predictive capability of LES depends on the success of its SGS closures and the numerical scheme. Nevertheless, it can be argued that LES provides the most optimum means of capturing the unsteady physical features of the flow. It is also easy to predict that LES will continue to replace RAS in many applications. In addition to these three general simulation methodologies there are some schemes which are labeled as “hybrid”. That is, they are based on the combination of some of these methods [1-2]. The success of these methods is dependent on that of the original methods. The analytical methods that have been developed for two-phase turbulent flows are generally based on three different approaches: (i) Eulerian-Eulerian approach [3-10], (ii) Eulerian-Lagrangian approach [7,11-20], and (iii) Lagrangian-Lagrangian approach [21-23]. In the first approach, the continuum transport equations for both phases are solved. These equations are somewhat similar and are often obtained by some sort of volume averaging which is conceptually different than the ensemble averaging in RAS or space averaging in LES. In the second approach, the continuum carrier fluid equations 99 are solved in its instantaneous form in DNS or in its “averaged” form in RAS and LES over a fixed Eulerian grid system. However, the “dispersed” phase (particles, droplets, micro bubbles, ...) are described by a set of modeled Lagrangian equations which determine the position, velocity, temperature, and other properties [7,13-14,16] of the dispersed phase. In the third method, both phases are described in the Lagrangian context. The Eulerian-Eulerian and Eulerian—Lagrangian models have been employed for RAS, LES and DNS of a wide range of two-phase or two-fluid turbulent flows including free surface flows, or dispersed turbulent flows [3-4,6-8,11-16,24-27]. In comparison, there are limited numbers of models based on Lagrangian-Lagrangian approach. An example is the two-phase random vortex model of Salman and Soteriou 21]. The two- phase Lattice-Boltzmann and Molecular Dynamics Simulation models [22-23] may also be considered as examples of Lagrangian-Lagrangian models. These models have not been seriously employed for LES and are not discussed here. Neither we will discuss the RAS models that have been developed based on Eulerian—Eulerian and Eulerian- Lagrangian approaches. We refer the interested reader to review articles and books by Faeth [11], Crowe [13], Sirignano [.14] and others on this subject. Furthermore, the extent of the previous DNS studies of two-phase flows is too broad to be discussed here and is not the focus of this chapter. In Ref. [7], a comprehensive collection of the DNS studies is provided by categorizing them based on flow configuration and coupling between the phases. As indicated, homogeneous two-phase turbulent flows (reacting or non-reacting) have been the subject of a large number of studies [28-39] due to their relative simplicity. Inhomogeneous two-phase turbulent flows have also received a great deal of attention because of their extensive technological applications [36,40-56]. 100 From modeling and computational point of view there are many similarities among liquid-gas, solid-gas, and liquid-solid two-phase (or two-fluid) turbulent systems. For example, turbulent flows laden with droplets or solid particles all involve a large number of small elements with densities very different than the carrier fluid density [16]. These flows are often computed with similar Eulerian-Lagrangian or Eulerian-Eulerian models, although different models are usually needed in different regions of the flow, depending on the concentration, size, and shape of the “particles”. The two-phase/two- fluid systems are often classified by the flow regimes: dilute, moderately dense, and very dense. All of these regimes are present in sprays [11,13-14,57—59]. In dilute flow regime, the particle-particle interactions are rare and the volume fraction of the particle to carrier fluid is very small, making the continuum assumption and Eulerian—Eulerian formulation essentially invalid. In these regimes, the particles are often treated as point masses (or point volumes for micro bubbles) which are affected by the carrier gas through mass, momentum, and energy interactions. The coupling could be two—way, allowing the particles to also modify the carrier fluid. In the most early studies on the LES of dilute two-phase turbulent flows with one-way coupling [60-61], particles were tracked by using only the Eulerian filtered velocity field and the effects of subgrid scale velocity field on the particles were not taken into account. However, these effects are often considered in later studies. For example, Armenio et al. [62] studied the SGS carrier gas velocity effects on dispersed phase by comparing the DNS and LES results. The two-way coupling effects were incorporated in LES by Boivin et al. [63] through grid-averaged source terms in the filtered Navier-Stokes equations. Yuu et a1. [64] consider the two- way coupling between the carrier fluid and the particle phases at SGS level by modeling 101 the subgrid-scale turbulent mass flux of particles with a gradient transport closure similar to molecular transport Closures and have achieved a good comparison with experiment for a planar turbulent jet flow. The SGS effects of particles on the fluid (and vice versa) are also modeled with a more rigorous PDF-based closure [5]. The Eulerian-Lagrangian dilute two-phase LES models have been applied to turbulent flows in complex geometries. For example, Apte et al. [65] consider solid particle dispersion in a non-reacting swirling coaxial jet combustor. In this work, the incompressible Eulerian carrier fluid equations are solved over an unstructured grid system and the cell-averaged effects of particles on the turbulence are incorporated via a source term in the momentum equation. Also in this work, as other works discussed above, a modified nonlinear version of the Lagrangian Stokes equation is used for the particle velocity calculations [l3-14,16]. With the point mass approximation, the entire mass and heat transfer and surface stress around the particle surface are represented as integrated local variables in the physical space such that the continuous-fluid domain is not affected by the presence of particles [16]. Subsequently, the average effect of heat and mass transfer and surface stress is obtained from a set of semi-empirical relations. Sirignano [l4] and Crowe et al. [13] provide a comprehensive review of available momentum, heat and mass transfer models for small particles and droplets. Also, in a comprehensive review article, Loth [16] compares the point mass models for interphase momentum transfer. The models provide the drag and lift forces for solid particles or droplets in the particle momentum equation. The effects of turbulence, compressibility, particle spin, particle deformation, virtual mass, Basset history, stress gradient, and wall on the particle forces are discussed in details. By comparing their numerical results with 102 experiment, Apte et al. [65] indicated that the LES predictions are significantly more accurate than those of RAS, despite the fact that the effects of SGS fluctuations of the carrier fluid on the particles are ignored in their LES calculations. Consistently, by comparing various multiphase turbulent flow models for environmental and hydrodynamic applications, Lakehal [8] concluded that the RAS is inherently inadequate for these applications and should be transcended by the LES. In the dense regime of any particulate turbulent flow, the particle-particle collisions are significant, the particle number density is large and the volume displaced by the particles often have noticeable influence on the carrier fluid. Normally, a significant portion of the physical activities associated with the fluid-particle and particle-particle interactions take place at SGS which are not calculated in LES and have to be modeled. This imposes serious challenges on LES of dense two-phase turbulent flows. The particle tracking LES models discussed above are based on point mass (or point volume) approximation which is only suitable for dilute systems with insignificant particle-particle collisions and small particle to fluid volume fraction. In theory, it is possible to simulate the flow domain outside the particle in the boundary layer and wake region. However, these simulations, referred to as “resolved-volume simulations” [16] are very expensive and are usually limited to one or handful number of particles in simple flows and do not seem to be practical for LES. Also, as discussed below, Eulerian-Eulerian mathematical/computational models may be used for dense two-phase flow simulations. However, these models may not be viable for complex turbulent flows that involve preferential concentration of the particles and the simultaneous presence of both dense and dilute regimes. 103 For LES of dense two-phase particulate flows, efficient particle-particle and droplet-droplet collision models are available in the literatures [66-68] that can be and have been incorporated in the Lagrangian particle simulations. For example, by using their Monte Carlo model of inter-particle collision, Yamamoto et al. [69] have simulated a gas-particle turbulent flow in a vertical channel and have shown that by adding the particle-particle collision models, the LES results compare better with those of experiment. They also show that the turbulence attenuation by particles agree well with experiment when particle collision models are included in the LES. Gorokhovski and Chtab [68] have suggested a similar collision model for LES of dense particulate flows. Their model is computationally more efficient than that of Yamamoto et al. [69] but generate almost the same results. Also, to incorporate the finite particle volume and particle collisions effects in an incompressible suspending fluid, Apte et al. [70] have developed a mathematical/computational model that is a modification of that by Dukowicz [71] for spray computation. The model has been shown to produce accurate results in simple laminar flows and has been incorporated in an unstructured LES code [70] for simulations of spray in complex flow configurations. While the interface tracking LES models described in previous chapters appear to be promising for primary breakup and free surface simulations, they do not seem to be appropriate for secondary breakup and dispersed regions of spray, where there are extensive numbers of individual droplets. In these regions, the volume fraction occupied by the droplets or ligaments are high and there is considerable collisions between droplets, which can lead to coalescence (forms bigger drops), shattering (forms smaller droplets), and bouncing or glancing (leaves droplet sizes unchanged). For LES of 104 atomizing spray, Apte et al. [17] have also proposed a Lagrangian stochastic SGS model for droplet transport and break up for dense part of the spray. The breakup model is based on Kolmogorov concept of viewing solid particle —breakup as discrete random process and considers the droplet breakup process as series of uncorrelated events that are independent of initial size of the droplets. The LES results obtained with this Lagrangian stochastic SGS model was shown to be more accurate in comparison with the standard parcel-approach usually employed in spray computations. There are some ongoing efforts by us and others on the development of hybrid numerical schemes for simultaneous and coupled simulation of dense and dilute spray flow parts. The LES models of this type seem to be practical for simulations of sprays in turbulent flow configurations. Most of LES models are developed for isothermal two-phase turbulent flows; the application of LES to non-isothermal flows with evaporating/reacting dr0plets is somewhat limited [7,15] and it has been only recently that LES models for realistic combustion systems are developed and tested. From the limited number of LES studies on two-phase turbulent flows that involves evaporating and reacting droplets, we refer here to the works by Okong'o, N. and Bellan [72]; Leboissetier et al. [73]; Sankaran and Menon [74]; Patel et al. [19]; Ham et al. [75]; Mahesh [20]; Cuenot et al. [10]; Afshari et al. [76] and Li &Jaberi [77]. In the first two papers [72-73], the three-dimensional DNS data on a temporal mixing layer laden with evaporating droplets are used to assess, a priori and posteriori several different SGS models for carrier gas, droplet phase and evaporated vapor. It is shown that with the scale-similarity models for the carrier gas velocity, the predicted droplet distribution by LES compare reasonably well with that 105 obtained from DNS data. However, the mass, momentum and heat transfer between evaporated droplets and carrier gas are not accurately represented by the proposed deterministic closures. The two-phase reacting LES model in Refs. [19,74] is based on the Euelrian-Lagrangian approach. In this model, the canier gas equations are solved over fixed grid points. The mixing and reaction are implemented in one-dimensional (1D) domain via linear eddy model (LEM) which is coupled with 3D LES flow solver. The spray is based on a Lagrangian droplet model which includes an empirical vaporation model and a secondary break-up model for the droplets. A global multi-step reaction model is employed. The model is applied to realistic (complex) gas turbine combustors. Similarly, the two-phase reacting LES closures in Refs. [20,75] are based on the hybrid Eulerian-Lagrangian approach. In the Eulerian part, the zero Mach number Navier-Stokes equations are solved on the unstructured grid system.The spray part also employs similar Lagrangian closures with some empirical evaporation and break up submodels. The combustion is included via a flamelet/progress variable reaction submodel. This model is also applied to practical gas turbine combustors. The two-phase reacting LES model of Cuenot et al. [10] is fundamentally different than the above models as it solves the Eulerian equations for both phases. Our two-phase reacting LES/FMDF model [76-77] is also different than the other models mentioned above. It is based on an Eulerian-Lagrangian-Lagrangian method that is described in greater detail below. 3.2 LES/FMDF for Two-Phase Turbulent Reacting Flows The two-phase LES/FMDF model is based on an Eulerian-Lagrangian-Lagrangian mathematical/numerical methodology and can handle the two-way interactions between 106 particle and fluid phases and particle—particle interactions for low particle to fluid volume fraction ratio. Figure 3-1 shows various elements of the model and its computational flow solver in a block diagram. The gas-phase part of the model is based on a high-order compact finite-difference numerical scheme. The subgrid gas-liquid combustion is modeled with the two-phase scalar FMDF [26,78-82]. The spray is simulated with a non- equilibrium Lagrangian model and stochastic SGS closures. The two-way mass, momentum, and energy coupling between phases is implemented through series of source/sink terms. The LES/FMDF employs a variety of different fuels based on two reaction 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 the scalars (mass fractions and enthalpy). Finite-rate chemistry effects are explicitly and “exactly” included in this procedure since the chemistry is closed in the FMDF formulation. In (2), the LES/FMDF is employed in conjunction with an equilibrium fuel- oxidation model. This model is enacted via “flamelet” simulations which normally consider a laminar counterflow (opposed jet) flame configuration. For two-phase flow calculations only non-equilibrium models based on (1) are used. The two-phase LES/FMDF model has been used for detailed and large-scale simulations of various droplet-laden turbulent systems with and without droplet evaporation and combustion. Figure 3-2 shows for example the vorticity and mass fraction contours and droplet distribution in one of the experiments we have simulated. In this experiment, a lean premixed preheated air-decane flame is controlled by injection of a relatively small amount of liquid n-decane fuel. Also shown in this figure are the 107 interactions between Eulerian grid points, Lagrangian droplets and Monte Carlo Particles. There are basically three interacting fields: (i) the Eulerian finite difference field for the gas dynamic variables, (ii) the grid-free Lagrangian Monte Carlo field for gaseous species and temperature as obtained by FMDF, and (iii) another Lagrangian field for liquid-fuel droplets and spray. In “conventional” LES methods, the “resolved” field is obtained by solving the filtered form of the compressible Navier-Stokes, energy and scalar equations with the filtered equations being closed by appropriate SGS stress and scalar flux models. In reacting flows, additional models are normally required for source/sink terms. Here, we use the FMDF which has been implemented in two ways: (Formulation I) to consider only the SGS scalar quantities, and (Formulation II) to consider the SGS velocity—scalar- pressure quantities. Formulation I is more manageable computationally, and formulation 11 is more rigorous from the statistical standpoint. Most of our previous contributions are based on formulation I which considers the joint scalar (species and energy) FMDF. For two-phase flows, the modified version of the scalar FMDF is employed. 3.2.1 Governing Equations for Carrier Gas Field For a two-phase reacting system, the filtered Navier-Stokes, energy and scalar equations are in the following forms: at + ax, 2 (3'1) a

(. > a

(...) lu -> am at.) an- [at I L + l a): I L =— (3le + a); l — ax:- +18j +l (3.2) 108 + = — (it at, dxi 3x, a a , a a

1<¢a>L +

1<¢a>L<”I>L :_ l I _dM,- + at Bar,- dxi ax,- Ns <¢a>L

1 z

1R0L Z (3.3) (3.5) where (f(x,r)), and (f (x,t)) L =(pf),/(p), represent the filtered and the Favre— filtered values of the transport variable f(x,t), and ,0 , ui, P , TandH are the fluid density, velocity, pressure, temperature, and total enthalpy, respectively. The species’ mass fractions represented in a common form with ¢a 5 Ya,“ = 1,2, ...... N3. In equation (3.5) R0 is the universal gas constant and Wa is the molecular weight of species a. Note that, in equation (3.3), the viscous dissipation and pressure gradient terms are neglected. The viscous stress tensor TI-j, heat flux 511' and mass flux 1,“!(0' =1,2...,NS) are defined as: <2I'J>Iz<fl>ld ax} at, _3 dxk ,_ : _<#>,. 600. [I Pr ax,- 109 (3.6) (3.7) (3.8) (3.9) where it, K, Cp. Pr and Le are dynamic viscosity coefficient, the thermal conductivity coefficient, the specific heat, the Prandtl number and the Lewis number, respectively. The SGS closures appear in the above filtered equations include the SGS stress sz =

l i<”i”j>L _L<"I>LJ’ the SGS enthalpy flux N,- =

llLLJ and the SOS scalar flux M,“ =

[[L —L<¢a)L] [56,78-79, 82]. In multiphase reacting flows, additional models are required for the filtered source/sink terms. Here, modeling of these filtered source/sink terms is the subject of the probability formulation described below. For hydrodynamic SGS quantities, a diffusivity-type closure model is used [56,78-79,82]. The variable-density form of the model for SGS stress can be expressed as: 1 2 r0- : _2<'0>I ”t 8517),, ‘giskk >250 1 + 3%, K5055.)- (3.10) where <50- >1. is the resolved strain rate tensor and defined as: (3.11) 1 a BI“) (50);; ail-H a)? In the Smagorinsky model, the subgrid eddy viscosity v, and coefficient K SGS are given by the following formulations, Vt : CS (3)2 (5,7)Ll (3.12) 2 (3.13) — 2 K505 =C1(A) |L While in the modified kinetic energy viscosity (MKEV) closure [78], they are expressed BS: 110 v, = CR(Z)(8)“ 2 (3.14) K505 =C18 (3.15) <~i*>,<~z-*>,.-<.>., lit-1.1.. reference velocity in the x,- -direction. The subscript L2 denotes the filter at the :1: where E = , u,- =u,— —U,-,,ef and Ui,ref is the secondary level of size with (3L2 )> (A). A similar model is used for the closures of SGS enthalpy flux and SGS scalar flux N,- =-—y, at- L (3.16) 'l a M.“ =-r. <::_>L (3.17) where y, = (@117, = (,0) I J:— and Sc, is the subgrid Schmidt number, assumed to be Sc, constant and equal to the subgrid Prandtl number. Note that, the subgrid eddy viscosity v, is the same one used in modeling of SGS stress Fij. It must be indicated here that this model is not used directly in the FMDF but the modeled FMDF transport equation is constructed to be consistent with it. The effects of droplets and combustion on carrier gas are expressed through the mass, momentum, energy, and species source/sink terms Sp, Sui, SH , Sa, and 55. In conventional LES methods, the filtered equations for the scalars (i.e. equation (3.4)) are solved together with other equations. In these equations, the filtered chemical source/sink terms are not closed and need modeling. Here, the subgrid combustion model is based on 111 the FMDF methodology and the temperature and mass-fraction fields are obtained from the FMDF. The chemical source/sink terms are determined exactly with the knowledge of FMDF. 3.2.2 Two-Phase FMDF Equations for Scalar field The scalar filtered mass density function (FMDF) is a joint probability density function of the scalars at the subgrid-level and is defined as: FL (3.1-,0 = [:p<:.t)t<$,ecc_’.r)>0<: acid: (3.18) 0' s‘f‘JiQQJD = 60L -93(2,t)) = 116011.. -— $530260) (3.19) a=1 where C denotes the filter function and the term {f(fiLQQJ» is the “fine-grained” density [83-84]. The scalar field, QQJ):(¢l(£,t),¢2(§,t),¢3(g,t)---~,¢0(g,t)) , represents the mass fractions of the chemical species and the specific enthalpy (pma = N5 +1), and is obtained from the joint scalar FMDF. E = (T135313 ----,‘I’0) is the composition domain of scalar array and 6 denotes the delta function. The equation (3.18) implies that the FIVIDF is the mass weighted spatially filtered value of the fine- grained density. The integral property of the FMDF is such that KFL‘EWWE = KPQCIWA’ 106% = (Mari), (3.20) For further development, the mass-weighted conditional filtered mean of the variable Q(._r,t) is defined as: I: QQL'.t)p(£’,t)§(1P..Q(L',t))G(L’-20d; 1711315,!) (Miami), = (3.21) 112 which implies 0') For Q(~_\‘.I)=c, (QQJHE), =c (3.22) (1.) For Q(-_\3t) = Q‘tgw». (012.2131), = 0(3) (3.23) (iii) Integral property: E:l FL (E;£,t)d_‘lj = l L (3.24) where c is a constant, and Q(Q(gg,t)) E Q(5,t) denotes the case where the variable Q can be completely described by the compositional variable QQJ). From these properties, it follows that the filtered value of any function of the scalar variables (such as pzthQgJfl for low Mach number flows and Sa, =Sa@(§,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 [85] when the density can be completely expressed in terms of the compositional variables. By applying the method proposed by Lundgren [86] and Pope [87] to the original(unfiltered) two-phase scalar and energy equations, the transport equation for fine-grained density is obtained apart — 9m» + 87214.6(: — e090) 2 at 8x,- ajl.“ aag — 9091)) _ amaifl - 9%,!» (3 25) ea, awa as", ' as “601' -(x,t)) 8;» S 501’ —(x.t)) . _ p —— — —- a p -— — — _ a‘Pa + aI-Pa, + Sp§(g Q(£’I)) The transport equation for F L (feet) is obtained by multiplying equation (3.25) by the filter function 0(5' — g) and integrating over 1' space. The final result after some algebraic manipulation is 113 BFL a a (< 1 01,-“ __ _ \{I F = at +0.\',- (1"1(l’)|_>) L) d‘I’a 062) (ix,- $> FL ‘53—(502 (30171.) I \ a (3.26) _ a Willa + a lair—Wm , 0.12m 8%. 73(1) 8‘12. 2(2) [2(3) This is an exact transport equation for the two-phase FMDF which can be solved by some models for unclosed terms and standard numerical methods like finite- difference method. However, standard methods are very expensive due to added dimensions. The second term on the right-handside of equation (3.26) is the chemical reaction term and is in a closed form. The last three terms on the right-handside of this equation represent the effects of droplets on the carrier gas and are unclosed. The unclosed nature of SGS convection (second term on the left-handside) and mixing (first term on the right-handside) 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 first decomposed via <14i|1>IFL = <1I1>LFL +1051!” FL —LFL) (3.27) where the second term on the right-handside of equation (3.27) represents the influence of SGS convective flux and can be modeled as a F / (0.1a, exert-042% (3.28) The advantage of the decomposition [equation (3.27)] and the subsequent model [equation (3.28)] is that they yield results similar to that in conventional LES [88-90]. For example, the first Favre moments corresponding to equations (3.27) and (3.28) are (we), = (40,08), + (0.42.), int-Mia). ) (3.29) 114 in). (Outta). - (u.>L<¢a>.)= -r. a<::>L (3.30) The term within brackets in equation (3.29) is the generalized scalar flux. This makes equation (3.30) identical to equation (3.17). The closure adopted for the SGS mixing is based on the linear mean square estimation (LMSE) model [83,91], also known as the IBM (interaction by exchange with the mean) [92] a [< 1 _a_(_y§_¢g g> FL]:i(yaL)FLi (3.31) where $2,” (5,!)15 the “frequency of mixing within the subgrid”, which is not known a priori. This frequency is modeled as S2," = CQ(}/+ y, )/(<,0)1AZG) here but other models can be used too. For the first term on the right-hand side of equation (3.31), an additional minor assumption is made “ a F / __a_[ya(FL/p(s>_>>)z_a_[y < L (9).] (3.32) 8X1: 8x,- dxi ax,- This assumption is not necessary for the treatment of FMDF and is only adopted to establish consistency between the FMDF and conventional LES. Also, for the droplet effect terms, the following approximations are made a isglillFL ~ 6 153%“ d‘f’a A!) ~ a‘ya (10)] (3.33) a (Sp|:r__)[\1’aFL ~ 8 <5p>l‘1’a,FL . ~ (3.34) d‘l’a Mi) 3‘11“ (10)] 115 (SpigliFL 2.. (5);)le (335) 80!.) (p), ' With these approximations, the modeled two-phase FMDF transport equation becomes BFL a _ a .. .. 31F]. 48),) a _ V+5ziLFLi-5;[(7 +71) 6x, + 8‘1", Qm(\ya <¢a>L)FLi a a . a 1FL —— S (3)17 — (3.36) B‘I’a [ a Li a?“ (,0), a (Salim (5,)1F, awe! (I0), (,0), 3.2.3 Lagrangian Equations for Droplet field The droplet field is solved with a Lagrangian mathematical/computational method [93—94]. In this method, the evolutions of the droplet displacement vector (X ,- ), the droplet velocity vector (v,- ), the droplet temperature (Tp ), and the droplet mass (m p) are governed by the following equations ex,- __ = V. 3.37 dr , ( ) fl 2 11mg“- v.)+-l-n. (3.38) dt 7,: Fr dT 22?. 2 112m“ 4p ) + LV de (3.39) (1! TP "IPCL (1t “mp = —————’"”f3 ln(1+ BM) (3.40) (1! Tp 116 where the asterisk refers to the local fluid variables which are interpolated to the droplet positions. z'P : pp (dp )2 /l8,u* is the normalized droplet time constant and pp , dp are droplet density and droplet diameter (the droplet is assumed to be spherical), respectively. The Stokes drag, fl is modeled as [95] f 1+0.0545Rep+0.1Relp/2(1—0.03Rep) l: b (3.41) 1+a|Reb| * lli —Vl' rip/If and Reb=p*Ubdp/,u* are particle Reynolds where Rep = ,0* numbers based on slip velocity and blowing velocity U =—(dm /dt)/(7r,0*d2), b p 1) respectively. The coefficients “a” and “b” are functions of droplet Reynolds number as a = 0.09 + 0.077 exp(—0.4 Rep) and b = 0.4 + 0.77 exp(—0.04 Rep) . In the energy equation (3.39), 77 = fl/[exp(,6)—1] is an analytical evaporative heat transfer correction, Imp is evaluated at where the non-dimensional evaporation parameter ,8 = —1.5P,Tpn'1 p previous time step. Also, the parameters f2 and f3 in these equations are corrections to the equilibrium heat and mass transfer coefficients and are calculated by the following relations: _ Sh f NuCp f '7 3—3Sc = , 3.42 - 3P,CL ( ) where Pr and Sc are the Prandtl and Schmidt number of gas mixture at droplet location, respectively. A mass averaging [94] is used to calculate the gas mixture heat capacity, N? C I, = ECHO/$6! and C L is the heat capacity of the liquid droplet. The Nusselt (Nu) a=l 117 and Sherwood (511) numbers are modified for finite Reynolds number effects on the semi-empirical Ranz-Marshall correlations: Nu = 2 +0.552Re‘p’2 P)”, Sh = 2 +0.552Rej,’2 52/3 (3.43) The evaporation rate is dependent on the mass transfer number emu/m —Yv )/(1 —YV,5) , where rm =1...q../ new +(1—Z,,,,,,,,)MWC IMWVJ and YV are the non-equilibrium vapor surface mass fraction and vapor mass fraction away from droplet surface. The non-equilibrium vapor molar fraction at droplet surface, 2’11qu is given by Langmuir-Knudsen evaporation law [93] as 2L Ina/,5 : let/,5 —[—d—E_]fl (344) P P L 2.4,. = “if." exp 0 V l —i (3.45) P R IMWV TB,L Td where T3914 is the liquid droplet’s boiling temperature and Rois the universal gas constant. Lv 2128 —(C L —Cp,v 17'], is the droplet latent heat of evaporation and I28 is the vapor reference enthalpy. MWV and MWC are molecular weights of vapor and carrier gas, respectively. While the Knudsen layer thickness is given by .. 0 /2 ,u (ZnTpR IMWV)l LK = * (3.46) ScP The volumetric phase coupling source/sink terms appearing in the filtered carrier gas equations and two-phase FMDF equations are evaluated based on the volumetric averaging of the Lagrangian field as 118 Sp = 1 {— "'Pf3 ln(1+ BM )} (3.47) TP 1 dm Sm. 2757215 +vi—d7’3) (3.48) s~=—-.:7:l’""°‘:“2"fz (ii p (3.49) 31‘ ’1‘ 1 vv +u-u. >1: dmp 25122{(hm ”Lil—zi—l—ut 12,-)7} where the summation is taken over all droplets in a volume centered at each Eulerian grid point and It” 2 C p.va +h8 is the evaporated vapor enthalpy at droplet surface. a]: 17,-2m £01,- — p TP v,) is the drag force and a2 = C L/C p is the ratio of liquid droplet heat capacity to the gas capacity. 3.3 Numerical Solution As mentioned before, the two-phase LES/FMDF model is based on an Eulerian— Lagrangian-Lagrangian mathematical/computational methodology, namely (i) the conventional Favre-filtered Eulerian LES equations are solved with high-order finite- difference scheme for the velocity field, (ii) the Lagrangian non-equilibrium droplet equations are solved with standard time-differencing methods, and (iii) the FMDF equation is solved with Lagrangian Monte Carlo method for the scalar field. The numerical methods used for each part are solved very different and they are discussed in detail below. 119 3.3.1 An Eulerian Finite Difference Method for Velocity Field The filtered Eulerian carrier-gas equations are solved together with diffusivity- type closures for the SGS stress and scalar flux terms via “standard” finite-difference schemes. The time differencing is based on the classical third-order TVD Runge-Kutta (R-K) scheme, which belongs to the strong stability-preserving (SSP) R-K family proposed by Shu and Osher [96] and can be expressed as (7 “I = U " + mud") (3.50) *(2) 3~n1~<0 1 (1) U =—U" +—u +— AtL(U ) (3.51) 4 4 4 0”“ =éu" +§U(2) +§A1L(U(2)) (3.52) —. Here, the vector U z{

l’<'0>I<”i>L’

IL’

l<¢a>L} and the operator L(l7) represent all the source terms and spatial derivatives in the filtered equations. In the equations (3.50)-(3.52), l7" and U’Mdenote the primary variables at present and next time steps, respectively. While (7 m and (7 ‘2’ are solutions at intermediate time steps. For any flow variable, the spatial derivative of the variable is computed by the sixth-order central compact differencing scheme (COMP6) proposed by Lele [97]. Given a flow variable f , its spatial derivative inxdirection f 'can be obtained by solving the following tridiagonal system 14 f_i_________+1- fi—l +lfi+2-fi_2 (3.53) 9 2Ax +9 4Ax l I I :fi—1+fi'+':31fi+l:_ with using the Thomas algorithm. The fifth-order one-sided and lop-sided compact boundary schemes of Cook and Riley [98] are implemented for boundary points and points next to boundary can be written as 5 1 '+4 ' =— (3.54) f1 f2 4x1?“ fl] , , , 1 5 afl +f2 +0)"3 =—— be (3.55) Ax H: where a1 = —37/12,(13 = 2/3,a3= 3,a 4=—2/3, a5 =l/12 for the equation (3. 54) and a = 0.2142857143, bl = —0.6785714286, b2 = —0.1 19047619,b3 = 0.8571428571, b. = —0.07142857l43 and b5 = 0.0119047619 for the equation (3.55). The application of non-dissipative spatial schemes, such as the central compact schemes causes a pileup of energy at the smallest scales of the flow. To avoid numerical instability and to eliminate numerical oscillations at very small scales due to mesh non- uniformities and boundary conditions [99], the low-pass, hi gh-order implicit spatial filter [100-101] is applied to interior and near boundary points (no filter is needed for boundary points). This filtering operation is very different than the standard LES filtering operation but is an important component of the solution algorithm. The eighth-order central and sixth-order one-sided filtering are used in this study, 4 A A x a . affr—1+ft+a’fft+1 = 271((fnk1’fi-k) 5515N‘4 (3-56) k=0 A A A 3 aff,-_l +f,- +aff,,1— _ Z—(f,,k +f,-_ k) i=4and N—3 (3.57) k: 02 . ,. . 3 aff,_l +f, +o.ff,-+1 = gum-f, i=2and3 (3.58) k=l where f is the filtered value of the original variable f . The filtering parameter af lies in the range of (-0.5, 0.5), with higher values representing less filter. The coefficients ak and bk in equations (3.56)-(3.58) are the same used by others [100-101]. 3.3.2 Lagrangian Model for Dilute Spray Once the fluid velocity, density and temperature are obtained from LES and two- phase FMDF solvers, the droplet transport equations are integrated with a second-order time differencing scheme. The evaluation of the fluid velocity at the droplet locations is based upon a fourth-order accurate Lagrangian interpolation method. Also for the droplet phase, a stochastic velocity model is considered by which the residual or subgrid velocity of the carrier gas at the droplet location is constructed as 2 * 0P 1 a <“i>L (a: ) du = ———+ +G--u- — u- dt+ C EdW. 3.59 3x1 R60 axjaxj ’1 ‘ <‘>L 0 I ( I where Gij = —w[% AC0 )ng , E = CEkZ/3/Af anda) = E/k. The combined large- and small-scale fluid velocity is then used for calculations of droplet location and velocity. No model is needed for subgrid energy (or temperature), density and mass fraction, since they can be obtained exactly from the subgrid two-phase FMDF equation. 3.3.3 Lagrangian Monte Carlo Solution of FMDF From an operational standpoint, PDF or FMDF methods are implemented via stochastic differential equations. The full two—phase FMDF equation (3.36) is a partial differential equation in physical and composition space. The connection to stochastic differential equations is via the Fokker-Planck equation: the Fokker-Planck equation of the stochastic model is a transport equation for the FMDF. Rather than solving that 122 partial differential equation, and then computing statistics, the statistics can be obtained far more economically by averaging of the stochastic process. Hence, in FMDF modeling one actually solves a set of stochastic differential equations (SDEs), with the partial differential equations being a point of reference. The most convenient means of solving the SDE equations is via the “Lagrangian Monte Carlo” procedure [78-79]. With the Lagrangian procedure, the two-phase 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 (filtered flow velocity) and diffusion (molecular and subgrid). In addition, changes in the “composition space” occurs due to chemical reaction, SGS mixing, and droplet effects, such as evaporation. Again, all of these are implemented via a stochastic process described by a set a SDEs. These SDEs are fully consistent with the original two-phase FMDF transport equation (3.36). The two-phase FMDF represents the gaseous scalar and energy fields and is used to evaluate the local values of the temperature, density and species mass fractions at droplet location. The droplets in turn modify the species concentration and temperature values of the Monte Carlo particles or FMDF due to mass and energy coupling. Hence, the three-way coupling between the carrier gas velocity field, FMDF scalar field and droplet field are included in the computations. In physical space, the spatial transport of the two-phase FMDF is represented by the general diffusion process governed by the following stochastic differential equation (SDE) [102-103] dzmc (t) : QINC (EINC (t),t)dt + E(—X”1C (t),t)d_W_”lc , (3.60) 123 where 5"“ is the Lagrangian position of a Monte Carlo particle, ch and E are the “drift” and “diffusion” coefficients respectively, and Wm“ denote the Wiener process [104]. The drift and diffusion coefficients are obtained by comparing the Fokker-Planck equation corresponding to equation (3.60) with the spatial derivative terms in the two- phase FMDF transport equation (3.36), III(.‘ III(‘ 8 ~ ~ 2 (2S (z),z)=[(g)L +< 1) (75”? (3.61) p I E xnlC(t) E(_>g”"(1),1) = [(281102 + )7)/

l )“ijmcm (3.62) I In composition space, the subgrid mixing, reaction and droplet effects are implemented by altering the compositional makeup of the particles via the following equafion: . S a + S We) “(Dam=—Q:z[¢;(I)_<¢;(I)> 1+sal+ a ’ (3'63) dt L 70(9 (1)) 0(9 (0) where (I); =¢a(£(r),t)denotes the scalar value of the Monte Carlo particle at the + Lagrangian position vector 5(1) . The same notations are used for <83) , I <5 p>l+ and 9+0). The solutions of equations (3.60) and (3.63) should generate fully consistent statistics with those obtained by the original two-phase FMDF transport equation (3.36), according to the principle of equivalent systems [84,105]. In the Monte Carlo method, each particle is initially assigned a weight which represent the mass of carrier gas and satisfies the following equation in each cell, 124 2W," z (,0), AV (3.64) we cell where Wm and AV are the weight of particle m and the volume of the ensemble domain used for particle averaging, respectively. In an attempt to reduce the computational cost, a procedure involving the use of non—uniform weights is also considered. This procedure allows a smaller number of particles to be imposed in regions where a low degree of variability is expected. Conversely, in regions of high varying character, a larger number of particles are allowed. This is akin to grid compression in the finite-difference (or finite-volume) schemes. At locations far away from chemical reaction zone and spray, no Monte Carlo particles are deployed. For the case with droplet evaporation, the weights of particles are modified according to the following equations to include the added mass (density) effect on FMDF. flat. = +AV (3.65) <5e>+ "’ WASP), I_ Zwm me cell (3.66) Equation (3.65) is consistent with the filtered continuity equation (3.1). Statistical information, such as the first filtered moment, at any point is obtained by considering an ensemble of N E computational particles residing within an ensemble domain characterized by the length scaleA E- This is necessary because, with probability one, no particle will coincide with the point considered [105]. From a numerical standpoint, specification of the size of the ensemble domain is important. Ideally, it is desired to obtain the statistics from the Monte Carlo solution when the size of sample 125 domain is infinitely small (AE ——> 0) and the number of particles within this domain is infinitely large. With a finite number of particles, if the ensemble domain is too small there may not be enough particles to construct reliable statistics. A larger ensemble domain decreases the statistical error, but increases the spatial error which manifests 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 [84- 85]. This limitation does not diminish the capability of the scheme, but exemplifies the importance of the parameters which govern the statistics. Fourth-order interpolation is implemented for the transfer of information from the Eulerian grid points to the Monte Carlo particles. While the transfer of information from Monte Carlo particles to the Eulerian grid points is accomplished by using ensemble averaging as described above. The ensemble or Favre-filtered average value of a transport quantity like Q(Q) is constructed from the weighted average, ZWmQ(Q—rrt) L = "M5 (3.67) 2 Wm mEAE Again, the approximations in equation (3.64) and (3.67) improve when A E ——> 0 and the number of particles within the ensemble domain increases [84]. 3.4 Consistency of Eulerian LES and Lagrangian FMDF Solutions An important issue in a hybrid Eulerian-Lagrangian mathematical/computational methodology like LES/FMDF, is the consistency of its sub-components. Figure 3-3 126 shows some of important attributes of conventional LES and FMDF parts of LES/FMDF. It is shown that some of the variables are calculated by finite difference method (FD), some by Monte Carlo method (MC), and some by both. That is, there is a “redundancy” in the determination of some of the quantities. This seems to be a disadvantage of LES/FMDF as extra efforts are needed to calculate the redundant variables. However, the redundant variables provide an opportunity for testing the numerical accuracy of the model. We have examined the consistency of the two-phase FMDF and its Lagrangian Monte Carlo solver with the finite-difference and conventional LES models in our LES/FMDF calculations for all simulated cases. This is done for (1) non-reacting flows without spray, (2) non—reacting flows with spray and droplet evaporation, (3) reacting flows without spray, and (4) reacting flows with fuel spray and evaporation. In the discussion below, the scalar (temperature and species mass fraction) values generate by the finite-difference method on the Eulerian grid points are denoted by FD and the results obtained by ensemble averaging of the Monte Carlo particles are denoted by MC. For the density, MCl denotes the density calculated through equation (3.67) with replacing A Q(_) by 1/,0(§3) and MC2 denotes the density calculated via equation (3.64). The consistency between the FD and MC indicates the accuracy and reliability of the computed values for both of them. We would like to emphasize that it is possible to compare the FD results with the MC results in our simulation because the corresponding reaction/droplet source/sink terms in FD are computed from the MC. Such terms are not available in “standard” LES/FD methods. Sample results obtained by our two-phase LES/FMDF method for the dump combustor experiment [106] are presented in Figure 3-2, where the contours of the 127 carrier gas vorticity and droplets are shown. In the dump combustor experiment, a lean premixed air-fuel flame is controlled by injection of a relatively small amount of liquid fuel. We have simulated this experiment with and without (Heptane) fuel droplet evaporation and combustion and with realistic flow and droplet parameters. Contours of the instantaneous temperature as obtained by LES-FD and FMDF-MC are compared in Figure 3-4 for the case with combustion but no fuel spray. There seems to be a good agreement in the instantaneous results, even though the values obtained from the MC particle seems to be a little noisy, which is due to the limited number of MC particles. The corresponding fuel mass fraction contours (not shown) also indicate a nearly perfect match between FD and MC results. Figure 3-5 shows the radial variations of the filtered temperature and the fuel mass fraction at three locations from the inlet for the Eulerian FD and Lagrangian MC fields in the dump combustor for the case that there is combustion but liquid-fuel spray is off. Evidently, the MC predictions are in good agreement with the FD results in all cases considered, which again demonstrate the good consistency of the LES/FMDF method for single-phase reacting flows. Similarly, the instantaneous contours of fuel mass fraction (Figure 3—6), the instantaneous contours of temperature (Figure 3-7) and the radial plots of both fuel mass fraction and temperature at different axial locations (Figure 3-8) for the case with combustion and liquid-fuel spray, indicate a good consistency between FD and MC values. With the evaporation of liquid-fuel, the cooling effect on the temperature and added effect of sprayed fuel on the gas field are shown to be important in the above contour and radial plots. To further demonstrate the consistency between LES-FD and FMDF-MC and to establish the numerical accuracy of LES/FMDF model, the contours of instantaneous 128 values of the filtered temperature and the radial variations of filtered temperature, fuel mass fraction and density in a double swirl spray burner [107] are compared in Figure 3- 9 and Figure 3-10, respectively. In the double-swirl burner the combustion is controlled by liquid fuel kerosene spray. There is no fuel in the swirling air streams. Evidently, the FMDF-MC predictions of the instantaneous temperature are very close to those of LES- FD. The radial variations of the filtered temperature, fuel mass fraction and density as obtained from the LES-FD data at different axial locations are shown to be very close to those obtained from FMDF data. The discrepancy between MC2 and FD values of the density at larger radial locations is due to that fact that no Monte Carlo particles are deployed there in the original simulations. This can be and is easily fixed by adding particles. The above results demonstrate the consistency of the two-phase LES/FMDF methodology and its numerical solution in the presence of spray and combustion. Based on these results, we can conclude that the MC and FD formulations in the two-phase LES/FMDF method are both mathematically correct and numerically accurate. Again, we note that in the reacting two-phase flows, the reaction and droplet source/sink terms in the FD equations are obtained from the MC particles. This is for the testing of the numerical methods and for showing the consistency. Such information is not available in “standard” LES-FD methods. 3.5 Summary In this chapter, a new Eulerian-Lagrangian-Lagrangian two-phase LES/FMDF mathematical/numerical methodology is presented for two-phase turbulent reacting flows which include with two-way mass, momentum and energy coupling between phases. In 129 this two-phase LES/FMDF methodology, the filtered Eulerian equations for the carrier gas are calculated by a high-order compact finite-difference numerical scheme. The subgrid species energy transport, mixing and combustion are modeled with the two- phase scalar FMDF transport equation, which is solved by the Lagrangian Monte Carlo method. The liquid droplet/spray is simulated with a non-equilibrium Lagrangian model and stochastic SGS closures. The two-way mass, momentum, and energy coupling between phases is implemented through series of source/sink terms. The main features of the two-phase LES/FMDF model are (1) Large scale, unsteady, non-universal, geometry-depended quantities are explicitly computed in LES/FMDF. (2) FMDF accounts for the effects of chemical reactions and droplet source/sink in an exact manner and maybe be used for various types of chemical reactions. (3) Two-phase LES/FMDF can be implemented via complex chemical kinetics models and is applicable to 3D simulations of hydrocarbon flames in complex geometries. (4) Two-phase FMDF contains high-order information on sub-grid (or small- scale) fluctuations. (5) The Lagrangian Monte Carlo solution of the two-phase FMDF is free of artificial (diffusion) numerical errors, which is typical of FD simulations on fixed grid points. The demonstrated good consistency between LES—FD and FMDF-MC results for various turbulent flows with and without fuel spray and combustion indicates that the MC and FD formulations in the two-phase LES/FMDF are both mathematically correct and numerically accurate. The results shown in this chapter (and those not shown) also indicate that the two-phase LES/FMDF model is a reliable and affordable methodology for detailed calculations of turbulent spray combustion. The model is based on a robust 130 mathematical/computational framework and can be continuously improved and applied to increasingly more complex systems without having to change its basic framework. 131 3.6 Figures LES of Two-Phase Turbulent Reacting Flows A New Lagrangian-Eulerian-Lagrangian Methodology K Filtered continuity and momentum equations \ via a generalized multi-block high-order finite . . , difference Eulerian scheme for high Reynolds LGasdynamIcs He'd j - number turbulent flows in complex geometries carious closures for sub grid stresses j ‘ / Scalar Field Filtered Mass Density Function (FMDF) (mass fractions - equation via Lagrangian Monte Carlo method - i and temperature) L Ito Eq. for convection, diffusion 81 reaction J , Lagrangian model for droplet equations with Droplet Field 1 full mass, momentum and energy couplings l l (spray) 1 between phases and a stochastic sub grid 1 velocity model Kinetics: (1) reduced kinetics schemes with 1 direct ODE or ISAT solvers, and (II) flamelet Chemistry library with detailed mechanisms or complex reduced schemes. _ Fuels: methane, propane, decane, kerosene, heptane, JP-10 Figure 3-1. A block diagram showing different components of the LES/FMDF and its Lagrangian-Eulerian-Lagrangian flow solver. 132 Spray-Controlled ‘ Dump Combustor I I I 285 Wall 14 § - Eulerian'Grld e3 8 M M t 5w 0 Monte Carlo Particles 1. 838, omen um. If . g “3 Terms from Droplets quuld Fuel Droplets u s ' ° / . e e Eulenan Cell , , e . -~e~ - Solver *— , : ‘— e . Eulorlan Flnlte Interptlonl a: M 111 ”WWW. GM Favre Fllter garflclu 0 Figure 3-2. Contours of the instantaneous vorticity magnitude, droplets, Monte Cano particles and grid layout in a spray-controlled dump combustor as obtained by LES/FMDF. 133 (,0): Conventional LES (11,-) L Solver (LES-FD) A = 48—23;,Izgiaja (T). l <¢a>L [conserved scalar I 1 j (Th A (1050,) (15.1) I L J Figure 3-3. The attributes of LES/FMDF methodology and its LES-FD and FMDF subcomponents. (b) Figure 3-4. Contours of instantaneous filtered temperature in the premixed dump combustor without fuel spray: (a) FD values and (b) MC values. 134 \l Temperature C») J?- 01 O) N «h 01 1'! VIII IT Temperature 0) 01 A YYT fi‘Tffi'IIVT T T1 1 Temperature GD '9 2 l A g2_2 L2 4‘ A2 Figure 3-5. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations in the premixed dump combustor without fuel 0.5 r/D 1 1.5 spray as obtained by LES/FMDF. Fuel Mass Fraction Fuel Mass Fraction IrTY III! I I 0.04 . (E 0.03 f 0.02 ’ Fuel Mass Fraction 0.01 oflTVITIIY 135 (a) (b) Figure 3-6. Contours of instantaneous filtered fuel mass fraction in the premixed dump combustor with fuel spray: (a) FD values and (b) MC values. (b) Figure 3-7. Contours of instantaneous filtered temperature in the premixed dump combustor with fuel spray as obtained by LES/FMDF: (a) FD values and (b) MC values. 136 Temperature _. Am 00 a. 01 o: \r \l O) I VIT Fuel Mass Fraction 0.2 0.15" c . .9 ~ : ‘5 0.15 e 5» e 3 4 IL 9 ill a (U 0.1% E 3 E o _ '5' y '- 2» E 0.05‘ 0.1 x/D 2 FD MC 7. 0.2 6f. C so I .9 MC 2 5; '3 0.15» .—. 1 ,, ‘15 4“- a, 2 E 3 e 33‘ 2 m b - " 2: ii 1 : X/D=6 ' ~ 0 1 l . . . 2_l O’- 1 l A 2 . 0 0.5 rID 1 1.5 o 0.5 rID Figure 3-8. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations as predicted by LES/FMDF in the premixed dump combustor with fuel spray. 137 (a) (b) Figure 3-9. Contours of Instantaneous filtered temperature in the double swirl spray burner as obtained by LES/FMDF: (a) FD values and (b) MC values. 138 \JCD A I Temperature N (A) -h 01 O) \l (I) Temperature 10 co #‘ 01 a) _L O Fuel Mass Fraction Fuel Mass Fraction 0 0.5 rlR 1 1.5 2 XIR=2 l . . . . I . I . I #J . . A 1 0 . . . A l r . . 1 . 1.5 2 O 0.5 r/R 1 1.5 2 O . 1 L 2 I 1 . l I I o 0.5 rlR 1 Figure 3-10. Radial variations of instantaneous filtered temperature and fuel mass fraction at different axial locations in the double swirl spray burner as predicted by LES/FMDF. 139 3.7 [1] [2] l3] l4] [5] [6] [7] [8] [9] [10] [11] [12] References Basu, D., Hamed, A. and Das, K., “DES and hybrid RANS/LES models for unsteady separated turbulent flow predictions,” AIAA paper 2005-0503, Jan.2005. Kenjeres, S. and Hanjalic, K., “Transient analysis of Rayleigh-Benard convection with a RANS model,” International Journal of Heat and Fluid Flow, Vol. 20, 1999, pp.329-340. Scardovelli, R. and Zaleski, 8., “Direct numerical simulation of free-surface and interfacial flow,” Annual Review of Fluid Mechanics, Vol. 31, 1999, pp.567-603. Shyy, W. and Narayanan, R., Fluid Dynamics at Interfaces, Cambridge University Press, 1999. Pandya, R. and Mashayek, F., “Two-fluid large-eddy simulation approach for particle-laden turbulent flows,” International Journal of Heat and Mass Transfer, Vol. 45, 2002, pp.4753-4759. Sethian, J. A. and Smereka, P., “Level set methods for fluid interfaces,” Annual Review of Fluid Mechanics, Vol. 35, 2003, pp.341-372. Mashayek, F. and Pandya, R., “Analytical description particle/droplet- laden turbulent flows,” Progress in Energy Combustion Science, Vol. 29, 2003. PP.329-379. Lakehal, D. “On the modeling of multiphase turbulent flows for environmental and hydrodynamic applications,” International Journal of Multiphase Flow, Vol. 28, 2002, pp.823-863. Herrmann, M, “Modeling primary breakup: A three-dimensional Eulerian level set/vortex sheet method for two-phase interface dynamics,” Annual Research Briefs, Center for turbulence research: Stanford, 2003, pp.185- 196. Cuenot, B., Boileau, M., Pascaud, S., Mossa, J.-B., and Poinsot, T., “Large eddy simulation of two-phase reacting flows,” ECCOMAS CFD 2006, Egmond Aan Zee, The Netherlands, 2006. ~ Faeth, G. M., “Mixing, transport and combustion in sprays,” Progress in Energy Combustion Science, Vol. 13, 1987, pp.293-345. Shirolkar, J. S., Coimbra, C. F. M. and McQuay, M. Q., “Fundamental aspects of modeling turbulent particle dispersion in dilute flows,” Progress in Energy Combustion Science, Vol. 22, 1996, pp.363-399. 140 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] Crowe, C., Sommerfeld, M. and Tsuji, Y., Multiphase Flows with Droplets and Particles, CRC Press, 1998. Sirignano, W. A., Fluid Dynamics and Transport of Droplets and Sprays, Cambridge University Press, 1999. Bellan, J., “Perspectives on Large Eddy Simulation for sprays: Issues and solutions,” Atomization and Sprays, Vol. 10, 2000, pp.409-425. Loth, E., “Numerical approaches for motion of dispersed particles, droplets and bubbles,” Progress in Energy Combustion Science, Vol. 26, 2000, pp.161-223. Apte, S. V., Gorokhovski, M. and Moin, P., “LES of atomizing spray with stochastic modeling of secondary breakup,” International Journal of Multiphase Flow, Vol. 29, 2003, pp.1503-1522. Menon, S. and Patel, N., “Subgrid modeling for simulation of spray combustion in large-scale combustors,” AIAA Journal, Vol. 44, 2006, pp.709-723. Patel, N., Kirtas, M., Sankaran, V. and Menon, S., “Simulation of spray combustion in a lean-direct injection combustor,” Proceedings of the Combustion Institute, Vol. 31, 2007, pp.2327-2334. Mahesh, K., Constantinescu, G., Apte, S., laccarino, G., Ham, F., and Moin, P., “Large-eddy simulation of reacting turbulent flows in complex geometries,” Journal of Applied Mechanics, Vol. 73, 2006, pp.374-381. Salman, H. and Soteriou, M., “Lagrangian simulation of evaporating droplet sprays,” Physics of Fluids, Vol. 16, 2004, pp.4601-4622. Seta, T., Kono, K. and Chen, S., “Lattice Boltzmann method for two-phase flows,” International Journal of Modern Physics B, Vol. 17, 2003, pp.169- 172. Chen, S. and Doolen, G. D., “Lattice Boltzmann method for fluid flows,” Annual Review of Fluid Mechanics, Vol. 30, 1998, pp.329-364. Cressman, J. R., Davoudi, J., Goldburg, W. I. and Schumacher, J., “Eulerian and Lagrangian studies in surface flow turbulence,” New Journal of Physics, Vol. 6, 2004, pp.53-77. Menon, S., “Subgrid combustion modeling for LES of single and two- phase reacting flows,” Advances in LES of Complex Flows, Friedrich, R. and Rodi, W. (eds), Proceedings of the Euromech Colloquium 412, Munich, Germany, 2000, pp.329—351. 141 [251 [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] Madnia, C. K., Jaberi, F. A. And Givi, P., “Large eddy simulation of heat and mass transport in turbulent flows,” Handbook of Numerical Heat Transfer 2““, Minkowycz, W. J., Sparrow, E. M. and Murthy, J. Y. (eds), Wiley: New York, 2006, pp.167-189. Moin, p., “Large eddy simulation of multi-phase turbulent flows in realistic combustors,” Progress in Computational Fluid Dynamics, An International Journal, Vol. 14, 2004, pp.237-240. Squires, K. D. and Eaton, J. K., “Particle response and turbulence modification in isotropic turbulence,” Physics of Fluids, Vol. 2, 1990, pp.1191-1203. Squires, K. D. and Eaton, J. K.,”Measurements of particle dispersion obtained from direct numerical simulations of isotropic turbulence,” Journal of Fluid Mechanics, Vol. 226, 1991 , pp.1-35. Squires, K. D. and Eaton, J. K., “Preferential concentration of particles by turbulence,” Physics of Fluids, Vol. 3, 1991 , pp.1 169-1 178. Elghobashi, S. and Truesdell, G. 0., “Direct simulation of particle dispersion in a decaying isotropic turbulence,” Journal of Fluid Mechanics, Vol. 242, 1992, pp.655-700. Elghobashi, S. and Truesdell, G. C., “On the two-way interaction between homogeneous turbulence and dispersed solid particles (I): turbulence modification,” Physics of Fluids, Vol. 5, 1993, pp.1790-1801. Truesdell, G. C. and Elghobashi, 8., “On the two-way interaction between homogeneous turbulence and dispersed solid particles (ll): particle dispersion,” Physics of Fluids, Vol. 6, 1994, pp.1405-1407. Wang, LP and Maxey, M. R, “Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence,” Journal of Fluid Mechanics, Vol. 256, 1993, pp.27-68. Mashayek, F., Jaberi, F. A., Miller, R. S. and Givi, P., “Dispersion and polydispersity of droplets in stationary isotropic turbulence,” International Journal of Multiphase Flow, Vol. 23, 1997, pp.337-355. Mashayek, F., “Direct numerical simulation of evaporating droplet dispersion in forced low Mach number turbulence,” International Journal of Heat and Mass Transfer, Vol.41, 1998, pp.2601-2617. Jaberi, F. A. and Mashayek, F., “Temperature decay in two-phase turbulent flows,” International Journal of Heat and Mass Transfer, Vol. 43, 2000, pp.993-1005. 142 [38] [39] [40] [41] [42] [43] [44] [45] [45] [47] [48] [49] Mashayek, F. and Jaberi, F. A., “Particle dispersion in forced isotropic low Mach number turbulence,” International Journal of Heat and Mass Transfer, Vol. 42, 1999, pp.2823-2836. Maxey, M. R., Patel, B. K., Chang, E. J. and Wang, L.-P., “Simulations of dispersed turbulent multiphase flow,” Fluid Dynamics Research, Vol. 20, 1997, pp.143-156. Mashayek, F., “Droplet-turbulence interactions in Iow-Mach-number homogeneous shear two-phase flows,” Journal of Fluid Mechanics, Vol. 376, 1998. PP.163-203. Mclaughlin, J. B., “Aerosol particle deposition in numerically simulated channel flow,” Physics of Fluids, Vol. 1, 1989, pp.1211-1224. Ounis, H., Ahmadi, G. and Mclaughlin, J. B., “Brownian particle deposition in a directly simulated turbulent channel flow,” Physics of Fluids, Vol. 5, 1993. pp. 1 427-1 432‘. Brooke, J. W., Kontomaris, K., Hanratty, T. J. and Mclaughlin, J. B., “Turbulent deposition and trapping of aerosols at a wall,” Physics of Fluids, Vol. 4, 1992, pp.825-834. Brooke, J. W., Hanratty, T. J. and Mclaughlin, J. B., “Free-flight mixing and deposition of aerosols,” Physics of Fluids, Vol. 6, 1994, pp.3404-3415. Pedinotti, S., Mariotti, G. and Banerjee, 8., “Direct numerical simulation of particle behavior in the wall region of turbulent flows in horizontal channels,” International Journal of Multiphase Flow, Vol. 18, 1992, pp.927- 941. Chen, M. and Mclaughlin, J. B., “A new correlation for the aerosol deposition rate in vertical ducts,” Journal of Colloid and Interface Science, Vol. 169, 1995. 90437-455. Chen, M., Kontomaris, K. and Mclaughlin, J. B., “Dispersion, growth, and deposition of coalescing aerosols in a direct numerical simulation of turbulent channel flow,” Gas-Particle Flow, ASME-FED, Vol. 228, 1995, pp.27-32. Rouson, D. W. l., Eaton, J. K. and Abrahamson, S. D., “A direct numerical simulation of a particle-laden turbulent channel flow,” Department of Mechanical Engineering Report TSD-101, Stanford University , Stanford, CA, 1997. Liao, S., Mashayek, F. and Guo, D., “Numerical simulations of particle- laden axisymmetric turbulent flows,” Numerical Heat Transfer, Part A: Application, Vol. 39, 2001, pp.847-855. 143 [50] [51] [52] [53] [54] [55] [55] [57] [58] [59] [60] [61] [52] Ling, W., Chung, J. N. and Crowe, C. T., “Direction numerical simulation of the two-way coupled interaction between particles and mixing layer,” Proceedings of Royal Society, London A. Vol. 437, 2000, pp.2931-2955. Miller, R. S. and Bellan, J., “Direct numerical simulation and subgrid analysis of a transitional droplet laden mixing layers,” Physics of Fluids, Vol. 12, 2000, pp.650-671. Okong’o, N. and Bellan, J., “A priori subgrid analysis of temporal mixing layers with evaporating droplets,” Physics of Fluids, Vol. 12, 2000, pp.1573-1591. Mashayek, F., “Numerical investigation of reacting droplets in homogeneous shear turbulence,” Journal of Fluid Mechanics, Vol. 405, 2000, pp.1-36. Kuerten, J. G. M., “Subgrid modeling in particle-laden channel flow,” Physics of Fluids, Vol. 18, 2006, 025108. Almeida, T. and Jaberi, F. A., “Direct numerical simulations of a planar jet laden with evaporating droplets,” International Journal of Heat and Mass Transfer, Vol. 49, 2006, pp.2113-2123. Almeida, T. and Jaberi, F. A., “Large eddy simulation of a particle-laden turbulent jet flow with a stochastic two-phase subgrid closure,” International Journal of Heat and Mass Transfer, in press, 2007. Lefebrvre, A. H., Atomization and Sprays Hemisphere Publishing, New York, 1989. Lin, S. P. and Reitz, R. D., “Droplet and spray formation from a liquid jet,” Annual Review of Fluid Mechanics, Vol. 30, 1998, pp.85-105. Lasheras, J. C. and Hopfinger, E. J., “Liquid jet instability and atomization in a coaxial gas stream,” Annual Review of Fluid Mechanics, Vol. 32, 2000, pp.275-308. Yeh, F. and Lei, U., “On the motion of small particles in a homogeneous isotropic turbulent flow,” Physics of Fluids, Vol. 3, 1991, pp.2571-2586. Yeh, F. and Lei, U., “On the motion of small particles in a homogeneous turbulent shear flow,” Physics of Fluids, Vol. 3, 1991 , pp.2758-2776. Armenio, V., Piomelli, U. and Fiorotto, V., “Effect of the subgrid scales on particle motion,” Physics of Fluids, Vol. 11, 1999, pp.3030-3042. 144 [63] [54] [65] [65] [57] [68] [69] [70] [71] [721 [73] [74] Boivin, M., Simonin, O. and Squires, K. D., “On the prediction of gas-solid flows with two-way coupling using large eddy simulation,” Physics of Fluids, Vol. 12, 2000. 902080-2090. Yuu, S., Ueno, T. and Umekage, T., “Numerical simulation of the high Reynolds number slit nozzle gas-particle jet using subgrid-scale coupling,” Chemical Engineering Science, Vol. 56, 2001, pp.4293-4307. Apte, S.V., Mahesh, K., Moin, P. and Oefelein, J. 0., “Large eddy simulation of swirling particle-laden flows in a coaxial jet combustor,” International Journal of Multiphase Flow, Vol. 29, 2003, pp. 1 31 1-1331. Schmidt, D. P. and Rutland, C. J., “A new droplet collision algorithm,” Journal of Computational Physics, Vol. 164, 2000, pp.62-80. Post, S. L. and Abraham, J., “Modeling the outcome of drop-drop collisions in diesel sprays,” International Journal of Multiphase Flow, Vol. 28, 2002, pp.997-1019. Gorokhovski, M. and Chtab, A., “Hypothetical heavy particles dynamics in LES of turbulent dispersed two-phase channel flow,” Annual Research Briefs, Center for turbulence research: Stanford, 2003. Yamamoto, Y., Potthoff, M., Tanaka, T., Kajishima, T. and Tsuji, Y., “Large eddy simulation of turbulent gas-particle flow in a vertical channel: effect of considering inter-particle collisions,” Journal of Fluid Mechanics, Vol. 442, 2001, pp.303-334. Apte, S.V., Mahesh, K. and Lundgren, T., “An Eulerian-Lagrangian model to simulate two-phase/particulate flows,” Annual Research Briefs, Center for turbulence research: Stanford, 2003. Dukowicz, J. K., “Particle-fluid numerical model for liquid sprays,” Journal of Computational Physics, Vol. 35, 1980, pp.229-253. Okong’o, N. and Bellan, J., “Consistent large eddy simulation of a temporal mixing layer laden with evaporating drops. Part 1: direct numerical simulation, formulation, and a priori analysis,” Journal of Fluid Mechanics, Vol. 499, 2004, pp.1-47. Leboissetier, A., Okong’o, N. and Bellan, J., “Consistent large eddy simulation of a temporal mixing layer laden with evaporating drops. Part 2: A posteriori modeling,” Journal of Fluid Mechanics, Vol. 523, 2005, pp.37- 78. Sankaran, V. and Menon, S., “LES of spray combustion in swirling flows,” Journal of Turbulence, Vol. 3, 2002, pp.1-23. 145 [75] [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] [86] [87] Ham, F., Apte, S. V., laccarino, G., Wu, X., Herrmann, M., Constantinescu, G., Mahesh, K. and Moin, P., “Unstructure LES of reacting multiphase flows in realistic gas-turbine combustors,” Annual Research Briefs, Center for turbulence research: Stanford, 2003. Afshari, A., Almeida, T., Mehravan, K. and Jaberi, F. A., “Large scale simulation of turbulent combustion and propulsion systems,” Proceedings of the Seventeeth ONR Propulsion Meeting, MA, June 2004. Li, 2. And Jaberi, F. A., “A new two-phase stochastic subgrid model for large eddy simulation of spray combustion,” Journal of Fluid Mechanics, to be submitted. 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,” Journal of Fluid Mechanics, Vol. 401, 1999, pp.85-121. James, S. and Jaberi, F. A., “Large scale simulations of two-dimensional nonpremixed methane jet flames,” Combustion and Flame, Vol. 123, 2000, 90465-487. Jaberi, F. A., “Large eddy simulation of turbulent premixed flame via filtered mass density function,” AIAA paper 1999-0199, Jan. 1999. Givi, P., “Filtered density function for subgrid scale modeling of turbulent combustion,” AIAA Journal, Vol. 44, 2006, pp.16-23. Afshari, A., Jaberi, F. A. and Shih, T. l.-P., “Large-eddy simulation of turbulent flow in an axisymmetric dump combustor,” AIAA Journal, in press, 2007. O’Brien, E. E., “The probability density function (PDF) approach to reacting turbulent flow,” Turbulent Reacting Flow, edited by P. A. Libby and F. A. Williams, chapter 5, Springer-Verlag, Heidelberg, 1980, pp.185- 218. Pope, S. B., “PDF method for turbulent reacting flows,” Progress in Energy Combustion Science, Vol. 11, 1985, pp.119-192. Colucci, P. J., Jaberi, F. A., Givi, P., and Pope, S. B., “Filtered density function for large eddy simulation of turbulent reacting flows,” Physics of Fluids, Vol. 10, 1998, pp.499-515. Lundgren, T. 8., “Model equation for nonhomogeneous turbulence,” Physics of Fluids, Vol. 12, 1969, pp.485-497. Pope, S. B., “The probability approach to modeling of turbulent reacting flows,” Combustion and Flame, Vol. 27, 1976, pp.299-312. 146 [88] [89] [90] [91] [92] [93] [94] [95] [95] [97] [98] [99] [100] Germano, M., “Turbulence: the filtering approach,” Journal of Fluid Mechanics, Vol. 238, 1992, pp.325-336. Gao, F. and O’Brien, E. E., “A large-eddy simulation scheme for turbulent reacting flows,” Physics of Fluids, Vol. 5, 1993, pp.1282-1284. Salvetti, M. V. and Banerjee, S., “ A priori tests of a new dynamic subgrid- scale model for finite-difference large-eddy simulations,” Physics of Fluids, Vol. 7, 1995. 002831-2847. Dopazo, C. and O’Brien, E. E., “Statistical treatment of non-isothermal chemical reactions in turbulence,” Combustion Science and Technology, Vol. 13, 1976, 91199-122. Borghi, R., “Turbulent combustion modeling,” Progress in Energy Combustion Science, Vol. 14, 1988, pp.245-292. Miller, R. S., Harstad, K. and Bellan, J., “Evaluation of equilibrium and non- equilibrium evaporation models for many-droplet gas-liquid flow simulations,” International Journal of Multiphase Flow, Vol. 24, 1998, pp. 1 025-1 055. Khatumria, V. C. and Miller, R. 8., “Numerical simulation of a fuel droplet laden exothermic reacting mixing layer,” Intemational Journal of Multiphase Flow, Vol. 29, 2003, pp.771-794. Bellan, J. and Harstad, K., “The details of the convective evaporation of dense and dilute clusters of drops,” International Journal of Heat and Mass Transfer, Vol. 30, 1987, pp.1083-1093. Shu, C. W. and Osher, S., “Efficient implementation of essentially non- oscillatory shock-capturing schemes,” Journal of Computational Physics, Vol. 77, 1988, pp.439-471. Lele, S. K., “Compact finite different schemes with spectral-like resolution,” Journal of Computational Physics, Vol. 103, 1992, pp. 1 6-42. Cook, A. W. and Riley, J. J., “Direct numerical simulation of a turbulent reacting plume on a parallel computer,” Journal of Computational Physics Vol. 129, 1996. 99263-283. Carpenter, M. H, Gottlieb, D. and Abarbanel, S., “The stability of numerical boundary treatments for compact high-order finite-difference schemes,” Journal of Computational Physics, Vol. 108, 1993, pp.272-295. Gaitonde, D. V. and Visbal, M. R., “Pade-type high-order boundary filters for the Navier-Stokes equations,” AIAA Joumal, Vol. 38, 2000, pp.2103- 2112. 147 [101] [102] [103] [104] [105] [106] [107] Visbal, M. R. and Gaitonde, D. V., “Very high-order spatially implicit schemes for computational acoustics on curvilinear meshes,” Journal of Computational Acoustics, Vol. 9, 2001 , pp. 1 259-1 286. Risken, H., The Fokker-Planck Equation, Methods of Solution and Applications, 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, Academies Press, New York, NY, 1981. Pope, S. B., “Lagrangian PDF methods for turbulent flows,” Annual Review of Fluid Mechanics, Vol. 26, 1994, pp.23-63. 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, VA, Aug. 2002 Linck, M., Gupta, A. K., Bourhis, G. and Yu, K., “Combustion characteristics of pressurized swirling spray flame and unsteady two- phase exhaust jet,” AIAA paper 2006-0377, Jan. 2006. 148 CHAPTER 4 Applications of Two-phase LES/FMDF Model to Spray Combustion 4.1 Introduction With the advancements in computational power and numerical algorithms, the high fidelity models such as those based on large eddy simulation (LES) are becoming more popular and are more frequently employed for the development of low-emission, compact and/or efficient combustors. LES can provide detailed time-dependent spatial data for model assessment and improved understanding of turbulent combustion in realistic conditions. Despite great contributions that have been made in the LES of turbulent combustion in recent years, there remain several important challenges in the subgrdi-scale (SGS) modeling and numerical implementation of LES. In particular, modeling and numerical simulations of turbulent spray combustion has proven to be very difficult. A major challenge is to develop affordable models which can properly describe the complicated interactions among turbulence, combustion and spray at a wide range of time and length scales. This chapter discusses some of the results we have recently obtained with our new two-phase LES methodology (as described in chapter 3) for high Reynolds number compressible turbulent reacting flows in a spray—controlled lean premixed dump combustor [1] and non-premixed double-swirl spray burner [2]. The new LES model is implemented via an Eulerian-Lagrangian-Lagrangian numerical scheme [3]. The Eulerian gas-phase part of the scheme employs a generalized hi gh-order finite-difference 149 method and is applicable to turbulent flows in complex geometries. This is based on compact differencing and monotonicity-preserving methods and a multi-block strategy. The subgrid gas-liquid combustion is modeled with a two-phase PDF-based methodology, term the filtered mass density function (FMDF) [4-8]. The two-phase LES/FMDF model is capable of capturing the complex interactions among turbulence, combustion, and spray, and has shown to be applicable to a variety of turbulent flames (slow, fast, nonpremixed, premixed, and partially premixed). The spray is simulated with a Lagrangian model [9-10] which allows full two-way mass, momentum, and energy coupling between phases. For systematic assessment of the two-phase LES/FMDF methodology and its submodels, simulations of low-speed and high-speed, single-phase and two-phase turbulent reacting flows in various flow configurations have been conducted and the generated results have been compared with the experimental and DNS data whenever it is possible. Among the flows that have been simulated, there are double-swirl spray burner [2] and spray-controlled square dump combustors [,1]. These combustion systems are geometrically simple but are relevant to those used in air-breathing propulsion systems. 4.2 Results and Discussions The two-phase LES/FMDF is employed for systematic analysis of turbulent spray combustion in two different combustion systems mentioned above. The performance of the combustor in a typical liquid-fuel combustion system is dependent on the complicated and often coupled effects of various parameters such as input/output operating flow conditions, the geometry, the fuel spray type and the fuel/chemistry 150 characteristics. The effects of some of these parameters on the combustion and turbulence are discussed below. 4.2.1 Double swirl spray burner The schematics of “Enclosed” combustion chamber with co-annular air passages and central fuel injector nozzle for the double swirl spray burner [2] are shown in Figure 4-1. The swirl number,S which characterizes the degree of tangential momentum or swirl in the incoming air, is known to have significant effects in the stability and structure of the flame and is defined by Gupta [11] as: Ge S : ROGa (4.1) where R0 is the outer radius of the swirl. The terms G9 and Ga are the axial flux of angular momentum and the axial flux of axial momentum, respectively, and are defined as G = IROpu u rzdr (4.2) G _ 1R0 2 a — . (pug +P)rdr. (4.3) R1 where an, Mg and R,- are the axial velocity, tangential velocity and inner radius of the swirl, respectively. For the vane angle of 450 /500(inner/outer annulus) used here, the corresponding swirl number, S is 0.8/ 1.0 for the inner/outer annulus. There are obviously many important parameters that affect the combustion in this burner. However, here we only consider the effects of enclosures (wall and endplate with nozzle), droplet to gas mass loading ratio, spray angle and injected droplet size distribution on the combustion. 151 4.2.1.1 Enclosed and Unenclosed Burners Numerical simulations of double-swirl spray burner with and without enclosure are canied out to examine the effects of wall and flame confinement on the flame and turbulent structures. Figure 4-2 shows that the flame in the enclosed case is considerably wider than that in unenclosed case. This is reflected in the temperature distribution and other variables’ behavior. Also, the flame seems to be longer in the enclosed case in comparison to that in the unenclosed case. These observations are quantitatively consistent with the experimental results [1] and the radial plots of temperature obtained from LES/FMDF data as shown in Figure 4-3. Additionally, it is observed that the strength and the width size of the “vorticity layer” in the enclosed case are larger than those in the unenclosed case. The difference between enclosed and unenclosed flames is due to the entrained cold air in the unenclosed case, which cools the flame and restrains its expansion. The end-plate or nozzle also plays an important role in enhancing the combustion, as expected. Again, the numerical results are found to be consistent with the experimental results [2]. Interestingly, Figure 4-4 shows a significant number of survived fuel droplets accumulated in the region close to the exit nozzle, where the turbulence is considerably suppressed by the droplets. The accelerated axial velocity of the outgoing gas in that region also contributes to the damping of turbulence. These are all due to side and end walls which significantly affect the near-field combustion and turbulence. The end-wall or nozzle effects seem to be much more significant. 152 4.2.1.2 Effects of Mass Loading Ratio In the spray combustion, the evaporation of liquid fuel droplets has two major effects on the flow and flame. First, it provides gaseous fuel for chemical reaction and sustains or promotes combustion. Second, the drag, added mass and the cooling effects of droplets on hot carrier gas significantly change the temperature, turbulence and species distributions. These two effects often counter each other in the overall process of spray combustion. In this respect, the mass loading ratio (MLR), which is defined as the ratio of fuel flow rate to air flow ratio, is a critical parameter that not only characterizes the degree of richness/leanness of the flame but also indicates how much the physical structure of the turbulence and the flame are affected by mass, heat and momentum couphngs Figure 4.5 shows that when the temperature of the air at the inlet section (exits of inner/outer swirler) is Tm 2 500K , the flame can be maintained for both MLR=0.04 and MLR=0.1. The gas temperature in the central region is lower in the case with MLR=0.1 as shown in Figure 4-6. When the inlet temperature is decrease to ambient temperature (Tm =300K ), the flame may no longer be maintained and is lifted for larger mass loading ratios as shown in Figure 4-7. This is primarily due to cooling effects of the droplets and also the “richness” of the fuel-air mixture. Note that the momentum (drag) effect of the larger number of droplets causes considerable damping of the turbulence, more so when MLR is higher. The droplet damping (and cooling) effects on the turbulence is important in both reacting and nonreacting flows as the results in the Figure 4-5 and 4-6 suggest. 153 4.2.1.3 Effects of Droplet Size Distribution In this study, two types of droplet size distributions are considered both with the same mean diameter. One is the uniform distribution, and the other is a nearly Gaussian distribution. Figure 4-8 shows that the profiles of droplet temperature in different axial locations are very different for cases with uniform and Gaussian inlet droplet size distribution. As the fuel droplets pass through the flame and evaporate, their temperatures increase and their mass decrease in a rate that is dependent on many factors including the turbulent velocity and temperature structure and the droplet size distribution. In the cases with inlet Gaussian droplet size distribution, the droplet temperature becomes more uniform in the radial direction as they move downstream. Figure 4-9 also indicates that at locations close to inlet, the droplet mass or diameter is almost constant in radial direction for the case with uniform drOplet size distribution. However, further downstream the droplet mass is decreased in the radial direction with a distribution that is very different for uniform and Gaussian inlet droplet size distributions. The profiles of droplet axial velocities (not shown) are qualitatively similar to that of droplet mass or diameter. The axial velocity of the injected droplet with smaller diameter is decreased more by the momentum interaction with the carrier gas. 4.2.1.4 Effects of Spray Angle The instantaneous contours of the gas temperature and vorticity in Figure 4-10 indicate that with smaller spray angle, the flame is longer and narrower in the region close to the inlet. This is confirmed in the Figure 4-11, where the radial variations of the gas temperature for different spray angles are plotted. It is also shown in Figure 4-11 that the flame temperature for the case with higher spray angle is higher that that for the case 154 with smaller spray angle. The effect of spray angle on flow and turbulence is also significant. For example, with the increase of spray angle, the vorticity field grows more and covers a wider region at locations close to inlet. 4.2.2 Spray-Controlled Dump Combustor In the dump combustor experiment [1], a lean premixed air-fuel flame is controlled by injection of a relatively small amount of liquid fuel. We have simulated this experiment with and without (Heptane) fuel droplet evaporation and combustion and with realistic flow and droplet parameters. Sample results are presented in Figure 3-2 in chapter 3, where the contours of the cartier gas vorticity and droplets are shown. The two-phase LES/FMDF results are consistent with the experiment and indicate that the spray tend to decrease the pressure oscillations in the combustor. Expectedly, the cooling of gas due to (fuel) droplet evaporation decreases the local temperature at the location close to droplet injection point. However, the maximum fuel droplet evaporation may or may not occur at high temperature flame zones. Depending on the spray parameters (droplet size, spray and injection angles, injection frequency & duty cycle, etc.) the characteristics of the combustion is noticeably different for different spray parameters. In general, the spray may have “negative” or “positive” influence on the combustion in terms of flame temperature, product species concentration, pressure, etc. The effect on turbulence is also complicated. 4.2.2.1 Effect of Droplet Size The effects of average droplet size or Sauter mean diameter (SMD) of the injected fuel are shown in Figure 4-12 and 4-13, where contours and radial variations of temperature and fuel mass fraction of gas mixture are plotted, respectively. Evidently and 155 expectedly, the overall number of evaporated droplets is increased by employing smaller droplet size. However, the cooling and added fuel mass effects of droplets appear very differently at different locations for different droplet SMD. For small and medium size droplets, the cooling and added fuel mass occur primary in the location close to the inlet. On the other hand, the cooling and added fuel mass of larger droplets take effect at downstream locations, which is due to the fact that evaporation time scale is increased with droplet diameter and decreased with droplet temperature. In the later case, the droplets will have more effects on the turbulence. Note that, the droplet’s evaporation cannot be always described by just the cooling effect. The evaporated fuel does also increase the local equivalence ratio of the lean premixed mixture, leading to higher reaction rate and temperature, depending on the local flow conditions and amount of added evaporated fuel. 4.2.2.2 Effects of Injection Angle and Spray Angle Figure 4-14 shows the effects of injection angle B on the radial variations of temperature, fuel mass fraction and pressure fluctuations in the gas phase at axial locations of x/D=2 and x/D=6. Evidently, for the case with larger injection angle, the peak droplet concentration moves close to inlet, where the cooling and added fuel mass effects have significant effect on combustion and turbulence. However, at the downstream locations, the pressure oscillation is larger for the case with large injection angle. Our results (not shown) also indicate that, as the spray angle or in gradually increased from 30 degree, the droplet distribution becomes more homogeneous at downstream locations, and the overall dr0plet evaporation increases. The evaporated fuel 156 increase the equivalence ratio of the incoming lean premixed mixture, leading to higher reaction rates and temperature at downstream locations. 4.2.2.3 Effects of Injection Frequency and Pulsed Duty Cycle In Figure 4-15, the effects of pulsed spray duty cycle for a fixed injection frequency of 38Hz and the effects of pulsed spray injection frequency for a fixed duty cycle of 50% on the carrier gas pressure fluctuations are shown. The duty cycle is the percentage of the time that the spray is on in each cycle. These results are compared with those obtained. for the case with continuous spray. Due to dampening effect of the droplets on the pressure, the oscillations are less significant when spray is on, which is consistent with the experiment. However, the effect of spray on the pressure field is dependent on the distribution of added fuel. In the cases considered here, the added fuel significantly changes the intensity and the physical structure of the pressure, depending on the injection frequency and the duty cycle. This is clearly shown in Figure 4-15. 4.3 Summary Large eddy simulations of two-phase turbulent reacting flows are conducted via a PDF-based two-phase subgrid combustion model, termed the filtered mass density function (FMDF). The numerical solution of LES/FMDF is based on an Eulerian- Lagrangian-Lagrangian scheme. A high-order finite difference method is used for the simulation of velocity field and the scalars are obtained by the solution of FMDF equation with a Lagrangian Monte Carlo scheme. The spray is simulated with a Lagrangian model which allows two-way mass, momentum and energy coupling between phases. The LES/FMDF can be employed in conjunction with non-equilibrium and equilibrium reaction models and reduced and detailed chemical kinetics 157 mechanisms. LES of two-phase turbulent flames in various flow configurations are conducted and the generated results are employed for detailed analysis of these flames under various operating conditions. The numerical results indicate that the new two- phase LES/FMDF methodology is a reliable and affordable methodology for detailed calculations of turbulent combustion. The effects of droplet size distribution, spray angle, fuel-air equivalent ratio and several other parameters on combustion in two different burners are studied: (1) the double swirl spray burner, (2) the spray-controlled lean premixed combustor. The results clearly indicate the significance of these parameters and the possibility of modifying the combustion by changing them. 158 4.4 Figures Quart]. 4.71 ID 3 .9 l l) Swirlcrs" D - 62 min i 4.161) 1 : H l E Atomization 22 J er Air Flow Outer Air Figure 4-1. A schematic of double swirl-stabilized spray burner and its fuel injector. 159 Temperature Vorticity (a) (b) Figure 4-2. Contours of the instantaneous the filtered temperature and vorticity in the double swirl spray burner: (a) with enclosure and (b) without enclosure. 3 3_ 7E Enclosed 7: cf 63 .552 55; E 3 g : 4_ 4- E E g E “,3: ”3: 1- - 1- 2 2_ 15 12 c: I I I C: I I I 0 0.5 rIR 1 1.5 2 0 0.5 rIR 1 1.5 2 Figure 4-3. Radial variations of the instantaneous filtered temperature at different axial locations for both enclosed and unenclosed double swirl spray burners. (a) (b) Figure 4-4. The 3D instantaneous vorticity lso-surface and fuel droplets' positions in the double swirl spray burner: (a) with enclosure and (b) without enclosure. 161 Temperature Vorticity (a) (b) Figure 4-5. Contours of the instantaneous filtered temperature and vorticity in an unenclosed double-swirl spray burner with different mass loading ratios, MLR and inlet air temperature ofT,-,, = 500K , (a) MLR=0.04 and (b) MLR=0.1. 162 8_ 8. 7:- MLR=O.1 7: g -—o—MLR=0.04 ; 6, 6- 2 . E 35: 5,5: 111 ~ g : 54’ 4_ 3 31 o ~ . 1- ; 11’ g 2_ 2 1: 1: 0:2. 2 . 1 1 I . . 1 1 0: .4. . . I 1 1 . . I_. . . . 0 1 r/R 2 3 0 1 r/R 2 3 Figure 4-6. Radial variations of instantaneous filtered temperature at different axial locations for an unenclosed double swirl spray burner with different mass loading ratios, MLR and inlet air temperature ome = 500K . 163 Temperature Vorticity Figure 4-7. Contours of the instantaneous filtered temperature and vorticity in an unenclosed double-swirl spray burner with different mass loading ratios, MLR and inlet air temperature of T," =300K , (a) MLR=0.04; (b) MLR=0.1; and (c) MLR=0.3. 164 1.7 _ 1.7 I x/R=l.0 _ 1. 1.6 : 6 : 1.5 1.5 : . 1.4 g 1-3 i 1.2 1'3 r 1.1 1.211111LIIIIIJL1L 1.071111lIIIIIIIILlJILI 0 0,2 ”R 0,4 0,6 0,3 0 0.2 rm 0.4 0.6 0.8 1.7 2 g! 1.6 , 1.5 D. t— 1.4 1.3 ller rlllllllllrlllrf looilllllljllIJLIIIIll 1.2111111111141111 0 0,2 r/R 0,4 0.6 0,3 0 0.2 r/R 0.4 0.6 0.8 1'7 - 1.7 : 1.6 E : 12.; 3‘ 1-6 j. Z . ‘. ":31, I" - . 1.5 53". . t .1? 3...? 141” 1.5 - 12* " . , r ' 'thfiz'ff. a. ' : a W“ 9113 :_ 1.4 C XIR=4.0 é XIR=490 I 1.2 :- 1'35 1.15- 102:11LL1L11111111L 1.0:IIIIIIIIIIIILLITI 0 0,2 r/R 0,4 0.6 0,3 0 0.2 r/R 0.4 0.6 0.8 (a) (b) Figure 4-8. Radial variations of the droplet temperature in an unenclosed double swirl spray burner for (a) uniform droplet size distribution and (b) nearly Gaussian droplet size distribution. 165 513-05 _ 4E 05 E x/R=l.0 - i ‘ 313-05 : 2 c Q.‘ _ 213-05 : 113-05 0E.05 :1 L l I l l l l l l l l l l l 0 0.2 ”R 0.4 0.6 0.8 5E-05 _ 315-05 413-05 E =2” ‘. 6E-05 313-05 E E a, E 413-05 213-05 : 113-05 E 2505 OEOOS : l l l 1 l l l l l l l J l 1“: 0E-05 . ‘ .' ‘ ' . . 0 0.2 ”R 0.4 0.6 0.8 0 0.2 0.4 ”110.6 0.8 1.0 513-05 813-05 3: . : x/R=4.0 .1. ,. ‘ , x/R=4.0 415'“ E 6E-05 3.5-.1.“ - 313-05 : h. i . E _ a. 3 .d, 3* 213-05 : .,,-,,..:=.r.1 : . alt-xii: :5: . 113-0 E 1% 9:3 3% 5 : fa“. _ ‘3'": 2. fl" Jug. (a) (b) Figure 4-9. Radial variations of the droplet mass in an unenclosed double swirl spray burner for (a) uniform droplet size distribution and (b) nearly Gaussian droplet size distribution. 166 Vorticity Temperature (a) (b) Figure 4-10. Contours of the instantaneous filtered values of temperature and vorticity in an unenclosed double swirl spray burner with (a) spray angle of a = 200 and (b) spray angle ofa = 60°. 167 Temperature 0 J. 1 l 1 I l 1 l l l l 1 1 l (a) 3 Temperature 'TilfjiiTrT 1 1 l L I l l l l l l t l l 0 10,22 3 (b) Figure 4-11. Radial variations of the instantaneous filtered temperature in an unenclosed double swirl spray burner with different spray anglesa at (a) axial location 0fo R = 2 and (b) axial location 0fo R = 4 . 168 Temperature Fuel mass fraction Figure 4-12. Contours of the instantaneous filtered temperature and fuel mass fraction in a spray-controlled dump combustor for different droplet sizes: (a) SMD=30pm; (b) SMD=45pm; (c) SMD=60pm; and (d) SMD=90pm. 169 10 10 e e 8 5 E - a g 6 g, . E E 4 b 0 Q E- H j 2 . 0.3_ 03 :1 x/D=2.0 8 “9:” .2 a *5 0.2 g 0.2 a . a: . h m V‘ 3 a 0.1 E 0.1 "a; 1 8 1 Ln . E“ ' 0.0,-...111...1.J..1.. 0001 ..111...1. .1 .. 0 0.5 m, 1 1.5 0 0.5 W 1 1.5 Figure 4-13. Radial variations of the filtered temperature and fuel mass fraction of the gas mixture in a spray-controlled dump combustor for different droplet sizes (SMD=30,45,60,90) at two axial locations of x/D=2, 4. 170 Temperature Temperature 0 ”015‘ +201“ 115““ 0.3 _ 0.3 ; x/D=2.0 P 0.2 '. j 0.1 1 Fuel mass fraction Fuel mass fraction L+A P 0 0.5 r/D l 1.5 Figure 4-14. Radial variations of the filtered temperature, the filtered fuel mass fraction and the RMS pressure fluctuations of the gas mixture in a spray- controlled dump combustor for different injection angles ,6 at two axial locations of x/D=2, 6. 171 Heqmancy 38 112 Figure 4-15. Instantaneous lso-surface of the filtered pressure of the gas mixture in the spray-controlled dump combustor for different pulsed spray injection frequencies (11 and duty cycles D. 172 4.5 [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] References Hsu, 0., Pang, B. and Yu, K., “Active control studies for advanced propulsion temperature dependence of critical fuel flux,” Proceedings of the Fourteenth ONR Propulsion Meeting, VA, Aug. 2002 Linck, M., Gupta, A. K., Bourhis, G. and Yu, K., “Combustion characteristics of pressurized swirling spray flame and unsteady two- phase exhaust jet,” AIAA paper 2006-0377, Jan. 2006. Afshari, A., Almeida, T., Mehravan, K. and Jaberi, F. A., “Large scale simulation of turbulent combustion and propulsion systems,” Proceedings of the Seventeeth ONR Propulsion Meeting, MA, June 2004. 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,” Journal of Fluid Mechanics, Vol. 401 , 1999, pp.85-121. James, S. and Jaberi, F. A., “Large scale simulations of two-dimensional nonpremixed methane jet flames,” Combustion and Flame, Vol. 123, 2000, pp.465-487. Jaberi, F. A., “Large eddy simulation of turbulent premixed flame via filtered mass density function,” AIAA paper 1999-0199, Jan. 1999. Sheikhi, M. R. H., Drozda, T. G., Givi, P., Jaberi, F. A. and Pope, S. 8., “Large eddy simulation of turbulent non-premixed piloted methane jet flame (Sandia Flame 0),” Proceedings of the Combustion Institute, Vol. 30, 2005, pp.549-556. Afshari, A., Jaberi, F. A. and Shih, T. |.-P., “Large-eddy simulation of turbulent flow in an axisymmetric dump combustor,” AIAA Journal, in press, 2007. Miller, R. S., Harstad, K. and Bellan, J., “Evaluation of equilibrium and non- equilibrium evaporation models for many-droplet gas-liquid flow simulations,” International Journal of Multiphase Flow, Vol. 24, 1998, pp.1025-1055. Apte, S. V., Gorokhovski, M. and Moin, P., “LES of atomizing spray with stochastic modeling of secondary breakup,” International Journal of Multiphase Flow, Vol. 29, 2003, pp.1503-1522. Gupta, A. K., Lilley, D. G. and Syred, D., Swirl Flows, Abacus Press, Kent, England, 1984. 173 Wlllllljllllljjlllfllllfllllljlljljllllllll“