AIAA-2000-0337 An Overviewof the NCCSpray/Monte-Carlo-PDF Computations M.S.Raju DynacsEngineeringCompany,Inc. Brook Park,OH 38th Aerospace Sciences Meeting & Exhibit January 10-13, 2000 / Reno,NV For permission to copy or republish, contact the American Institute of Aeronautics and Astronautics 1801 Alexander Bell Drive, Suite 500, Reston, VA 20191-4344 1 ! E W An Overview of the NCC Spray/Monte-Carlo-PDF Computations M.S. Raju* Dynacs Engineering Co., Inc., NASA Glenn Research Center 2001 Aerospace Parkway Brook Park, Ohio-44142 1 ABSTRACT by summarizing the results of a 3D test case with periodic boundary conditions. For the 3D case, the This paper advances the state-of-the-art in parallel performance of all the three solvers (CFD, spray computations with some of our recent contri- PDF, and spray) has been found to be good when butions involving scalar Monte Carlo PDF (Proba- the computations were performed on a 24-processor bility Density Function), unstructured grids and par- SGI Origin work-station. allel computing. It provides a complete overview of 2 NOMENCLATURE the scalar Monte Carlo PDF and Lagrangian spray computer codes developed for application with un- A pre-exponent of an Arrhenius structured grids and parallel computing. Detailed reaction rate term comparisons for the case of a reacting non-swirling a non-unity exponent of an Arrh- spray clearly highlight the important role that chem- enius reaction rate term istry/turbulence interactions play in the modeling of an outward area normal vector reacting sprays. The results from the PDF and non- of the nth surface, ms PDF methods were found to be markedly different Bk Spalding transfer number and the PDF solution is closer to the reported exper- b non-unity exponent of an imental data. The PDF computations predict that Arrhenius reaction rate term some of the combustion occurs in a predominantly CD drag coefficient premixed-ftame environment and the rest in a pre- Cv specific heat, J/(kg K) dominantly diffusion-flame environment. However, C¢ a constant in Eq. (39) the non-PDF solution predicts wrongly for the com- cn convection/diffusion coefficient bustion to occur in a vaporization-controlled regime. of the nth face, kg/s Near the premixed flame, the Monte Carlo particle D turbulent diffusion coefficient, m2/s temperature distribution shows two distinct peaks: d droplet diameter, m one centered around the flame temperature and the Ea activation energy of an other around the surrounding-gas temperature. Near Arrhenius reaction rate term the diffusion flame, the Monte Carlo particle tem- h specific enthalpy, J/kg perature distribution shows a single peak. In both Jia diffusive mass flux vector, kg/ms cases, the computed PDF's shape and strength are k turbulence kinetic energy, m2/s 2 found to vary substantially depending upon the prox- I_ latent heat of evaporation, J/kg imity to the flame surface. The results bring to the I_,_]1 effective latent heat of evaporation, fore some of the deficiencies associated with the use J/kg (defined in Eq. (11)) of assumed-shape PDF methods in spray computa- Mi molecular weight of ith tions. Finally, we end the paper by demonstrating species, kg/kg-mole the computational viability of the present solution mk droplet vaporization rate, kg/s procedure for its use in 3D combustor calculations m_o initial mass flow rate associated with kth droplet group, kg/s *Engineering Specialist, Associate fellow AIAA. Copyright Na_ number of time steps employed in (c) 1999 by the author. Published by the AIAA with permission. the PDF time-averaging scheme ¢ represents a set of scalars of the joint PDF Nf number of surfaces contained in _¢ independent composition space a given computational cell p density, kg/m 3 Nffq total number of Monte Carlo o" dimensionality of C-space particles per grid cell r stress tensor term, kg/ms 2 total number of computational cells 0 void fraction nk number of droplets in kth group P pressure, N/m 2 Subscripts Prandtl number P joint scalar PDF f represents conditions associated with fuel gas constant, J/(kg K) g global or gas-phase Re Reynolds number i index for the coordinate or species components rk droplet radius, m j index for the species component rko initial droplet radial location, m k droplet group or liquid-phase 8k droplet radius squaredl r_, m2 l liquid-phase 8talc liquid source contribution of the m conditions associated with Nrn gas-phase continuity equation n nth-face of the computational cell Srnle liquid source contribution of the o initial conditions or oxidizer gas-phase energy equation p conditions associated with the liquid source contribution of the Smlm properties of a grid cell gas-phase momentum equations s represents conditions at the droplet liquid source contribution of the Srnl$ surface or adjacent computational cell gas-phase species equations t conditions associated with time liquid source contribution of the a variable Sc_ index for the scalar component of T temperature, K the joint PDF equation t time, s , partial differentiation with respect ith velocity component, m/s to the variable followed by it ith velocity component of kth droplet group, m/s Superscripts ½ volume of the computational cell, m3 Wo chemical reaction rate, 1/s " Favre averaging gas-phase chemical reaction rate, 1/s - time averaging or average based on the Monte Cartesian coordinate in the ith direction m Carlo particles present in a given cell mass fraction of jth species II fluctuations spatial vector X. mole fraction At local time step used in the 3 INTRODUCTION PDF computationsl _-.......... Spray combustion=is =:of interest in a wide' At/ local time step in the flow solver, s variety of applications: g_turbine combustors, in-: At91 global time step in the spray solver, s At. fuel injection time step, s ternal combustion engines (diesel and spark-ignition), Atml allowable time step in the spray solver, s liquid-rocket motors, and industrial furnaces. A large AV computational cell volume, m3 fraction of the world's energy needs axe met by the 6 Dirac-delta function combustion of the liquid fuels. Because of its abun- dant use, the need for developing better predictive rate of turbulence dissipation, m2/s 3 tools for aiding in the design of ever more efficient, ej species mass fraction at the droplet surface pollution-free, and stable combustion devices has re- species mass fraction at the droplet surface ceived considerable attention. 1-i° re turbulent diffusion coefficient, kg/ms A thermal conductivity, J/(ms K) The physical modeling of turbulent spray flames requires consideration of various complex and rate- dynamic viscosity, kg/ms controlling processes associated with turbulent trans- turbulence frequency, 1/s port, mixing, and chemical kinetics, fluid-dynamic characteristicosffuelinjectionandsprayformation,siderations but also depends on the computational transporat ndvaporizatiocnharacteristicosfindivid- efficiency considerations as determined by the com- ualdropletsa,ndtheinteractionofturbulencewith puter memory and turnaround times afforded by the chemicaklinetics,amongothers.Theinteractionof present-day computers. With the aim of developing turbulencewith chemicaklineticsmayoccurover an efficient solution procedure for use in multidimen- a wide-rangeofdisparatetime andlengthscales.sional combustor modeling, we extended the scalar Turbulencpelaysanimportantrolein determiningMonte Carlo PDF approach to the modeling of sprays theratesofmassandheattransferc,hemicarleac- with parallel computing in order to facilitate large- tions,andliquid-phaseevaporatioinnmanypracti- scale combustor applications/In this approach, the calcombustiodnevicesT.heinfluenceofturbulence mean gas-phase velocity and turbulence fields are de- in agaseoudsiffusionflamemanifestsitselfin sev- termined from the solution of a conventional CFD eralforms,rangingfromtheso-calledwrinkledor method, the scalar fields of species and enthalpy from stretchefdlameletrsegimetothedistributedcombus-a modeled PDF transport equation using a Monte tionregime,dependinugponhowturbulenceinter- Carlo method, and a Lagrangian-based dilute spray actswithvariousflamescale1s.1-12Althoughmost model is used for the liquid-phase representation. oftheturbulentspraycombustiomnodels are based The application of this method showed reasonable on either diffusion or premixed gaseous flame theo- agreement when detailed comparisons were made for ries, combustion in a spray flame is more complex two different cases involving an open and a confined and seldom occurs in a single mode. swirl-stabilized spray flames/ There are several approaches used in the study It is well known that considerable effort usually of turbulent spray flames. 1-1° E1 Banhawy and goes into generating structured-grid meshes for grid- Whitelaw 2 in their study on the prediction of the ding up practical combustor geometries which tend local flow properties in a spray flame employed an to be very complex in shape and configuration. The assumed shape PDF model with values for its mean grid generation time could be reduced considerably and variance obtained from the solution of additional by making use of existing, automated unstructured transport equations. Raju and Sirignano 4 in their grid generators. 1_-1s With the aim of advancing the study on multicomponent spray computations in a current multi-dimensional computational tools used swirl-stabilized center-body combustor employed the in the design of advanced technology combustors, two eddy break-up model of Spalding 13 to take some ef- new computer codes - LSPRAY s and EUPDF 19 - fect of turbuler/ce on combustion. The eddy break-up were developed here, thereby extending our previ- model might be applicable in cases where chemical ous work 7 on the Monte Carlo PDF and sprays to reaction rates are either very fast or slow compared unstructured grids as a part of the National Com- with the turbulent time scales. bustion Code (NCC) activity. NCC is being devel- Most of the turbulence closure models for reac- oped in the form of a collaborative effort between tive flows have difficulty in treating nonlinear reaction NASA GRC, aircraft engine manufacturers, and sev- rates. 11-12 The use of assumed shape PDF methods eral other government agencies and contractors. 17 was found to provide reasonable predictions for pat- The unstructured 3D solver is designed to be mas- tern factors and NOx emissions at the combustor sively parallel and accommodates the use of an un- exit. 14However, their extension to multi-scalar chem- structured mesh with mixed elements comprising of istry becomes quite intractable. The solution proce- either triangular, quadrilateral, and/or tetrahedral dure based on the modeled joint composition PDF type: The ability to perform the computations on transport equation has an advantage in that it treats unstructured meshes allows representation of com- the nonlinear reaction rates without any approxima- plex geometries with relative ease. The application of tion. This approach holds the promise of modeling the unstructured grid extension to a confined swirl- various important combustion phenomena relevant stabilized spray flame provided reasonable agreement to practical combustion devices such as flame extinc- with the available droplet velocity measurements. 9 tion and blow-off limits, and unburnt hydrocarbons A current status of the the use of the paral- (UHC), CO, and NOx predictions. 11-12'14 lel computing in turbulent reacting flows involving The success of any numerical tools used in the sprays, scalar Monte Carlo PDF and unstructured multidimensional combustor modeling depends not grids was described in Ref. 10. It also outlines sev- only on the modeling and numerical accuracy con- eral numerical techniques developed for overcoming someofthehighcomputetrime-and-storalgimeita- the gas-phase in Eulerian coordinates derived for the tionsplacedbytheuseofMonteCarlosolutionmeth- multicontinua approach. 21 This is done for the pur- ods.Theparallepl erformancoefboththePDFand pose of identifying the interphase source terms arising CFDcomputationwsasfoundtobeexcellenbtutthe from the exchanges of mass, momentum, and energy resultsweremixedforthespraymoduleshowingrea- with the liquid-phase. sonablpeerformancoenmassivelpyaralleclomputers likeCrayT3D;butitsperformancweaspooronthe The conservation of the mass leads to: workstationclusters.In ordertoimprovethepar- allel performance of the spray module, two different b yo],,+ = s ,o = nk (1) domain decomposition strategies were developed and k the results from both strategies were summarized. I° The main objective of our present work is to For the conservation of the jth species, we have: to investigate the importance of considering chem- istry/turbulence interactions in the calculation of a reacting spray. This was done by making de- tailed comparisons for the case of a reacting non- swirling spray for which experimental data was re- sm ,= (2) ported by McDonell and Samuelson. 2° The compar- k isons involved predictions from two different sets of where computations, one in which the solution for the tem- perature and species fields is obtained from the use of the scalar Monte Carlo PDF method and in the other they are obtained from the solution of a conventional Z (vj = 0 and Z ej = 1 CFD solution. The second objective is to demon- J J strate the computational viability of the present so- For the momentum conservation, we have: lution procedure for its use in 3D computations. This was done by summarizing the results of a 3D test case which was designed primarily to examine the appli- + + [pv4,.,- cability of the spray particle search algorithm in 3D computations and of the newly-implemented periodic [OV_rij],=j - [(1 - O)Vc_ij],=j = s,,_trn= boundary conditions. For clarity, the overall solution procedure from Ref. 10 is repeated here. However, 4_ for a detailed account of the parallel performance to- k k gether with the development and implementation of the parallel method, the interested reader is referred where _ = the void fraction of the gas which is defined to Ref. 10. as the ratio of the equivalent volume of gas to a given First, complete details of the overall solution volume of a gas and liquid mixture. For dilute sprays, procedure with a particular emphasis on the PDF the void fraction is assumed to be equal to one. The and spray algorithms are presented along with sev- shear stress rq in Eq. (3) is given by: eral other numerical issues related to the coupling between the CFD, spray, and PDF solvers. It is fol- lowed by the results-and-discussion section where the 2 vii = I.t[ui,rj + u.i,x,] -- "_6ijui,x_ application of the method to predict the local prop- erties of two different cases are presented. Finally, we For the energy conservation, we have: conclude the paper by presenting a brief summary of important results along with the parallel performance of the three solvers (CFD, PDF, and spray). [_Vch],,+_SVcuih],x,-[OVc)_T,_,],_,-[(1-8)VcAtT,:,],_, 4 GOVERNING EQUATIONS FOR THE GAS PHASE = (4) k Here, we summarize the conservation equations for 5 GAS-PHASE SCALAR JOINT PDF The second term on the right-hand side repre- EQUATION sents transport in the scalar space due to molecu- lar mixing. A mathematical description of the mix- The transport equation for the density-weighted joint ing process is rather complicated, and the interested PDF of the compositions, P, is: reader is referred to Ref. 12. Molecular mixing is accounted for by making use of the relaxation to the ensemble mean submodel. 11 _],, + [paip'],., {Transient} {Mean convection} < pjai,., I¢_ >= -c#_(¢_ - $_) (7) + [p_(c__)p'],¢o= -[p < u71¢_> _,., {Chemical reactions} {Turbulent convection} where w = e/k, and C¢ is a constant. For a conserved scalar in a homogeneous turbulence, this I -[P<_ 1jai,=,l¢->- p-],,o-[_< ps_lC->p-],¢. (5) model preserves the PDF shape during its decay, but there is no relaxation to a Gaussian distribution} 2 {Molecular mixing} {Liquid-phase contribution} However, the results of Ref. 14 indicate that the choice between the different widely-used mixing mod- where els is not critical in the distributed reaction regime of = chemical source term Wa premixed combustion as long as the turbulent mixing for the c_-th frequencies are above I000 Hz. Most of the practical composition variable, combustors seem to operate at in-flame mixing fre- < u_II ]¢_.> = conditional average of quencies of 1000 Hz and above. The application of Favre velocity this mixing model seemed to provide some satisfac- fluctuations, tory results when applied to flows representative of 1 ct < ;J_,_, [¢ > = conditional average those encountered in the gas-turbine combustion} 4 of scalar dissipation, The third term on the right-hand side represents and the contribution from the spray source terms: 1 < ;s_ ]¢ > = conditional average of spray source terms. The terms on the left-hand side of the above 1 1 < -so 1¢ >= -- _'_-kmk(¢o. - ¢_) (s) equation could be evaluated without any approxima- p - pAY tion, but the terms on the right-hand side of the equa- where ¢o = y_,a= 1,2,...,s= o'- 1 tion require modeling. The first term on the right represents transport in physical space due to turbu- 1 1 lent convection. 12 Since the joint PDF, /3, contains < -s_ I¢ >= _ _ n_m_(-l_,_jj + h_. - ¢_) no information on velocity, the conditional expecta- p -- pAy (9) tion of < uIiI J¢ > needs to be modeled. It is mod- eled based on a gradient-diffusion model with infor- where ¢_ = h and is defined by: mation supplied on the turbulent flow field from the a-1 flow solver} 2 h = _,ihi (10) i=l - < u7I¢ > P= r_p,_, (6) where The fact that the turbulent convection is modeled as a gradient-diffusion makes the turbulent model no better than the k - e model. The uncertainties asso- hi = h°li + CpiyidT, rtJ ciated the use of a standard k- cturbulence model to swirling flows are well known32 Some of the modeling }:_tt uncertainties associated with the use of the standard Cpi = "_i(Ati + A2iT + A3iT 2+ A4iT 3+ AhiT4), k - cmodel would be addressed in our future studies with the implementation of a non-linear k - ¢ devel- h°li is the heat of formation of ith species, R_ is the oped for the modeling of swirling flows} 2 universal gas constant, cas is a mass fraction of the evaporatinsgpecieastthedroplet surface, and lk,,]: is the effective latent heat of vaporization as modified c9 = _ 1+ (15) by the heat loss to the droplet interior: For droplet size, the droplet regression rate is (11) determined from one of three different correlations depending upon the droplet-Reynolds-number range. When Rek > 20, the regression rate is determined Here we assumed that the spray source terms based on a gas-phase boundary-layer analysis valid could be evaluated independent of the fluctuations in for Reynolds numbers in the intermediate range. 23 the gas-phase compositions of species and enthalpy. The other two correlations, valid when Rek < 20, are Eqs. (8)-(10) represent the modeled representation taken from Clift et al. _4 for the conditional averages of the spray contribution to the PDF transport equation. ds_ -2 I__[[[2;ne,]j1/2 f(B,) ;fR. >20 -E= Pk 6 LIQUID-PHASE EQUATIONS ds___£ The spray model is based on the multicontinua ap- dt = _Pp._kL [1+ (1+ Rek) '13] Re°£°771n(1 + Bk) proach which allows for resolution on a scale greater if 1 < Rek <_20 (16) than the average spacing between two neighboring droplets. A Lagrangian scheme is used for the liquid- -d-sk = ---#_ [1+ (1+ Re_) 1/3] ln(l+Bk) if Rek < phase equations as it eliminates errors associated dt pk with numerical diffusion. The vaporization modelofa where Bk is the Spalding transfer number defined in polydisperse spray takes into account the transient ef- Eq. (22). The function f(Bk) is obtained from the fects associated with the droplet internal heating, the solution of Emmon's problem. 2s The range of validity forced convection effects associated with droplet in- of this function was extended in l_ju and Sirignano 4 ternal circulation and the phenomena associated with to consider the effects of droplet condensation. boundary layers and wakes formed in the intermedi- The internal droplet temperature is determined ate droplet Reynolds number range. 4The present for- based on a vortex model. 23 The governing equation mulation is based on a deterministic particle-tracking for the internal droplet temperature is given by: method and on a dilute spray approximation which is applicable for flows where the droplet loading is >,, low. Not considered in the present formulation are + the effects associated with the droplet breakup, the "3F = 17_c--,-,=_--p__;_,-,[_ _ (1 + C(Oa ) TEJ droplet/shock interaction, the multi-component na- (17) ture of liquid spray and the phenomena associated with dense spray effects and super-criticM conditions. where The spray method provided some favorable results (18) when applied to several unsteady and steady-state rk-_-- calculations.4- lo For the particle position of the kth droplet where a represents the coordinate normal to the group, we have: streamsurface of a Hill's Vortex in the circulating fluid, and C(t) represents a nondimensional form of dxik the droplet regression rate. The initial and boundary = uik (12) dt conditions for Eq. (17) are given by: For the droplet velocity: "-- tinjection, T;, = T,Lo (19) dui_ 3 CDpa,Rek (13) 1 [c,,p,] dt - 16 pkr_ [ui9 -- Uik] 0__ 0_ = l"ff [--_--_ J r_ at where OT_ O:=1, O_ Rek =2_--_-p[(u9- _,k).(u9- uk)]1/2 (14) (21) I-tas where a = 0 refers to the vortex center, and c_ = 1 grid spacing and average droplet spacing and velocity. refers to the droplet surface. However, our experience showed that for the case of a The Spalding transfer number is given by: steady-state solution, a time step based on the aver- age droplet lifetime yields better convergence. 4-s Its - - (Y:"- Y:) (22) value typically ranges between 1 and 2 milli-seconds Ik,,:: (1- y:,) for the case of reacting flows. The spray computations facilitate fuel injection y-:: ,= l+Ma( x:-1,- :) through the use of a single fuel injector comprised of different holes, s-s However, multiple fuel injection in where Ma is the molecular weight of the gas excluding a steady-state calculation could be simulated by sim- fuel vapor. ply assigning different initial conditions for the spatial Based on the assumption that phase equilibrium locations of the droplet groups associated with each exists at the droplet surface, the Clausius-Clapeyron one of the different holes. For a polydisperse spray, relationship yields the spray computations require inputs for the num- ber of droplet groups in a given stream and for the Pn :h initial droplet locations and velocities. However, the number of droplets in a given group and their sizes In Eq. (14), the molecular viscosity is evaluated at a could be either input directly or computed from a reference temperature using Sutherland's equation properly chosen function for the droplet size distribu- tion. The specified initial inputs should be represen- tative of the integrated averages of the experimental T3/2 conditions:- :0 pg,(T_]) = 1.4637 10-6 "'_1 (25) T,_: + 120 One correlation typical of those used for the droplet size distribution is taken from Ref. 2: where T,,:= [Tg+ _T_, (26) The droplets may evaporate, move along the [ d ]3.5 , _ ,o.4dd wall surfaces, and/or reflect with reduced momen- dnn = 4.21 IOs [d-'_32J e-:6"gs(s_7) --d32 (27) tum upon droplet impingement with the combustor walls. In our present computations, subsequent to the droplet impingement with the walls, the droplets where n is the total number of droplets and dn is are assumed to flow along the wall surfaces with a the number of droplets in the size range between d velocity equal to that of the surrounding gas. and d-4-dd. The Santer mean diameter, d32, could be either specified or estimated from the following 7 DETAILS OF DROPLET FUEL correlation: 26 INJECTION The success of any spray model depends a great deal on the specification of the appropriate injector 2_roz_. (28) d32 -" Bd pg-'-_ m exit conditions. However, a discussion involving the physics of liquid atomization is beyond the scope of this subject matter. In our present computations, the where Bu is a constant, VT is the average relative liquid fuel injection is simulated by introducing a dis- velocity between the liquid interface and the ambi- cretized parcel of liquid mass in the form of spherical ent gas, and _* is a function of the Taylor number, droplets at the beginning of every fuel-injection time step. For certain cases, the fuel-injection time step, A typical droplet size distribution obtained Atiz, needs to be determined based on the resolution from the above correlation in terms of the cumulative permitted by the length and time scales associated percentage of droplet number and mass as a function with several governing parameters such as average of the droplet diameter is shown in Fig. 1.7 22 where subscript n refers to the nth-face of the com- 1.e putational cell. The coefficient c,_ represents the _g _ _ Number / c .9 o transport by convection and diffusion through the 18 nth-face of the computational cell, p. The convec- _J 16 .o tion/diffusion coefficients in the above equation are .7 '[ determined by one of the following two expressions: "- 12 2a n ._ .5 _ E In -o B .4 E 2a_n.an . 6 .3 c, = rnaz[10.5fia_ntU_n I,I'¢(A_-p _V ) ] - 0.5_an.U_a 4 .2 O" and ,-1 2 .1 J Cp "-- Z Cn J e | In both the above expressions for cn, a cell-centered finite-volume derivative is used to describe the vis- Drop|et diameter, mlc_-ons cous fluxes; but an upwind differencing scheme is used for the convective fluxes in the first expression and a hybrid differencing scheme in the second. Figure 1 Droplet-size distribution. 9.1 Numerical Method Based on 8 CFD SOLUTION ALGORITHM Approximate Factorization The gas-phase mass and momentum conservation The transport equation is solved by making use of an equations together with the standard k-e turbulence approximate factorization scheme) 2 Eq. (29) can be equations with wall functions are solved by making recast as: of a finite-volume solver with an explicit fourth-stage Runge-Kutta scheme. Further details of the code can _,,(¢__,+ At)= be found in Refs. 16 & 18. (I+AtR)( (I+AtS)( I+AtM)( I+AtT)_p (¢, _)+O(At 2) 9 PDF SOLUTION ALGORITHM (30) In order to facilitate the integration of the Monte where I represents the unity operator and T, M, S, Carlo PDF method in a finite-volume context, the and R denote the operators associated with spatial volume integrals of convection and diffusion in Eq. transport, molecular mixing, spray, and chemical re- (5) were first recast into surface integrals by means of actions, respectively. The operator is further split a Gauss's theorem. -_7Partial integration of the PDF into a sequence of intermediate steps: transport equation would yield: _(¢_,t) = (I + AtT)fiv(C_.,t ) (31) cpAt ifip(__, t + At) = (I - #--_)ifip (__,t) _v (C__,t) ----(I + AtM)_(C_,t) (32) c_ At + _ _(_¢,_)- A@o(c_)v'],¢o Tt p_(¢_., t) = (I + AtS)W (¢__,t) (33) - A_[<1S,_,=#,1_¢> p-],_:- A_[<}_ I_¢> p-3,¢_ (29) ,_p(¢_., t + At) = (I + AtR)p"_(¢_., t) (34)