## Abstract

Predicting the spray evolution using simulations requires accurate modeling of the turbulent gas-phase flow field. In this study, the high-order spectral-element method (SEM), implemented in the code Nek5000, was used to provide highly resolved solutions to the turbulent flow field. Spray modeling capabilities were implemented into the Nek5000 code. The spray is modeled in a Lagrangian–Eulerian (LE) framework, where the liquid is represented by discrete parcels of droplets. The method for coupling liquid and gas in the context of SEM is described, which allows for very fine meshes to be used without affecting the stability of the solution. Large-eddy simulations (LES) of the eight-hole ECN Spray G gasoline injector were conducted. Numerical results are compared against experimental data for liquid penetration, droplet size and gas velocity. The morphology of the multiplume spray is compared against experimental data. The effect of different spray injection inputs is analyzed. It was found that using a plume direction of 33 deg and an injection cone angle of 30 deg produced the best results overall. This work shows the applicability of SEM for spray modeling applications, where use of a high-order flow solver can help us understand the multiplume spray aerodynamics and how it leads to plume collapse under certain conditions. Results also highlight the need for tuning spray input parameters in the LE framework, even when high-fidelity gas flow solutions are possible.

## 1 Introduction

Gasoline direct injection (GDI) systems have become a growing trend in modern gasoline engines. Improved efficiency and cold-start performance has been obtained with direct injection [1]. However, GDI engines tend to have higher particulate mass and particulate number emissions compared with port fuel injection (PFI) engines [2]. For this reason, it is important to understand and control the fuel-air mixing processes in GDI engine designs.

Computational fluid dynamics (CFD) may be used to accelerate the discovery of knowledge and exploration of new injection strategies. A typical approach to model fuel injection sprays in engines is the Lagrangian–Eulerian (LE) approach. The liquid droplets are modeled as Lagrangian particles, while the gas phase is solved in an Eulerian reference frame, typically solved by second-order finite-volume methods.

Most GDI injectors are multihole, where interplume interactions and aerodynamics play an important role in the spray evolution, and may lead to plume collapse under certain conditions, causing the merged plumes to be ill-directed and uncontrollable [3–5]. Therefore, interplume aerodynamics must be well captured by simulations of GDI sprays [6], which requires numerically accurate solution of the gaseous flow field. To achieve this, high-order methods may be used to solve the turbulent flow field with minimal numerical errors [7].

This work introduces the use of the high-order spectral element method (SEM) [8] for the solution of the gas phase in GDI sprays using the LE approach (SEM-LE). The high-order solution in SEM achieves lower numerical error (higher accuracy) at similar levels of refinement than second-order methods (finite volume, or else). The SEM-LE approach was used previously by our group to perform high-fidelity large eddy simulations (LES) of turbulent particle-laden jets [9] and nonevaporative diesel sprays [10]. The objective of this work is to perform high-fidelity LES of GDI sprays to gain insight into the flow aerodynamics, where the high-resolution solution reduces uncertainty due to numerical accuracy.

This paper begins by presenting the governing equations used to model GDI sprays in an LE framework. Then, an overview of the numerical method is presented, followed by a description of the computational setup. Numerical results are then compared against experiments and conclusions are drawn based on macro-and local quantities such as liquid penetration, droplet sizes and gas velocity.

## 2 Governing Equations

In this section, we present the governing equations used to model the gas-phase flow field in an Eulerian reference frame. The Lagrangian droplet dynamics equations are also presented.

### 2.1 Gas-Phase Transport Equations.

where *p*_{1} is the hydrodynamic pressure, *p*_{0} is the thermodynamic pressure, and $\epsilon =\gamma Ma2$, where *γ* is the ratio of specific heat capacities and Ma is the Mach number. The thermodynamic pressure, *p*_{0}, is only a function of time. In this work, it was assumed constant, i.e., $dp0/dt=0$.

*N*-component gaseous mixtures are

_{g}Here, $q\u02d9$, *h _{i}*, $m\u02d9i$, and

*Y*are the heat transfer rate due to conduction/convection, enthalpy, mass flowrate due to phase change, and mass fraction of species

_{i}*i*, respectively. The term $\Delta hi=hi(Td)\u2212hi(T)$ is the difference in enthalpy calculated at droplet temperature,

*T*, and at local gas temperature,

_{d}*T*. This represents the heat loss due to the mass transfer from the evaporating liquid. The gas density is $\rho =\rho g\theta g$, where

*ρ*mixture density, $\theta g=1\u2212\theta l$ is the void fraction and

_{g}*θ*is the liquid volume fraction. Terms including void fraction,

_{l}*θ*, in Eqs. (2)–(5) account for volume displacement of the gas-phase due to the presence of liquid droplets. The source term, $fd$, represents the momentum transferred from the liquid phase to the gas [15]. All other body forces (gravity, electromagnetic, etc.) acting on the gas phase are neglected in this study.

_{g}*D*is the diffusion coefficient. Gas mixture heat capacity is computed as $cp=\u2211i=1Ngcp,iYi$, where $cp,i$ is the heat capacity of gaseous species

*i*. Density of the mixture is given by the ideal gas law

where $W=(\u2211i=1NgYi/Wi)\u22121$ is the molecular weight of the mixture and *W _{i}* is the molecular weight of each species, and

*R*is the universal gas constant. To calculate the mixture's thermal conductivity (

*k*) and viscosity (

*μ*), Wilke's mixing rule is used [16,17].

### 2.2 Liquid Droplet Equations.

*T*and

*T*are the gas and droplet temperatures, respectively;

_{d}*h*is the latent heat at

_{L}*T*,

_{d}*d*is the droplet diameter, and

_{d}*ρ*is liquid density. The Spalding mass transfer number is defined as

_{l}*Y*

_{Fs}is the fuel vapor mass fraction at the surface of the droplet,

*Y*is the fuel vapor mass fraction in the gaseous mixture,

_{g}*W*is the fuel molecular weight,

_{F}*W*is the molecular weight of the ambient gas, and

_{g}*p*is the vapor pressure of the fuel at droplet temperature. The Sherwood (Sh

_{v}*) and Nusselt (Nu*

_{d}*) numbers are calculated using the Abramzon and Sirignano evaporation model [19]. Gas properties at the film surrounding the droplet are used to compute the nondimensional numbers. Film properties are calculated according to the 1/3-rule [20]*

_{d}where the subscript *v* indicates quantity taken at film, *s* indicates the surface of the droplet, and *g* is the freestream gas mixture.

The forces $mdg,\u2009FD,d,\u2009FC,d$ and $FB,d$ are gravity, drag, interparticle and buoyancy forces, respectively. The mass of the droplet is $md=\rho l\pi dd3/6$. For this study, interparticle interactions and gravity forces were neglected: $FC,d=FB,d=g=0$.

*C*

_{D}*A*is the projected frontal area of the droplet, which is assumed to be spherical. The drag coefficient of a sphere is given by [21]

_{f}*y*, is solved for each droplet based on the Taylor analogy breakup model [22]. The droplet distortion parameter is used to calculate the dynamic drag coefficient [21], which varies linearly between the coefficient of a sphere and that of a disk, and is given by

_{d}### 2.3 Liquid Injection and Droplet Breakup.

*d*

_{0}, are introduced at the injection location. Droplet diameter is corrected to account for area contraction coefficient,

*C*, of the liquid jet passing through the nozzle hole

_{a}Rate of injection is specified from experimental measurements, and is controlled in the simulation by managing the blob injection frequency and the number of droplets per parcel. At each time-step, the mass from ROI is accumulated until it reaches $md,0=\rho lnd\pi dd,03/6$, at which point a blob is injected.

*V*

_{0}, is added to the blob injection velocity to account for the spreading of the spray near the nozzle [18]. The magnitude of the normal velocity component is calculated as

where *θ* is a uniform random number between 0 and Θ, and Θ is a user-defined cone angle. The direction of the normal component is random and orthogonal to the plume axis.

Atomization and breakup in this work are modeled using the Kelvin-Helmholtz/Rayleigh-Taylor (KH–RT) breakup models. Parent droplets (blobs) undergo primary breakup under the KH (wave) mechanism only [18], while child droplets may undergo primary or secondary breakup. In the latter case, the droplets are evaluated for secondary breakup (RT model) [23], and if no breakup occurs, they are evaluated for primary breakup (KH model).

## 3 Numerical Method

In this section, we present an overview of the numerical methods used to solve the Eulerian gas-phase fields, the Lagrangian liquid-phase field, and how the two solutions are coupled.

### 3.1 Gas Phase.

*Nek5000*code [24]. The equations are solved in their weak form, where one must find $u,T,Yi\u2208XbN$ and $p\u2208YN$ such that

*N*th-order spectral element basis functions, $XbN$ is the subset of

*X*satisfying the Dirichlet conditions on the boundary, $\u2202\Omega $, and $X0N$ is the subset of

^{N}*X*satisfying the homogeneous Dirichlet conditions on $\u2202\Omega $. Furthermore, the inner product is defined as

^{N}*E*is the total number of elements. The solutions for $u,\u2009p,T$ and

*Y*are represented by piecewise-defined

_{i}*N*th-order tensor-product Lagrange polynomial interpolants

where $r=(r,s,t)=(r1,r2,r3)$ are local coordinates in the element $\Omega \u0302=[\u22121,1]3$, whose image is mapped isoparametrically to each element, $xe(r)$, and $hi(\xi )$ are *N*th-order Lagrange polynomials having nodes at Gauss-Legendre-Lobatto (GLL) quadrature points, $\xi jN\u2208[\u22121,1]$, such that $hiN(\xi jN)=\delta ij$, with *δ _{ij}* being the Kronecker delta. In each element, there are $(N+1)3$ GLL quadrature or grid points, where the coefficients $uijke$ represent the numerical solution at the points. The test functions take the form$v=\u2211i=1N\u2211j=1N\u2211k=1N[hi(r)hj(s)hk(t)]$.

where *w _{i}* is the quadrature weight associated with GLL point

*ξ*. The pointwise Jacobian $Jijke=|\u2202xpe/rqe|ijk$ is associated with the mapping from local to global coordinates, $xe(r)$. $Me$ is the diagonal local mass matrix, and $u\xafe$ is the vector of $(N+1)3$ entries, $uijke$, ordered by lexicographical ordering $l=i+(N+1)(j\u22121)+(N+1)2(k\u22121)$. A fully assembled system for the global degrees-of-freedom can be assembled by gathering the local contributions at each nodal point, resulting in block diagonal system, which operates on the concatenated vector of solution coefficients $u\xafL={u\xafe}e=1E$. As a result, the work is effectively cast as highly vectorizable matrix–matrix products.

_{i}The method shows exponential convergence of the numerical solution as a function of polynomial order, *N* [8], making it a very accurate and efficient method. For more details on the method, the reader is referred to Ref. [7].

*BDF*2), and nonlinear terms are solved by a characteristics-based explicit scheme

**u**at the foot of the characteristics, and

*β*s are the standard coefficients for

_{j}*BDF2*derived from polynomial interpolation through timepoints. The values for $u\u0303n\u2212j$ are computed by solving the convective initial-boundary value problem

where *t ^{n}* is the time-step, and $u\u0303$ is treated as a passive vector field which is convected by velocity

**u**. The pure convection problem has the effect of propagating the initial condition forward, along the characteristics of the convecting field

**u**. The subproblem is solved explicitly using fourth-order Runge–Kutta (RK4) integrator with a small time-step $\Delta s<\Delta t$ to maintain stability. The OIFS method allows for larger time-step, $\Delta t$, than would be allowed with regular explicit/extrapolation schemes [7], which offsets the additional cost of having to solve the IBV subproblem, Eq. (40). In this study, a $CFL=|u|\Delta t/\delta x=2$ is used, where $\Delta t$ is the numerical time-step and $\Delta x$ is the grid spacing based on GLL points. More details on the OIFS temporal scheme may be found in Refs. [7,25].

To solve the momentum equation, a fractional step method is used where an initial guess for the velocity field is obtained from the OIFS scheme. The divergence operation is applied to the equation to arrive at the Poisson equation for pressure, which is solved using projection techniques. The computed pressure gradient field is used to correct the velocity field to satisfy the divergence condition, Eq. (28). The time derivatives in Eqs. (30) and (31) are solved implicitly using the variable-step *k*th-order time integrator, CVODE [26,27].

#### 3.1.1 Sub-Grid Scale Model.

In Eq. (42), $\chi >0$ sets the strength of the source term, typically set to satisfy $0<\chi \Delta t<1$ for stability reasons, where $\Delta t$ is the time-step for the gas-phase solution [31]. The choice of Lagrange polynomials as basis function for SEM makes $H$ semipositive definite, making $\u2212\chi H(u)$ purely dissipative. For LES, the filter cutoff is usually set to 65–90% of the highest-frequency modes resolved by the basis functions. The effect of this filter diminishes as polynomial order is increased, allowing for grid-independent results [29–31]. The relaxation term was also applied to the heat and species transport equations. No turbulence or turbulent dispersion model were used to account for SGS flow. Therefore, the dissipation from the HPF-RT model was used to approximate the energy dissipation occurring at subgrid scales.

### 3.2 Liquid Phase.

The liquid phase is handled by a Lagrangian particle solver *ppiclF* [32], which has been coupled with the *Nek5000* solver. The equations for droplet dynamics, droplet evaporation, liquid injection and breakup models (see Secs. 2.2 and 2.3) were implemented in the *Nek5000*-*ppiclF* framework for this work.

Parcel positions and velocities at each instant are calculated using Eqs. (13) and (14) for a single droplet. The system of ODEs resulting from Eqs. (13)–(16) is solved using an explicit strong stability-preserving third-order Runge–Kutta method (*SSP-RK3*).

Subcycling was used for the Lagrangian phase. The integration time-step for parcels was set to $\Delta ts=\Delta t/ns$, where $\Delta t$ is the time-step used for the gas-phase numerical solution and *n _{s}* is the number of subcycles. In this study,

*n*= 10 was used, which helped with the stability of the explicit method. In the current approach, the particle solution is advanced through time using a frozen gas velocity field over the duration of $\Delta t$. Interpolation of gaseous fields at parcel locations, $xd$, is done only for the first RK stage of each subcycle using spectrally accurate and efficient interpolation procedures [33]. These values are used to compute the right-hand-sides of Eqs. (14)–(16), and are kept constant throughout the RK solve for a single subcycle. After each subcycle, droplets are evaluated for breakup and distortion parameter values are updated.

_{s}### 3.3 Liquid-Gas Coupling.

*n*

_{d}where *δ _{f}* is the filter half-width, such that when $|r|=\delta f/2$, the kernel has decayed to half of its $|r|=0$ value. The value for

*δ*must be determined based on accuracy and computational cost [15]. For computational efficiency, the kernel is cut off at a specific radius,

_{f}*r*, to limit the influence of parcels to near-by grid points. In the current study, the cutoff is at $rc\u223c\delta f$, where the value of the kernel is near zero.

_{c}## 4 Computational Setup

The gasoline spray simulations with SEM are described in this section. Operating conditions and numerical parameters are presented.

### 4.1 Operating Conditions.

Standard conditions (nonflash boiling, nonreacting) from ECN Spray G [34] were simulated in this study for a single-component fuel (iso-octane). Ambient gas and liquid conditions are summarized in Table 1. In simulations, only nitrogen (*N*_{2}) was included as ambient gas. The injector has eight symmetrical holes, each with a drill angle of 37 deg relative to the injector axis. A schematic describing the nozzle geometry and reference angles is included in Fig. 1, and the relative coordinates of reference points T, N and H are listed in Table 2, where *z* is the axial direction. The injector has a nominal nozzle diameter of $d0=0.165$ mm and a counterbore diameter of $D0=0.388$ mm. Further details about the injector geometry or ambient conditions may be found on the ECN website [34].

Ambient gas temperature | 573 K |

Ambient gas pressure | 6.0 bar |

Ambient gas density | 3.5 kg/m^{3} |

Ambient gas | Nitrogen |

Fuel injection pressure | 200 bar |

Fuel | Iso-octane |

Fuel temperature at nozzle | 363 K |

Ambient gas temperature | 573 K |

Ambient gas pressure | 6.0 bar |

Ambient gas density | 3.5 kg/m^{3} |

Ambient gas | Nitrogen |

Fuel injection pressure | 200 bar |

Fuel | Iso-octane |

Fuel temperature at nozzle | 363 K |

Point | x (mm) | y (mm) | z (mm) |
---|---|---|---|

T | 0 | 0 | 0 |

N | 0.790 | 0 | −0.258 |

H | 0.513 | 0 | −0.625 |

Point | x (mm) | y (mm) | z (mm) |
---|---|---|---|

T | 0 | 0 | 0 |

N | 0.790 | 0 | −0.258 |

H | 0.513 | 0 | −0.625 |

Experimentally measured ROI [6] was used. Equal mass flowrate was applied across the holes. Parcels were injected at the center of the counterbore exit (point *N* in Fig. 1). Experiments had an electronic end of injection at 0.68 ms ASOI and an actual end of injection at 0.76 ms ASOI. In this work, simulations ended at 0.6 ms ASOI. This was done to avoid additional difficulties related to modeling end of injection, including incomplete atomization and injector tip wetting [35,36], which are outside of the scope of this paper. The simulated time was long enough to determine whether the simulations could capture the quasi-steady behavior in axial gas velocity observed in experiments [3,6].

### 4.2 Domain and Mesh.

The domain was cylindrical with diameter *D *=* *100 mm and length *L *=* *60 mm. The nozzle geometry was not included in the simulations. The injection point for blobs (point N) was located 0.270 mm away from the head of the cylinder, while the injector tip is 0.528 mm away from the wall. Wall boundary conditions were imposed everywhere except at the bottom face (opposing the nozzle), where a robust outflow condition was used [37]. The gas flow field was initially quiescent ($u=0$) with a uniform temperature of 573 K, which are appropriate assumptions for conditions in the experimental vessel [6,38].

A conforming, unstructured mesh consisting of hexahedral elements was used, as required by the solver [24]. Polynomial order of *N *=* *5 was used in this study, which was previously found to be appropriate for predicting spray penetration in LES of nonevaporative sprays with SEM [10]. Different levels of refinements were applied throughout the domain, as illustrated in Fig. 2. Near the nozzle, elements have approximately the same size as the nozzle diameter ($\Delta xe\u2248d0=0.165\u2009mm$), which for the given polynomial order (*N *=* *5) results in an effective grid spacing of $\Delta xe/N=\Delta x=0.033$ mm. Minimum grid spacing for LES of Spray G found in the literature is typically between 0.125 mm [6] and 0.375 mm [39]. Therefore, the current level of resolution is about 4–11 times finer than in those studies, where finite-volume method was used.

The mesh becomes coarser outside of the near-nozzle region, with a ratio $\Delta xcoarse/\Delta xfine\u223c3$ between adjacent refinement regions. The mesh had a total of 270k elements and 33.9 million unique grid points. Each simulation ran on 720 Intel Broadwell (Intel Xeon E5-2695v4) cores for about 3.8 days.

The Gaussian filter half-width was set to $\delta f=2d0=2\Delta xe=0.330$ mm in this work. The ratio of $\delta f/\Delta xe=2$ was found to accurately predict momentum transfer in nonevaporative sprays, while preserving the stability of the solution [10].

### 4.3 Spray Breakup and Injection Parameters.

Five different setups were used for the Spray G simulations. All simulations had the same area contraction coefficient (*C _{a}*), blob mass ($md,0$), and breakup parameters (except for

*B*

_{1}), as listed in Table 3. The value for

*C*was obtained from experimental measurements in Ref. [40], where $Ca=0.635\u22120.687$ was reported for the current injector at injection pressures close to the simulated conditions. In this study, $nd\u223c0.1$ droplets per blob were used. This resulted in about $\u223c91,000$ injected blobs ($\u223c11,400$/hole) for a simulation time of 0.6 ms ASOI. After breakup and evaporation, approximately 450k parcels were present in the domain at 0.6 ms ASOI. The number of injected blobs is relatively low compared to what is typically done with finite-volume codes [41]. That is because the current method did not require increasing injected blobs to maintain stability of the solution at the current level of grid refinement. The effect of this parameter was not analyzed in this study, and will be analyzed in future work.

_{a}Area contraction coefficient, C_{a} | 0.65 |

Blob mass, $md,0$ | 8.43 × 10^{−5} mg |

KH model size constant, B_{0} | 0.61 |

RT model time constant, $C\tau $ | 1.0 |

RT model size constant, C_{RT} | 0.6 |

Area contraction coefficient, C_{a} | 0.65 |

Blob mass, $md,0$ | 8.43 × 10^{−5} mg |

KH model size constant, B_{0} | 0.61 |

RT model time constant, $C\tau $ | 1.0 |

RT model size constant, C_{RT} | 0.6 |

For this study, the KH model time constant, *B*_{1}, and injection cone angle (Θ) were varied to examine their effects on jet penetration, spray width and gas velocity. The effect of using fixed versus variable plume direction was also analyzed, using experimentally measured profile [38] shown in Fig. 3. Finally, the effect of using the measured nozzle diameter ($d0=0.170$ mm [34]) versus the nominal nozzle diameter ($d0=0.165$ mm) was also examined. In this case, mass flowrate was kept the same, resulting in a lower injection velocity and larger droplets when the real nozzle diameter was used. The different cases are listed in Table 4. The effect of each parameter is evaluated by comparing against the baseline case, case 1. The value of $B1=80$ was found to accurately predict spray penetration in a previous study on nonevaporating sprays using the current method [10], and therefore was used for the baseline case. The cone angle $\Theta =25\u2009deg$ was chosen from experimental values calculated using tomographic data and schlieren boundary, at 50% intermittency [34,38].

## 5 Results

In this section, numerical results are compared against experimental data in terms of liquid penetration, jet width, and axial gas velocity. All simulation results are obtained from a single injection event for each case.

### 5.1 Liquid Jet.

*θ*, over the line of sight

_{l}In this study, the ECN-recommended threshold of PLV = 0.2 × 10^{−3} mm^{3}/mm^{2} was used to determine the liquid boundary.

Figure 4 shows that all cases over-predict the liquid penetration during early injection (<0.2 ms ASOI), but case 3 (30 deg cone angle) produced very good agreement with experimental data at later times. Cases 2 and 5, with variable plume and larger nozzle diameter, respectively, showed improved penetration with respect to the baseline case. A smaller *B*_{1} value produced lower penetration early in the injection (<0.1 ms ASOI), but overall the effect of this parameter on liquid penetration was not significant.

The width of the spray at an axial location of 15 mm downstream of the injector tip is shown in Fig. 5. This width is measured as the distance between the liquid boundaries in the radial direction. Case 1 showed an early rise in width, which is due to the faster liquid penetration in this case. Cases 1, 3, and 4 peaked at around the same value ($\u223c27.5$ mm), while case 2 resulted in a larger peak ($\u223c30$ mm) and closer agreement with experiments. The increased radial penetration is due to the higher plume direction angle early during injection, which varies from 35.7 deg to 33.5 deg during the first 0.2 ms ASOI, compared to 33 deg in the fixed plume direction cases. Width in case 5 also rose to a similar peak value as case 2, but the behavior was delayed. The larger width in case 5 might be due to the presence of larger droplets at the edge of the spray compared to other cases, as will be seen in Sec. 5.2. This increases the PLV value near the edge of the spray, thus enlarging the area enclosed by the PLV threshold.

A comparison of projected liquid volume fraction between Cases 1-3 at 0.62 ms ASOI is shown in Fig. 6, where Cases 2 and 3 were the best-performing setups in terms of spray width and liquid penetration, respectively. The contours from LES and experiments correspond to the ECN-recommended threshold of PLV = 0.2 × 10^{−3} mm^{3}/mm^{2}. The figure confirms the behavior in Fig. 4, i.e., liquid penetration is better predicted in Cases 2 and 3 than in case 1. Case 2 produces a wider spray, showing better agreement with experiments at the location of *z* = 15 than the other two cases. However, it overpredicts the radial penetration at the tip of the jet. The liquid distribution between the plumes seems to be better predicted in case 3, since the contour between plumes is in better agreement with experiments, indicated by the red line penetrating less toward the nozzle between plumes in case 3 compared to Cases 1 and 2.

### 5.2 Droplet Diameter.

Numerical results are compared against experimental data obtained using phase-Doppler interferometry (PDI) [42]. Droplet sizes were measured at a fixed axial location of *z* = 15 mm, for different radial distances, *r*, relative to the injector axis. The droplet diameters were sampled along a line that crosses the center of the plume, as shown in Fig. 7. Experimental measurements were done for a single plume, over multiple injection events. Simulation results were an average of the eight plumes from a single injection event.

The SMDs were obtained over a grid of 1 mm resolution in the radial direction. For numerical results, a box of 1 mm^{3} was used at each sampling station to increase sampling. In experiments, the temporal resolution is of 0.1 ms, while in simulations the temporal resolution is 0.025 ms. Results for SMD in simulations and experiment as a function of radial location and time are shown in Fig. 8. It should be mentioned that the experimental setup struggled to obtain reliable measurements near the center of the plume ($r\u223c10$ mm), resulting in a low droplet count near this region [42].

All cases underpredict SMD near the center of the plume, $r\u223c10$ mm. SMD increases toward the edge of the plume (*r *<* *7 mm and *r *>* *13 mm) in most cases, as well as toward the spray tip (*t *=* *0.2 ms ASOI), where there is closer agreement with experimental values. As expected, case 5, with the larger nozzle diameter, predicts larger droplets in these regions, overpredicting droplet size near the spray tip and toward the injector axis (*r *<* *7 mm). Using a smaller *B*_{1} value (case 4) underpredicts SMDs more drastically than other cases. Furthermore, case 2 fails to capture the presence of droplets at $5<r<6$ mm at times earlier than 0.4 ms ASOI. Therefore, using the variable plume direction from experiments (case 4) did not produce accurate distribution of liquid droplets. Cases 1 and 3 yield very similar results, with case 3 having slightly better prediction of SMD toward the edge of the plume.

### 5.3 Gas Velocity.

Numerical results for gas velocity at a location of (*x*,*y*,*z*) = (0,0,15) mm downstream of the nozzle are compared against experimental PIV data [6] in Fig. 9. The negative value of axial velocity means that the gas is moving toward the nozzle due to recirculating flow. While spray width was improved by using experimental data to impose plume direction, axial velocity in case 2 does not follow the same trend as in the experiment, having a slower decay in time. The closest agreement with experiment was achieved in case 1, with Cases 3 and 4 showing similar trends but stabilizing at a lower-magnitude velocity. Using a larger nozzle diameter did not improve reduction of axial velocity at this location, compared to the baseline case. All cases showed delayed start of recirculation compared to the experiment. While the difference between numerical and experimental results is not negligible, current results are an improvement compared to results in previous studies using the finite-volume method [6], since current results are able to capture the axial recirculation while matching experimental liquid penetration.

The effect of injection parameters is more evident when axial velocity is evaluated closer to the nozzle, as depicted in Fig. 10. The case with the larger injection cone angle (case 3) produced a velocity magnitude of almost double that of the other cases near the nozzle (*z *<* *11 mm). All cases tend to converge at similar axial velocity for *z *>* *11 mm. This highlights that the aerodynamics around the plumes are more sensitive to the different injection parameters in the near-nozzle region, where it is difficult to measure gas velocity experimentally.

The reason for the discrepancy between case 1 and case 3 may be explained in terms of the proximity between the plumes, as illustrated by the velocity field in Fig. 11. The white lines in the figure indicate the contour of the liquid jet, using a threshold of liquid volume fraction of $\theta l=0.001$. The flow patterns in both simulations are very similar, but the case with 30 deg cone angle produced a narrower area between the plumes in the near-nozzle region, with the edges of the wider plumes coming closer to each other. The gas is thus entrained through a smaller area and accelerated as a result. In contrast, further away from the nozzle the liquid core becomes more diffused in case 3, which could be due to the fact that parcels are more dispersed in the direction normal to the plume direction, as imposed by the definition of the injection cone angle (see Eq. (27)). Therefore, the gap between the plumes become larger in case 3, which contributes to the decrease in recirculation velocity magnitude at this location.

## 6 Conclusions

A novel approach for simulating gasoline sprays was presented, where the spectral element method was used for the solution of the gaseous field in an Eulerian reference frame. The liquid droplets were modeled using the well-known stochastic Lagrangian parcels approach. The current method allowed for grid resolutions which are much finer (*x*4–*x*11) than what is typically used with second-order methods, such as the finite volume method.

The ECN Spray G under standard conditions (nonreacting, nonflash-boiling) [34,43] was simulated in this study. Five different configurations were used, by changing injection cone angle, different spray breakup model parameters, nozzle diameter and using experimentally measured plume direction.

Using the experimental plume direction profile and a larger cone-angle degree were both strategies that reduced liquid penetration in simulations compared to the baseline case (33 deg plume direction angle, 25 deg cone angle), and improved agreement with experiments. Using a smaller *B*_{1} value did not seem to affect numerical results significantly. The case with time-varying plume direction (case 2) also predicted a wider jet, which agreed well with experimental data. However, case 2 had the worst performance amongst all cases in terms of predicting axial velocity at a point 15 mm downstream of the injector tip, since it did not capture the trend observed in PIV data, showing that this strategy might not be beneficial for the overall prediction of the spray flow dynamics. At this location (*z* = 15), the baseline case was closer to capturing the behavior in experiments than other cases, since it produced lower (higher-magnitude) axial velocity. Other cases (with the exception of case 2) followed the same trend, but stabilized at a slower axial velocity. On the other hand, the gas behavior was very different between cases closer to the nozzle, where the case with larger cone angle predicted more vigorous recirculation, indicated by a higher magnitude of axial velocity. This suggests that a wider cone angle could potentially trigger faster plume collapse due to increased axial recirculation [3,4].

The current modeling approach allowed a highly resolved solution of the gaseous field in the recirculation region, which reduces uncertainty in results due to numerical accuracy in the gas-phase solution. Current results show improved prediction of recirculation compared to LES studies done with finite volume method [6]. However, this study shows that, even with improved accuracy for the gas phase, numerical results are highly sensitive to input parameters for the Lagrangian liquid phase, which agrees with observations from previous LES studies [6]. This highlights the need for models and techniques to predict appropriate input parameters for the Lagrangian phase and reduce the amount of tuning necessary for validating spray simulations.

Droplet sizes determine the breakup, momentum, heat and mass transfer processes in sprays, therefore, simulations should be able to accurately capture droplet diameters. Numerically predicted droplet sizes in Cases 1 and 3 match experimental droplet size at the edges and at the tip of the plume. Toward the center of the plume, however, droplet sizes in Cases 1 and 3 were underpredicted compared to PDI data.

It should be noted that current simulations did not include interdroplet collision/coalescence models. It is possible that including such models in future simulations could improve predicted SMDs, since smaller droplets could coalesce to form larger droplets, especially toward the center of the plume where liquid volume fraction is higher.

Overall, case 3 (33 deg cone angle) performed better than other cases, since it was able to accurately predict liquid penetration, the spray morphology (see Fig. 6) and the trend in axial velocity between the plumes, leading to believe that the gas velocity near the nozzle would behave as in case 3. Droplet sizes were also close to experimental values, especially near the edges and tip of the spray.

Future work includes varying other spray parameters, like area contraction coefficient and combining some of the cases studied here, for example, using a variable plume direction with a higher cone angle. The effect of increasing the number of injected blobs on the predicted spray will also be studied. Furthermore, a one-way coupled approach will be implemented into the *Nek5000* solver, which will use maps from internal nozzle flow simulations as inputs for the Lagrangian phase [3]. Finally, droplet collision/coalescence models will be implemented in the current modeling approach to examine its effect in predicted SMDs and liquid penetration.

## Acknowledgment

The submitted paper has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. This research was conducted as part of the Partnership to Advance Combustion Engines (PACE) sponsored by the U.S. Department of Energy (DOE) Vehicle Technologies Office (VTO). A special thanks to DOE VTO program managers Mike Weismiller and Gurpreet Singh. The authors would like to thank the Laboratory Computing Resource Center (LCRC) for providing the computational resources to run the simulations for this work. We would also like to thank the multiphase flow modeling team at Argonne, especially Lorenzo Nocivelli, Roberto Torelli, and Gina Magnotti, whose input was key for the advancement of this work.

## Funding Data

U.S. Department of Energy (DOE) Vehicle Technologies Office (VTO) (No. 100006177; Funder ID: 10.13039/100011884).

## References

**290**, p.

**10**(