A Model for Fuel Spray Formation with Atomizing Air

The formation of a liquid spray emanating from a nozzle in the presence of atomizing air was studied using a computational model approach that accounted for the deformation and break up of droplets. Particular attention was given to the formation of sprays under non-swirling flow conditions. The instantaneous fluctuating fluid velocity and velocity gradient components were evaluated with the use of a probability density function (PDF)-based Langevin equation. Motions of atomized fuel droplets were analyzed, and ensemble and time averaging were used for evaluating the statistical properties of the spray. Effects of shape change of droplets, and their breakup, as well as evaporation, were included in the model. The simulation results showed that the mean-square fluctuation velocities of the droplets vary significantly with their size and shape. Furthermore, the mean-square fluctuation velocities of the evaporating droplet differed somewhat from non-evaporating droplets. Droplet turbulence diffusivities, however, were found to be close to the diffusivity of fluid point particles. The droplet velocity, concentration, and size of the simulated spray were compared with the experimental data and reasonable agreement was found.


Introduction
Understanding the formation and subsequent dynamics of turbulent liquid sprays are of considerable interest to the design and operation of jet and rocket engines.However, the large range of length and time scales over which these liquid sprays form and evolve lead to considerable computational complexities.Therefore, there have been numerous studies in the literature on various aspects of spray processes.
Extensive reviews on turbulent spray combustion were reported by Faeth [1] and Law [2] and more recently by Jiang et al. [3] and Gutheil [4] and Jenny et al. [5].Shang et al. [6] and Berlemont et al. [7] presented the procedure for Lagrangian simulation of droplet evaporation in turbulent flows.Kvasnak et al. [8] included the effect of droplet deformation and breakup to their models.Recent advances for using Eulerian-Lagrangian approach for spray modelling was reported by Kolakaluri et al. [9], Keller et al. [10], and Doisneau et al. [11], among others.Hoyas et al. [12] performed two-dimensional Eulerian-Lagrangian spray simulations in the near-nozzle region.
Kvasnak et al. [8] examined the formation of liquid sprays in the absence of atomizing air.Their analysis included effects of primary breakup for liquid ellipsoidal droplets based on the Taylor Analogy Breakup (TAB) model developed by O'Rourke and Amsden [23], and later modified by Ibrahim [24], and Kvasnak and Ahmadi [25].Recent reviews of the spray simulation and breakup methods were presented by Jenny et al. [5].
This work was concerned with understanding the effects of droplet deformation, breakup, and evaporation on the dispersion of deforming ellipsoidal spray droplets in turbulent fuel injector flow fields with atomizing air, with a focus towards practical applications.The other goal was to develop a computational tool for the simulation of practical liquid spray fuel injectors.Kvasnak et al. [8] described a detailed numerical procedure for generating ensembles of simultaneous sample particle trajectories for estimating particle velocity and dispersion statistics.The presented study showed that, in many cases of practical importance to fuel injectors with atomizing air, the deformation time scale of the fuel droplets was comparable with their evaporation time scale.The droplet nonsphericity was also found to be an important factor for modeling the fuel injector performance.The simulation results were compared with the experimental data of McDonell and Samuelsen [22] for a simplex atomizer with atomizing air, where qualitative agreement was found.

Equations of Motion for Ellipsoidal Particles
In this section, the equations of motion for an ellipsoid of revolution suspended in a turbulent velocity field as described by Fan and Ahmadi [26] and Zhang et al. [27] are outlined.Figure 1 shows the schematic of an ellipsoidal particle moving in a general flow field and the corresponding inertial (x, y, z) and particle ( x, ŷ, ẑ) coordinate systems.The co-motion frame, ( x, ŷ, ẑ), with its origin being at the particle centroid and having axes parallel to the inertial frame, is also shown in this figure .Fluids 2018, 3, x FOR PEER REVIEW 2 of 18 solvers.Earlier, McDonell and Samuelsen [22] provided experimental data on droplet size and velocity distributions in a simplex atomizer.Kvasnak et al. [8] examined the formation of liquid sprays in the absence of atomizing air.Their analysis included effects of primary breakup for liquid ellipsoidal droplets based on the Taylor Analogy Breakup (TAB) model developed by O'Rourke and Amsden [23], and later modified by Ibrahim [24], and Kvasnak and Ahmadi [25].Recent reviews of the spray simulation and breakup methods were presented by Jenny et al. [5].
This work was concerned with understanding the effects of droplet deformation, breakup, and evaporation on the dispersion of deforming ellipsoidal spray droplets in turbulent fuel injector flow fields with atomizing air, with a focus towards practical applications.The other goal was to develop a computational tool for the simulation of practical liquid spray fuel injectors.Kvasnak et al. [8] described a detailed numerical procedure for generating ensembles of simultaneous sample particle trajectories for estimating particle velocity and dispersion statistics.The presented study showed that, in many cases of practical importance to fuel injectors with atomizing air, the deformation time scale of the fuel droplets was comparable with their evaporation time scale.The droplet nonsphericity was also found to be an important factor for modeling the fuel injector performance.The simulation results were compared with the experimental data of McDonell and Samuelsen [22] for a simplex atomizer with atomizing air, where qualitative agreement was found.

Equations of Motion for Ellipsoidal Particles
In this section, the equations of motion for an ellipsoid of revolution suspended in a turbulent velocity field as described by Fan and Ahmadi [26] and Zhang et al. [27] are outlined.Figure 1 shows the schematic of an ellipsoidal particle moving in a general flow field and the corresponding inertial (, , ) and particle ( , ,) coordinate systems.The co-motion frame, ( , ,̂ ), with its origin being at the particle centroid and having axes parallel to the inertial frame, is also shown in this figure.The equations of motion of a non-spherical rigid particle moving in an arbitrary flow field, in the presence of hydrodynamic and gravitational forces are given as Goldstein [28], and Gallily and Cohen [29].The equations of motion of a non-spherical rigid particle moving in an arbitrary flow field, in the presence of hydrodynamic and gravitational forces are given as Goldstein [28], and Gallily and Cohen [29].

Rotational Motion
Here, m p is the mass of the particle, t is the time, v is the translational velocity vector acting at the particle center of gravity, f h is the hydrodynamic drag force acting on the particle, f L is the lift force, and g is the acceleration of body force (i.e., gravity).In Equations ( 2)-( 4), I x, I ŷ, I ẑ are the moments of inertia of the particle with respect to the particle principal axes, ω x, ω ŷ, ω ẑ are the angular velocities of the particle with respect to the principal particle axes, and T h x , T h ŷ , T h ẑ are the components of the hydrodynamic torque acting on the particle.Equation (1) expresses the translational motion in the inertial frame, (x, y, z), while Equations ( 2)-( 4) for the rotational motion are stated in the particle frame, ( x, ŷ, ẑ).It should be emphasized that Equations ( 2)-( 4) are for rigid bodies.Here, it is assumed that these equations are approximately valid for droplets with high viscosity.
The particle position may be easily obtained by integrating where x = (x, y, z) is the position vector.Evaluation of the particle orientation, however, is more complicated.Recognizing that the rotational motion is written in particle coordinates, a change of frame is required.The transformation between the co-motion and the particle frame of reference was given by Goldstein [28] as where x = ( x, ŷ, ẑ) and x = x, ŷ, ẑ , and A is the transformation matrix.Fan and Ahmadi [26] showed that the Euler four parameters or quaternions may be used to avoid the singularity produced by the use of Euler's angles for full rotational motion of a rigid body.Here, the Euler four parameters, shown in Figure 2, are defined as where e is the unit vector along the axis of rotation, φ is the angle of rotation, and the superscript T denotes a matrix transpose.Because the most general rotation of a rigid body has only three degrees of freedom, Euler's four parameters are subject to a constraint given as The transformation matrix A in Equation ( 6) may be stated in terms of quaternions as Hughes [30], Fluids 2019, 4, 20

of 17
The time rates of change of quaternions are given by The initial conditions of the quaternions are obtained from the initial direction cosines, a ij , by the following relations: The non-uniqueness given by Equations ( 12) and ( 13) does not introduce any difficulty because all solutions correspond to the same transformation matrix and result in identical particle trajectories. . ( The non-uniqueness given by Equations ( 12) and ( 13) does not introduce any difficulty because all solutions correspond to the same transformation matrix and result in identical particle trajectories.

Ellipsoids of Revolution
Ellipsoids of revolution cover a wide variety of particle shapes (e.g., spheres, prolate spheroids, and oblate spheroids).In fact, all quadratic surfaces of finite extent are, by definition, ellipsoids.A particle entrained in a viscous fluid is subject to various hydrodynamic forces and torques such as drag, lift, Faxen, Bassett history, etc.In this work, only the leading order hydrodynamic forces of drag and lift are included in the equations of motion for small ellipsoidal particles.Explicit expressions for the hydrodynamic forces and torques for an ellipsoid of revolution are described in this section.Throughout this paper, the y-axis is assumed to be along the particle major axis as shown in Figure 1.

Nonlinear Hydrodynamic Drag
The general hydrodynamic forces and torques acting on an ellipsoidal particle in an arbitrary creeping flow field were described by Happel and Brenner [31].For larger particles moving outside of the creeping flow regime, certain Reynolds number corrections are needed.Here, for evaporating ellipsoidal droplets is used.In this equation, μ is the viscosity of the fluid, B is the Spalding transfer number for evaporating droplets, and Re a is the particle Reynolds number based

Ellipsoids of Revolution
Ellipsoids of revolution cover a wide variety of particle shapes (e.g., spheres, prolate spheroids, and oblate spheroids).In fact, all quadratic surfaces of finite extent are, by definition, ellipsoids.A particle entrained in a viscous fluid is subject to various hydrodynamic forces and torques such as drag, lift, Faxen, Bassett history, etc.In this work, only the leading order hydrodynamic forces of drag and lift are included in the equations of motion for small ellipsoidal particles.Explicit expressions for the hydrodynamic forces and torques for an ellipsoid of revolution are described in this section.Throughout this paper, the y-axis is assumed to be along the particle major axis as shown in Figure 1.

Nonlinear Hydrodynamic Drag
The general hydrodynamic forces and torques acting on an ellipsoidal particle in an arbitrary creeping flow field were described by Happel and Brenner [31].For larger particles moving outside of the creeping flow regime, certain Reynolds number corrections are needed.Here, for evaporating ellipsoidal droplets is used.In this equation, µ is the viscosity of the fluid, B is the Spalding transfer number for evaporating droplets, and Re a is the particle Reynolds number based on the minimum radius of the particle and the flow/particle velocity difference.This equation is a simple generalization of the nonlinear drag for a sphere to ellipsoids.The Spalding transfer number correction to the coefficient of drag is also a simple extension from the spherical case.The particle Reynolds number is defined as where u = (u x ,u y ,u z ) is the flow velocity vector at the centroid of the particle in the absence of the particle, and ν is the kinematic viscosity of the fluid.
In Equation ( 14), a double hat denotes a quantity expressed in the co-motion frame.The transformation of the translation matrix is given by The particle-frame translation matrix K for an ellipsoid of revolution is a diagonal matrix, i.e., with and where λ is the droplet to gas viscosity ratio, and β = b/a is the particle aspect ratio.
It should be noted that the correction factors for evaporation and nonlinear drag used in Equation ( 14) involving B and Re a , and the correction factors involving λ, introduced in Equations ( 18) and (19) to account for the internal circulation of the ellipsoidal droplets, are identical to those used for a spherical particles.In the absence of a rigorous analysis for nonlinear drag or liquid ellipsoids, the expressions given in Equations ( 14), ( 18) and ( 19) are expected to be reasonable approximations for high values of λ and/or slightly deformed spheres.The low aspect ratio asymptotes of Equations (18) and (19) reduce to the Stokes drag for a spherical droplet.The effective drag coefficient vector with components along the ellipsoid three principal directions are given as where C D is the coefficient of drag for a spherical particle.That is, Figure 3 shows the variation of the effective coefficient of drag in the x −, ẑ −, and ŷ-directions as a function of Reynolds number and aspect ratio.Figure 3a shows the effective drag coefficient in x− and ẑ-directions.Figure 3b shows the drag coefficient in the ŷ-direction.

Shear Induced Lift
The inertial lift force acting on an arbitrary body of symmetry in a shear field as obtained by Harper and Chang [32] is given as In the limit of spherical particle,  approaches 1 and the translational matrix becomes diagonal.

Shear Induced Lift
The inertial lift force acting on an arbitrary body of symmetry in a shear field as obtained by Harper and Chang [32] is given as Here, Re G = | Ĝ|a 2 ν is the Reynolds number based on the flow shear rate where | Ĝ| is the magnitude of the principal shear rates.The lift tensor L is given as In the limit of spherical particle, β approaches 1 and the translational matrix becomes diagonal.That is, K = K = 6I and Equation ( 22) reduced to Fluids 2019, 4, 20 which is in agreement with the Saffman expression for the lift force for a small spherical particle in a shear field.

Hydrodynamic Torque
In the absence of a rigorous theoretical model for the torque acting on a deformable ellipsoid, the components of hydrodynamic torque acting on a solid ellipsoid of revolution as obtained by Jeffery [33] are used as a first approximation.These are given by where are the elements of the deformation-rate and spin tensors.Here, The instantaneous carrier-fluid velocity gradients (evaluated at the ellipsoidal droplet centroid) in the particle frame needed in Equations ( 22), ( 28) and ( 29) may be obtained by using the transformation, where Ĝ and Ĝ stand for a dyadic expressed in the particle and the co-motion frames, respectively.Since the viscosity of the fuel droplets are more than several hundred times that of the atomizing gas, the use of Equations ( 25)-( 27) for ellipsoidal droplets may be, in part, justified.In an earlier study, Stone [34] reported that droplets (with high liquid-to-gas viscosity ratio) undergo roughly rigid body rotation in shear flows.

Equations of Droplet Motion
The governing equations used for the liquid droplets and the surrounding turbulent airflow field were described by Kvasnak et al. [8].Therefore, only a brief outline is given here.The equations of motion of an ellipsoid of revolution in an arbitrary fluid flow field may be written as Fluids 2019, 4, 20 8 of 17 where the transformation matrix K was defined by Kvasnak et al. [8].
In these equations, the m p is the mass of an ellipsoid of revolution, which is given as Here, ρ p denotes the density of the particle.The particle mass moments of inertia in the particle frame are given by For a given flow field, the equations of motion given by Equations ( 33)-( 38) may be numerically integrated for determining the droplet position and orientation as functions of time.

Langevin Model for the Instantaneous Velocity and Velocity Gradient
The instantaneous velocity and velocity gradient fields are governed by the joint velocity-velocity gradient PDF described by Haworth and Pope [35,36], and Girimaji and Pope [37].The corresponding associated Langevin equations are: and Here, u i is the instantaneous fluid velocity vector, h ij = ∂u i / ∂x j is the instantaneous velocity gradient, and W i and W ij are the vector and tensor valued zero mean incremental Weiner processes with a variance of dt.Additional details of the Langevin equation models for simulating the velocity-velocity gradient of the turbulent flow field were described by Kvasnak et al. [8].
In the present study, the mean flow and root mean square properties of the flow were generated by the steady state solution of the Reynolds Averaged Navier-Stokes (RANS) equations using the Reynolds stress transport turbulence model.The ANSYS-FLUENT code was used for this purpose.The results were then implemented into the Langevin equation models given by Equations ( 39) and (40) to generate the instantaneous flow velocity and velocity gradient fields.

Computational Scheme
The full spray simulation program is a modified version of the program used by Kvasnak et al. [8] with the addition of the atomizing air and evaporation model.It is a standalone post-processor code for a given turbulent velocity field.The code uses a Lagrangian approach and solves the equations of motion for solid or liquid ellipsoidal particles undergoing translation, rotation, deformation, breakup and evaporation.The program can also handle the dispersion of a cloud of randomly sized and shaped particles.A semi-implicit backward Euler method is used for the integration of the translational, rotational, and deformation equations for the ellipsoid.
At each time-step, the instantaneous fluid velocity, its first derivative, and the fluid temperature at the location of the particle are evaluated via a Newton's iteration scheme along with a trilinear interpolation in generalized coordinates in order to provide communication in their block structured mesh algorithm, allowing for the efficient evaluation of the fluid velocity and its derivatives at the centroid of the particle.First, the Newton iteration method is used to find the particle position in grid coordinates.Then, a trilinear interpolation is used to find the fluid velocity, temperature, and velocity gradients at the particle location.Additional details of the trilinear interpolation were discussed by Kvasnak and Ahmadi [25] and Kvasnak et al. [8].

Results
In this section, simulation results for a turbulent liquid spray assisted by non-swirling atomizing air are presented and discussed.

Simplex Atomizer
The fuel injector used in the present simulation is identical to the one used in the experimental study of McDonell and Samuelsen [22].In the experiment, a simplex injector tip was mounted centrally within a passage allowing airflow.The geometry of the injector is shown in Figure 4.The simplex tip has a diameter of 0.5 mm and was angled at 41 • to the primary direction of flow.In the experiment, the atomizing airflows through a cylindrical cone gap of varying thickness with the gap size of 2.5 mm at the orifice.The fuel spray flowed through a 4.9 mm diameter orifice before issuing into the ambient air.Methanol with a mass flow rate of 1.26 g/s was injected vertically downward as a cone annular sheet.The spray pressure drop was 420 kPa.The atomizing airflow rate was 1.32 g/s with a pressure drop of 3.73 kPa and 13.79 kPa for cases without and with swirling flows, respectively.The experiment was performed at room temperature.

Results
In this section, simulation results for a turbulent liquid spray assisted by non-swirling atomizing air are presented and discussed.

Simplex Atomizer
The fuel injector used in the present simulation is identical to the one used in the experimental study of McDonell and Samuelsen [22].In the experiment, a simplex injector tip was mounted centrally within a passage allowing airflow.The geometry of the injector is shown in Figure 4.The simplex tip has a diameter of 0.5 mm and was angled at 41° to the primary direction of flow.In the experiment, the atomizing airflows through a cylindrical cone gap of varying thickness with the gap size of 2.5 mm at the orifice.The fuel spray flowed through a 4.9 mm diameter orifice before issuing into the ambient air.Methanol with a mass flow rate of 1.26 g/s was injected vertically downward as a cone annular sheet.The spray pressure drop was 420 kPa.The atomizing airflow rate was 1.32 g/s with a pressure drop of 3.73 kPa and 13.79 kPa for cases without and with swirling flows, respectively.The experiment was performed at room temperature.In this study, the experimental conditions of McDonell and Samuelsen [22] with idealized geometry are used and the spray formation from the atomizer is simulated numerically.Liquid methanol is injected as a string of spheres with diameters of 100 μm.The spheres are oriented around the center of the nozzle tip.As they exit through the small orifice at the center of the injector, the particle radial velocity induces the formation of a conical annular sheet.Atomizing air is supplied through the ring surrounding the orifice as an inlet boundary condition.
The fuel injector is simulated using the following steps.The steady state airflow field is first generated using the Reynolds stress transport turbulence model of the ANSYS-FLUENT code.The aspect ratio of the droplet is set to 1.001 and the number of droplets injected is determined by the volume flow rate of liquid through the nozzle.After injection, the droplets are allowed to translate, rotate, deform, breakup, and evaporate in the turbulent gas flow field.Evolution of this process leads to an atomized spray of droplets.As the droplets leave the nozzle, they form a conical annular sheet similar to that observed in the experiment.One major feature of the present simulation is that In this study, the experimental conditions of McDonell and Samuelsen [22] with idealized geometry are used and the spray formation from the atomizer is simulated numerically.Liquid methanol is injected as a string of spheres with diameters of 100 µm.The spheres are oriented around the center of the nozzle tip.As they exit through the small orifice at the center of the injector, the particle radial velocity induces the formation of a conical annular sheet.Atomizing air is supplied through the ring surrounding the orifice as an inlet boundary condition.
The fuel injector is simulated using the following steps.The steady state airflow field is first generated using the Reynolds stress transport turbulence model of the ANSYS-FLUENT code.The aspect ratio of the droplet is set to 1.001 and the number of droplets injected is determined by the volume flow rate of liquid through the nozzle.After injection, the droplets are allowed to translate, Fluids 2019, 4, 20 10 of 17 rotate, deform, breakup, and evaporate in the turbulent gas flow field.Evolution of this process leads to an atomized spray of droplets.As the droplets leave the nozzle, they form a conical annular sheet similar to that observed in the experiment.One major feature of the present simulation is that all of the particles formed by the spray are tracked for their lifetime (i.e., either they evaporate or leave the computational space through an outer boundary surface).
Figure 5 shows a close up of the grid used in the simulation of the simplex atomizer near the injector orifice.The computational space is axisymmetric with a radius of 50 mm and a length of 200 mm.Body-fitted coordinates are used in the computation so that the detailed feature of the atomizer nozzle can be accounted for.A grid sensitivity study was performed using different size meshes.Finally, a grid of 201 by 51 cells that is uniformly expanding in both the streamwise and radial directions is selected and used for the present simulations.A computational grid near the tip of the nozzle is relatively dense to improve the resolution in this critical region.Nevertheless, since the nozzle is very small, only a couple of cells cover the interior of the nozzle.The boundary conditions in the radial direction are symmetry at r = 0 and continuous outflow at r = 50 mm.The boundary conditions in the streamwise direction is no slip wall at x = 0, except for the nozzle inlet which is at a constant velocity and continuous outflow at x = 200 mm.The inlet boundary conditions are uniform velocity equal to that specified in the experiment.For the injector tip here, the inflow boundary condition is 8.35 m/s.all of the particles formed by the spray are tracked for their lifetime (i.e., either they evaporate or leave the computational space through an outer boundary surface).
Figure 5 shows a close up of the grid used in the simulation of the simplex atomizer near the injector orifice.The computational space is axisymmetric with a radius of 50 mm and a length of 200 mm.Body-fitted coordinates are used in the computation so that the detailed feature of the atomizer nozzle can be accounted for.A grid sensitivity study was performed using different size meshes.Finally, a grid of 201 by 51 cells that is uniformly expanding in both the streamwise and radial directions is selected and used for the present simulations.A computational grid near the tip of the nozzle is relatively dense to improve the resolution in this critical region.Nevertheless, since the nozzle is very small, only a couple of cells cover the interior of the nozzle.The boundary conditions in the radial direction are symmetry at r = 0 and continuous outflow at r = 50 mm.The boundary conditions in the streamwise direction is no slip wall at x = 0, except for the nozzle inlet which is at a constant velocity and continuous outflow at x = 200 mm.The inlet boundary conditions are uniform velocity equal to that specified in the experiment.For the injector tip here, the inflow boundary condition is 8.35 m/s.

Non-Swirling Atomizing Air
This section describes the steady state flow field of the injector with non-swirling atomizing air as evaluated by the ANSYS-FLUENT code assuming an inlet velocity of 26.8 m/s.
Figure 6 shows the streamwise velocity contours downstream of the nozzle.This figure clearly shows the development of the shear layer at the edge of the atomizing air jet and the suppression of the shear layers generated by the nozzle flow.As the flow develops the velocity decays in amplitude and spreads radially.A small recirculation region created at the tip of the atomizing air orifice can also be seen in this figure.

Non-Swirling Atomizing Air
This section describes the steady state flow field of the injector with non-swirling atomizing air as evaluated by the ANSYS-FLUENT code assuming an inlet velocity of 26.8 m/s.
Figure 6 shows the streamwise velocity contours downstream of the nozzle.This figure clearly shows the development of the shear layer at the edge of the atomizing air jet and the suppression of the shear layers generated by the nozzle flow.As the flow develops the velocity decays in amplitude and spreads radially.A small recirculation region created at the tip of the atomizing air orifice can also be seen in this figure .Computed contours of turbulence kinetic energy that are shown in Figure 7 indicate the presence of two shear layers which are associated with the nozzle jet or atomizing air.The nozzle jet shear layer diminishes after short distances, while the atomizing air shear layer expands across the region.The peak in the turbulence kinetic energy occurs near the center of the shear layer; however, the kinetic energy is stronger at the point where the annular jet reaches the axis of the injector.
as evaluated by the ANSYS-FLUENT code assuming an inlet velocity of 26.8 m/s.
Figure 6 shows the streamwise velocity contours downstream of the nozzle.This figure clearly shows the development of the shear layer at the edge of the atomizing air jet and the suppression of the shear layers generated by the nozzle flow.As the flow develops the velocity decays in amplitude and spreads radially.A small recirculation region created at the tip of the atomizing air orifice can also be seen in this figure.Computed contours of turbulence kinetic energy that are shown in Figure 7 indicate the presence of two shear layers which are associated with the nozzle jet or atomizing air.The nozzle jet shear layer diminishes after short distances, while the atomizing air shear layer expands across the region.The peak in the turbulence kinetic energy occurs near the center of the shear layer; however, the kinetic energy is stronger at the point where the annular jet reaches the axis of the injector.Figures 8 and 9 show the computed contours of components of the Reynolds stress tensor  and   .From Figure 8, it is seen that the peak value of the normal stress  occurs near the middle of the shear layer.Figure 9 shows the turbulence shear stress contours, and indicates that the shear stress   has large values in the downstream vicinity of the nozzle and reaches its peak near the center of the shear layer.These fluctuation flow statistics are used in the Langevin equation model for generating the instantaneous fluid velocity at the location of droplets.The spray formation is described in the subsequent section.Figures 8 and 9 show the computed contours of components of the Reynolds stress tensor u 2 y and u x u y .From Figure 8, it is seen that the peak value of the normal stress u 2 y occurs near the middle of the shear layer.Figure 9 shows the turbulence shear stress contours, and indicates that the shear stress u x u y has large values in the downstream vicinity of the nozzle and reaches its peak near the center of the shear layer.These fluctuation flow statistics are used in the Langevin equation model for generating the instantaneous fluid velocity at the location of droplets.The spray formation is described in the subsequent section.Computed contours of turbulence kinetic energy that are shown in Figure 7 indicate the presence of two shear layers which are associated with the nozzle jet or atomizing air.The nozzle jet shear layer diminishes after short distances, while the atomizing air shear layer expands across the region.The peak in the turbulence kinetic energy occurs near the center of the shear layer; however, the kinetic energy is stronger at the point where the annular jet reaches the axis of the injector.Figures 8 and 9 show the computed contours of components of the Reynolds stress tensor  and   .From Figure 8, it is seen that the peak value of the normal stress  occurs near the middle of the shear layer.Figure 9 shows the turbulence shear stress contours, and indicates that the shear stress   has large values in the downstream vicinity of the nozzle and reaches its peak near the center of the shear layer.These fluctuation flow statistics are used in the Langevin equation model for generating the instantaneous fluid velocity at the location of droplets.The spray formation is described in the subsequent section.

Spray Simulation
The qualitative comparison of the predicted spray angle in the presence of atomizing air with the experimental data reported by McDonell and Samuelsen [22] is presented in Figure 10.The left-hand side of the figure corresponds to the simulated spray while a picture of the experimental spray is shown on the right side.In this figure, the location and sizes of all droplets in the flow 0.02 s after the spray is turned on are shown.As was noted before, the computational domain is 50 mm wide and 200 mm long.The spray figure is truncated at an axial distance of 100 mm, which is far from the downstream boundary condition.Each initial droplet breaks on the average into four daughter droplets in this flow field through the primary breakup mechanisms.A total of 25,000 droplets are formed in this spray during the period of 0.02 s.Qualitative agreement is seen between the computational spray and the experimental spray with a difference in spray angles of a few percent.In Figure 11, the radial distance of droplet position downstream of the nozzle exit at axial locations of 25, 50, and 75 mm are shown.At 25 mm downstream, the spray is an annular cone with a thickness of approximately 2 mm, and the radial distribution function is symmetric with a very small standard deviation.Further downstream at 50 mm, the spray is still an annular cone and has a thickness of approximately 5 mm with no droplets in the core region.The distribution function is also roughly symmetric across the spray region.Further downstream at 75 mm, the spray is still an

Spray Simulation
The qualitative comparison of the predicted spray angle in the presence of atomizing air with the experimental data reported by McDonell and Samuelsen [22] is presented in Figure 10.The left-hand side of the figure corresponds to the simulated spray while a picture of the experimental spray is shown on the right side.In this figure, the location and sizes of all droplets in the flow 0.02 s after the spray is turned on are shown.As was noted before, the computational domain is 50 mm wide and 200 mm long.The spray figure is truncated at an axial distance of 100 mm, which is far from the downstream boundary condition.Each initial droplet breaks on the average into four daughter droplets in this flow field through the primary breakup mechanisms.A total of 25,000 droplets are formed in this spray during the period of 0.02 s.Qualitative agreement is seen between the computational spray and the experimental spray with a difference in spray angles of a few percent.

Spray Simulation
The qualitative comparison of the predicted spray angle in the presence of atomizing air with the experimental data reported by McDonell and Samuelsen [22] is presented in Figure 10.The left-hand side of the figure corresponds to the simulated spray while a picture of the experimental spray is shown on the right side.In this figure, the location and sizes of all droplets in the flow 0.02 s after the spray is turned on are shown.As was noted before, the computational domain is 50 mm wide and 200 mm long.The spray figure is truncated at an axial distance of 100 mm, which is far from the downstream boundary condition.Each initial droplet breaks on the average into four daughter droplets in this flow field through the primary breakup mechanisms.A total of 25,000 droplets are formed in this spray during the period of 0.02 s.Qualitative agreement is seen between the computational spray and the experimental spray with a difference in spray angles of a few percent.In Figure 11, the radial distance of droplet position downstream of the nozzle exit at axial locations of 25, 50, and 75 mm are shown.At 25 mm downstream, the spray is an annular cone with a thickness of approximately 2 mm, and the radial distribution function is symmetric with a very small standard deviation.Further downstream at 50 mm, the spray is still an annular cone and has a thickness of approximately 5 mm with no droplets in the core region.The distribution function is also roughly symmetric across the spray region.Further downstream at 75 mm, the spray is still an In Figure 11, the radial distance of droplet position downstream of the nozzle exit at axial locations of 25, 50, and 75 mm are shown.At 25 mm downstream, the spray is an annular cone with a thickness of approximately 2 mm, and the radial distribution function is symmetric with a very small standard deviation.Further downstream at 50 mm, the spray is still an annular cone and has a thickness of approximately 5 mm with no droplets in the core region.The distribution function is also roughly symmetric across the spray region.Further downstream at 75 mm, the spray is still an annular cone with a thickness of approximately 15 mm.The distribution function is skewed with its peak moving toward the outer edge of the spray cone and a longer tail toward the inner cone.Figure 12 shows the variation of the droplet minimum diameter at axial distances of 25, 50, and 75 mm downstream of the nozzle exit.At a downstream position of 25 mm, the droplet minimum diameter is in the range of 30 μm to 40 μm, with the central part of the spray cone containing slightly larger sizes compared to the edges.At 50 mm downstream, the droplet minimum diameter is scattered in range of 40 μm to 60 μm with the mean of about 50 μm.While at 75 mm downstream, the droplet minimum diameter varies more narrowly between 70 μm to 80 μm with a mean of about 75 μm.Figure 13 shows the predicted droplet aspect ratio at downstream axial positions of 25, 50, and 75 mm.It is observed that the droplets at the distance of 25mm are highly elongated having aspect ratios observed between 11 and 13 with the highest values appearing on the edges of the spray.At a downstream position of 50 mm, the droplets are still quite elongated with aspect ratios ranging from 6 to 11.Further downstream at axial distance of 75 mm from the injector, the droplets become nearly spherical with some having aspect ratios about 1.5 to 2. It can be seen that the droplets deform significantly between the axial distances of 25 and 50 mm and particularly from 50 to 75 mm.The reason is that the velocity gradient decrease with distance from the spray nozzle and becomes quite small at the axial distance of about 75 mm.In the absence of shear induced deformation, the droplets tend to become spherical under the action of surface tension.Figure 12 shows the variation of the droplet minimum diameter at axial distances of 25, 50, and 75 mm downstream of the nozzle exit.At a downstream position of 25 mm, the droplet minimum diameter is in the range of 30 µm to 40 µm, with the central part of the spray cone containing slightly larger sizes compared to the edges.At 50 mm downstream, the droplet minimum diameter is scattered in range of 40 µm to 60 µm with the mean of about 50 µm.While at 75 mm downstream, the droplet minimum diameter varies more narrowly between 70 µm to 80 µm with a mean of about 75 µm.
annular cone with a thickness of approximately 15 mm.The distribution function is skewed with its peak moving toward the outer edge of the spray cone and a longer tail toward the inner cone.Figure 12 shows the variation of the droplet minimum diameter at axial distances of 25, 50, and 75 mm downstream of the nozzle exit.At a downstream position of 25 mm, the droplet minimum diameter is in the range of 30 μm to 40 μm, with the central part of the spray cone containing slightly larger sizes compared to the edges.At 50 mm downstream, the droplet minimum diameter is scattered in range of 40 μm to 60 μm with the mean of about 50 μm.While at 75 mm downstream, the droplet minimum diameter varies more narrowly between 70 μm to 80 μm with a mean of about 75 μm.Figure 13 shows the predicted droplet aspect ratio at downstream axial positions of 25, 50, and 75 mm.It is observed that the droplets at the distance of 25mm are highly elongated having aspect ratios observed between 11 and 13 with the highest values appearing on the edges of the spray.At a downstream position of 50 mm, the droplets are still quite elongated with aspect ratios ranging from 6 to 11.Further downstream at axial distance of 75 mm from the injector, the droplets become nearly spherical with some having aspect ratios about 1.5 to 2. It can be seen that the droplets deform significantly between the axial distances of 25 and 50 mm and particularly from 50 to 75 mm.The reason is that the velocity gradient decrease with distance from the spray nozzle and becomes quite small at the axial distance of about 75 mm.In the absence of shear induced deformation, the droplets tend to become spherical under the action of surface tension.Figure 13 shows the predicted droplet aspect ratio at downstream axial positions of 25, 50, and 75 mm.It is observed that the droplets at the distance of 25mm are highly elongated having aspect ratios observed between 11 and 13 with the highest values appearing on the edges of the spray.At a downstream position of 50 mm, the droplets are still quite elongated with aspect ratios ranging from 6 to 11.Further downstream at axial distance of 75 mm from the injector, the droplets become nearly spherical with some having aspect ratios about 1.5 to 2. It can be seen that the droplets deform significantly between the axial distances of 25 and 50 mm and particularly from 50 to 75 mm.The reason is that the velocity gradient decrease with distance from the spray nozzle and becomes quite small at the axial distance of about 75 mm.In the absence of shear induced deformation, the droplets tend to become spherical under the action of surface tension.The experimental data show a continuous decrease of diameter toward the centerline, and the presence of small droplets in the core region.The simulation, however, predicts no droplets in the core region.It is conjectured that the smaller droplet size in the experiment and the presence of very small droplets in the core is due to the secondary breakup mechanisms.In the present study, the secondary break was not included in the computational model as the bulk of the mass of the spray is accounted for by the primary breakup mechanism.The absence of the secondary breakup in the model is perhaps the main cause for the discrepancy seen in Figure 14 for smaller particles.In addition, it should be pointed out that the present model prediction suggests that the droplets are highly elongated at short distances from the nozzle and this could introduce error in the experimental Phase Doppler Anemometry (PDA) measurements.Further downstream, the variations of the Sauter diameter are similar to that shown for 25 mm downstream except for the spreading of the spray to larger radial segment.As noted before, the simulated D 32 overestimates the experimental data by approximately 15-20%.In addition, the   [22].Here, for the model predictions, the Sauter mean diameter of the equivalent volume spherical droplet is used for non-spherical droplets.The simulated droplets have a roughly constant Sauter diameter of about 75 µm to 80 µm, while the experimental data show a Sauter diameter of about 50 µm to 65 µm at the axial distance of 25 mm, and a diameter in the range of 50 µm to 70 µm at a distance of 50 and 75 mm.The experimental data show a continuous decrease of diameter toward the centerline, and the presence of small droplets in the core region.The simulation, however, predicts no droplets in the core region.It is conjectured that the smaller droplet size in the experiment and the presence of very small droplets in the core is due to the secondary breakup mechanisms.In the present study, the secondary break was not included in the computational model as the bulk of the mass of the spray is accounted for by the primary breakup mechanism.The absence of the secondary breakup in the model is perhaps the main cause for the discrepancy seen in Figure 14 for smaller particles.In addition, it should be pointed out that the present model prediction suggests that the droplets are highly elongated at short distances from the nozzle and this could introduce error in the experimental Phase Doppler Anemometry (PDA) measurements.Figure 14 compares the predicted Sauter mean diameters at distances of 25, 50, and 75 mm downstream from the injector with the experimental data of McDonell and Samuelsen [22].Here, for the model predictions, the Sauter mean diameter of the equivalent volume spherical droplet is used for non-spherical droplets.The simulated droplets have a roughly constant Sauter diameter of about 75 μm to 80 μm, while the experimental data show a Sauter diameter of about 50 μm to 65 μm at the axial distance of 25 mm, and a diameter in the range of 50 μm to 70 μm at a distance of 50 and 75 mm.The experimental data show a continuous decrease of diameter toward the centerline, and the presence of small droplets in the core region.The simulation, however, predicts no droplets in the core region.It is conjectured that the smaller droplet size in the experiment and the presence of very small droplets in the core is due to the secondary breakup mechanisms.In the present study, the secondary break was not included in the computational model as the bulk of the mass of the spray is accounted for by the primary breakup mechanism.The absence of the secondary breakup in the model is perhaps the main cause for the discrepancy seen in Figure 14 for smaller particles.In addition, it should be pointed out that the present model prediction suggests that the droplets are highly elongated at short distances from the nozzle and this could introduce error in the experimental Phase Doppler Anemometry (PDA) measurements.Further downstream, the variations of the Sauter diameter are similar to that shown for 25 mm downstream except for the spreading of the spray to larger radial segment.As noted before, the simulated D 32 overestimates the experimental data by approximately 15-20%.In addition, the Further downstream, the variations of the Sauter diameter are similar to that shown for 25 mm downstream except for the spreading of the spray to larger radial segment.As noted before, the simulated D 32 overestimates the experimental data by approximately 15-20%.In addition, the simulation also did not capture the smaller droplets in the central region of the annular cone that are generated by the secondary breakup.
Figure 15 compares the predicted streamwise droplet velocity, u x , with the experimental data of McDonell and Samuelsen [22] at 25, 50, and 75 mm downstream of the nozzle.The flow velocity profiles as predicted by the ANSYS-FLUENT code are shown by solid lines in the figure for comparison.The predicted and experimental droplet velocities are larger than the flow velocities.The experimental data of McDonell and Samuelsen [22], however, show a sharp increase in velocity toward the centerline.As noted before, in this region, the droplets are smaller and are generated by the secondary breakup mechanism.The predicted particle velocities are shown only for the range of spray layer as there are no particles in the core region.As was mentioned before, the present computational model for spray formation lacks the secondary breakup mechanisms and, therefore, the formation small particles are missed.
Figure 15 also shows that the simulated droplet velocity across the spray in the region 10 ≤ r ≤ 15 mm at axial distance of 25 mm is in reasonable agreement with the experimental data.Similarly, the predicted velocity at distance of 50 and 75 mm are comparable with the experimental data for the range that the simulations results are available.The experimental data of McDonell and Samuelsen [22], however, show a sharp increase in velocity toward the centerline.As noted before, in this region, the droplets are smaller and are generated by the secondary breakup mechanism.The predicted particle velocities are shown only for the range of spray layer as there are no particles in the core region.As was mentioned before, the present computational model for spray formation lacks the secondary breakup mechanisms and, therefore, the formation small particles are missed.Figure 15 also shows that the simulated droplet velocity across the spray in the region 10  15 mm at axial distance of 25 mm is in reasonable agreement with the experimental data.
Similarly, the predicted velocity at distance of 50 and 75 mm are comparable with the experimental data for the range that the simulations results are available.

Conclusions
In this work, the formation of a spray in a fuel injector with atomizing airflow field is studied.The equations of motion for the translation and rotation of an ellipsoidal droplet are presented.The evaporation rate of a non-spherical droplet is included.Physical models for the deformation and breakup of liquid droplets in the presence of shear flow and pressure gradient fields are described.A Langevin equation for generating instantaneous air velocities that satisfy the Navier-Stokes approximately (via a closed PDF formulation) is used for evaluating the droplet trajectories.In addition, a novel application of the Newton iteration and trilinear interpolation in generalized coordinates is used as the basis for the spray simulation algorithm, in which every droplet is traced through the flow field.The general conclusions of the study are: • The presented engineering model for simulation of spray leads to results that are in qualitative agreement with the experimental data for the studied cases.
• The average droplet size due to primary breakup mechanisms is in good agreement with the experimental data in regions of the flow where the droplets are nearly spherical.
• In regions of the spray where the droplets are elongated, the predicted droplet diameter overestimates the experimental data.
• The droplet velocity is in reasonable agreement with the experimental data.

Conclusions
In this work, the formation of a spray in a fuel injector with atomizing airflow field is studied.The equations of motion for the translation and rotation of an ellipsoidal droplet are presented.The evaporation rate of a non-spherical droplet is included.Physical models for the deformation and breakup of liquid droplets in the presence of shear flow and pressure gradient fields are described.A Langevin equation for generating instantaneous air velocities that satisfy the Navier-Stokes approximately (via a closed PDF formulation) is used for evaluating the droplet trajectories.In addition, a novel application of the Newton iteration and trilinear interpolation in generalized coordinates is used as the basis for the spray simulation algorithm, in which every droplet is traced through the flow field.The general conclusions of the study are:

•
The presented engineering model for simulation of spray leads to results that are in qualitative agreement with the experimental data for the studied cases.

•
The average droplet size due to primary breakup mechanisms is in good agreement with the experimental data in regions of the flow where the droplets are nearly spherical.

•
In regions of the spray where the droplets are elongated, the predicted droplet diameter overestimates the experimental data.

•
The droplet velocity is in reasonable agreement with the experimental data.

•
Droplet aspect ratio affects the evaporation rate of droplets.

•
Evaporation affects the fluctuation velocity of droplets, but it has little effect on the droplet diffusivity.

•
The presented model and simulation algorithm provides a viable tool for the computational prediction of spray simulations.

•
Secondary breakup mechanisms may play a major role in the atomization of sprays with complex atomizing airflow patterns and must be included in future studies.

( 22 )
Here, Re = | | is the Reynolds number based on the flow shear rate where | | is the magnitude of the principal shear rates.The lift tensor L is given as

Figure 5 .
Figure 5. Close-up view of the computational grid of the Simplex atomizer.

Figure 5 .
Figure 5. Close-up view of the computational grid of the Simplex atomizer.

Figure 7 .
Figure 7. Turbulent kinetic energy of the flow with non-swirling atomizing air.

' 2 yu
of the flow with non-swirling atomizing air.

Figure 7 .
Figure 7. Turbulent kinetic energy of the flow with non-swirling atomizing air.

Figure 7 .
Figure 7. Turbulent kinetic energy of the flow with non-swirling atomizing air.

' 2 yu
of the flow with non-swirling atomizing air.

Figure 8 .
Figure 8. u 2 y of the flow with non-swirling atomizing air.

Figure 9
Figure 9. ' ' x y u u of the flow with non-swirling atomizing air.

Figure 10 .
Figure 10.Qualitative comparison of spray angle with non-swirling atomizing air.

Figure 9 .
Figure 9. u x u y of the flow with non-swirling atomizing air.

Figure 10 .
Figure 10.Qualitative comparison of spray angle with non-swirling atomizing air.

Figure 10 .
Figure 10.Qualitative comparison of spray angle with non-swirling atomizing air.

FluidsFluids
2018, 3, x FOR PEER REVIEW 14 of 18 annular cone with a thickness of approximately 15 mm.The distribution function is skewed with its peak moving toward the outer edge of the spray cone and a longer tail toward the inner cone.

Figure 14
Figure 14 compares the predicted Sauter mean diameters at distances of 25, 50, and 75 mm downstream from the injector with the experimental data of McDonell and Samuelsen [22].Here, for the model predictions, the Sauter mean diameter of the equivalent volume spherical droplet is used for non-spherical droplets.The simulated droplets have a roughly constant Sauter diameter of about 75 μm to 80 μm, while the experimental data show a Sauter diameter of about 50 μm to 65 μm at the axial distance of 25 mm, and a diameter in the range of 50 μm to 70 μm at a distance of 50 and 75 mm.The experimental data show a continuous decrease of diameter toward the centerline, and the presence of small droplets in the core region.The simulation, however, predicts no droplets in the core region.It is conjectured that the smaller droplet size in the experiment and the presence of very small droplets in the core is due to the secondary breakup mechanisms.In the present study, the secondary break was not included in the computational model as the bulk of the mass of the spray is accounted for by the primary breakup mechanism.The absence of the secondary breakup in the model is perhaps the main cause for the discrepancy seen in Figure14for smaller particles.In addition, it should be pointed out that the present model prediction suggests that the droplets are highly elongated at short distances from the nozzle and this could introduce error in the experimental Phase Doppler Anemometry (PDA) measurements.

Figure 14 .
Figure 14.Sauter mean diameter of droplets at axial distances of 25, 50, and 75 mm with non-swirling atomizing air.

Figure 14
Figure14compares the predicted Sauter mean diameters at distances of 25, 50, and 75 mm downstream from the injector with the experimental data of McDonell and Samuelsen[22].Here, for the model predictions, the Sauter mean diameter of the equivalent volume spherical droplet is used for non-spherical droplets.The simulated droplets have a roughly constant Sauter diameter of about 75 µm to 80 µm, while the experimental data show a Sauter diameter of about 50 µm to 65 µm at the axial distance of 25 mm, and a diameter in the range of 50 µm to 70 µm at a distance of 50 and 75 mm.The experimental data show a continuous decrease of diameter toward the centerline, and the presence of small droplets in the core region.The simulation, however, predicts no droplets in the core region.It is conjectured that the smaller droplet size in the experiment and the presence of very small droplets in the core is due to the secondary breakup mechanisms.In the present study, the secondary break was not included in the computational model as the bulk of the mass of the spray is accounted for by the primary breakup mechanism.The absence of the secondary breakup in the model is perhaps the main cause for the discrepancy seen in Figure14for smaller particles.In addition, it should be pointed out that the present model prediction suggests that the droplets are highly elongated at short distances from the nozzle and this could introduce error in the experimental Phase Doppler Anemometry (PDA) measurements.

Figure 14 .
Figure 14.Sauter mean diameter of droplets at axial distances of 25, 50, and 75 mm with non-swirling atomizing air.

Figure 14 .
Figure 14.Sauter mean diameter of droplets at axial distances of 25, 50, and 75 mm with non-swirling atomizing air.

Fluids 2018, 3 ,
x FOR PEER REVIEW 16 of 18 simulation also did not capture the smaller droplets in the central region of the annular cone that are generated by the secondary breakup.Figure 15 compares the predicted streamwise droplet velocity, u x, with the experimental data of McDonell and Samuelsen [22] at 25, 50, and 75 mm downstream of the nozzle.The flow velocity profiles as predicted by the ANSYS-FLUENT code are shown by solid lines in the figure for comparison.The predicted and experimental droplet velocities are larger than the flow velocities.

Figure 15 .
Figure 15.Axial velocity of droplets at axial distances of 25, 50, and 75 mm with non-swirling atomizing air.

Figure 15 .
Figure 15.Axial velocity of droplets at axial distances of 25, 50, and 75 mm with non-swirling atomizing air.