ANSYS Forte Theory Manual 18.2 [PDF]

Citation preview

Forte Theory Manual



ANSYS, Inc. Southpointe 2600 ANSYS Drive Canonsburg, PA 15317 [email protected] http://www.ansys.com (T) 724-746-3304 (F) 724-514-9494



Release 18.2 August 2017 ANSYS, Inc. and ANSYS Europe, Ltd. are UL registered ISO 9001: 2008 companies.



Copyright and Trademark Information © 2017 ANSYS, Inc. Unauthorized use, distribution or duplication is prohibited. ANSYS, ANSYS Workbench, AUTODYN, CFX, FLUENT and any and all ANSYS, Inc. brand, product, service and feature names, logos and slogans are registered trademarks or trademarks of ANSYS, Inc. or its subsidiaries located in the United States or other countries. ICEM CFD is a trademark used by ANSYS, Inc. under license. CFX is a trademark of Sony Corporation in Japan. All other brand, product, service and feature names or trademarks are the property of their respective owners. FLEXlm and FLEXnet are trademarks of Flexera Software LLC.



Disclaimer Notice THIS ANSYS SOFTWARE PRODUCT AND PROGRAM DOCUMENTATION INCLUDE TRADE SECRETS AND ARE CONFIDENTIAL AND PROPRIETARY PRODUCTS OF ANSYS, INC., ITS SUBSIDIARIES, OR LICENSORS. The software products and documentation are furnished by ANSYS, Inc., its subsidiaries, or affiliates under a software license agreement that contains provisions concerning non-disclosure, copying, length and nature of use, compliance with exporting laws, warranties, disclaimers, limitations of liability, and remedies, and other provisions. The software products and documentation may be used, disclosed, transferred, or copied only in accordance with the terms and conditions of that software license agreement. ANSYS, Inc. and ANSYS Europe, Ltd. are UL registered ISO 9001: 2008 companies.



U.S. Government Rights For U.S. Government users, except as specifically granted by the ANSYS, Inc. software license agreement, the use, duplication, or disclosure by the United States Government is subject to restrictions stated in the ANSYS, Inc. software license agreement and FAR 12.212 (for non-DOD licenses).



Third-Party Software See the legal information in the product help files for the complete Legal Notice for ANSYS proprietary software and third-party software. If you are unable to access the Legal Notice, contact ANSYS, Inc. Published in the U.S.A.



Table of Contents 1. Forte Theory Manual Introduction ......................................................................................................... 1 2. Basic Governing Equations ..................................................................................................................... 3 2.1. Conservation Equations for Turbulent Reacting Flow ......................................................................... 3 2.1.1. Species Conservation Equation ................................................................................................. 4 2.1.2. Fluid Continuity Equation ......................................................................................................... 4 2.1.3. Momentum Conservation Equation .......................................................................................... 4 2.1.4. Energy Conservation Equation ................................................................................................. 4 2.1.5. Gas-phase Mixture Equation of State ........................................................................................ 5 2.2. Turbulence Models ........................................................................................................................... 5 2.2.1. Reynolds-Averaged-Navier-Stokes (RANS) Approach ................................................................. 5 2.2.2. Large-Eddy Simulation Approach ............................................................................................. 7 2.3. Chemical Kinetics Formulation .......................................................................................................... 9 2.3.1. Optional Semi-empirical Soot Model ...................................................................................... 10 2.4. Turbulence-Kinetics Interaction Model ............................................................................................ 10 3. Boundary Conditions ............................................................................................................................ 13 3.1. Wall Boundaries .............................................................................................................................. 13 3.1.1. Wall Conditions for the Gas-phase Species and Fluid Continuity Equations .............................. 13 3.1.2. Wall Conditions for the Momentum Equation .......................................................................... 13 3.1.3. Wall Conditions for the Energy Equation ................................................................................. 14 3.1.4. Turbulence Model Wall Conditions .......................................................................................... 15 3.2. Inflow and Outflow Boundaries ....................................................................................................... 16 3.2.1. Velocity Inflow Boundaries ..................................................................................................... 16 3.2.2. Outflow Boundaries ............................................................................................................... 16 3.2.3. Pressure Inflow Boundaries ..................................................................................................... 16 3.2.4. Pressure Outflow Boundaries .................................................................................................. 16 3.3. Periodic Boundaries ........................................................................................................................ 17 3.4. Axis-of-Symmetry Boundaries ......................................................................................................... 17 4. Initial Conditions ................................................................................................................................... 19 4.1. Initialization of the Fluid Properties ................................................................................................. 19 4.2. Swirl Profiles ................................................................................................................................... 19 5. Discretization and Solution Methods ................................................................................................... 21 5.1. Discretization of the Governing Equations ....................................................................................... 21 5.1.1. Temporal Differencing Method ............................................................................................... 21 5.1.2. Spatial Differencing Method ................................................................................................... 21 5.2. SIMPLE Method .............................................................................................................................. 23 5.2.1. Convective Flux Discretization ................................................................................................ 24 5.3. Chemistry Solver Options ................................................................................................................ 25 5.3.1. Operator Splitting Method and Parallel Implementation ......................................................... 25 5.3.2. Dynamic Adaptive Chemistry ................................................................................................. 26 5.3.3. Dynamic Cell Clustering ......................................................................................................... 28 6. Spray Models ......................................................................................................................................... 29 6.1. Solid-Cone Spray Models ................................................................................................................ 29 6.1.1. Nozzle Flow Model ................................................................................................................. 29 6.1.1.1. Discharge coefficient ..................................................................................................... 30 6.1.1.2. Effective Injection Velocity and Effective Flow Exit Area .................................................. 31 6.1.1.3. Spray Angle ................................................................................................................... 32 6.1.1.4. Empirical Nozzle Discharge Coefficient and Spray Angle ................................................. 32 6.1.2. Kelvin-Helmholtz / Rayleigh-Taylor Breakup Model ................................................................. 32 6.1.2.1. Kelvin-Helmholtz Breakup ............................................................................................. 33 6.1.2.2. Rayleigh-Taylor Breakup ................................................................................................ 35 Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



iii



Theory Manual 6.1.3. Unsteady Gas-Jet Model ......................................................................................................... 36 6.2. Hollow-Cone Spray Models ............................................................................................................. 38 6.2.1. Linearized Instability Sheet Atomization (LISA) Model ............................................................. 39 6.2.1.1. Film Formation .............................................................................................................. 39 6.2.1.2. Sheet Breakup ............................................................................................................... 40 6.2.1.3. Atomization .................................................................................................................. 41 6.2.2. Taylor-Analog-Breakup Model ................................................................................................ 41 6.3. Droplet Collision and Coalescence Model ........................................................................................ 42 6.3.1. Radius of Influence (ROI) Collision Model ................................................................................ 42 6.3.2. Collision Mesh Model ............................................................................................................. 43 6.4. Multi-Component Droplet Vaporization Model ................................................................................ 43 6.4.1. Liquid Phase Balance Equation ............................................................................................... 44 6.4.2. Governing Equations for Gas Phase ........................................................................................ 44 6.4.3. Vapor-Liquid Equilibrium ........................................................................................................ 45 6.4.4. Determination of Surface Temperature ................................................................................... 45 6.4.5. Modeling Boiling, Including Flash Boiling ................................................................................ 46 6.5. Wall Impingement Model ................................................................................................................ 46 6.5.1. Mass Fraction of Secondary Droplets ...................................................................................... 47 6.5.2. Size of Secondary Droplets ..................................................................................................... 47 6.5.3. Velocity of Secondary Droplets ............................................................................................... 48 6.6. Wall Film Model .............................................................................................................................. 49 6.6.1. Wall Film Governing Equations ............................................................................................... 49 6.6.2. Wall Separation Criterion ........................................................................................................ 50 7. Turbulent Flame Propagation Model .................................................................................................... 51 7.1. Discrete Particle Ignition Kernel Flame Model .................................................................................. 51 7.2. Autoignition-Induced Flame Propagation Model ............................................................................. 52 7.3. G-equation Model ........................................................................................................................... 52 7.3.1. General Formulation .............................................................................................................. 52 7.3.2. Laminar and Turbulent Flame Speed Correlations ................................................................... 53 7.3.2.1. Laminar Flame Speeds ................................................................................................... 53 7.3.2.1.1. Power-law Formulation ......................................................................................... 53 7.3.2.1.2. Laminar Flame Speed at Reference State ............................................................... 54 7.3.2.1.3. Temperature and Pressure Dependencies .............................................................. 54 7.3.2.1.4. Diluent Effect ........................................................................................................ 55 7.3.2.1.5. Laminar Flame Speed Through Table Look-up ........................................................ 55 7.3.2.2. Turbulent Flame Speeds ................................................................................................ 55 7.3.2.3. Turbulent Flame Brush Thickness ................................................................................... 56 7.3.3. Flame-front Heat Release Calculation ...................................................................................... 56 7.3.4. Post-flame and End Gas Kinetics ............................................................................................. 57 7.3.5. Flame Quench Model ............................................................................................................. 58 8. Method of Moments .............................................................................................................................. 59 8.1. Description and Properties of a Particle Population .......................................................................... 59 8.1.1. Moments of Particle-Size Distribution Functions ..................................................................... 59 8.1.2. Total Particle Number of a Particle Population ......................................................................... 60 8.1.3. Total and Average Particle Mass .............................................................................................. 60 8.1.4. Total and Average Geometric Properties of a Particle Population ............................................. 61 9. Engine Crevice Model ............................................................................................................................ 63 9.1. Diagram and Assumptions .............................................................................................................. 63 9.1.1. Isentropic ............................................................................................................................... 64 9.1.2. Isothermal ............................................................................................................................. 64 9.2. Input Parameters ............................................................................................................................ 65 10. References ........................................................................................................................................... 67



iv



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 1: Forte Theory Manual Introduction The ANSYS Forte CFD Package is designed for internal combustion engine design applications. The modeling approach described in this Theory manual is therefore tailored to provide the most accurate solutions possible for these applications, using computing resources that are practical for every-day design activities. ANSYS Forte takes advantage of well-established theoretical representations of 3-D fluid flow, spray dynamics and combustion behavior. The dynamics of spray combustion in diesel engines are controlled by both turbulent mixing dynamics and fuel combustion kinetics. Spray dynamics and fuel vaporization are typically the dominant contributors to the creation of stratified fuel/air mixtures. Subsequently, the ignition and combustion chemical kinetics are controlled by a complex network of reactions between fuel and air species under these stratified conditions. Both spray dynamics and chemical kinetics contribute to source terms in the reacting-flow transport equations. For flame propagation in spark-ignition engine, combustion kinetics in the end-gas also controls important phenomena, such as engine knocking and emissions production. For advanced-concept engines that are based on compression-ignition strategies or dual-fuel combustion, chemical kinetics becomes even more important. The ANSYS Forte CFD Package introduces important breakthroughs in chemistry-solution techniques that greatly enhance the accuracy achievable by engine simulation within commercial design time lines. These techniques reduce simulation time by as much as two orders of magnitude when compared to conventional CFD. Chemistry models that were previously thought of as only practical for 0-D simulations now become practical for full 3-D engine simulations with moving pistons and valves. Better handling of chemistry with multi-component fuel representation makes predictive simulation possible within the schedule constraints of the concept phase of design. For direct-injection engines, ANSYS Forte also sets a new standard for accuracy in the representation of fuel-spray droplet breakup and vaporization. For example, true multi-component fuel vaporization models can now work hand-in-hand with multi-component chemistry models. In addition, new approaches to representing the gas-transport part of the spray reduce the sensitivity of the spray model to the grid. These advances together with better chemistry allow more accurate simulations at reduced computational cost, and without the intervention of expert calibration. ANSYS Forte builds on models and sub-models that have been well validated against experimental data over a broad range of conditions and over many years by engine-simulation experts. This manual describes the model assumptions and solution techniques employed. The ANSYS Forte User's Guide provides end-user instructions for model inputs and for the mechanics of using the ANSYS Forte Simulate and Visualizer interfaces.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



1



2



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 2: Basic Governing Equations 2.1. Conservation Equations for Turbulent Reacting Flow In internal combustion engines, the fuel-air mixture before combustion and the burned products after combustion are the working fluids. In the turbulent reacting flow representation of ANSYS Forte, the basic fluid dynamics are governed by the Navier-Stokes equations. Model transport equations of mass, momentum, and energy conservation laws are formulated for the compressible, gas-phase flows, and represent the turbulent nature of the flow. When liquid sprays are injected into the flow, exchange functions are used to account for the interaction between the gas phase and the liquid droplets. Beyond these models, the main assumptions made in the derivation of the governing equations are the use of the ideal gas law for the gas-phase equation of state, the use of Fick’s law for mass diffusion, the assumption of Newtonian fluid, and the use of Fourier’s law for thermal diffusion. In this chapter, we present the governing equations for the gas phase of the working fluids in the engine. The model equations for the liquid sprays are presented in Spray Models (p. 29). Turbulent flow is characterized by a wide range of flow length scales as well as significant and irregular variations of the flow field. Forte offers two options of turbulence modeling. The first option is the Reynolds-Averaged-Navier-Stokes (RANS) approach, which aims at capturing the ensemble average of the flow field from many realizations of flows under equivalently set conditions. Since an important effect of turbulence is more effective transport and mixing of fluid compared to a laminar flow, the ensemble average of the turbulent transport and mixing is analogous to a large-scale diffusion. The RANS approach removes the necessity of resolving small-scale structures and fluctuations seen in individual flow realizations, while retaining the main effects of turbulence on the averaged flow and combustion characteristics. To accomplish this, the Favre average is employed to represent an instantaneous quantity, such as the flow velocity vector , into an ensemble average and a fluctuating part , as . In this approach, the average part is defined as a conventional density-weighted average by , while the fluctuation is defined to satisfy , where the over-bar represents an averaging operator. The second option is the Large-eddy Simulation (LES) approach, which simulates individual flow realizations instead of the ensemble average of the flows. LES is expected to capture a substantial degree of unsteadiness and smaller flow length scales, which are omitted in the ensemble average. Due to the restriction on the computational cost, only the larger three-dimensional unsteady turbulent motions are directly resolved by the mesh resolution, whereas the effects of the smaller-scale motions are modelled as “sub-grid scale" (SGS) effects. In this approach, the resolved flow field is regarded as a result of filtering the actual flow field, with the filter size being the mesh size. The filtering eliminates the need to resolve the flow scales smaller than the filter size, which is accounted for by the SGS models. In general, the smaller the filter (mesh) size is, the better the scale-resolution of the turbulent flow field that is expected. In the LES approach, the instantaneous quatity of the flow field (e.g., ), is decomposed into a filtered component and a resicual (or SGS) component, also denoted by and , as . In this approach, Again, the filtered part is defined by the Favre averaging while the fluctuation is defined to satisfy , where



, where is the filtered density field. is the filtered density field.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



3



Basic Governing Equations The governing equations in ANSYS Forte are formulated to solve the ensemble-averaged flow field in the RANS approach and to solve the filtered flow field in the LES approach. The basic form of the equations for the flow field is presented in a unified way as follows, while the turbulence models themselves are described in Turbulence Models .



2.1.1. Species Conservation Equation The gas-phase working fluids in the combustion engine are modeled as a mixture of individual gas components, or species, and this composition changes during the engine cycle due to flow convection, molecular diffusion, turbulent transport, interactions with fuel sprays, and combustion. The conservation equation for the mass of species is: (2.1) where is the density, subscript is the species index, is the total number of species, is the flow velocity vector. Application of Fick’s Law of diffusion results in a mixture-averaged molecular diffusion coefficient . The Φ term accounts for the effects of ensemble-averaging or filtering of the convection term, that is,



, which has to be modeled.



and



are source terms due to chemical reac-



tions and spray evaporation, respectively.



2.1.2. Fluid Continuity Equation The summation of Equation 2.1 over all species gives the continuity equation for the total gas-phase fluid: (2.2)



2.1.3. Momentum Conservation Equation The momentum equation for the fluid considers the effects of convection, pressure force, viscous stress, and turbulent transport, as well as the impact from liquid sprays and body force: (2.3) where p is the pressure, is the rate of momentum gain per unit volume due to the spray, g is the specific body force, is the viscous shear stress given by (2.4) in which is the laminar kinematic viscosity, I is the identity tensor, and superscript T means transpose of a tensor. The stress    accounts for the effects of ensemble-averaging or filtering of the nonlinear convection term, that is, . In the RANS approach, it is called the Reynolds stress; in the LES approach, it is called the SGS stress. In both cases, turbulence models are needed to provide closure.



2.1.4. Energy Conservation Equation Based on the First Law of Thermodynamics, the change of internal energy has to be balanced by the pressure work and heat transfer. For the flow problems relevant to internal combustion engines, the effects of convection, turbulent transport, turbulent dissipation, sprays, chemical reactions, and enthalpy diffusion of a multi-component flow should also be considered. The internal energy transport equation is: 4



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Turbulence Models (2.5) where I is the specific internal energy, J is the heat flux vector accounting for contributions due to heat conduction and enthalpy diffusion: (2.6) λ is the thermal conductivity, which is related to the thermal diffusivity α and heat capacity by , T is the fluid temperature, and h k is the specific enthalpy of species k. is the dissipation rate of the turbulent kinetic energy, which will be defined in Turbulence Models . and are source terms due to chemical heat release and spray interactions, respectively. The term accounts for the effects of ensemble-averaging or filtering of the convection term, that is, . Again, it needs to be properly modeled from the turbulence approach.



2.1.5. Gas-phase Mixture Equation of State The thermodynamic state relations for the gas-phase mixture are assumed to obey the ideal gas law. The mixing of gas components is assumed to follow the Dalton model, that is, each component behaves as an ideal gas as if it were alone at the temperature and the volume of the mixture, (2.7) where R



u



is the universal gas constant, and W k is the molecular weight of species k.



2.2. Turbulence Models We now describe the turbulence models needed to provide closure to the undetermined terms,   ,   , and H, in Equation 2.1, Equation 2.3, and Equation 2.5 , respectively.



2.2.1. Reynolds-Averaged-Navier-Stokes (RANS) Approach As introduced at the beginning of this chapter, the RANS approach aims aims to simulate the ensembleaveraged flow field. The most widely used approach is to model the turbulent transport processes with gradient-diffusion assumptions. For the momentum equation, the deviatoric components of the Reynolds stress are assumed to be proportional to the mean deviatoric rate of strain. The Reynolds stress tensor   , is defined as (2.8) in which



is the turbulent kinematic viscosity, and



is the turbulent kinetic energy, defined by: (2.9)



The turbulent viscosity



is related to the turbulent kinetic energy



and its dissipation rate



by: (2.10)



where is a model constant that varies in different turbulence model formulations, shown in Table 2.1: Constants in the standard and RNG k - ε models . Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



5



Basic Governing Equations The turbulent flux term    in the species transport Equation 2.1 is modeled as: (2.11) in which is the turbulent diffusivity. Similarly, the turbulent flux term H in the energy Equation 2.5 is modeled as: (2.12) in which is the turbulent thermal conductivity and is related to the turbulent thermal diffusivity and heat capacity by . The turbulent mass and thermal diffusivity are related to the turbulent viscosity by: (2.13) (2.14) where and are the turbulent Schmidt and Prandtl numbers, respectively. As seen in Equation 2.10 , the calculation of turbulent viscosity requires that the turbulent kinetic energy and its dissipation rate to be modeled. In ANSYS Forte, both the standard and the advanced (based on Re-Normalized Group Theory) k-ε model formulations are available. These consider velocity dilatation in the ε- equation and spray-induced source terms for both k and ε equations. The standard Favre-averaged equations for k and E are given in Equation 2.15 and Equation 2.16 : (2.15)



(2.16)



In these equations, , , , , are model constants, which are listed and described in Table 2.1: Constants in the standard and RNG k - ε models . The source terms involving



are calculated based on the droplet probability distribution function



(cf. Ref.Amsden 1997 5 ). Physically, is the negative of the rate at which the turbulent eddies are doing work in dispersing the spray droplets. was suggested by Amsden 5 based on the postulate of length scale conservation in spray/turbulence interactions. The advanced (and recommended) version of the k - ε model is derived from Re-Normalized Group (RNG) theory, as first proposed by Yakhot and Orszag 78 . The k equation in the RNG version of the model is the same as the standard version, but the ε equation is based on rigorous mathematical derivation rather than on empirically derived constants. The RNG ε equation is written as (2.17)



where the



6



in the last term of the right-hand side of the equation is defined as



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Turbulence Models



(2.18) with (2.19) (2.20) and



is the mean strain rate tensor, (2.21)



Compared to the standard ε equation, the RNG model has one extra term, which accounts for nonisotropic turbulence, as described by Yakhot and Orszag 78 . Values of the model constants , , , and used in the RNG version are also listed in Table 2.1: Constants in the standard and RNG k - ε models . In the ANSYS Forte implementation, the RNG value for the variable is based on the work of Han and Reitz 21 , who modified the constant to take the compressibility effect into account. According to Han and Reitz 21 , (2.22) where m =0.5, n =1.4 for an ideal gas, and (2.23) with (2.24) Using this approach, the value of varies in the range of -0.9 to 1.726 21 , and in ANSYS Forte is determined automatically, based on the flow conditions and specification of other model constants, η 0 and β . Han and Reitz 21 applied their version of the RNG k- ε model to engine simulations and observed improvements in the results compared to the standard k- ε model. For this reason, the RNG k- ε model is the default and recommended turbulence model in ANSYS Forte. Table 2.1: Constants in the standard and RNG k - ε models68



Standard



0.09



1.44



1.92



-1.0



1.0



RNG



0.0845



1.42



1.68



Equation 2.22 1.39



0.769 1.39



4.38



0.012



2.2.2. Large-Eddy Simulation Approach In the LES approach, the filtered flow quantities are solved by the transport equations on the CFD mesh. For example,



now refers to the filtered flow velocity. The undetermined stress term Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



in Equation 2.3 7



Basic Governing Equations is called the sub-grid-scale (SGS) stress tensor, which is essential in LES modeling. ANSYS Forte offers two types of LES models, the Smagorinsky model and the dynamic structure model. The Smagorinsky model 65 is based on a viscosity assumption and accounts for the dissipative nature of turbulent flows, i.e., kinetic energy is dissipated from the large-scale to the small-scale motions. It relates the SGS stress to the strain rate of the filtered flow field: (2.25) The SGS viscosity,



, is modeled as: (2.26)



in which



is the local CFD mesh size,



is the magnitude of the filtered strain rate tensor (see Equa-



tion 2.20 , but note that the overbar now means filtering), and grid kinetic energy, , is modeled as:



is a model constant, 0.17. The sub(2.27)



where .



is a model constant, 0.101. It is noted that both



The sub-grid flux term



and



depend on the CFD mesh size,



in the species transport Equation 2.1 is modeled as: (2.28)



in which is the SGS turbulent diffusivity. Similarly, the sub-grid flux term H in the energy equation (Equation 2.5 ) is modeled as: (2.29)



in which capacity



is the SGS thermal conductivity and is related to the SGS thermal diffusivity and heat by . The SGS mass and thermal diffusivity are related to the SGS viscosity by: (2.30) (2.31)



where and are the turbulent Schmidt and Prandtl numbers, respectively. The viscosity-based assumption makes Smagorinsky model effective in dissipating the kinetic energy from large-scale to small-scale motions. However, such dissipation can be excessive, and energy-containing flow structures may not be adequately resolved unless a very fine mesh is used. Although, its numerical stability in complex flows is an advantage. The dynamic structure model 52 is a non-viscosity and similarity-based model, relating the SGS stress to the “Leonard stress” derived from a larger filter: (2.32)



8



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chemical Kinetics Formulation where L is the Leonard stress tensor,



, and



. An explicit filtering operation (denoted



as ) needs to be performed to compute the Leonard stress tensor, and the explicit filter size is chosen as twice of the local CFD mesh size. Essentially, the first term on the right-hand-side of Equation 2.32 can be viewed as a similarity model relating the SGS stress to the Leonard stress, with a dynamic scaling factor of representing the ratio of sub-grid kinetic energy to sub-filter kinetic energy. The second term is an optional added viscosity-based model 71 . It is suggested to be activated for fuel injection and spray simulations for numerical stability purposes. The viscosity is modeled as (2.33)



where sprays.



and



. The symbols



and



are injection velocity and nozzle diameter in



To use the dynamic structure model, the sub-grid flux terms and can still adopt the viscosity-based models shown in Equation 2.28 and Equation 2.29 . The SGS viscosity is calculated as: (2.34) where closure



. A transport equation for the sub-grid kinetic energy,



, needs to be solved to provide



(2.35)



where



. The dissipation rate of the sub-grid kinetic energy has been modeled as



.



The dynamic structure model offers good prediction of the sub-grid stresses 59 , and has been validated in a number of benchmark flow problems 71 , 72 . It is recommended as the default option for the LES approach in ANSYS Forte.



2.3. Chemical Kinetics Formulation The chemical reactions that occur in combustion simulations can be described by chemical kinetic mechanisms that define the reaction pathways and the associated reaction rates leading to the change in species concentrations. In detailed chemical kinetic mechanisms, the reversible (or irreversible) reactions involving K chemical species can be represented in the general form 11 (2.36) The production rate of the k th species in the i th reaction can be written as (2.37) where q



i



is the rate of progress of reaction i.



ANSYS Forte uses the Chemkin-Pro chemistry solver. Details of the various types of reactions available in Chemkin-Pro can be found in the Chemkin-Pro Theory Manual , in these chapters: • Chapter 3: Gas-phase Chemical Rate Expressions:



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



9



Basic Governing Equations • Chapter 4: Surface Chemical Rate Expressions: The summation of



over all the reactions gives the chemical source term



in the species continuity



equation (Equation 2.1 ) as (2.38) Correspondingly, the chemical heat release term in the energy equation is given by (2.39) where Q



i



is the heat of reaction of reaction i at absolute zero, (2.40)



and



is the heat of formation of species



at absolute zero.



2.3.1. Optional Semi-empirical Soot Model The soot model used in ANSYS Forte is the two-step soot model consisting of competing formation and oxidation steps. The Hiroyasu soot formation 26 and Nagle and Strickland-Constable oxidation 44 models were used. The governing equations for the soot model are as follows: (2.41) (2.42) (2.43) (2.44) where p is pressure, n is a constant, M s is soot mass, with an additional subscript f for “formation” or o for “oxidation”. M pre is the mass of the soot precursor, K f is the soot formation rate, A sf is the pre-exponential factor for the global soot-formation reaction, E f is the activation energy for soot formation, MW c is the molecular weight of carbon, ρ s is soot density, D s is the assumed soot particle diameter, and R total is the Nagle and Strickland-Constable oxidation rate. More details about the soot model are given by Vishwanathan and Reitz 73 . The default values of model constants used in FORTÉ are: n=0.5, Asf = 40 cm3/mol-s, Ef =12,500 cal/mol, ρ s=2 g/cm3, and Ds=25 nm.



2.4. Turbulence-Kinetics Interaction Model ANSYS Forte includes a generalized turbulence-chemistry interaction model that is adapted from the work of Kong et al. [30 , 31 ], which was previously applied to simulations of Diesel engines. This model is included in ANSYS Forte as an option and is intended for simulating turbulence effects on combustion kinetics. This mixing time-scale model considers that the combustion chemistry should be partly controlled by the breakup of turbulent eddies due to the imperfect mixing of fuel and oxidizer in an actual engine



10



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Turbulence-Kinetics Interaction Model process. The model assumes every species moves towards its local equilibrium values with a time scale of τ eff and thus the effective (or actual) production rate of species k are expressed as (2.45) where the effective time scale is related to the chemical time scale τ mixing time scale τ mix as



chem



and the turbulent scalar (2.46)



The local chemical time scale is defined as the time it would take the mixture to reach equilibrium under the conditions within the computational cell. The turbulent scalar mixing time scale is obtained from the local turbulent kinetic energy and dissipation rate: (2.47) A relationship between the effective species production rate tion rate



and the kinetic-only species produc-



is derived from Equation 2.49 as (2.48)



The effective species production rates that are directly used in kinetics integration are then (2.49) In this way the gas composition and temperature remain consistent. The model constant C tki is one of the turbulence model constants and is provided as the Mixing Time Coefficient user input parameter in ANSYS Forte when the Turbulence Kinetics Interaction option is turned on. A value of 1.5 is found to be appropriate, based on engine simulation studies. To ensure that the effective species production rates do not go to zero when the turbulent mixing time scale is very large, the Turbulence Kinetics Interaction model is turned OFF when the τ mix value is above 1 millisecond.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



11



12



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 3: Boundary Conditions There are several types of boundary conditions that can be specified in ANSYS Forte. These include inflow, outflow, rigid (moving or stationary) walls, periodic, and symmetry boundaries. There are several types of rigid-wall conditions for momentum and energy equations. Wall boundary conditions for the momentum equation include free-slip velocity, no-slip velocity, and turbulent law-of-the-wall conditions. Wall boundary condition options for the energy equation include adiabatic walls and fixed-temperature walls. For in-cylinder engine simulations, the turbulent law-of-the-wall velocity condition and fixedtemperature walls are usually employed. The following sections describe in more detail each of the boundary-condition options available.



3.1. Wall Boundaries Heat transfer between the working fluid and the cylinder wall in an engine can significantly affect engine performance, efficiency and emissions. The heat flux through the chamber walls is mainly due to gasphase convection, fuel-film conduction and radiation. Under typical engine operating conditions, gasphase convective heat transfer is the dominant factor 23 . At rigid walls, boundary conditions must be specified or assumed for each of the gas-phase governing equations. Typically, the boundary layer of an engine in-cylinder flow is thin relative to the practical computational grid size. For this reason, velocity and temperature wall functions are often employed to solve for the near-wall shear stress and heat transfer.



3.1.1.Wall Conditions for the Gas-phase Species and Fluid Continuity Equations A flux balance at the surface is required for each gas-phase species in the system. The surface boundary condition for the gas-phase species equation is: (3.1) Where K g is the number of gas-phase species in the system, u is the convective velocity normal to the surface, V k is the diffusive velocity of the species and is the sum of any net external flux imposed at the boundary. ANSYS Forte does not currently allow for the possibility of mass deposition onto or surface chemistry at wall boundaries, so such net flux would be related to specified flow velocities only. For a no-slip velocity condition, then, the net flux (including diffusive and convective flux) of each species to the wall must be zero, and the sum of all species fluxes (or the net mass flux of gas) to the surface is also zero.



3.1.2. Wall Conditions for the Momentum Equation Fluid momentum boundary conditions on rigid walls are introduced by imposing a value of the velocity at the wall. On no-slip walls, the gas velocity is set equal to the wall velocity: (3.2) The wall stress is then determined implicitly through Equation 2.3 . On free-slip and turbulent law-ofthe-wall boundaries the normal gas velocity is set equal to the normal wall velocity, (3.3)



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



13



Boundary Conditions and the two tangential components of the wall stress are explicitly specified. For free-slip walls the tangential components of the wall stress are zero. For turbulent law-of-the-wall conditions the tangential components are determined by matching to a logarithmic profile: (3.4)



where



is the Reynolds number based on the gas velocity relative to the wall,



is evaluated a distance y from the wall, and u stress by



*



, which



is the shear speed, which is related to the wall shear (3.5)



In Equation 3.4 and Equation 3.5 it is assumed that y is small enough to be in the logarithmic region or the laminar sublayer region of the turbulent boundary layer. The Reynolds number R c defines the boundary between these two regions. The constants κ , c lw , R c , and B in Equation 3.4 are related to the k- ε model constants. Commonly used values are found in the KIVA-II documentation 4 .



3.1.3. Wall Conditions for the Energy Equation In ANSYS Forte, the wall heat transfer model of Han and Reitz 22 is used to calculate the gas-phase wall heat transfer. In the near wall region of wall-confined in-cylinder engine flows, the following assumptions hold: 1. Gradients normal to the wall are much greater than those parallel to the wall; 2. The flow velocity is directed parallel to the flat wall; 3. Pressure gradients are neglected; 4. Viscous dissipation and enthalpy diffusion effects on energy flux are neglected; 5. Radiation heat transfer is neglected. With these assumptions, a one-dimensional energy conservation equation in the near wall region can be written as (3.6) with (3.7) where J is heat flux, y is distance to the wall, ρ is the gas density, c



p



mixture, p is pressure, and is volumetric heat release, and k and k thermal conductivity, respectively.



is the specific heat of the gas t



are laminar and turbulent



Integration of Equation 3.6 over the distance y from the wall ( y =0) with consideration of the gasdensity variation and the increase of the turbulent Prandtl number in the boundary layer, the temperature profile equation (wall function) is formulated as (3.8)



14



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Wall Boundaries where



, y



+



and G



+



are defined as (3.9)



and (3.10)



Here, u* is the friction velocity defined in terms of the turbulence kinetic energy, k , as



, y is



the normal distance to the wall, J w is the heat flux through the wall, is the average volumetric chemical heat release, and ν is the laminar kinematic viscosity. The corresponding formulation for the wall heat flux is given as (3.11) where T and T



w



are the gas temperature and wall temperature, respectively. Following the implement-



ation of Han and Reitz 22 , we neglect the source term



, such that Equation 3.11 becomes (3.12)



Equation 3.12 states that due to the density variation, the wall heat flux is proportional to the logarithm of the ratio of the gas temperature to the wall temperature. Other energy-equation boundary conditions on rigid walls are introduced by directly specifying either the wall temperature or the wall heat flux . For adiabatic walls, we set J w equal to zero. For fixed-temperature walls that are also either free slip or no slip; the wall temperature is prescribed and J w is determined implicitly from Equation 2.5 . When turbulent law-of-the-wall velocity conditions are used, frictional heating serves as a source to the internal energy and it has the form (3.13) where f



w



is the heating per unit area of wall, and v is the tangential gas velocity relative to the wall.



3.1.4. Turbulence Model Wall Conditions In calculations of turbulent flow, boundary conditions are also needed for the turbulent kinetic energy k and its dissipation rate ε . These are taken to be



and (3.14) where k and ε are evaluated a distance y from the wall and



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



15



Boundary Conditions



(3.15)



3.2. Inflow and Outflow Boundaries The varieties of inflow and outflow boundaries in ANSYS Forte are: • Velocity inflow • Continuative outflow • Pressure inflow • Pressure outflow The usage and input requirements for each of these are described below.



3.2.1. Velocity Inflow Boundaries To specify a velocity inflow boundary requires a non-zero value for the normal velocity. The actual u, v, and w velocity components to be applied will be computed based on the boundary orientation. The information required at an inflow boundary includes the specific turbulent kinetic energy and turbulence length scale and species densities. The turbulence parameters are assumed to be constant; these values can either be specified or borrowed from the initial values of the neighboring region. The initial species densities outside the boundary are assumed to be at the reference ambient pressure, which is equal to the initial pressure in the neighboring region. Alternatively, an inlet mass flow rate can be specified instead of inlet velocity. The mass flow rate will be converted to inlet velocity based on the flow density and the inlet boundary area.



3.2.2. Outflow Boundaries For a continuative outflow boundary, the u, v, and w velocity components at the boundary vertex are set equal to those of the logical inside neighbor vertex. If reverse flow is allowed at such outflow boundaries, the species composition outside the boundary must also be specified. In ANSYS Forte, the initial species composition in the neighboring region is used for the reverse flow, and the initial pressure of the neighboring region is used as the reference pressure.



3.2.3. Pressure Inflow Boundaries For pressure inflow boundaries, a pressure vs. crank angle profile can be specified at the boundary. The ambient turbulence parameters and mixture composition are specified in a similar way as for the velocity inflow boundaries (Velocity Inflow Boundaries ).



3.2.4. Pressure Outflow Boundaries At a pressure outflow boundary, the pressure is imposed either at the boundary or at some distance beyond the boundary. Using a distance beyond the boundary enables the boundary to absorb acoustic waves, and reduce any tendency they may have to be reflected back into the system. The turbulence parameters and mixture composition for reverse flow are specified in a similar way as for continuative outflow boundaries (Outflow Boundaries ).



16



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Axis-of-Symmetry Boundaries



3.3. Periodic Boundaries Periodic boundaries are only used when the flow field is assumed to have an N -fold periodicity about the z-axis. This is the common case, for example, when the in-cylinder geometry of a diesel engine is represented by a sector mesh. For example, a 45-degree sector mesh would correspond to an eightfold ( N =8) periodicity. When periodic boundaries are used, the computational region is composed of points in a pie-shaped sector , where θ satisfies and . The periodic boundaries are those for which θ = 0 and θ = 2 π / N . The conditions imposed on these boundaries can be inferred from the assumed N -fold periodicity. For a scalar quantity q a constraint sets q(r, θ ,z) = q(r, θ + 2 π /N, z) , where r = g . For a vector v the constraint is that v (r, θ +2 π /N, z) = R ⋅ v (r, θ , z) , where R is the rotation matrix corresponding to the angle 2 π /N.



3.4. Axis-of-Symmetry Boundaries Also for periodic geometries (sectors), there will be an axis of symmetry, at which conditions must be specified. For sectors, an axis-of-symmetry boundary is handled as a degenerate cell face with zero area that lies on an axis.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



17



18



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 4: Initial Conditions Most Forte simulations are of transient engine conditions. For transient simulations, initial conditions for all variables and fluid properties are required. This section describes the initial conditions applied to the governing equations as well as the user input required to establish the initial conditions.



4.1. Initialization of the Fluid Properties Fluid properties are used to determine the initial values of the state variables and to initialize other terms in the governing equations. The fluid therefore must be initialized for any transient simulation. Determination of initial fluid properties requires: • Initial pressure for each computational region. • Initial temperature for each computational region. • Initial species composition for each computational region, in terms of mole or mass fractions. For in-cylinder regions, this may be the composition of pure air (e.g., 79% nitrogen, 21% oxygen), a mixture of exhaust gas (for exhaust-gas recirculation) and air, or a pre-mixed charge of pre-vaporized fuel, air, and/or exhaust gas. • Initial turbulence kinetic-energy density for each computational region. In an internal combustion engine, the initial turbulence kinetic-energy density (TKEI) is specified as a fraction of the total kinetic energy that is turbulent kinetic energy, where the total kinetic energy is determined based on the mean piston speed. In non-engine applications TKEI is the actual turbulent kinetic energy with a value typically about 10% of the mean flow kinetic energy 0.5* u 2. • Initial turbulence length scale (TLS) for each computation region. If TLS > 0.0, its value is used to determine an initial (uniform) value of ε , which will be proportional to the TLS value. If TLS = 0.0, the initial ε value will instead be proportional to the distance from the cell center to the nearest solid wall.



4.2. Swirl Profiles In addition to the fluid properties described in Initialization of the Fluid Properties, initial momentum resulting from swirl flow may be important for some types of simulations. Internal combustion engines in particular are designed to impart a significant amount of swirl motion into the incoming air. This aids in turbulent mixing and enhances combustion efficiency. From experimental observation, modelers have determined that a Bessel function profile for the swirl component of the flow more accurately represents the flow under IC engine conditions. Figure 4.1: Bessel-function swirl velocity profile illustrates the Bessel function velocity profile provided in ANSYS Forte and compares it with a simpler “wheel flow” model that is commonly used in other CFD representations. The comparison is shown for the same swirl number. The quantity α is a dimensionless constant that defines the initial azimuthal velocity profile and lies between 0.0 (the wheel flow limit) and 3.83 (zero velocity at the wall). In Figure 4.1: Bessel-function swirl velocity profile, the Bessel profile is illustrated for different values of α . The Bessel function profile is defined to give the same angular momentum as wheel flow with the same swirl number. Thus, higher values of α correspond to the higher initial slope of the Bessel curve. A value suggested by Wahiduzzaman and Ferguson 74 for typRelease 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



19



Initial Conditions ical engine applications is about 3.11. The Bessel function profile in ANSYS Forte is defined such that it gives the same angular momentum as a wheel-flow profile that has the same swirl number. Thus the initial slope of the α = 3.11 curve is necessarily higher than the corresponding slope for wheel flow. A second input quantity is the initial swirl ratio of the air rotation rate to crankshaft rotation rate (i.e., rpm). When viewed from the positive z direction (the top of the engine cylinder), the swirl is clockwise if SWIRL < 0 and counterclockwise if SWIRL > 0. Figure 4.1: Bessel-function swirl velocity profile



20



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 5: Discretization and Solution Methods 5.1. Discretization of the Governing Equations The governing equations are discretized with respect to the spatial coordinates of the system, on a computational grid based on a control volume approach. In addition, to provide time-accurate solutions, the equations are further discretized with respect to time following the operator-splitting method.



5.1.1. Temporal Differencing Method To integrate the equations in time, a temporal differencing of the equations is performed. For a sequence of discrete times t n



n



( n =0, 1, 2,...), the time interval for the next step is given as Δ t



n



= t



n+ 1



- t



is the time step, and the integer n is referred to as the time “cycle” number. Throughout this docu-



mentation, we use a superscript n on a variable to represent “the value at time t when we use the abbreviated “ Δ t ” without a superscript, we mean Δ t



n



n



”. Additionally,



. Using this notation, the



difference approximation to the derivative δ Q/ δ t is the first-order expression: (Q t.



n+ 1



-Q



n



)/ Δ



During time integration, ANSYS Forte employs three stages of solution for each time step. The time stepping employs the operator-splitting method to separate the chemistry and spray source terms and the flow transport. The flow transport solution is based on the Arbitrary-Lagrangian-Eulerian (ALE) method 4 . The hydrodynamic time step is adaptively controlled to maximize the solution efficiency, accuracy and stability. The three stages considered are: 1. Stage 1: This stage solves the chemistry and spray source terms in the species and energy transport equations. The calculation is based on the Lagrangian coordinate in which cells move with the fluid and spray droplets are followed through collision, oscillation, and break-up processes, along with mass and energy terms due to the chemistry, gas and spray interactions. 2. Stage 2: A coupled, implicit solution of the acoustic-mode terms, the spray momentum source term, and the terms due to diffusion of mass, momentum, and energy. This stage of calculation also includes the remaining source terms in the turbulence equations. The acoustic-mode terms include the pressure gradient in the momentum equation and velocity dilation terms in the mass and energy equations. 3. Stage 3: The third stage is a rezone stage, in which the flow field is frozen and then remapped onto a revised computational mesh after wall motion is accounted for. The convective transport is calculated by moving the Lagrangian mesh result from stage 2 to the revised mesh relative to the fluid.



5.1.2. Spatial Differencing Method ANSYS Forte uses the Arbitrary Lagrangian-Euler (ALE) method [27 , 51 ] for spatial differencing of the governing equations. For 3-D geometries, this considers a mesh consisting of arbitrary hexahedrons. A control-volume or integral-balance approach 27 is employed, which preserves the local conservation properties of the differential equations. The computational mesh consists of spatial regions, which are composed of cells, the corners of which are the vertices. The vertices are allowed to move in an arbitrarily prescribed manner, enabling piston and valve motion, for example. In addition, when differencing Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



21



Discretization and Solution Methods the momentum equations, a different control volume, called a momentum cell, is used. The momentum cells center about the vertices of regular computational cells. An important advantage of locating the velocities at cell vertices in the ALE method is that this requires no interpolation when determining vertex motion in the Lagrangian stage of the calculation. However, classical ALE method solutions are susceptible to parasitic modes in the velocity field. This problem has been alleviated in FORTÉ by the centering of velocities on cell faces, using the method described by Amsden 3 . Vertex velocities are retained, and momentum associated with the vertices is conserved, but normal velocity components on cell faces are used to compute cell volume changes in Stage 2 and the volume of fluid transported across a cell face in Stage 3. Accelerations of the cell-face velocities due to pressure gradients are calculated by constructing momentum control volumes centered about the cell faces. The volume of any momentum control volume are calculated once the volumes of the main computational cells are determined. In the finite-difference approximations of ANSYS Forte, primary velocities are located at the vertices, such that (5.1) Thermodynamic quantities, on the other hand, are primarily located at cell centers: (5.2) where Q = p, ρ , T, I, or ρ k , as well as k and ε . Quantities needed at points where they are not primarily located are obtained by averaging neighboring values. The discretized equations are formed by integrating the differential term in question over the volume of each cell (or momentum cell). Volume integrals of gradient or divergence terms are usually converted into surface area integrals using the divergence theorem. The volume integral of a time derivative is related to the derivative of the integral using Reynolds’ transport theorem 70 . Volume and surface area integrals are performed under the assumption that the integrands are uniform within cells or on cell faces. In this way, area integrals over surfaces of cells are represented as summations over cell faces (or sub-faces): (5.3) When differencing diffusion terms for a cell-centered quantity Q , it is necessary to evaluate (∇ Q ) α ⋅ A α . Explanation of this calculation requires reference to Figure 5.1: The six points used to define the gradient of cell-centered quantity Q on cell face α . , where the points x l and x r refer to the centers of the cells on either side of face α , and x t , x b , x f and x d are the centers of the four edges bounding face α . We first solve for coefficients a l r , a tb , and a fd such that (5.4) Note that since the mesh may be non-orthogonal, the vector x l - x r is not necessarily parallel to A α , and thus a tb and a fd may be nonzero. The finite-difference approximation to (∇ Q ) α ⋅ A α is obtained by applying a dot product with (∇ Q ) α on both sides of both sides of Equation 5.4 with (∇ Q ) α and neglecting terms that are second- and higher-order in the cell dimensions: (5.5) In Equation 5.5 Q t is a simple average of the values of Q in the four cells surrounding cell edge “ t, ” and Q b , Q f , and Q d are defined analogously, relative to edge “b”, “f”, and “d”.



22



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



SIMPLE Method Area integrals over momentum cell faces are ordinarily converted into area integrals over regular cell faces. This procedure is explained by considering a quantity Q that is uniform within the regular cells, and then consider the volume of overlap between a regular cell ( i, j, k ) and the momentum cell associated with one of its vertices. Three faces of this overlap volume (for example, a,b,c ) are faces of the momentum cell in question, with outward area vectors A ’ α , while the other three (e.g., d,e,f ) are surfaces of the regular cell ( i, j, k ), with outward area vectors ¼ A α . The divergence theorem shows that the integral



over the entire surface of this overlap volume is zero, such that: (5.6)



In this way, the integral



over the three momentum cell faces in question is represented by (5.7)



and the area vectors A ’ a do not need to be explicitly evaluated. A similar procedure is used to express the outward normal areas A ” a of faces of the cell-face control volumes in terms of the regular cell face areas A a . Figure 5.1: The six points used to define the gradient of cell-centered quantity Q on cell face α .



5.2. SIMPLE Method ANSYS Forte employs implicit methods in solving the algebraic finite-volume equations that result from the differencing methods described in Discretization of the Governing Equations . The advantages of implicit methods are: 1. In transient (unsteady) flow calculations, the time-step size is limited by the temporal accuracy desired, not by stability constraints that are typical of explicit methods. 2. When temporal accuracy is not desired, steady states can be achieved much more quickly by using large time steps.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



23



Discretization and Solution Methods Specifically, ANSYS Forte uses a modified version of the SIMPLE method 49 , which is a two-step iterative procedure used to solve for the flow field variables. Velocities at each time step in the numerical simulation need to be computed from time-advanced pressure gradients. This requires an iterative solution procedure, because the time-advanced pressures depend on accelerations and the velocities computed from the pressures. The SIMPLE method extrapolates the pressure, iteratively solves for velocities, then temperature, and finally the pressure. Within the ANSYS Forte 3-stage time-step advance procedure, Stage 2 involves selecting a predicted value of the pressure p b . The pressure field is then frozen during solution of other flow quantities where diffusion terms are handled fully implicitly. Then the resulting values of the diffusion terms are frozen and we solve for the corrected pressure field using equations that difference the pressure terms implicitly. In the second step of this iteration ANSYS Forte simultaneously solves the cell-face velocity equations, the volume change equations, and a linearized form of the equation of state. By algebraically eliminating the volumes and cell-face velocities from these equations in favor of the pressures, we are essentially solving a Poisson equation for the pressure in this step. Next, the predicted and corrected pressures are compared. If they agree to within a specified convergence tolerance, the equations have been solved, and we proceed to Stage 3. If the difference between the pressure fields exceeds the convergence tolerance, the corrected pressure field becomes the new predicted pressure field, and we return to the first step of the iteration and repeat the process. Each pass through the two steps is referred to as an outer iteration. The mass fractions are used in the calculation of the Stage 2 pressure, but the values of the Stage 2 pressures and velocities do not influence the solution of the mass fractions. Thus the species diffusion terms are solved to update the species mass fractions before beginning the outer iteration This results in a considerable computational time savings over other schemes, such as those that calculate implicit convection during the SIMPLE iteration and which include the mass fraction equations in the outer iteration. We often have many chemical species (10s to 100s) in our combustion applications, such that solving species and mass equations in the outer iteration would greatly increase computational times. For the Stage 2 calculation of k and ε , the flow field influences their values through the turbulence production terms after completion of the outer iteration. The finite-difference equations have been derived in such a way as to assure this one-way coupling. Mathematically, the values of k and ε influence the flow through the turbulent diffusivity and the turbulent pressure



. This coupling could be ac-



counted for by using Stage 2 values of k and ε to evaluate the turbulent diffusivity and but this would greatly increase computational times and is not necessary for stability. Furthermore, it is usually not necessary for accuracy when time steps are used that satisfy the stability constraints. For the ANSYS Forte implementation of the SIMPLE method, then, the only equations in the outer iteration are the momentum equation, internal energy equation, and the pressure equation.



5.2.1. Convective Flux Discretization Convection terms are solved using the quasi-second-order upwind (QSOU) method, as described by Amsden, et al. 4 . In ANSYS Forte, convection is calculated in Stage 3 in a sub-cycled explicit fashion that offers some significant advantages over implicit methods. Furthermore, because the equations for the species mass fractions, k and ε are weakly coupled to the flow-field solution, these equations are not included in the outer iteration.



24



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chemistry Solver Options



5.3. Chemistry Solver Options One of the major strengths of ANSYS Forte is its ability to handle realistic fuel-combustion mechanisms that include large numbers of species and reactions, within simulation times that are practical for design. This is accomplished through a number of unique and proprietary solver components within the Forte chemistry solver. The use of large detailed chemical kinetics to model fuel combustion gives rise to a large system of stiff nonlinear ordinary differential equations (ODEs). In traditional CFD solution approaches, this can result in very long CPU times. To overcome this CPU time barrier, ANSYS Forte employs several advanced solution strategies that drastically improve the chemistry solution efficiency by two to three orders of magnitude without compromise of accuracy. The solution techniques include a dynamic adaptive chemistry (DAC) method [33 , 34 ], a dynamic cell clustering (DCC) method 35 and an advanced proprietary sparse-matrix solution methodology developed at Reaction Design. The DAC and DCC methods are enabled by the use of an operator splitting approach to solving the species-conservation and energy equations, which is described in Operator Splitting Method and Parallel Implementation . The DAC and DCC methods are then described in Dynamic Adaptive Chemistry and Dynamic Cell Clustering , respectively.



5.3.1. Operator Splitting Method and Parallel Implementation ANSYS Forte employs an advanced operator-splitting method to solve the species-conservation and energy-conservation equations for time-accurate transient simulations. This method splits the transport equation into two sub-equations and solves the sub-equations with overlapping time-steps; the first step handles the kinetics and the second step handles the transport (e.g., diffusion or convection). The chemistry solver described in this section handles the first step. For this step, the calculation is performed on a cell-by-cell basis, allowing coupling of all of the species and energy together, since the transport terms are handled only in the second step. The transport of the species (i.e., the second step) is then handled using the methods described in SIMPLE Method . With the source term being handled by an advanced chemistry solver, the CFD step is well behaved and can proceed using reasonable time-step values that are limited only by the desired transport accuracy. To illustrate the method using equations, we first define the full species equation that we want to solve in each cell in the computational domain as: (5.8)



where is the net production rate of species k due to all chemical reactions that occur in that cell, V cell is the volume of the cell, ρ is the mass density of the fluid, Y k is the species mass fraction, C j represents the convective mass flow of the species into a cell from neighboring cell j , and D j is the diffusive mass flow of the species into a cell from the j th neighboring cell. Note that this equation is simplified just to illustrate the operator-splitting methodology. We can approximate the time derivative on the left-hand-side of Equation 5.8 , for a discrete time step, t , so that the equation becomes: (5.9)



where the superscript “ n ” represents evaluation of the terms at the new time and “ o ” represents evaluation at the old time. With the convective terms and the production terms all evaluated at the Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



25



Discretization and Solution Methods new time, the equation is fully implicit and fully coupled within the cell. Equation 5.9 can be rearranged to solve for the new value of the species mass fraction as: (5.10)



Now, since it is extremely difficult to solve this equation exactly, we choose to split the equation by letting: (5.11) (5.12)



In this way, we lose some implicitness in the chemistry production term, but we gain the fact that Equation 5.11 is now solved for each cell independently with no influence of the neighboring cells. This allows efficient solution of all of the species simultaneously instead of equation by equation. With operator splitting, the first step balances production and destruction of chemical kinetics against the timerate-of-change of the species in a cell, neglecting transport into or out of the cell. In this step the solution algorithm couples all species calculations together within each cell. In this calculation, there is no transport into or out of the cell, such that each cell can be considered independently of the others. As a result of this “half” step, the chemistry solver returns a new map of species concentrations over all cells, before transport occurs. The second step allows species to transport into/out of the cell, solving each species one at a time, but over all cells simultaneously. In this “half” step there is no further production or destruction of the species (i.e., there are no chemical source terms in this sub-step); only transport terms are considered. The key advantage to the operator splitting is that in the first “half-step”, the solution of the net species production rates due to chemical reaction is achieved on a cell-by-cell basis. On this basis, we solve for all of the species simultaneously with temperature, which means that the equations are fully coupled. The second advantage is that we can apply a state-of-the-art ordinary-differential equation (ODE) solver to the task of integrating the stiff simultaneous equations within each cell. The operator-splitting approach allows each algorithm to do what it does best. The transport algorithm is applied to transport in time; the stiff kinetic solver is applied to net production/destruction due to kinetics and resulting change of species in a “closed” system with time. The “sub-cycling” or integration of the species equation over the specified transport time-step is done very efficiently using ANSYS Forte's advanced, proprietary ODE solver that performs adaptive time-stepping with strict error control to assure a highly accurate solution in the most efficient time. The adaptive (sub) time steps are adjusted based on the largest rate of change for any species or temperature. The key advantage of this method is the adaptive selection of time step and sophistication in resolving both the trace and major species in the chemistry step. Furthermore the proprietary ODE solution strategy takes advantage of the sparsity of the species-to-species reaction matrix in a typical combustion mechanism to further gain major efficiency benefits. ANSYS Forte also takes advantage of the operator-splitting method through implementation of parallel solution for the cell-by cell calculations on a parallel platform. In such cases, the chemistry calculation is easily load balanced over any number of processors.



5.3.2. Dynamic Adaptive Chemistry The Dynamic Adaptive Chemistry method also takes advantage of the operator-splitting method described above in Operator Splitting Method and Parallel Implementation , by focusing on the local conditions of each cell. This method performs on-the-fly mechanism reduction for each local condition prior to solving the chemistry terms in the transient chemistry integration step of the operator-splitting method. 26



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chemistry Solver Options The reduction is dynamically applied cell-by-cell (or cluster-by-cluster if the Dynamic Cell Clustering method is used) at each time step. In this way the dynamic chemistry solution often uses locally valid, smaller mechanisms instead of the full mechanism, causing significant time savings for the overall transient integration. To ensure validity over a wide range of thermochemical conditions, comprehensive kinetic mechanisms for the combustion of realistic fuels typically include hundreds of species and thousands of elementary reactions. However, much smaller subsets of species and reactions often are adequate to capture the dominant reaction pathways for specific, local conditions over a short time span (typically taken to be the hydrodynamic time step in CFD calculations). To make use of this fact, the dynamic adaptive chemistry (DAC) method reduces the comprehensive detailed mechanism to locally valid smaller mechanisms. This operation is performed on the fly (i.e., during the dynamic simulation, at every time step). It is based on a skeletal mechanism-reduction method called Directed-Relation-Graph with Error Propagation (DRGEP) 38 method, which offers very efficient and very accurate reduction. Given a detailed kinetic mechanism and a specific thermochemical state X( T, p, yk ), where k is the species index ( k =1, , K ) and yk is the mass fraction of species k , the directed relation graphs in the DRGEP method are constructed such that one vertex (species) is connected to all others by directed edges. These edges are weighted by the immediate dependence of one species on another. This dependence is quantified by the normalized contribution of species B to A , defined by (5.13)



where i is the reaction index ( i =1, …, I ), υ Ai is the stoichiometric coefficient of species A in the i th reaction, and ω i is the progress variable of reaction i . r AB is a measure of the error introduced to the production rate of A due to elimination of all the reactions that contain B . Certain species deemed of primary importance are selected as initial species in the reduced mechanism. Then starting from each of these pre-selected initial species, a breadth-first search (BFS) is performed to identify the species on which the initial species depends to form a subsidiary set. Consequently, the union of the subsidiary sets of all the initial species forms the active species set of the reduced mechanism. Thus the mechanism reduction is equivalent to identifying vertices to which there exist “strong” paths that connect them to a vertex in the initial set. The strength of the connection between the species being visited and the initial species diminishes as we proceed along a path. This diminution can be used to control the search depth. To quantify the decreasing dependence, an “ R -value” is defined at each vertex V with reference to the initial vertex V 0: (5.14) where Ω is the set of all possible paths leading from V 0 to V , and Π r ij is the chain product of the weights of the edges along the given path. Based on this definition, vertex V will be marked as “reachable” if R V 0(V) is larger than a user-defined threshold value ε R . Thus, all the reachable vertices starting from V 0 constitute the subsidiary set of V 0 , and the union of such sets gives the species that are active in the reduced mechanism. This method has been tested on both diesel and gasoline surrogate fuels in Liang et al. [33 , 34 ]. The studies proved that {fuel, CO and HO2} is an effective choice for the search-initiating species set, and suggested that ε R be set in the range of 10 -5 to 10 -4 . Both the search-initiating species set and the search threshold are provided as user inputs in ANSYS Forte (Initial Species and Search Tolerance, respectively).



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



27



Discretization and Solution Methods The DRGEP method extracts a group of active species for the existing local thermochemical conditions. Consequently, a reaction is included in the reduced mechanism only if all the reactants and products are active species (third bodies are not counted as participants). Species not in the active set are treated as inactive, with their mass fractions kept fixed. Given a system that involves m active and n inactive species , where the superscripts “ a ” and “ i ” denote active and inactive species, respectively, the formulation can be expressed as



(5.15)



In Equation 5.15 Equation 5-15 , ordinary differential equations are formulated with respect to only active species. When the rate functions are evaluated, however, all the species are considered, thus eliminating the need to include the third-body species in the reduced mechanisms.



5.3.3. Dynamic Cell Clustering In the operator-splitting method for transient simulation, the chemistry-solution portion of the species equations are solved on a cell-by-cell basis. As a result, the equations solved are independent of the mass and volume of a specific cell. In this way, cells that have the same temperature, pressure, and initial species mass fractions will yield the same result. To take advantage of this fact, ANSYS Forte uses a Dynamic Cell Clustering (DCC) method to group computational cells of high similarity into clusters using an efficient data-clustering method and that requires solution of the kinetic equations only once for each cluster. The optimal number of clusters is dynamically determined for each CFD time step. The clustering algorithm is solely based on the cells' thermochemical states and is independent of their locations in the CFD mesh. In addition, the clustering algorithm is highly automated and requires minimal user-provided inputs. For typical hydrocarbon combustion computation, cell temperature and equivalence ratio are employed as the clustering indices. The only user-provided control parameters are the maximum dispersions of temperature and equivalence ratio allowable in each cluster. Smaller values will result in a larger number of clusters. The DCC method includes three major steps: (1) grouping cells into clusters using an evolutionary dataclustering algorithm; (2) solving chemical kinetic equations based on cluster averaged state variables; and (3) mapping the cluster averaged solution back to the individual cells while preserving the initial temperature and species stratification. In ANSYS Forte, if only one clustering feature is selected, temperature will be the feature. If two features are selected, both temperature and equivalence ratio will be used by default.



28



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 6: Spray Models ANSYS Forte includes advanced models to simulate multi-component fuel-spray dynamics and spray interactions with flowing multi-component gases. The modeled sub-processes include: nozzle flow, spray atomization, droplet breakup, droplet collision and coalescence, droplet vaporization, and wall impingement. Different sub-models are applied to represent the spray atomization and droplet breakup processes of solid-cone sprays and hollow-cone sprays, respectively. For solid-cone sprays, the initial spray conditions at a nozzle exit are determined through either a discharge coefficient or a nozzle-flow model, and the Kelvin-Helmholtz / Rayleigh-Taylor (KH/RT) model 9 is used for droplet breakup. For hollow-cone sprays, both the spray initialization and breakup processes are represented by the Linearized Instability Sheet Atomization (LISA) model 60 . Other spray sub-models are shared by both solid-cone and hollow-cone spray simulations. The spray-atomization, droplet-breakup and droplet-collision models in ANSYS Forte employ several advanced techniques that minimize the dependency of results on mesh size and time-step size. For solid-cone sprays, an unsteady gas-jet breakup model is used by default. For both solid-cone and hollowcone models, additional options include a radius-of-influence (ROI) collision model and a collision mesh method. Reducing the mesh-size and time-step dependency enables accurate solutions with less computational time compared to conventional spray models, by allowing a coarser mesh near the nozzle and larger time steps during spray injection. This chapter is organized to provide information on the main sub-models used by the various spray options. Solid-Cone Spray Models describes the solid-cone spray model, which is most often used for diesel-spray injection, while Hollow-Cone Spray Models describes the hollow-cone spray model. Submodels specific to these types of sprays are included within these sections. Droplet Collision and Coalescence Model provides information regarding additional options that are available for both types of sprays. Multi-Component Droplet Vaporization Model describes the multi-component vaporization model and Wall Impingement Model the wall-impingement model.



6.1. Solid-Cone Spray Models To initialize the spray for solid-cone injections, there are two options: use of an empirical nozzle discharge coefficient and use of a nozzle-flow model to determine nozzle discharge characteristics. These options are described in Nozzle Flow Model . Following this discussion are details on the Kelvin-Helmholtz / Rayleigh-Taylor droplet-breakup models in Kelvin-Helmholtz / Rayleigh-Taylor Breakup Model , and the unsteady gas-jet model in Unsteady Gas-Jet Model .



6.1.1. Nozzle Flow Model The nozzle-flow model describes the instantaneous flow conditions inside the nozzle. The purpose of the nozzle-flow model is to provide initial conditions for the spray model. The input parameters to the nozzle-flow model are: 1. mass flow rate 2. ambient gas pressure



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



29



Spray Models 3. physical properties of the liquid fuel 4. geometrical hole diameter 5. L/D ratio (see Figure 6.1: Flow through nozzle passage ) 6. R/D ratio, where R is the radius of curvature of the injector entrance region (see Figure 6.1: Flow through nozzle passage ) Using these inputs, the nozzle-flow model determines the instantaneous discharge coefficient (Cd), spray angle, effective injection velocity and effective flow exit area. The flow exit area is then used to determine the initial injected liquid droplet size. Figure 6.1: Flow through nozzle passage



6.1.1.1. Discharge coefficient The volumetric mean flow velocity of the liquid fuel at the inlet of the nozzle passage, calculated as



, can be (6.1)



where



is fuel mass flow rate,



is liquid fuel density, A is nozzle cross-sectional area, and D is nozzle



diameter. The mean mass flow through the nozzle exit is always lower than that predicted by the Bernoulli equation, because of flow losses. The losses are caused by the acceleration and formation of the velocity profile at the inlet, the expansion after the vena contracta, and the wall friction. To quantify this difference, the discharge coefficient, , is defined as: (6.2)



where p 1 and p 2 are pressures at location 1 and 2, respectively, in Figure 6.1: Flow through nozzle passage .



30



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Solid-Cone Spray Models By using tabulated inlet loss coefficients ( K inlet) and the Blasius or the laminar equation for wall friction, the discharge coefficient is given by: (6.3) with f defined as (6.4) and Re



D



is the Reynolds number based on nozzle diameter.



A first estimation of the inlet pressure, p



1



, for a turbulent flow yields: (6.5)



Next ANSYS Forte checks whether the flow under this condition is already cavitating. Assuming a flat velocity profile, and using Nurick’s expression for the size of the contraction (Figure 6.1: Flow through nozzle passage ), continuity gives the velocity at the smallest flow area: (6.6) (6.7) In the case of cavitation the potential flow theory allows the application of the Bernoulli equation from point 1 to c without any losses: (6.8) If p vena is lower than vapor pressure p vapor, it is assumed that the flow must be fully cavitating and a new inlet pressure and discharge coefficient are calculated by: (6.9) (6.10)



6.1.1.2. Effective Injection Velocity and Effective Flow Exit Area In the cases of non-cavitating turbulent flows, the exit velocity of the droplets is set equal to the mean velocity, U mean, and the diameter of the jet at the nozzle exit is equal to the nozzle diameter. When the nozzle is cavitating, those values are calculated as: (6.11) (6.12) (6.13)



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



31



Spray Models



6.1.1.3. Spray Angle The spray angle in the nozzle flow model is estimated using an aerodynamic model. This approach is based on Taylor’s analysis of high-speed liquid breakup due to the unstable growth of surface waves and the resulting mass shedding. In this approach, the spray angle, , is expressed as: (6.14) Where A is 3+0.28( L/D) and L/D is the length-to-diameter ratio of the nozzle. The function f(T) is approximated by (6.15)



6.1.1.4. Empirical Nozzle Discharge Coefficient and Spray Angle Instead of applying the nozzle-flow model, the user can also specify empirical constant values for discharge coefficient and spray angle. The empirical discharge coefficient is then used to calculate the effective injection velocity and initial drop size at the nozzle exit.



6.1.2. Kelvin-Helmholtz / Rayleigh-Taylor Breakup Model The spray atomization and droplet breakup of solid-cone sprays are modeled by the Kelvin-Helmholtz / Rayleigh-Taylor (KH /RT) hybrid breakup model [9 , 66 ]. Figure 6.2: KH/RT breakup model for solidcone sprays delineates the way the KH/RT breakup models are applied. The KH breakup model is applied within a specified Breakup Length, L , from the nozzle exit (region A). It strips small droplets off the jet (represented by parent parcels, or “blobs”), but the jet still maintains itself as a dense liquid core. Beyond the Breakup Length (region B), the RT model is used in conjunction with the KH model to predict secondary break-up [75 , 9 ].



32



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Solid-Cone Spray Models Figure 6.2: KH/RT breakup model for solid-cone sprays



6.1.2.1. Kelvin-Helmholtz Breakup The Kelvin-Helmholtz (KH) model is based on linear stability analysis of a liquid jet 32 and is used to model the jet’s primary breakup 56 region. From the linear stability analysis any perturbation applied to a liquid-gas interface can be expanded as a Fourier series, with its fastest growing mode contributing to the eventual breakup and generation of new droplets. The growth rate and wave length of the fastest growing mode was found numerically 56 to be: (6.16)



(6.17)



Here is the wave length of the fastest growing wave, is its growth rate, r is surface tension, the dimensionless gas Weber number We g is defined as:



p



is the jet radius,



(6.18) with U



rel



being the magnitude of the liquid-gas relative velocity: (6.19)



Here



is the CFD gas-phase velocity. Similarly, the liquid Weber number is defined as: Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



33



Spray Models



(6.20) and the dimensionless Ohnesorge number is: (6.21) The Reynolds number is defined as: (6.22) where



is the dynamic viscosity of the liquid. The last dimensionless number is the Taylor number: (6.23)



The radius of new droplets (denoted as r c ), which are created during the primary breakup process, is related to the wave length as follows: (6.24) where



is the Size Constant of the KH breakup model.



ANSYS Forte employs the concept of “ blob injection” 57 , in which the liquid-jet injection is represented with parcels of blobs that have their initial characteristic sizes equal to the effective nozzle diameter. Using this concept, primary breakup is modeled as a process where new droplets (with radius r c ) are formed from parent droplets (with radius r p ). The radius change of the parent droplet due to mass lost to its “child” droplets is described by the rate equation 57 : (6.25) where r



c



is less than or equal to r



p



. Here, the breakup time scale



is calculated as: (6.26)



In Equation 6.26 , is the Time Constant of KH breakup. This is a user-controlled input for which a general-purpose default value is provided in ANSYS Forte. In general it is impractical to track each droplet formed from the parent droplet. Using the “droplets in parcel” concept, new droplets cannot leave their parent parcel until their mass reaches or exceeds 3% of the average injected parcel mass. The new droplets are then grouped into a “shed parcel” (or child parcel) that is separated from its parent. This parcel-shedding process helps to predict a realistic droplet size distribution in the primary breakup region. Also, as the number of parcels is increased, the local statistical resolution is improved. The value 3% is provided as a default model constant in ANSYS Forte, which may be adjusted to change the size-distribution resolution. In summary, the KH breakup model involves two steps: the first step gradually reduces the parent droplet size by the rate equation (Equation 6.25 ); the second step creates a new parcel from the parent. During the second breakup step, two different options can be used to determine the droplet sizes and the number of droplets in the parent-child parcels.



34



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Solid-Cone Spray Models • In the first option, the number of droplets in the parent parcel will remain the same before and after the breakup. The size of the droplets in the parent parcel will simply be adjusted based on the mass taken away by the child parcel. This is the default option used in ANSYS Forte. • In the second option, it is argued that the second breakup step involves a numerical regrouping of droplets and is not an actual breakup event. Therefore, in addition to mass conservation of parent-child parcels, the original Sauter Mean Diameter (SMD) is also conserved. Imposing this SMD conservation constraint typically leads to a slower breakup rate.



6.1.2.2. Rayleigh-Taylor Breakup Beyond the Breakup Length from the nozzle exit, the Rayleigh-Taylor (RT) model is used together with the KH model to predict secondary breakup of spray droplets [66 , 9 ]. The “ breakup length ” is predicted by Levich's theory 32 as: (6.27) in which C b is called the RT Distance constant in ANSYS Forte. The RT model is applied beyond the RT Distance constant times the Breakup Length from the nozzle exit. The RT model considers the instability generated when two fluids with different densities accelerate in a direction that is perpendicular to their interface. The frequency and wavelength of the fastest-growing wave are given by Bellman 10 as: (6.28)



(6.29) In the case of a high-speed droplet moving in air, a is the deceleration due to drag. As in the KH model, this fastest growing wave is assumed to lead to the breakup of droplets. The radius of the newlyformed droplet and the time of breakup are predicted as: (6.30) (6.31) in which B RT , C RT are two constants. B RT is the Size Constant of RT Breakup. C RT is the Time Constant of RT Breakup. Unlike the KH model, the standard implementation of the RT model does not strip small child droplets off from their parent; instead, a catastrophic breakup of the parent droplet into small droplets is assumed, and the liquid column is eventually disintegrated and dispersed into the ambient gas. To reduce the time step dependency of the RT model, the RT breakup process is modeled by a rate equation similar to that of the KH model, Equation 6.25 : (6.32) in which r c is predicted by Equation 6.30 . The parent droplet breaks up continuously at each time step, and the time-step dependency is avoided. Assuming r c and does not change with time, Equation 6.32 can be solved analytically as: Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



35



Spray Models (6.33) Here r p,0 is the radius of parent droplet at the start of break-up. Thus, required to break the parent droplet.



is the characteristic time



In practice, neither rate equation (Equation 6.25 ) nor (Equation 6.32 ) can be solved analytically, since both r c and and (or ) change with time, r p and U rel can change significantly. Therefore, a sub-cycle approach is used to solve these rate equations to eliminate time-step dependency.



6.1.3. Unsteady Gas-Jet Model Any mesh-dependency of the KH-RT breakup model is mainly due to the calculation of the liquid-gas relative velocity U rel in Equation 6.19 , in which is taken to be the CFD cell gas velocity. In ANSYS Forte, the unsteady gas jet model [2 , 56 , 79 ] is used to remove this mesh-size dependency for the liquid droplet-ambient gas coupling. The gas-jet model is based on unsteady gas-jet theory 1 , in which the axial droplet-gas relative velocity is modeled without use of discretization on the CFD mesh (Figure 6.3: Unsteady gas-jet model ). Figure 6.3: Unsteady gas-jet model



In the gas-jet model, the jet tip develops with respect to time according to: (6.34) where x is jet tip penetration, y is the local spray-axial location of the particle, K is an entrainment constant, U inf,eff is the " effective injection velocity", d eq is the equivalent diameter which is related to nozzle diameter d noz and liquid-gas density ratio by:



36



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Solid-Cone Spray Models



(6.35) The downstream spray-axial location x



0



marks the start of jet-velocity decay and is calculated as: (6.36)



The effective injection velocity is determined as an integral of the response to changes of the injection speed from the start of injection t 0 to the current time t : (6.37) The response function R takes the form: (6.38) where



is a response time scale, which is related to the local flow time scale



by a Stokes number: (6.39)



This assumption is justified by the fact that a fluid particle responds to a change in the surrounding gas velocity exponentially 14 . The time scale of the response is determined by how long the local fluid particle has resided in the flow and by the local spray-axial location of the particle (denoted as y ). The local gas-jet speed along the spray-axis direction is correspondingly calculated as: (6.40) Assuming axis symmetry, the gas-jet velocity at any radial location r within the jet cross-section can be calculated as 2 : (6.41)



Using the gas-jet velocity in Equation 6.41 , the droplet-gas relative motion is modeled as: where



is the droplet velocity vector, C D is the drag coefficient,



is the local gas-phase turbulent



fluctuating velocity vector, and is acceleration due to gravity. is the gas-phase mean flow velocity with the axial component corrected by Equation 6.41 . The radial and azimuthal components are retained, so that other flow velocity components are not influenced. The gas-jet-corrected relative velocity U



rel



: (6.42)



only applies beyond a distance x 0 from the nozzle. In the near-nozzle region ( y < x 0 ), where the droplet velocity is very close to the injection velocity, a parabolic profile is used, which merges with Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



37



Spray Models the profile for y > x 0 (see Figure 6.4: Effective injection velocity in the unsteady gas-jet model. ), such that the jet velocity is:



(6.43)



in which is a constant with a fixed value of 0.6. Applying the parabolic entrained gas-jet velocity profile along the spray axis is very effective in reducing the mesh-size dependency of the breakup models. In ANSYS Forte, the Gas Entrainment Constant K in Equation 6.34 is a user-defined input. A larger value of K promotes gas entrainment and thus reduces spray penetration length. In addition, the Effective Distance Factor is a non-dimensional distance factor used to control the effective range of the gas jet model. The actual effective distance is this factor times the breakup length and is measured from nozzle exit to the spray droplets. Figure 6.4: Effective injection velocity in the unsteady gas-jet model.



6.2. Hollow-Cone Spray Models Hollow-cone sprays are typically created by pressure-swirl injectors. The injector contains internal swirl vanes that produce rotational motion in the liquid. The liquid forms a film along the inside walls of the injector with an air core in the center. The liquid film becomes a free sheet when exiting the injector, and the tangential velocity of the liquid becomes radial. As the radial distance from the centerline to the sheet increases, the rotation decreases due to conservation of angular momentum. Additionally, because of the conservation of mass, the sheet thins as it progresses, and subsequently disintegrates into droplets, forming a hollow cone spray. Once droplets are formed, their behavior is governed by secondary breakup, drag, collision, coalescence and vaporization.



38



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Hollow-Cone Spray Models The modeled processes in ANSYS Forte are sketched in Figure 6.5: Modeled processes in a hollow-cone spray. . Both inwardly-opening (with an air core but no pintle) and outwardly opening (with a pintle but no air core) are considered. These two sub-models governing the behavior of these sprays are described in the following sub-sections. Figure 6.5: Modeled processes in a hollow-cone spray.



6.2.1. Linearized Instability Sheet Atomization (LISA) Model The transition from internal injector flow to a fully developed hollow-cone spray is modeled using the linearized instability sheet atomization (LISA) model 60 . In the LISA model, the process is divided into three stages: 1. Film formation 2. Sheet breakup 3. Atomization The modeling of each of these stages is described in the following sub-sections.



6.2.1.1. Film Formation The centrifugal motion of the liquid within the injector creates a liquid film surrounding an air core. The thickness of the film, t f , is related to the mass flow rate, , by (6.44) where



is the liquid density, u the axial component of velocity at the injector exit, d



0



the injector



hole diameter. u is related to the total velocity U by (6.45) where the cone half-angle is a user-specified input parameter for a specific injector. The total velocity U is related to the pressure drop across the injector exit by



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



39



Spray Models



(6.46) Based on similarity considerations between swirl ports and nozzles, the discharge coefficient k v is set to a fixed value, 0.7, but to guarantee that the size of the air core is non-negative, the following expression is used for k v , (6.47)



6.2.1.2. Sheet Breakup The sheet-breakup model assumes that a two-dimensional, viscous, incompressible liquid sheet of thickness 2h moves with velocity U through a quiescent, invisid, incompressible gas medium. A spectrum of infinitesimal disturbances of the form (6.48) is imposed on the initially steady motion and produces fluctuating velocities and pressures for both the liquid and the gas, where is the initial wave amplitude, is the wave number, and is the complex growth rate of the surface disturbances. The most unstable disturbance has the largest value of , denoted by , and is assumed to be responsible for sheet breakup. Thus, it is desired to obtain a dispersion relation from which the most unstable disturbance can be deduced. Two solutions that satisfy the liquid governing equations subject to the boundary conditions at the upper and lower interfaces 62 . For the first solution, called the sinuous mode, the waves at both surfaces are exactly in phase. For the second solution, the so-called varicose mode, the waves are radians out of phase. Typically the sinuous mode is sufficient for representing engine conditions, and thus a simplified form of the dispersion relation for pressure-swirl atomizers is used: (6.49) Here



is the liquid kinematic viscosity, Q is the gas/liquid density ratio



, and



is the surface



tension. Once the unstable waves on the sheet surface grow to a critical amplitude, ligaments are formed due to the sheet breakup. The breakup time for this process is formulated based on an analogy with the breakup length of cylindrical liquid jets, i.e., (6.50) where is the critical amplitude at breakup, Ω is found by numerically maximizing Equation 6.49 as a function of k . The corresponding breakup length L is estimated by (6.51)



where the quantity



is given a constant value 12. Based on a mass balance, the resulting



ligament diameter at the point of breakup is derived as (6.52)



40



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Hollow-Cone Spray Models where K s is the wave number corresponding to the maximum growth rate, Ω. Based on the assumption that the sheet is in the form of a cone with its vertex at a point behind the injector orifice, the sheet half-thickness h at the breakup position L , is approximately (6.53) where (6.54) Assuming that breakup occurs when the amplitude of the most unstable wave is equal to the radius of the ligament d L , then a mass balance gives the drop size d D (6.55) where the most unstable wavelength K



L



is given by (6.56)



based on an analogy to Weber’s result for growing waves on cylindrical, viscous liquid columns.



6.2.1.3. Atomization As a consequence of the sheet breakup process described above, fuel droplets are introduced into the computational domain with certain initial conditions. Subsequently, the droplets are subject to secondary breakup (modeled using the Taylor-Analog-Breakup model), collision and coalescence, aerodynamic drag, vaporization, and wall impingement. The models used to describe the atomization physics are described in the following sections.



6.2.2. Taylor-Analog-Breakup Model The Taylor-Analog-Breakup (TAB) model 48 is used to model the secondary breakup of the droplets in hollow-cone sprays. The TAB model exploits an analogy between a distorted droplet and an oscillating spring-mass-system. The external forces acting on the mass, the restoring force of the spring, and the damping force are analogous to the gas aerodynamic force, the liquid surface tension force, and the liquid viscosity force, respectively. The force balance on the droplet gives (6.57) where t is time, y is the normalized (by the drop radius) drop distortion parameter, σ is the surface tension coefficient, and the subscripts g and l denote the gas and liquid phase, respectively. It is assumed that breakup occurs if and only if y > 1 48 . When this condition is satisfied, the droplet breaks up into smaller children droplets with sizes determined by an energy balance taken before (subscript 1) and after (subscript 2) the breakup as (6.58)



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



41



Spray Models In ANSYS Forte’s TAB model used in the present study, the droplet size after breakup, r to follow a Rosin-Rammler distribution 4 .



2



, is assumed



6.3. Droplet Collision and Coalescence Model Droplet collision and coalescence are important physical phenomena in dense sprays. The model of O’Rourke 48 identified three primary collision outcomes: 1. Permanent coalescence 2. “Grazing” (drops separate after collision, possibly forming satellite droplets) 3. Shattering Munnannur and Reitz 41 suggested several additional collision outcomes: 1. Bounce 2. Stretching 3. Reflexive separation Which outcome occurs depends on forces acting on the coalesced pair of droplets. At low Weber numbers, surface forces dominate over liquid inertia forces, and the droplets coalesce permanently. At higher Weber numbers, the liquid inertia forces become more important and the grazing collisions occur. With further increase of Weber number, the dominant liquid inertia forces cause shattering of the colliding droplets, forming a group of small droplets. All of the outcomes listed above are in ANSYS Forte to provide accurate predictions of the results of a successful collision. The probable number of collisions of the two colliding droplet parcels can be written as: (6.59) in which N radii, and V



and N 2 are the numbers of droplets within the two parcels, r 1 and r col is the control volume in which the parcels are allowed to collide.



1



2



are the droplets



6.3.1. Radius of Influence (ROI) Collision Model In the collision model of O’Rourke, 48 spray particles are allowed to collide only if they reside in the same computational cell. In a cylindrical mesh, the cell size around the spray axis can be very small, preventing collisions from occurring. The Radius-of-Influence (ROI) collision model 41 is used in ANSYS Forte to remove both mesh-size dependency and time-step dependency for the droplet collision process. In the ROI method, one particle is allowed to collide with another only if this particle resides within in the radius of influence of the other. Collision partners of a certain parcel are searched for within a “sphere of influence” centered at the parcel's location. All parcels within this sphere are possible collision partners for collide with the centered parcel. This approach removes the dependency on the CFD mesh size (as seen in O'Rourke's model 48 ), and the dependency on collision mesh size if a separate collision mesh is used 63 . The issue of time-step dependency was also investigated in Munnannur 2007 42 in light of the observation that the time step must be properly chosen such that the average distance traveled by a particle 42



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Multi-Component Droplet Vaporization Model within the time step is a fraction of the average cell length. Using the ROI approach, the “ mean collision time” (MCT) for each cell is estimated as: (6.60) in which u d,max, cell is the maximum droplet velocity magnitude in the cell, R inf is the Radius of Influence, N cell is number of parcels in the cell, and C col is a constant. The global MCT is then the minimum of the cell MCTs. In ANSYS Forte, the MCT is calculated for each parcel based on the ROI concept: (6.61) in which u d is the average droplet velocity magnitude, and N ROI is number of parcels in the sphere of influence. The global MCT is taken as the minimum of all parcels' MCTs. The time step is dynamically controlled, such that if Δ t col is smaller than the computational time step, Δ t col is chosen as the computational time step; if Δ t col is greater than the time step, collisions occur only when the accumulated time of Δ t col is reached. To avoid imposing a small collision time step Δ t col on the entire engine simulation, the collision calculations are sub-cycled within the CFD time step Δ t if Δ t col < Δ t . After each sub-cycle, all the droplet properties except their spatial locations are updated to reflect the collision outcomes, and these updated properties are used to predict the collision process in the next sub-cycle. In the ROI model, the radius-of-influence is an input parameter in ANSYS Forte that assumes the suggested value of Munnannur 2007 42 by default.



6.3.2. Collision Mesh Model In addition to the ROI model, ANSYS Forte includes an alternative approach to mitigating the mesh dependency of the collision model 67 . This is the collision mesh model, which employs a pseudo cylindrical collision mesh separate from the CFD computational mesh. The collision mesh is based on a cylindrical coordinate system that is coaxial with the spray. To further reduce non-axisymmetry of the predicted spray structure, the collision mesh is rotated at a random angle around the spray axis at every time step. In this method, only parcels in the same collision mesh cell, instead of a CFD cell, are allowed to collide. Use of the collision mesh model requires specification of the diameter, height, and number of grids along the azimuthal direction of the collision mesh.



6.4. Multi-Component Droplet Vaporization Model ANSYS Forte uses a discrete multi-component (DMC) fuel-vaporization model 54 model to represent the vaporization of spray droplets. The DMC vaporization model tracks the individual components (molecules) of an actual surrogate fuel during the evaporation process and allows coupling with the reaction kinetics of the individual fuel components. In the DMC model, an explicit form of the equation that determines the heat flux from the surrounding gas mixture to the droplet-gas interface is obtained from an approximate solution of the quasi-steady energy equation. The model is formulated to track each component of the fuel regardless of the direction of the process, i.e., whether evaporation from the droplet surface or condensation into the droplet is occurring.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



43



Spray Models The DMC vaporization model considers a spherical liquid droplet that consists of a finite number of components vaporizing without chemical reactions in a gaseous environment. Radiation and secondorder effects such as the Soret and Dufour effects are assumed to be negligible.



6.4.1. Liquid Phase Balance Equation With no absorption of ambient gas into the spherical liquid droplet, a general form of the governing equation for the change in the liquid fuel distribution is (6.62) where y i is the mass fraction of component i in the liquid droplet, ρ i is the mass density of the liquid fuel, R is the droplet radius, and is the vaporization rate per unit area of species i. The change of liquid droplet energy is obtained from the conservation equation of energy for the two phase system consisting of the droplet and the surrounding gas mixture as (6.63)



where c v,l is the specific heat of the liquid fuel, q i is the heat transfer rate from the droplet surface to the interior per unit area, and T d and T s are the average droplet temperature and surface temperature, respectively.



6.4.2. Governing Equations for Gas Phase The transport portion of the conservation equation of species in the gas phase is (6.64) where v and ρ are the velocity and density of the gas mixture, respectively, yi and Di are the mass fraction and diffusion coefficient of species i, and s g,i is the source term from vaporization. Summation of Equation 6.64 gives the species conservation equation for the two-species system (fuel and air) as (6.65) where y F is the total mass fraction of fuel species, D is the average diffusion coefficient of the fuel species, and Sg is the total vaporization source term. The energy conservation equation for the gas phase is (6.66) where T is the temperature, k is the thermal conductivity,



is the mixture specific heat, C



PA



is the



specific heat of air, and is the average value of the product of specific heat and the diffusion coefficient of the fuel species. The last term in Equation 6.66 represents energy transport due to interdiffusion of species.



44



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Multi-Component Droplet Vaporization Model



6.4.3. Vapor-Liquid Equilibrium The equilibrium at the interface between the liquid droplet and the surrounding gas is based on the assumption that the chemical potentialfor the liquid phase, l , and the vapor phase, v , are equal for each species, i . Assuming an ideal solution, the surface mass fraction of fuel vapor can be determined using Raoult’s law. For a mixture of discrete components, Raoult’s law is (6.67) where p i is the partial pressure of species i in the vapor phase at the droplet surface, Psat,i(T) is the vapor pressure of species i at temperature T, x is the mole fraction and the subscripts v and l denote the vapor phase and liquid phases, respectively. The species vapor pressure is given by the ClausiusClapeyron equation.



6.4.4. Determination of Surface Temperature The surface temperature of the droplet is determined from a heat and mass transfer balance at the interface between the droplet and the surrounding gas. There are two regimes of heat transfer, i.e., heat transfer occurring from the inside of the droplet to the surface, qi , and heat transfer occurring from the outer gas to the surface, qo . The rate of heat transfer balances the required heat for vaporization at the surface (6.68) where L(Ts) is the latent heat of the fuel at the surface temperature, Ts, and rate.



is the mass vaporization



The heat transfer from inside the droplet is modeled as a convective heat-transfer process with internal circulation taken into account. The effective heat transfer coefficient for the outer flux is determined from an approximate solution of the energy equation for the vapor phase with the effects of inter-diffusion and Stefan flow considered. An explicit equation that relates the vaporization rate, , to the temperatures of the droplet and the surrounding gas mixture can be derived as 54 (6.69)



where h i,eff is the heat transfer coefficient inside the droplet, which is determined from the thermal conductivity, λ , and the unsteady equivalent thickness of the thermal boundary layer, r 0 is the droplet radius, Sh is the Sherwood number, Nu is the Nusselt number, C p is the average specific heat of the gas mixture including fuel vapor, K is a correlation factor defined by Ra and Reitz 2003 53 , [ CA ] is the inter-diffusional difference of energy flux between fuel and air, , is the average diffusion coefficient of fuel species, y F0 and y Fsur are the mass fractions of fuel at the interface and far away, respectively, and T sur is the surrounding gas temperature. The rate of mass transport at the droplet surface is calculated using the high mass transfer rate equation with Spalding’s transfer number 64 (6.70) where g



m



is the mass transfer coefficient determined from



transfer number,



, and B



M



is Spalding’s



.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



45



Spray Models



6.4.5. Modeling Boiling, Including Flash Boiling Assuming that 1) the droplet surface distortion can be neglected so that the droplet maintains its spherical shape, 2) there is no sudden break-up of the droplet due to internal phenomena such as microexplosions, and 3) the droplet surface temperature remains at the boiling temperature as long as the assumption of equilibrium is valid, a similar formulation to Equation 6.69 is obtained for the boiling process by setting the surface temperature equal to the boiling temperature, T b , and the surface mass fraction to unity. The boiling formulation is written as (6.71)



where h i,eff is the coefficient for the contribution of heat transfer by internal circulation at the saturation temperature, α sh is the heat transfer enhancement through the effect of nucleation. The vaporization rate calculated based on Spalding’s mass transfer number (Equation 6.70 ) is no longer valid under these conditions.



6.5. Wall Impingement Model When an airborne spray droplet hits a wall surface, the wall impingement model in ANSYS Forte determines the outcome of the collision between the droplet and the wall, depending on the Reynolds number and Weber number of the incident droplet and the surface condition. Four impingement regimes are considered, including stick, rebound, spread and splash (Figure 6.6: Wall impingement regimes. ). For a wetted surface, the stick regime occurs as the impact energy carried by the incident droplet is very low and the droplet adheres to the film in nearly spherical form. As the impact energy increases, the air layer trapped between the droplet and the surface causes low energy loss and the droplet rebounds. When the impact energy is further increased to a certain level, the impinging droplet spreads and merges with the film. Finally, splash occurs at very high impact energy, where the incident droplet flies away from the impinging site and breaks up into many smaller secondary droplets. Figure 6.6: Wall impingement regimes.



The regime transition criteria for a wetted wall as used in 8 are employed in ANSYS Forte: 1. Stick: 46



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Wall Impingement Model 2. Rebound: 3. Spread:



and



4. Splash: Re n and We n are the Reynolds number and Weber number of the incident droplet, respectively. They are defined as: (6.72) (6.73) where ρ i is the liquid density, w0 is the droplet’s velocity normal to the impact surface, D is the droplet diameter, μ is the liquid viscosity, σ is the surface tension. Among these impingement regimes, splash is the most important and complex one. It is commonly observed in engineering applications, such as port fuel injection engines. The splash threshold H cr is proposed by Han et al. 20 in the form of: (6.74) where β is non-dimensional surface roughness (roughness height/incident droplet diameter), δ is nondimensional film thickness (film thickness/incident droplet diameter). The absolute value of surface roughness is exposed as a user input in ANSYS Forte. Typical roughness values range from 0.1 micron (for very smooth surfaces, such as glass surfaces) to 5 micron (for rough surfaces, such as metal surfaces). Splash results in the rebounding of many smaller secondary droplets from the impinging location. The properties of the secondary droplets are modeled statistically, following the formulations derived by Han et al. 20 . These properties are described below.



6.5.1. Mass Fraction of Secondary Droplets Experiments show that a fraction of the incident droplet will remain on the surface after splash. Based on experimental observation, the ratio of the total mass of the secondary droplets m and the incident droplet mass M increases rapidly from zero to a value around 0.75, and then remains constant. The correlation for the mass fraction is given as (6.75) where



6.5.2. Size of Secondary Droplets The size of the secondary droplets d is modeled using a Nukiyama-Tanasawa distribution: (6.76) The mean size d



m



is derived as (6.77)



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



47



Spray Models where D is the incident droplet diameter, ρ



l



is the liquid density, ρ



g



is the ambient gas density.



6.5.3. Velocity of Secondary Droplets The velocity of a secondary droplet is modeled statistically based on experimental observations. The velocity is given as (6.78) where w and v are the normal and tangential velocity components of the secondary droplet, is the unit normal vector to the impinged wall surface, is the unit vector tangent to the impinged surface and in the plane made by and the incident droplet velocity, and is defined as . The normal velocity component w is modeled using a Nukiyama-Tanasawa distribution (6.79) where the mean of the distribution w m is defined as a function of the incident angle α (the angle between the surface normal and the incident droplet velocity) and azimuthal angle ϕ as (6.80) The azimuthal angle is defined as the angle made by the tangential velocity of the secondary droplet and the vector , and it lies in the interval [- π, π ]. It is statistically chosen according to the distribution proposed by Naber and Reitz 43 (6.81) Where P is a random number in the interval [0, 1] and γ is a parameter related to the incident angle α by (6.82) The tangential velocity of the secondary droplet is described by a Gaussian distribution as (6.83)



Where the mean



and variance δ of the distribution are given as (6.84)



According to Han et al. 20 , the constants are set to be A = 0.7 , B = 0.03, c



1



= 0.1, c



2



= 0.02306.



The term ζ in Equation 6.78 is used to describe the radial distance from the edge of the film to the point of impact at any azimuthal angle ϕ . It is assumed to be proportional to the radial velocity (6.85)



48



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Wall Film Model



6.6. Wall Film Model Liquid film may form on a wall as a consequence of sticking, splashing or spreading impingement. Wall film dynamics is influenced by the impinging sprays, the wall conditions, and the gas flow near the wall. The impinging spray affects the film through mass, tangential momentum, and energy addition. The wall affects the film through the no-slip boundary condition and heat transfer. The gas flow affects the film dynamics through tangential shear stresses and the mass and energy exchange between the film and the gas. The wall film model by O’Rourke and Amsden [46 , 47 ] is employed in ANSYS Forte. This model uses a particle numerical method to represent the wall film. In the particle method, each spray particle is tracked as a spherical particle while it moves on the wall surface. In calculating the mass and energy exchange between the wall film and the ambient gas through vaporization, the wall film particles within the same computational cell are converted to a thin film based on the total volume of the liquid and the boundary surface area associated with the corresponding cell.



6.6.1. Wall Film Governing Equations The wall film mass equation is (6.86) where ρ



l



is the liquid density, h is the film thickness,



mean velocity of the film, is the wall velocity. re-entrainment, or vaporization.



is the surface gradient operator,



is the



is the mass source per unit area due to impingement,



The motion of the wall film particles on the wall is governed by the momentum equation (6.87)



where p f is the pressure due to impingement, δ p is the shear stress on the gas side of the wall film, μ temperature T



f



,



and



f l



is the pressure difference across the film, τ w is the film viscosity evaluated at the mean film



are the momentum and liquid mass source terms (per unit area), re-



spectively, due to impingement, is the unit vector normal to the wall, of , is acceleration due to gravity.



is the unit vector in the direction



Based on Han et al. 20 , the incident droplet velocity of a sticking or spreading particle is treated as (6.88) where w p,0 and v p,0 are the normal and tangential velocity of the incident droplet. , ϕ , defined in the wall impingement model section.



,



are



The wall film energy equation is (6.89)



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



49



Spray Models where C



vl



is the liquid specific heat, T



the film surface, T



w



l



is the mean film temperature, T



s



is the gas temperature at



is the wall temperature, λ is the heat conductivity of the liquid film,



are the energy and mass source terms due to impingement. Note that



and



is computed based



on the internal energy of the incident droplets.



6.6.2. Wall Separation Criterion When the wall film approaches an expanding corner, the effect of its momentum may drive the film to separate from the wall surface. The separation criterion proposed by Wegener et al. 76 is employed in ANSYS Forte. This criterion is based on a force balance analysis that considers the surface tension forces at the top and the bottom of the film, and the gravitational force. A critical force ratio is derived by solving the momentum equation of the film. For an expanding angle θ shown in Figure 6.7: Film separation from wall at an expanding corner. , the critical force ratio is given as (6.90) where ρ l is the liquid density, uf is the velocity magnitude of the film, h is film thickness, σ is surface tension, g is gravitational acceleration, and L b is a characteristic breakup length for the thin film breakup, which takes the form of (6.91) For this correlation, the Reynolds number of the film is defined as (6.92) And the Weber number is based on the relative velocity between the gas and liquid film (6.93) where ρ g is the gas density and ug is the gas velocity magnitude. When the critical force ratio becomes greater than one, the inertial force is considered to be great enough to trigger the onset of film separation. Figure 6.7: Film separation from wall at an expanding corner.



50



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 7: Turbulent Flame Propagation Model Spark-ignition engines are characterized by flame initiation near the spark location followed by flame propagation into the engine cylinder. ANSYS Forte employs the G-equation model to track the propagation of fully developed, premixed or partially premixed, turbulent flames for this application. When the flame is initiated by the spark, the ignition-kernel flame has a structure that is typically smaller than the average grid size in the computational mesh. During this time, then, the kernel flame front is first tracked by a group of discrete “particles”. The calculation switches from this kernel flame model to the G-equation model after the flame structure grows bigger than a characteristic flow length scale. The Kernel Flame model and the G-Equation model are described in this chapter.



7.1. Discrete Particle Ignition Kernel Flame Model ANSYS Forte tracks the growth of the ignition kernel by using the Discrete Particle Ignition Kernel Flame (DPIK) model by Fan, Tan, and Reitz [16 , 69 ]. By assuming a spherical-shaped kernel, the flame front position is marked by Lagrangian particles, and the flame surface density is obtained from the number density of these particles in each computational cell. Assuming the temperature inside the kernel to be uniform, the kernel growth rate is: (7.1) where r k is the kernel radius, ρ the kernel region. The plasma velocity S



plasma



u



is the local unburnt gas density, and ρ



k



is the gas density inside



is given as: (7.2)



where ρ



u



and h



u



are the density and enthalpy of the unburned mixture. ρ



density and internal energy of the mixture inside the kernel. η



eff



k



and u



k



are the



is the electrical energy discharge rate,



is the electrical energy transfer efficiency due to heat loss to the spark plug.



and η



eff



are



user-specified inputs (Energy Release Rate and Energy Transfer Efficiency, respectively) in ANSYS Forte. A typical value of η eff is 0.3, as suggested by Heywood 25 . The laminar flame speed in Equation 7.8 was multiplied by a stretch factor, I 0 , which accounts for strain and curvature effects, and the modified correlation is used as the turbulent flame speed, S T , in the kernel stage. C ms is the stretch factor input present in the user interface describing the stretching effect of turbulence on flame speed. A larger value will increase the flame strain cause by turbulence and reduce flame speed. Its effect can be large when flame propagation is weak, for example, under high EGR or lean burn conditions. A typical range is 0.5-1.0. I 0 follows the suggested expression of Herweg et al. 24 :



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



51



Turbulent Flame Propagation Model



(7.3)



Note that curvature effects are also considered in the combustion model by the last term of Equation 7.5 . The chemistry processes in the kernel-growth stage are treated in the same way as in the G-equation combustion model (described in G-equation Model ). Although the transport equation of is not solved here, the field is constructed based on the positions of the kernel particles, thus providing the necessary information for the chemical heat release calculations. The transition from the kernel model to the turbulent G-equation model is controlled by a comparison of the kernel radius with a critical size that is proportional to the locally averaged turbulence integral length scale, viz., (7.4) where Cm 1 is a model constant. Cm 1 is provided as a user input ( Kernel Flame to G-equation Switch Constant) in ANSYS Forte with typical value 2.0.



7.2. Autoignition-Induced Flame Propagation Model Detailed autoignition kinetics is used for ignition. Here we briefly describe the autoignition transition criteria 61 from ignition kernel to G-equation. ANSYS Forte checks two criteria: critical temperature and size of the ignition kernel. After the ignition kernel is formed, the flame propagation model is activated. Computational cells with temperatures greater than the critical temperature become ignition sites. For each of these cells, if the ignition kernel radius is greater than the TLS (turbulent length scale), a G ( x , t ) = 0 surface is initialized. This surface divides the areas where the gas temperature is lower or higher than the critical temperature. The area inside the surface is the area of ignition.



7.3. G-equation Model 7.3.1. General Formulation The G-equation combustion model is based on the turbulent premixed combustion flamelet theory of Peters (2000) 50 . This theory addresses two regimes of practical interest: 1. The corrugated flamelet regime where the entire reactive-diffusive flame structure is assumed to be embedded within eddies of the size of the Kolmogorov length scale η; and 2. The thin reaction zone regime where the Kolmogorov eddies can penetrate into the chemically inert preheat zone of the reactive-diffusive flame structure, but cannot enter the inner layer where the chemical reactions occur. For application of the G-equation model to IC engine applications, this theory was further developed by Tan et al. 69 and by Liang et al. [36 , 37 ]. The G-equation model consists of a set of Favre-averaged level-set equations. This includes the equations for the Favre mean, , and its variance,



52



, as well as a model equation for the turbulent/laminar flame



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



G-equation Model surface area ratio σ T . Application of the equation for the turbulent/laminar flame surface area ratio results in an explicit expression for the turbulent flame speed . Together with the Reynolds-averaged Navier-Stokes equations and the turbulence modeling equations, these provide a complete set of equations to describe premixed turbulent flame-front propagation. The equation set used in ANSYS Forte is: (7.5) (7.6) where ∇ || denotes the tangential gradient operator; is the fluid velocity; is the velocity of the moving vertex; ρ u and ρ b are the average densities of the unburned and burned mixtures, respectively; D T is the turbulent diffusivity; is the Favre mean flame front curvature; c s , a 4 , b 1 , and b 3 are modeling constants (cf. ref. 50 ] ); and are the Favre mean turbulent kinetic energy and its dissipation rate from the RNG k - ε model; u ′ is the turbulence intensity.



7.3.2. Laminar and Turbulent Flame Speed Correlations The prediction of the turbulent burning velocity plays a crucial role in the modeling of SI engine combustion. The laminar flame speed is one of the most important scaling factors in most of the published correlations for the turbulent flame speed.



7.3.2.1. Laminar Flame Speeds Laminar flame speed is an intrinsic property of a premixed fuel/air mixture for given unburned temperature, pressure, and composition (species mass fractions). The composition can be characterized using two derived parameters: equivalence ratio and diluent fraction. Equivalence ratio can be defined using the concentrations of C, H, and O atoms as: (7.7) Diluent fraction can be estimated by assuming the inert diluent mixture as complete combustion products. The diluent may consist of CO 2 , H 2 O, N 2 , as well as excess O 2 for lean mixtures. It is a function of equivalence ratio, mass fractions of CO 2 , H 2 O, and the C/H/O ratio in the fuel species. In ANSYS Forte, two options are available to specify laminar flame speeds. One is using the power-law correlations, the other is through table look-up. For details on using table look-up, see Laminar Flame Speed Through Table Look-up , Laminar Flame Speed Through Table Look-up , on Laminar Flame Speed Through Table Look-up .



7.3.2.1.1. Power-law Formulation The “power law” formula 40 can be written as: (7.8) where the subscript ref means the reference condition (typically at 298 K and 1 atm) and the superscript 0 means the flame is planar and unstretched. F dil is a factor accounting for the diluent’s effect.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



53



Turbulent Flame Propagation Model



7.3.2.1.2. Laminar Flame Speed at Reference State The reference flame speed can be specified using two different formulas. The first formula is: (7.9) Values for , and for selected fuels are listed inTable 7.1: Values for B the Metghalchi formula according to Metghalchi et al. 21 . Table 7.1: Values for B



m



Fuel



, B φ , and φ



m



m



, B φ , and φ



m



in



in the Metghalchi formula



(cm/sec)



(cm/sec)



methanol



36.9



-140.5



1.11



propane



34.2



-138.7



1.08



iso-octane



26.3



-84.7



1.13



gasoline



30.5



-54.9



1.21



Unfortunately, Equation 7.9 gives negative flame speeds for very lean or very rich mixtures. One practical solution to this problem is to follow the expression proposed by Gülder 14 in which the flame speed will never be driven negative, viz., (7.10) where , , , and are data-fitting coefficients. For iso-octane, a group of values for the coefficients in Equation 7.10 can be obtained by correlating the data of Metghalchi et al. within the range 0.65 < f < 1.6. For example, then, coefficients obtained for iso-octane's reference flame speed are: =26.9; = -0.134; =3.86; =1.146



7.3.2.1.3. Temperature and Pressure Dependencies The exponents and in Equation 7.8 describe the temperature and pressure dependencies. These exponents were treated as generic for all fuel types in Metghalchi et al. 21 and correlated as functions of equivalence ratio as: (7.11) (7.12) Values obtained using least-squares fits are given as = 2.98;



= -0.8;



= -0.38;



= 0.22;



=1 =1



An additional set of parameters for gasoline were reported by Rhodes et al. 35 : = 2.4; = -0.357;



54



= -0.271; = 0.14;



= 3.51 = 2.77



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



G-equation Model



7.3.2.1.4. Diluent Effect This effect by the diluent is usually accounted for by a term such as F sions for this factor take the form:



dil



in Equation 7.8 . The expres(7.13)



or (7.14) depending on whether mole fraction or mass fraction is selected, where is the mole fraction of diluent; is the mass fraction of diluent; are empirical constants. Fitted values from the literature are: (7.15) (7.16)



7.3.2.1.5. Laminar Flame Speed Through Table Look-up The laminar flame speed value returned by table look-up is the unstretched flame speed, which is comparable to the value calculated using Equation 7.8 . It does not contain the stretch factor I 0 described in Equation 7.3 . The details of working with flame speed table lookup is described in FlameSpeed Table Editor , Flame-Speed Table Editor , in the ANSYS Forte User's Guide .



7.3.2.2. Turbulent Flame Speeds The ratio between turbulent flame speed and laminar flame speed is given as (7.17)



in which the term



, called a progress variable, takes the form



is the laminar flame speed; l I and l flame thickness, respectively. The term I form



F



are the turbulence integral length scale and the laminar P in Equation 7.17 is a progress variable, which takes the (7.18)



in Equation 7.18 is the Flame Development Coefficient, a model constant relevant to spark ignition. Physically, the progress variable models the increasingly disturbing effect of the surrounding eddies on the flame front surface as the ignition kernel grows from the laminar flame stage into the fully developed turbulent stage. is needed only when a spark exists. The model constant is available in the User Interface under Models > Spark Ignition. It is provided with an appropriate default value that is applicable in most cases. Increasing the value of this variable will expedite the transition from the laminar kernel flame to the fully developed turbulent flame. See Spark Ignition in the ANSYS Forte User's Guide . The other model constants, b 1, b 3 and a 4 are generic for any turbulent flame (both kernel flame or fully-developed flame). They need to be provided as long as the G-equation model is used. Suggested Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



55



Turbulent Flame Propagation Model values were provided by Peters 2000 28 , and options for user definition are available in the User Interface. Here we present brief explanations of these constants: : Coefficient in equation , where is fully-developed turbulent flame speed , is turbulence velocity, defined as sqrt(2/3*k) with k being turbulent kinetic energy. The suggested value is 1.5, which was obtained by fitting experimental data. : Coefficient in equation , where and are turbulent diffusivity and laminar diffusivity, respectively. Suggested value is 1.0, which is based on Damkohler's publication in 1940 8 . A subsequent DNS study by Wenzel 77 suggested a similar value, 1.07. : Coefficient in equation , where is turbulent diffusivity bulence integral length scale. Suggested value is 0.78.



, is turbulence velocity, is tur-



7.3.2.3. Turbulent Flame Brush Thickness The turbulent flame brush thickness



is defined as the square root of the Favre variance of G: (7.19)



For fully developed flames, the turbulent flame brush thickness is correlated with turbulence integral length scale as [39], (7.20) where b



2



= 1.78.



In the context of spark ignition, the progress variable I P defined in Equation 7.18 is used to describe the growth of flame brush thickness from the laminar stage to the fully developed stage. Thus, the flame brush thickness in spark ignition engine is modeled as (7.21)



7.3.3. Flame-front Heat Release Calculation For turbulent premixed or partially premixed flames in SI engines, the heat release due to turbulent flame-front propagation constitutes a major portion of the total in-cylinder heat release. For this reason, we refer to this as the “ primary heat release.” The ANSYS Forte G-equation model assumes a flame structure shown schematically in Figure 7.1: Numerical description of the turbulent flame structure . For the heat-release calculation, the ANSYS Forte model builds on this assumption using the method described by Liang et al. 37 . This method uses the sub-grid scale unburnt/burned volumes of the flame-containing cells and assumes that the mean flamefront surface cuts every flame-containing cell into two parts, an unburned volume V u and a burned volume V b , as shown in Figure 7.1: Numerical description of the turbulent flame structure . As the mean flame front sweeps forward, the mixture within the sweeping volume tends to reach a local instantaneous thermodynamic equilibrium, following a constant pressure, constant enthalpy process. The pressure is assumed to be homogeneous across the flame within the cell, consistent with deflagration wave theory. The sub-grid scale volumes are tracked for every time step based on the coordinate information of the cell vertices and the flame surface piercing points. 56



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



G-equation Model Figure 7.1: Numerical description of the turbulent flame structure



The species conversion rate and the associated primary heat release rate are calculated by first determining the equilibrium concentration and adiabatic flame temperature under constant pressure and constant-enthalpy constraints. Then for a flame-containing cell, we assume the mean flame front sweeps throughout the unburned volume V u within time step dt . Taking into account the portion of the mixture within the swept volume V s that reaches the equilibrium state, the species conversion rate of species k for time step dt reads: (7.22) where ρ and ρ k are cell-averaged densities of the total mixture and species k and where the swept volume is approximated by (7.23) where A



f



is the mean flame-front area within the cell, and S



T



the turbulent flame speed.



Finally the conversion rate of species k can be written as (7.24) Numerical integration of Equation 7.24 using an explicit Euler method gives the updated species density .



7.3.4. Post-flame and End Gas Kinetics In ANSYS Forte, the species conversion and heat release in the regions outside the turbulent flame brush are modeled as chemical-kinetics-controlled processes. Specifically, the computational cells in the post-flame zone and the end-gas zone are simulated in the same manner as described in Chemistry Solver Options . In these zones, the chemical processes are assumed to be primarily controlled by the reaction pathways defined by the chemical kinetic mechanisms. The turbulence-kinetics interaction model can also be included for this portion in these regions.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



57



Turbulent Flame Propagation Model



7.3.5. Flame Quench Model For partially premixed turbulent flames, the unburned mixture at certain locations can be leaned out with excess air and therefore the flame may stop propagating due to charge stratification. This section briefly describes how the method in ANSYS Forte for modeling such local flame-quenching phenomena. For laminar flames, the flammability property of an unburned mixture is mainly governed by a balance between the heat transfer to the preheat zone due to chemical heat release and the heat loss from the preheat zone to the unburned mixture, while for turbulent flames, disturbances from turbulence have significant additional influence on the heat transfer balance. In the present combustion model, local flame quenching is modeled by examining whether the local flame condition crosses the border between the thin reaction-zone regime and the broken reaction-zone regime. Local laminar flame thickness will increase with a decrease of laminar flame speed , thus resulting in a proportional increase of the inner layer thickness. If the local laminar flame inner layer thickness is large enough so that the inner layer can be disturbed by the Kolmogorov eddies, the chemical reactions in the inner layer will cease due to excess heat and active species losses to the preheat zone, resulting in local flame quenching. Therefore, a comparison between the inner layer thickness and the Kolmogorov length scale is used as the criterion for local flame quenching, i.e., flame quenching occurs when the relation (7.25) is satisfied, where is determined based on the turbulence model, and a typical value of 1.0.



58



is a model constant with



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 8: Method of Moments The Particle Tracking option allows inclusion of dispersed, condensed-phase material in the form of particles. The Particle Tracking feature accounts for particle formation and destruction, and can be used to determine properties about the amount of particulates in the system. To describe some of the theory behind the Particle Tracking feature, we first need to define particles in this context. Like gaseous species, each type of particle corresponds to a symbolic representation in the chemistry set and has properties associated with it. To form particles from the gas phase, one or more nucleation events need to be defined to identify the particle properties at inception. Chemical composition and thermodynamic properties of the particle are defined in this way. The definition of the particle and its properties are described in Chapter 6 in the Chemkin-Pro Input Manual The nucleation “reaction” is used to define how the particle (or nucleus) is created from gas-phase species. The nucleation reaction is an irreversible reaction that provides the particle inception rate and defines the size and the surface coverage of the nucleus. After the nuclei are formed, they start to interact with each other as well as with the gas mixture around them. While particle-particle interactions such as coagulation are non-chemical processes, interactions between particles and the surrounding gas mixture can result in chemical processes taking place on the particle surface. These surface processes might result in mass growth or reduction of the particle or might just simply recondition the particle surface. To include the effects from all these surface processes, a surface mechanism is needed to describe all surface reactions and associated surface species on the particle. The Particle Tracking feature then determines the impact of individual surface reactions on the particle sizes from the expression of the surface reactions The sections below describe the concept and theory behind the Particle Tracking computations along with implementation and numerical solution considerations.



8.1. Description and Properties of a Particle Population 8.1.1. Moments of Particle-Size Distribution Functions The application of the method of moments to soot-particle formation was first reported by Frenklach and coworkers [17 , 6 , 18 ]. This method describes the average properties of a particle population. The method of moments tracks the evolution of an aerosol system by moments of its particle-size distribution function. The use of moments rather than the actual form of the particle-size distribution function overlooks variations among individual particles in the aerosol system. Also, the actual particle-size distribution function cannot be derived from its moments unless an assumption is made regarding the form of the distribution function. Since in many practical applications only the average properties of the aerosol system are sought, the history and properties of individual particles may not be important such that this approximation can still be quite useful. The loss in details of the particle-size distribution function due to the use of method of moments is compensated by computational speed with which the moments are calculated, which reduces the demand for computing resources. Detailed descriptions of the method of moments are reported elsewhere by Frenklach et al. [17 , 6 , 18 ]. Without making any assumptions about the form of the particle-size distribution function, the method of moments can provide overall properties of a particle system such as number density, total particle Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



59



Method of Moments volume fraction, total particle surface area density, and average particle size. To express these overall properties in terms of particle-size moments, we first define the particle-size moments. Given a particle-size distribution function , where represents a measure of particle size, e.g., particle mass or particle diameter, the r -th moment of this particle-size distribution function is defined in Equation 8.1 . (8.1)



The Particle Tracking module uses particle class, which is defined as the number of bulk species molecules in a particle core, as the measure of particle size. Both particle mass, and particle volume, are proportional to particle class. Because particle classes are discrete numbers, the number of class particles can be represented by a discrete function and Equation 8.1 is equivalent to Equation 8.2 . (8.2) In the following sections, the summation notation of Equation 8.2 will be used in formulation expressions and derivations.



8.1.2. Total Particle Number of a Particle Population From Equation 8.1 we see that the zero-th moment is the total particle number of the particle population, as shown in Equation 8.3 . (8.3)



8.1.3. Total and Average Particle Mass Since the mass of a particle is proportional to its class, the total mass of a particle population can be calculated as : (8.4) Note that is the mass of a bulk species molecule comprising the particle core and is constant. The total mass of a particle population can be expressed in a simple form as Equation 8.5 . (8.5)



Here, is the first moment of the discrete particle-size distribution function . The average particle mass of the population can be obtained by dividing the total particle mass by the number of particles and is given in Equation 8.6 , (8.6)



60



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Description and Properties of a Particle Population where



is the average particle class of the population.



8.1.4. Total and Average Geometric Properties of a Particle Population Similarly, the total particle volume is found to be (8.7) where [g/cm 3 ] is the bulk density of the particle core. By assuming that the particles are spherical, an average particle diameter can be obtained as (8.8) The total sphere-equivalent surface area of a particle population is given as (8.9) With these equations, all basic properties of a particle population are defined for the moment method. Several sections in the Chemkin-Pro Theory Manual, Method of Moments chapter, discuss the kinetics models used by the Particle Tracking feature to describe formation, growth, reduction, and interaction of the particles.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



61



62



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 9: Engine Crevice Model This model can be used to capture crevice effects when the crevice region around the piston is not explicitly represented in the computational mesh.



9.1. Diagram and Assumptions The schematic for the crevice model is given in Figure 9.1: Schematic diagram of the crevice volumes and planes used in the model , which reproduces Figure 15 of Namazian and Heywood 45 . The flow between regions 1 through 5 is pressure driven. In the model the pressure in region 1 is defined to be equal to the cylinder pressure, and zone 5 is at atmospheric pressure. The continuity equations, which give expressions for the pressure in zones 2, 3, and 4, are given in equation (6) of Namazian and Heywood 45 , and are reproduced here for completeness. Figure 9.1: Schematic diagram of the crevice volumes and planes used in the model



(9.1) (9.2) (9.3) where the subscript ‘o’ refers to a reference condition such as the one at intake valve closing, m is mass, p is pressure, t is time, and the double subscript means the flow from the first to the second



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



63



Engine Crevice Model number. For example, in Equation 9.1 , the 12 subscript refers to mass flow from region 1 above the top ring to region 2 behind the first ring. To determine the mass flow rates in equations, two types of flow conditions are used in the theory: isothermal and isentropic. They are used in the following situations:



9.1.1. Isentropic Isentropic flow is assumed when the flow is at high velocity through a small orifice. This holds for the flow through the two ring gaps at all times. It is also true for the flow over the rings when the ring lifts off of the seat. For example, if the top ring were “airborne”, there would be a small opening connecting zones 1 and 2 and another small opening between zones 2 and 3. The flow in this case will move quickly because of the lower restriction between zones 1 and 3, so it is assumed to be isentropic. The governing equations for this flow are as follows: (9.4) for the mass flow rate where C d is the discharge coefficient, A g is the cross-sectional area available to the flow, ρ is the local gas density, c is the local speed of sound and η is defined as (9.5)



where γ is the ratio of specific heats (usually 1.4) and the subscript, ‘i’, on the pressure terms denotes the region where the properties are calculated, with acceptable values running from 1 to 4. Equation 9.4 is equation (8) and Equation 9.5 is equation (2) in Namazian and Heywood 45 .



9.1.2. Isothermal This is the case when the flow is relatively slow. This occurs when the flow is into or out of regions 2 and 4, with the ring sealing one of the upper or lower surfaces of the groove in the piston. In essence, these two zones have one inlet unless the ring is not touching a surface, and therefore offer higher resistance to the flow. In ANSYS Forte, the isothermal flow is also considered to be laminar due to the small gaps in the rings. The governing equation in this case reduces to (9.6) where A is the area perpendicular to the mass flow, h is the distance between the ring and the ring groove, the subscripts u and d represent upstream and downstream pressure, W r is the width of the ring (the distance from the wall to where the ring touches zone 2), μ is the gas viscosity, R is the ideal gas constant, and T is the temperature of the gas. Equation 9.6 is reproduced from (9b) in Namazian and Heywood 45 .



64



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Input Parameters



9.2. Input Parameters Table 9.1: Crevice model input parameter labels in user interface. Label in Interface



Description



Crevice Height Top-to-Land



Height of the crevice region between the top of the crevice and the first compression ring.



Crevice Thickness



Distance between the Piston edge and the Cylinder Liner, or thickness of the crevice region.



Top Ring Gap Area



Cross-sectional flow area for crevice flow between the cylinder and the region beyond the first compression ring.



Second Ring Gap Area



Cross-sectional flow area for crevice flow between the first o-ring and the region beyond the second compression ring.



Volume Behind Top Ring



Volume behind the top compression ring.



Volume Between Ring 1 and 2



Volume of the region between the first (top) compression ring and the second compression ring.



Volume Behind Second Ring



Volume behind the second compression ring.



First Ring Thickness



Thickness of the first (top) compression ring.



Second Ring Thickness



Thickness of the second compression ring.



First Ring Width



Width (lateral) of the first compression ring.



Second Ring Width



Width (lateral) of the second compression ring.



Mass of the First Ring



Mass of the first (top) compression ring.



Mass of the Second Ring



Mass of the second compression ring.



Bottom Clearance - First Ring



Clearance distance on the bottom side of the first compression ring.



Bottom Clearance - Second Ring



Clearance distance on the bottom side of the second compression ring.



Top Clearance - First Ring



Clearance distance on the top side of the first compression ring.



Top Clearance - Second Ring



Clearance distance on the top side of the second compression ring.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



65



66



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



Chapter 10: References 1. Abani, N. and Reitz, R. D.,“Unsteady Turbulent Round Jets and Vortex Motion,” Physics of Fluids, 19(1): 125102, 2007, doi:10.1063/1.2821910. 2. Abani, N., Kokjohn, S., Park, S. W., Bergin, M., Munnannur, A., Ning, W., Sun, Y., and Reitz, R. D.,“An Improved Spray Model for Reducing Numerical Parameter Dependencies in Diesel Engine CFD Simulations, SAE Technical Paper 2008-01-0970,” SAE Technical Paper Series, 2008. 3. Amsden, A. A., J. D. Ramshaw, L. D. Cloutman, and P. J. O'Rourke,“Extensions to the KIVA Computer Program,” Los Alamos National Laboratory Report LA-10534-MS, 1985. 4. Amsden, A. A., P. J. O'Rourke, and T. D. Butler,“KIVA-II: A Computer Program for Chemically Reactive Flows with Sprays,” Los Alamos National Laboratory Report LA-11560-MS, Los Alamos National Laboratory, 1989. 5. Amsden, A. A.,“KIVA-3V: A Block-structured KIVA Program for Engines with Vertical or Canted Valves,” Los Alamos National Laboratory Report LA-13313-MS, Los Alamos National Laboratory, 1997. 6. Appel, J., H. Bockhorn, and M. Frenklach, Combustion and Flame, 121:122-136, 2000. 7. Artelt, C., H.-J. Schmid, and W. Peukert, "On the relevance of accounting for the evolution of the fractal dimension in aerosol process simulations," Journal of Aerosol Science, 34: 511-534, 2003. 8. Bai, C. and A. Gosman,“Mathematical Modeling of Wall Films Formed by Impinging Sprays, SAE Technical Paper 960626,” SAE Technical Paper Series 1996. 9. Beale, J.C., R.D. Reitz, Modeling spray atomization with the Kelvin-Helmholtz / Rayleigh-Taylor hybrid model. Atomization and Sprays, 9, 623-650, 1999. 10. Bellman R. and Pennington R. H.,“Effects of Surface Tension and Viscosity on Taylor Instability”, Quarterly of Applied Mathematics, 12(2): 151-162, 1954. 11. CHEMKIN-PRO, Reaction Design: San Diego, 2008. 12. CHEMKIN Theory Manual, Reaction Design, San Diego, CA, USA, 2009. 13. Chou, C.-P., P. Ho, and E. Meeks, 30th International Symposium on Combustion, Work-In-Progress Poster Section, 2F1-22, July 25-30, 2004. 14. Crowe C., Sommerfeld M., and Tsuji Y., Multiphase Flows with Droplets and Particles, CRC Press, Boca Raton, ISBN-10: 0849394694, 1998. 15. Damkohler, G., Z.,“Development of Ignition and Combustion Model for Spark-Ignition Engines, SAE Technical Paper 2000-01-2809,” SAE Technical Paper Series, 2000. 16. Fan, L., and Reitz, R. D., Elektrochem. Angew. Phys. Chem., 46: 601, 1940. 17. Frenklach, M. and Wang, H., in Soot Formation in Combustion: Mechanisms and Models, H. Bockhorn (Ed.), Springer-Verlag, pp. 165-192, 1994. Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



67



References 18. Frenklach, M., and S.J. Harris, Journal of Colloid and Interface Science, 118:252-261, 1987. 19. Gulder, O. L.,“Correlations of Laminar Combustion Data for Alternative S.I. Engine Fuels, SAE Technical Paper 841000,” SAE Technical Paper Series, 1984. 20. Han, Z., Z. Xu, and N. Trigui,“Spray/Wall Interaction Models for Multidimensional Engine Simulation,” International Journal of Engine Research, 2000. 1:1. 21. Han, Z., and R. D. Reitz,“Turbulence Modeling of Internal Combustion Engines Using RNG k - ε Models,” Combustion Science and Technology, 106, pp. 267-295, 1995. 22. Han, Z., and R. D. Reitz,“A Temperature Wall Function Formulation for Variable-density Turbulence Flows with Application to Engine Convective Heat Transfer Modeling,” International Journal of Heat Mass Transfer, 40, 3, pp. 613-625, 1997. 23. Haworth, D. C.,“A Review of Turbulent Combustion Modeling for Multidimensional In-Cylinder CFD, SAE Technical Paper 2005-01-0993,” SAE Technical Paper Series, 2005. 24. Herweg, R. and Maly, R. R.,“A Fundamental Model for Flame Kernel Formation in S.I. Engines, SAE Technical Paper 922243,” SAE Technical Paper Series, 1992. 25. Heywood, J. B., Internal Combustion Engine Fundamentals, McGraw-Hill, 1988. 26. Hiroyasu, H. and Kadota, T.,“Models for Combustion and Formation of Nitric Oxide and Soot in DI Diesel Engines, SAE Technical Paper 760129,” SAE Technical Paper Series, 1976. 27. Hirt, C. W., A. A. Amsden, and J. L. Cook,“An arbitrary Lagrangian-Eulerian computing method for all flow speeds,” Journal of Computational Physics, 14, 227-253, 1974. 28. Kazakov, A., and Frenklach, M.,“Dynamic Modeling of Soot Particle Coagulation and Aggregation: Implementation with the Method of Moments and Application to High-Pressure Laminar Premixed Flames”, Combustion and Flame, 114:484-501, 1998. 29. Koch, W., and Friedlander, S.K.,“The Effect of Particle Coalescence on the Surface Area of a Coagulating Aerosol”, Journal of Colloid and Interface Science, 140(2):419-427, 1990. 30. Kong, D., R.K. Eckhoff, and F. Alfert,“Auto-ignition of CH4/air, C3H8/air, CH4/C3H8/air and CH4/CO2/air using a 1 l ignition bomb,” Journal of Hazardous Materials, 40, pp. 69 - 84, 1995. 31. Kong, S.-C. and R. D. Reitz,“Use of Detailed Chemical Kinetics to Study HCCI Engine Combustion With Consideration of Turbulent Mixing Effects,” Transactions of the American Society of Mechanical Engineers, 124, pp. 702-707, 2002. 32. Levich, V. G., Physicochemical Hydrodynamics, ed. N.R. Amundson, Englewood Cliffs, N.J.: Prentice-Hall, 1962 33. Liang, L., Stevens, J.G., Farrell, J.T.“A Dynamic Adaptive Chemistry Scheme for Reactive Flow Computations.” Proceedings of the Combustion Institute, 32, pp. 527-534, 2009. 34. Liang, L., Stevens, J.G., Ram an, S., Farrell, J.T.,“The Use of Dynamic Adaptive Chemistry in Combustion Simulation of Gasoline Surrogate Fuels.” Combustion and Flame, 156:7, 1493-1502, 2009. 35. Liang, L., Stevens, J.G., Farrell, J.T.“A dynamic multi-zone partitioning scheme for solving detailed chemical kinetics in reactive flow computations,” Combustion Science and Technology, 181:11, 1345-1371, 2009.



68



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



36. Liang, L. and R. D. Reitz.“Spark Ignition Engine Combustion Modeling Using a Level Set Method with Detailed Chemistry, SAE TEchnical Paper 2006-01-0243,” SAE Technical Paper Series, 2006. 37. Liang, L., Reitz, R.D., Iyer, C.O. and Yi, J.,“Modeling Knock in Spark-Ignition Engines Using a G-equation Combustion Model Incorporating Detailed Chemical Kinetics, SAE Technical Paper 2007-01-0165,” SAE Technical Paper Series, 2007. 38. Lu, T. and C.K. Law,“A directed relation graph method for mechanism reduction,” Proceedings of the Combustion Institute, 30: 1333-1341, 2005. 39. Marr, J.A., PhD. Thesis, Dept. of Chemical Engineering, MIT, 1993. 40. Metghalchi, M., and Keck, J. C.,“Burning Velocities of Methanol, Ethanol and iso-Octane-Air Mixtures”, 19th Symposium (International) on Combustion/ The Combustion Institute, p. 275, 1983. 41. Munnannur, A. and Reitz, R.D.,“A Comprehensive Collision Model for Multi-Dimensional Engine Spray Computations,” Atomization and Sprays, 19:7 597-619, doi:10.1615/AtomizSpr.v19.i7, 2009. 42. Munnannur, A.,“Droplet Collision Modeling in Multi-Dimensional Engine Spray Computations,” Ph.D. Thesis, Mechanical Engineering Department, University of Wisconsin-Madison, 2007. 43. Naber, J., and R. D. Reitz,“Modeling Engine Spray/Wall Impingement, SAE Technical Paper 880107,” SAE Technical Paper Series, 1988. 44. Nagle, J., and Strickland-Constable, R. F.,“Oxidation of Carbon Between 1000-2000 °C,” Proceedings of the Fifth Carbon Conference, Vol. 1, Pergamon, New York, pp. 154-164, 1962. 45. Namazian, M., and J. B. Heywood,“Flow in the Piston-Cylinder-Ring Crevices of a Spark-Ignition Engine: Effect on Hydrocarbon Emissions, Efficiency and Power, SAE Technical Paper 820088,” SAE Technical Paper Series, 1982. 46. O'Rourke, P.J. and A.A. Amsden,“A Particle Numerical Model For Wall Film Dynamics in Port-Injected Engines, SAE Technical Paper 961961,” SAE Technical Paper Series, 1996. 47. O'Rourke, P.J. and A.A. Amsden, A Spray/Wall Interaction Submodel for the KIVA-3 Wall Film Model, SAE Technical Paper 2000-01-0271,” SAE Technical Paper Series, 2000. 48. O'Rourke, P. J., and A. A. Amsden,“The TAB Method for Numerical Simulations of Spray Droplet Breakup, SAE Technical Paper 872089,” SAE Technical Paper Series, 1987. 49. Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation: Washington D. C., 1980. 50. Peters, N., Turbulent Combustion, Cambridge University Press: Cambridge, UK, 2000. 51. Pracht, W. E.,“Calculating three-dimensional fluid flows at all speeds with an Eulerian-Lagrangian computing mesh,” Journal of Computational Physics, 17:2 132-159, 1975. 52. Promraning, E., Rutland, C. J.“Dynamic One-equation Nonviscosity Large-eddy Simulation Model,” AIAA J., 40:4 689-701, 2002. 53. Ra, Y., Reitz, R.D.,“The application of a multi-component vaporization model to gasoline direct injection engines,” International Journal of Engine Research, 4, 193-218, 2003.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



69



References 54. Ra, Y., Reitz, R. D.,“A vaporization model for discrete multi-component fuel sprays, International Journal of Multiphase Flow,” 35, 101-117, doi:10.1016/j.ijmultiphaseflow.2008.10.006, 2009. 55. Reid, R.C., J.M. Prausnitz, and B.E. Poling, The Properties of Gases and Liquids, 4th Edition, McGraw-Hill, New York, 1987. 56. Reitz, R. D.,“Modeling Atomization Processes in High-Pressure Vaporizing Sprays”, Atomization and Sprays, 3: 309-337, 1987. 57. Reitz, R. D. and Diwakar, R.,“Effect of Drop Breakup on Fuel Sprays, SAE Technical Paper 860469,” SAE Technical Paper Series, 1986. 58. Rhodes, D. B., and Keck, J.,“Laminar Burning Speed Measurements of Indolene-Air-Diluent Mixtures at High Pressures and Temperatures,” SAE Technical Paper 950047, SAE Tehcnical Paper Series, 1985. 59. Rutland, C. J.,“Large-eddy Simulations for Internal Combustion Engines – A Review,” International Journal of Engine Research, 12, pp. 1-31, 2011. 60. Schmidt, D. P., I. Nouar, P. K. Senecal, C. J. Rutland, J. K. Martin, R. D. Reitz, and J. A. Hoffman,“Pressure-Swirl Atomization in the Near Field, SAE Technical Paper 1999-01-0496,” SAE Technical Paper Series, 1999. 61. Singh, S., Liang, L., and Reitz, R. D.,“Development of a flame propagation model for dual-fuel partially premixed compression ignition engines,” Journal of Engine Research, 7, pp. 65-75, 2006. 62. Squire, H. B.,“Investigation of the Instability of a Moving Liquid Film,” British Journal of Applied Physics, 4, p. 167, 1953. 63. Schmidt, D. P. and Rutland, C. J.,“A New Droplet Collision Algorithm,” Journal of Computational Physics, 164:62-80, doi:10.1006/jcph.2000.6568, 2000. 64. Sirignano, W. A.,“Fuel droplet vaporization and spray combustion,” Progress in Energy and Combustion Science, 9, 291-322, 1983. 65. Smagorinsky, J.,“General Circulation Experiments with the Primitive Equations,” Mon. Weather Rev., vol. 99, no. 3, pp. 99-164, 1963. 66. Su, T. F., Patterson, M. A., Reitz, R. D., and Farrell P. V.,“Experimental and Numerical Studies of High Pressure Multiple Injection Sprays, SAE Technical Paper 960861,” SAE Technical Paper Series, 1996. 67. Sun, Y., Reitz, R. D.,“Modeling Low-Pressure Injections in Diesel HCCI Engines,” ILASS Americas, 20th Annual Conference on Liquid Atomization and Spray Systems, Chicago, 2007. 68. Tan, Z.,“Multi-Dimensional Modeling of Ignition and Combustion in Premixed and DIS/CI (Direct Injection Spark/Compression Ignition) Engines.” Ph.D. Thesis, University of Wisconsin-Madison, 2003. 69. Tan, Z.; Reitz, R.D.,“An Ignition and Combustion Model for Spark Ignition Engine Multi-dimensional Modeling,” Combustion and Flame, 145, pp. 1-15, 2006. 70. Thompson, P. A., Compressible-fluid Dynamics, McGraw-Hill: New York, 1972. 71. Tsang, C.-W., Trujillo, M.F., Rutland, C. J.,“Large-eddy Simulation of Shear Flows and High-speed Vaporizing Liquid Fuel Sprays,” Computational Fluids, 105, pp. 262-279, 2014.



70



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



72. Tsang, C.-W., Wang, Y., Wang, C., Shelburn, A., Liang, L., Puduppakkam, K., Modak, A., Naik, C., Meek, E., Rutland, C.,“Evaluation and Validation of Large-Eddy-Simulation (LES) for Gas Jet and Sprays, SAE Technical Paper 2018-01-0844,” SAE Technical Paper Series, 2017. 73. Vishwanathan, G., and Reitz, R. D.,“Numerical Predictions of Diesel Flame Lift-Off Length and Soot Distributions Under Low Temperature Combustion Conditions, SAE Technical Paper 2008-01-1331,” SAE Technical Paper Series, 2008. 74. Wahiduzzaman, S., and C. R. Ferguson,“Convective Heat Transfer from a Decaying Swirling Flow within a Cylinder,” 8th International Heat Transfer Conference, paper #86-IHTC-253, August 1986. 75. Wang, Y., Ge, H.-W., and Reitz, R.D.,“Validation of Mesh- and Timestep-Independent Spray Models for MultiDimensional Engine CFD Simulation, SAE Technical Paper 2010-01-0626,” SAE Technical Paper Series, 2010. 76. Wegener, J.L., et al.,“A Criterion for Predicting Shear-driven Film Separation at an Expanding Corner with Experimental Validation,” 21st Institute for Liquid Atomization and Spray Systems Annual Conference, 2008. 77. Wenzel, H., "Turbulent premixed combustion in the laminar flamelet and the thin reaction zones regime," Annual Research Briefs, Center for Turbulence Research, 1997: 237-252, 1997. 78. Yakhot, V., and S. A. Orszag,“Renormalization Group Analysis of Turbulence: I. Basic Theory,” Journal of Scientific Computing, 1:3, 1986. 79. Yang, X., Takamoto, Y., Okajima, A., Obokata, T. and Long, W.,“Comparison of Computed and Measured High-Pressure Conical Diesel Sprays, SAE Technical Paper 2000-01-0951,” SAE Technical Paper Series, 2000.



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.



71



72



Release 18.2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information of ANSYS, Inc. and its subsidiaries and affiliates.