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 [35]. 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.

In this work, a low-Mach formulation of the Navier–Stokes equations is used to account for compression effects due to thermal dilatation, while neglecting the effect of acoustic or shock waves [11]. Regular perturbation techniques can be used to decouple the waves from the governing equations [12], which leads to a decomposition of the pressure as:
p(x,t)=p0(t)+εp1(x,t)
(1)

where p1 is the hydrodynamic pressure, p0 is the thermodynamic pressure, and ε=γMa2, where γ is the ratio of specific heat capacities and Ma is the Mach number. The thermodynamic pressure, p0, is only a function of time. In this work, it was assumed constant, i.e., dp0/dt=0.

For multiphase flows, the gas-phase governing equations may be derived by applying volume filtering [1315]. The corresponding continuity, momentum, heat and species transport equations for Ng-component gaseous mixtures are
·u=1ρDρDt=1ρgDρgDt1θgDθgDt
(2)
ρg(ut+u·u)=p1+·τ+Sm
(3)
ρgcp(Tt+u·T)=·(kT)+Sh
(4)
ρg(Yit+u·Yi)=·Ji+Ss
(5)
τ=μ[u+(u)T2(·u)I/3]
(6)
Sm=1θgfd
(7)
Sh=1θg(q˙+i=1NgΔhim˙i)
(8)
Ss=1θgm˙i
(9)

Here, q˙, hi, m˙i, and Yi are the heat transfer rate due to conduction/convection, enthalpy, mass flowrate due to phase change, and mass fraction of species i, respectively. The term Δhi=hi(Td)hi(T) is the difference in enthalpy calculated at droplet temperature, Td, and at local gas temperature, T. This represents the heat loss due to the mass transfer from the evaporating liquid. The gas density is ρ=ρgθg, where ρg mixture density, θg=1θl is the void fraction and θl is the liquid volume fraction. Terms including void fraction, θg, 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.

The species mass diffusion flux Ji is approximated by Fick's law
Ji=ρgDYi
(10)
where D is the diffusion coefficient. Gas mixture heat capacity is computed as cp=i=1Ngcp,iYi, where cp,i is the heat capacity of gaseous species i. Density of the mixture is given by the ideal gas law
ρg=Wp0/RT
(11)

where W=(i=1NgYi/Wi)1 is the molecular weight of the mixture and Wi 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].

The first term on the right-hand side of Eq. (2) is calculated as follows:
1ρgDρgDt=1TDTDt+i=1NgWWiDYiDt=1ρgcpT[·kT+1θg(q˙+i=1NgΔhim˙i)]+1ρgi=1NgWWi(·Ji+1θgm˙i)
(12)

2.2 Liquid Droplet Equations.

The liquid phase is represented by groups of droplets, or parcels, that are tracked in a Lagrangian reference frame [18]. The droplets are assumed to have uniform properties within each parcel. The evolution of each droplet is described by
DxdDt=vd
(13)
DvdDt=1md(FD+FC+FB)d+g
(14)
DTdDt=1cp,lmd[πkdd(T(xd)Td)Nud+hLDmdDt]
(15)
DddDt=2ρgDρlddBMShd
(16)
DmdDt=12πρldd2DddDt
(17)
where xd is the position vector of the droplets and vd is their velocity, T and Td are the gas and droplet temperatures, respectively; hL is the latent heat at Td, dd is the droplet diameter, and ρl is liquid density. The Spalding mass transfer number is defined as
BM=YFsYg1YFs
(18)
YFs=WFWF+Wg(p0/pv1)
(19)
where YFs is the fuel vapor mass fraction at the surface of the droplet, Yg is the fuel vapor mass fraction in the gaseous mixture, WF is the fuel molecular weight, Wg is the molecular weight of the ambient gas, and pv is the vapor pressure of the fuel at droplet temperature. The Sherwood (Shd) and Nusselt (Nud) 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]
Tv=(2Ts+Tg)/3
(20)
Yv=(2Ys+Yg)/3
(21)

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,FD,d,FC,d and FB,d are gravity, drag, interparticle and buoyancy forces, respectively. The mass of the droplet is md=ρlπdd3/6. For this study, interparticle interactions and gravity forces were neglected: FC,d=FB,d=g=0.

The drag force acting on the droplets is calculated as in terms of the drag coefficient, CD
FD,d=12ρgCDAf|Ud|Ud
(22)
where Ud=u(xd)vd is the relative velocity between the droplet and the gas at the droplet position xd and Af is the projected frontal area of the droplet, which is assumed to be spherical. The drag coefficient of a sphere is given by [21]
CD,sphere={0.424,Re>100024Re(1+16Re2/3),Re1000
(23)
To account for droplet deformation, the distortion parameter, yd, 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
CD=CD,sphere(1+2.635yd)
(24)

2.3 Liquid Injection and Droplet Breakup.

For modeling liquid injection, the blob injection method was used. Parcels containing parent droplets (blobs) with diameter equal to the nozzle whole diameter, d0, are introduced at the injection location. Droplet diameter is corrected to account for area contraction coefficient, Ca, of the liquid jet passing through the nozzle hole
dd,0=Cad0
(25)
For cylindrical liquid jets, the droplet velocity is determined from the rate of injection (ROI), also referred to as mass flowrate, m˙l
U0=4m˙l/πρldd,02
(26)

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=ρlndπdd,03/6, at which point a blob is injected.

A velocity component normal to the plume axis, V0, 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
V0=U0tan(θ/2)
(27)

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.

The governing equations for the gas phase, Eqs. (2)(5), are solved in an Eulerian reference frame using the spectral element method (SEM) [8], implemented in the scalable MPI-based Nek5000 code [24]. The equations are solved in their weak form, where one must find u,T,YiXbN and pYN such that
(q,·u)=(1ρDρDt,q)
(28)
(v,ρgDuDt)=(·v,p1)(v,τ)+(v,Sm)
(29)
(q,ρgcpTt)=(q,ρgcpu·T)+(q,kT)+(q,Sh)
(30)
(q,ρgYit)=(q,ρgu·Yi)+(q,Ji)+(q,Ss)
(31)
for all test functions vX0N and qYN in the computational domain Ω. Density is removed from the equations using the equation of state, Eq. (11). Here, the PNPN formulation is used, where YN=XN is the set of continuous Nth-order spectral element basis functions, XbN is the subset of XN satisfying the Dirichlet conditions on the boundary, Ω, and X0N is the subset of XN satisfying the homogeneous Dirichlet conditions on Ω. Furthermore, the inner product is defined as
(p,q)=Ωpqdx
(32)
(u,v)=Ωu·vdx
(33)
The domain is decomposed into elements, Ω=e=1EΩe, where E is the total number of elements. The solutions for u,p,T and Yi are represented by piecewise-defined Nth-order tensor-product Lagrange polynomial interpolants
ue(r,s,t)=i=1Nj=1Nk=1Nhi(r)hj(s)hk(t)uijke
(34)

where r=(r,s,t)=(r1,r2,r3) are local coordinates in the element Ω̂=[1,1]3, whose image is mapped isoparametrically to each element, xe(r), and hi(ξ) are Nth-order Lagrange polynomials having nodes at Gauss-Legendre-Lobatto (GLL) quadrature points, ξjN[1,1], such that hiN(ξjN)=δ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 formv=i=1Nj=1Nk=1N[hi(r)hj(s)hk(t)].

The choice of Lagrangian interpolants and GLL quadrature points in SEM allows for the inner product to be accurately approximated by pointwise quadrature
(v,u)e=111111[hl(r)hm(s)hn(t)][hi(r)hj(s)hk(t)]uijkeJedrdsdti=1Nj=1Nk=1Nwiwjwk[hl(ξi)hm(ξj)hn(ξk)]
(35)
[hi(ξi)hj(ξj)hk(ξk)]uijkeJijke
(36)
=i=1Nj=1Nk=1NwiwjwkδilδjmδknJijkeuijke
(37)
:=Meu¯e
(38)

where wi is the quadrature weight associated with GLL point ξi. The pointwise Jacobian Jijke=|xpe/rqe|ijk is associated with the mapping from local to global coordinates, xe(r). Me is the diagonal local mass matrix, and u¯e is the vector of (N+1)3 entries, uijke, ordered by lexicographical ordering l=i+(N+1)(j1)+(N+1)2(k1). 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¯L={u¯e}e=1E. As a result, the work is effectively cast as highly vectorizable matrix–matrix products.

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].

The material derivative in Eq. (29) is solved using an operator-integration-factor splitting scheme (OIFS) for time integration, where linear terms are solved implicitly with second-order backward differentiation (BDF2), and nonlinear terms are solved by a characteristics-based explicit scheme
DuDtβ0Δtun+1Δtj=12βjũnj
(39)
The values ũnj represent the field u at the foot of the characteristics, and βjs are the standard coefficients for BDF2 derived from polynomial interpolation through timepoints. The values for ũnj are computed by solving the convective initial-boundary value problem
ũs+u·ũ=0,s[tnk,tn]
(40)
ũ(x,tnk)=u(x,tnk)
(41)

where tn is the time-step, and ũ 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 Δs<Δt to maintain stability. The OIFS method allows for larger time-step, Δ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|Δt/δx=2 is used, where Δt is the numerical time-step and Δ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 kth-order time integrator, CVODE [26,27].

3.1.1 Sub-Grid Scale Model.

For under-resolved flows at high Reynolds numbers, the solution can become unstable since the method has low numerical dissipation, making it susceptible to the growth of instabilities. Filter-based approaches may be used to remove energy due to high-frequency modes, having a stabilizing and dissipative effect similar to that of eddy viscosity-type SGS models [2831]. In the present study, the high-pass filter relaxation term (HPF-RT) method is used for large-eddy simulations simulations [30,31]. In this method, the high-pass filtered velocity, H(u), is applied as a source term to the momentum transport equation (Eq. (3)), which introduces dissipation to the system
Sm=1θgfdχH(u)
(42)

In Eq. (42), χ>0 sets the strength of the source term, typically set to satisfy 0<χΔt<1 for stability reasons, where Δ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 χ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 [2931]. 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 Δts=Δt/ns, where Δt is the time-step used for the gas-phase numerical solution and ns is the number of subcycles. In this study, ns = 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 Δ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.

3.3 Liquid-Gas Coupling.

The source terms for momentum, heat and mass transfer to the gas phase (Eqs. (7)(9)) are applied at grid points using a Gaussian projection filter [15] at each parcel location. The strength of the source term is determined as the product of the single-droplet momentum/heat/mass transfer times the total number of droplets in each parcel, nd
θl(x)=j=1NpG(r)ndVd,j
(43)
fd(x)=j=1NpG(r)ndFD,d,j
(44)
q˙(x)=i=jNpG(r)ndπkdd(T(xd,j)Td,j)Nud
(45)
Δhim˙i(x)=j=1NpG(r)ndm˙d,j[hi(Td,j)hi(T(xd,j))]
(46)
m˙i(x)=j=1NpG(r)ndm˙d,j
(47)
m˙d,j=Dmd,jDt
(48)
where r=xxd and Vd is the droplet volume. The kernel for the Gaussian filter is given by
G(r)=(πδf2/4ln2)3exp[|r|2δf2/4ln2]
(49)

where δf is the filter half-width, such that when |r|=δf/2, the kernel has decayed to half of its |r|=0 value. The value for δf must be determined based on accuracy and computational cost [15]. For computational efficiency, the kernel is cut off at a specific radius, rc, to limit the influence of parcels to near-by grid points. In the current study, the cutoff is at rcδf, where the value of the kernel is near zero.

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 (N2) 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].

Fig. 1
Schematic of Spray G geometry [34]
Fig. 1
Schematic of Spray G geometry [34]
Close modal
Table 1

Simulated conditions based on ECN Spray G operating conditions [34]

Ambient gas temperature573 K
Ambient gas pressure6.0 bar
Ambient gas density3.5 kg/m3
Ambient gasNitrogen
Fuel injection pressure200 bar
FuelIso-octane
Fuel temperature at nozzle363 K
Ambient gas temperature573 K
Ambient gas pressure6.0 bar
Ambient gas density3.5 kg/m3
Ambient gasNitrogen
Fuel injection pressure200 bar
FuelIso-octane
Fuel temperature at nozzle363 K
Table 2

Coordinates of reference points for hole #1 relative to nozzle tip [34]

Pointx (mm)y (mm)z (mm)
T000
N0.7900−0.258
H0.5130−0.625
Pointx (mm)y (mm)z (mm)
T000
N0.7900−0.258
H0.5130−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 (Δxed0=0.165mm), which for the given polynomial order (N =5) results in an effective grid spacing of Δxe/N=Δ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.

Fig. 2
Slice of the spectral-element mesh at the injector axis, colored by effective grid spacing, Δx (mm), using a logarithmic color scale. Grid points internal to the elements are not shown for clarity.
Fig. 2
Slice of the spectral-element mesh at the injector axis, colored by effective grid spacing, Δx (mm), using a logarithmic color scale. Grid points internal to the elements are not shown for clarity.
Close modal

The mesh becomes coarser outside of the near-nozzle region, with a ratio Δxcoarse/Δxfine3 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 δf=2d0=2Δxe=0.330 mm in this work. The ratio of δf/Δ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 (Ca), blob mass (md,0), and breakup parameters (except for B1), as listed in Table 3. The value for Ca was obtained from experimental measurements in Ref. [40], where Ca=0.6350.687 was reported for the current injector at injection pressures close to the simulated conditions. In this study, nd0.1 droplets per blob were used. This resulted in about 91,000 injected blobs (11,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.

Table 3

Fixed injection and spray breakup parameters

Area contraction coefficient, Ca0.65
Blob mass, md,08.43 × 10−5 mg
KH model size constant, B00.61
RT model time constant, Cτ1.0
RT model size constant, CRT0.6
Area contraction coefficient, Ca0.65
Blob mass, md,08.43 × 10−5 mg
KH model size constant, B00.61
RT model time constant, Cτ1.0
RT model size constant, CRT0.6

For this study, the KH model time constant, B1, 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 Θ=25deg was chosen from experimental values calculated using tomographic data and schlieren boundary, at 50% intermittency [34,38].

Fig. 3
Experimental profile of plume direction angle [38]
Fig. 3
Experimental profile of plume direction angle [38]
Close modal
Table 4

Simulation injection and spray parameters

Cased0 (μm)Plume directionCone angle, ΘB1
116533 deg25 deg80
2165Variable25 deg80
316533 deg30 deg80
416533 deg25 deg40
517033 deg25 deg80
Cased0 (μm)Plume directionCone angle, ΘB1
116533 deg25 deg80
2165Variable25 deg80
316533 deg30 deg80
416533 deg25 deg40
517033 deg25 deg80

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.

Liquid penetration was determined using the technique recommended by ECN. This consists of computing the projected liquid volume fraction (PLV) by integrating the liquid volume fraction, θl, over the line of sight
PLV(y,z)=θl(x,y,z)dx
(50)

In this study, the ECN-recommended threshold of PLV = 0.2 × 10−3 mm3/mm2 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 B1 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.

Fig. 4
Liquid penetration based on the ECN-recommended threshold—PLV = × 10−3 mm3/mm2. Experimental line from DBI data [38]. Shaded region represents experimental standard deviation.
Fig. 4
Liquid penetration based on the ECN-recommended threshold—PLV = × 10−3 mm3/mm2. Experimental line from DBI data [38]. Shaded region represents experimental standard deviation.
Close modal

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 (27.5 mm), while case 2 resulted in a larger peak (30 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.

Fig. 5
Jet width. Experimental line from DBI data in Ref. [38].
Fig. 5
Jet width. Experimental line from DBI data in Ref. [38].
Close modal

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 mm3/mm2. 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.

Fig. 6
Projected liquid volume (PLV) fraction (mm3/mm2). (a) Case 1—25 deg cone angle, (b) Case 2—variable plume direction, (c) Case 3—30 deg cone angle. Spray contours shown for PLV = 0.2 × 10−3 mm3/mm2. Single line—LES, double line—experiment [38].
Fig. 6
Projected liquid volume (PLV) fraction (mm3/mm2). (a) Case 1—25 deg cone angle, (b) Case 2—variable plume direction, (c) Case 3—30 deg cone angle. Spray contours shown for PLV = 0.2 × 10−3 mm3/mm2. Single line—LES, double line—experiment [38].
Close modal

5.2 Droplet Diameter.

Droplet size is presented here in terms of Sauter mean diameter, SMD, which is defined as
SMD=i=1Npnddd3i=1Npnddd2
(51)

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.

Fig. 7
Schematic of the multiplume spray, showing locations where droplet diameters were sampled
Fig. 7
Schematic of the multiplume spray, showing locations where droplet diameters were sampled
Close modal

The SMDs were obtained over a grid of 1 mm resolution in the radial direction. For numerical results, a box of 1 mm3 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 (r10 mm), resulting in a low droplet count near this region [42].

Fig. 8
Droplet size represented by SMD (μm). (a) Case 1, (b) case 2, (c) case 3, (d) case 4, (e) case 5, (f) experimental measurement [42].
Fig. 8
Droplet size represented by SMD (μm). (a) Case 1, (b) case 2, (c) case 3, (d) case 4, (e) case 5, (f) experimental measurement [42].
Close modal

All cases underpredict SMD near the center of the plume, r10 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 B1 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.

Fig. 9
Axial gas velocity 15 mm downstream of nozzle tip. The shaded area shows the experimental standard deviation [6].
Fig. 9
Axial gas velocity 15 mm downstream of nozzle tip. The shaded area shows the experimental standard deviation [6].
Close modal

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.

Fig. 10
Axial gas velocity along the injector axis at 0.58 ms ASOI. Experimental PIV data from Ref. [6].
Fig. 10
Axial gas velocity along the injector axis at 0.58 ms ASOI. Experimental PIV data from Ref. [6].
Close modal

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 θ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.

Fig. 11
Slice at yz-plane (x = 0) at 0.58 ms ASOI, colored by gas axial velocity. The vectors indicate the direction of gas velocity, but are not scaled by velocity magnitude. (a) Case 1–25 deg cone angle, (b) Case 3–30 deg cone angle. The solid line denotes the θl=0.001 level, indicating the location of the liquid jet.
Fig. 11
Slice at yz-plane (x = 0) at 0.58 ms ASOI, colored by gas axial velocity. The vectors indicate the direction of gas velocity, but are not scaled by velocity magnitude. (a) Case 1–25 deg cone angle, (b) Case 3–30 deg cone angle. The solid line denotes the θl=0.001 level, indicating the location of the liquid jet.
Close modal

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 (x4–x11) 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 B1 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

1.
Geiger
,
J.
,
Grigo
,
M.
,
Lang
,
O.
,
Wolters
,
P.
, and
Hupperich
,
P.
,
1999
, “
Direct Injection Gasoline Engines - Combustion and Design
,”
SAE
No. 1999-01-0170.10.4271/1999-01-0170
2.
Liu
,
H.
,
Wang
,
Z.
,
Long
,
Y.
,
Xiang
,
S.
,
Wang
,
J.
, and
Fatouraie
,
M.
,
2015
, “
Comparative Study on Alcohol-Gasoline and Gasoline-Alcohol Dual-Fuel Spark Ignition (DFSI) Combustion for Engine Particle Number (PN) Reduction
,”
Fuel
,
159
, pp.
250
258
.10.1016/j.fuel.2015.06.059
3.
Sphicas
,
P.
,
Pickett
,
L. M.
,
Skeen
,
S. A.
, and
Frank
,
J. H.
,
2018
, “
Inter-Plume Aerodynamics for Gasoline Spray Collapse
,”
Int. J. Engine Res.
,
19
(
10
), pp.
1048
1067
.10.1177/1468087417740306
4.
Rachakonda
,
S. K.
,
Paydarfar
,
A.
, and
Schmidt
,
D. P.
,
2019
, “
Prediction of Spray Collapse in Multi-Hole Gasoline Direct-Injection Fuel Injectors
,”
Int. J. Engine Res.
,
20
(
1
), pp.
18
33
.10.1177/1468087418819527
5.
Guo
,
H.
,
Nocivelli
,
L.
, and
Torelli
,
R.
,
2021
, “
Numerical Study on Spray Collapse Process of ECN Spray G Injector Under Flash Boiling Conditions
,”
Fuel
, 290, p.
119961
10.1016/j.fuel.2020.119961
6.
Sphicas
,
P.
,
Pickett
,
L. M.
,
Skeen
,
S.
,
Frank
,
J.
,
Lucchini
,
T.
,
Sinoir
,
D.
,
D'Errico
,
G.
,
Saha
,
K.
, and
Som
,
S.
,
2017
, “
A Comparison of Experimental and Modeled Velocity in Gasoline Direct-Injection Sprays With Plume Interaction and Collapse
,”
SAE Int. J. Fuels Lubricants
,
10
(
1
), pp.
184
201
.10.4271/2017-01-0837
7.
Deville
,
M. O.
,
Fischer
,
P. F.
,
Fischer
,
P. F.
, and
Mund
,
E. H.
,
2002
,
High-Order Methods for Incompressible Fluid Flow
, Cambridge Monographs on Applied and Computational Mathematics,
Cambridge University Press, Cambridge, UK
.
8.
Patera
,
A. T.
,
1984
, “
A Spectral Element Method for Fluid Dynamics: Laminar Flow in a Channel Expansion
,”
J. Comput. Phys.
,
54
(
3
), pp.
468
488
.10.1016/0021-9991(84)90128-1
9.
Colmenares
,
J. D.
,
Ameen
,
M. M.
,
Wu
,
S.
, and
Patel
,
S.
,
2021
, “
Large Eddy Simulation of Turbulent Particle-Laden Jets Using the Spectral Element Method
,”
AIAA
Paper No.
2021
0635
.10.2514/6.2021-0635
10.
Colmenares F
,
J. D.
,
Ameen
,
M. M.
, and
Patel
,
S. S.
,
2020
, “
Large-Eddy Simulation of Non-Vaporizing Sprays Using the Spectral-Element Method
,”
73rd Annual Meeting of the APS Division of Fluid Dynamics
, Virtual, Nov.
22
24
.10.1016/j.ijmultiphaseflow.2022.104155
11.
Tomboulides
,
A. G.
,
Lee
,
J. C. Y.
, and
Orszag
,
S. A.
,
1997
, “
Numerical Simulation of Low Mach Number Reactive Flows
,”
J. Sci. Comput.
,
12
(
2
), pp.
139
167
.10.1023/A:1025669715376
12.
Majda
,
A.
, and
Sethian
,
J.
,
1985
, “
The Derivation and Numerical Solution of the Equations for Zero Mach Number Combustion
,”
Combust. Sci. Technol.
,
42
(
3–4
), pp.
185
205
.10.1080/00102208508960376
13.
Capecelatro
,
J.
, and
Desjardins
,
O.
,
2013
, “
An Euler-Lagrange Strategy for Simulating Particle-Laden Flows
,”
J. Comput. Phys.
,
238
, pp.
1
31
.10.1016/j.jcp.2012.12.015
14.
Ling
,
Y.
,
Balachandar
,
S.
, and
Parmar
,
M.
,
2016
, “
Inter-Phase Heat Transfer and Energy Coupling in Turbulent Dispersed Multiphase Flows
,”
Phys. Fluids
,
28
(
3
), p.
033304
.10.1063/1.4942184
15.
Zwick
,
D.
, and
Balachandar
,
S.
,
2020
, “
A Scalable Euler-Lagrange Approach for Multiphase Flow Simulation on Spectral Elements
,”
Int. J. High Performance Comput. Appl.
,
34
(
3
), pp.
316
339
.10.1177/1094342019867756
16.
Wilke
,
C. R.
,
1950
, “
A Viscosity Equation for Gas Mixtures
,”
J. Chem. Phys.
,
18
(
4
), pp.
517
519
.10.1063/1.1747673
17.
Alkandry
,
H.
,
Boyd
,
I.
, and
Martinm
,
A.
,
2013
, “
Comparison of Models for Mixture Transport Properties for Numerical Simulations of Ablative Heat-Shields
,”
AIAA
Paper No. 2013-0303. 10.2514/6.2013-0303
18.
Reitz
,
R. D.
,
1987
, “
Mechanisms of Atomization Processes in High-Pressure Vaporizing Sprays
,”
Atomization Spray Technol.
,
3
, pp.
309
337
.https://uwmadison.app.box.com/v/AandS
19.
Abramzon
,
B.
, and
Sirignano
,
W. A.
,
1989
, “
Droplet Vaporization Model for Spray Combustion Calculations
,”
Int. J. Heat Mass Transfer
,
32
(
9
), pp.
1605
1618
.10.1016/0017-9310(89)90043-4
20.
Yuen
,
M. C.
, and
Chen
,
L. W.
,
1976
, “
On Drag of Evaporating Liquid Droplets
,”
Combust. Sci. Technol.
,
14
(
4–6
), pp.
147
154
.10.1080/00102207608547524
21.
Liu
,
A. B.
, and
Reitz
,
R. D.
,
1993
, “
Mechanisms of Air-Assisted Liquid Atomization
,”
Atomization Sprays
,
3
(
1
), pp.
55
75
.10.1615/AtomizSpr.v3.i1.30
22.
O'Rourke
,
P. J.
, and
Amsden
,
A. A.
,
1987
, “
The TAB Method for Numerical Calculation of Spray Droplet Breakup
,”
SAE
Paper
No. 872089
.10.4271/872089
23.
Xin
,
J.
,
Ricart
,
L.
, and
Reitz
,
R. D.
,
1998
, “
Computer Modeling of Diesel Spray Atomization and Combustion
,”
Combust. Sci. Technol.
,
137
(
1–6
), pp.
171
194
.10.1080/00102209808952050
24.
Fischer
,
P. F.
,
Lottes
,
J. W.
, and
Kerkemeier
,
S. G.
,
2019
, “
Nek5000 Version 19.0
,” accessed Apr. 1, 2021, https://nek5000.mcs.anl.gov
25.
Patel
,
S.
,
Fischer
,
P.
,
Min
,
M.
, and
Tomboulides
,
A.
,
2019
, “
A Characteristic-Based Spectral Element Method for Moving-Domain Problems
,”
J. Sci. Comput.
,
79
(
1
), pp.
564
592
.10.1007/s10915-018-0876-6
26.
Cohen
,
S. D.
,
Hindmarsh
,
A. C.
, and
Dubois
,
P. F.
,
1996
, “
CVODE, a Stiff/Nonstiff ODE Solver in C
,”
Comput. Phys.
,
10
(
2
), pp.
138
143
.10.1063/1.4822377
27.
Hindmarsh
,
A. C.
,
Brown
,
P. N.
,
Grant
,
K. E.
,
Lee
,
S. L.
,
Serban
,
R.
,
Shumaker
,
D. E.
, and
Woodward
,
C. S.
,
2005
, “
SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers
,”
ACM Trans. Math. Software (TOMS)
,
31
(
3
), pp.
363
396
.10.1145/1089014.1089020
28.
Fischer
,
P.
, and
Mullen
,
J.
,
2001
, “
Filter-Based Stabilization of Spectral Element Methods
,”
C. R. de L'Academie Des Sci. Ser. I Math.
,
332
(
3
), pp.
265
270
.10.1016/S0764-4442(00)01763-8
29.
Stolz
,
S.
,
Adams
,
N. A.
, and
Kleiser
,
L.
,
2001
, “
An Approximate Deconvolution Model for Large-Eddy Simulation With Application to Incompressible Wall-Bounded Flows
,”
Phys. Fluids
,
13
(
4
), pp.
997
1015
.10.1063/1.1350896
30.
Schlatter
,
P.
,
Stolz
,
S.
, and
Kleiser
,
L.
,
2006
, “
Analysis of the SGS Energy Budget for Deconvolution- and Relaxation-Based Models in Channel Flow
,”
Direct Large-Eddy Simul. VI
, 10(
1
), pp.
135
142
. 10.1007/978-1-4020-5152-2
31.
Negi
,
P. S.
,
Schlatter
,
P.
, and
Henningson
,
D. S.
,
2017
, “
A Re-Examination of Filter-Based Stabilization for Spectral-Element Methods
,” Report No. 1987.
32.
Zwick
,
D.
,
2019
, “
ppiclF: A Parallel Particle-In-Cell Library in Fortran
,”
J. Open Source Software
,
4
(
37
), p.
1400
.10.21105/joss.01400
33.
Offermans
,
N.
,
2017
, “
Gather-Scatter Library in Nek5000: Documentation of the GS Library Developed by James Lottes
,” Report No. 2.
34.
Sandia National Laboratories
,
2021
, “
Engine Combustion Network
,” Sandia National Laboratories, Albuquerque, NM, accessed Apr. 1, 2021, https://ecn.sandia.gov
35.
Arienti
,
M.
,
Wenzel
,
E. A.
,
Sforzo
,
B. A.
, and
Powell
,
C. F.
,
2021
, “
Effects of Detailed Geometry and Real Fluid Thermodynamics on Spray G Atomization
,”
Proc. Combust. Inst.
,
38
(
2
), pp.
3277
3285
.10.1016/j.proci.2020.06.039
36.
Wenzel
,
E. A.
, and
Arienti
,
M.
,
2021
, “
A New Approach for the Modeling and Simulation of Liquid/Vapor Phase Change at Engine-Relevant Conditions
,”
ILASS-Americas 31st Annual Conference on Liquid Atomization and Spray Systems
, Virtual, May 16–19, pp.
1
11
.https://www.osti.gov/servlets/purl/1862773
37.
Dong
,
S.
,
Karniadakis
,
G. E.
, and
Chryssostomidis
,
C.
,
2014
, “
A Robust and Accurate Outflow Boundary Condition for Incompressible Flow Simulations on Severely-Truncated Unbounded Domains
,”
J. Comput. Phys.
,
261
, pp.
83
105
.10.1016/j.jcp.2013.12.042
38.
Hwang
,
J.
,
Weiss
,
L.
,
Karathanassis
,
I. K.
,
Koukouvinis
,
P.
,
Pickett
,
L. M.
, and
Skeen
,
S. A.
,
2020
, “
Spatio-Temporal Identification of Plume Dynamics by 3D Computed Tomography Using Engine Combustion Network Spray G Injector and Various Fuels
,”
Fuel
,
280
, p.
118359
.10.1016/j.fuel.2020.118359
39.
Li
,
H.
,
Rutland
,
C. J.
,
Hernández Pérez
,
F. E.
, and
Im
,
H. G.
,
2021
, “
Large-Eddy Spray Simulation Under Direct-Injection Spark-Ignition Engine-Like Conditions With an Integrated Atomization/Breakup Model
,”
Int. J. Engine Res.
,
22
(
3
), pp.
731
754
.10.1177/1468087419881867
40.
Payri
,
R.
,
Gimeno
,
J.
,
Marti-Aldaravi
,
P.
, and
Vaquerizo
,
D.
,
2015
, “
Momentum Flux Measurements on an ECN GDi Injector
,”
SAE
Technical Paper No. 2015-01-1893.10.4271/2015-01-1893
41.
Senecal
,
P. K.
,
Pomraning
,
E.
,
Richards
,
K. J.
, and
Som
,
S.
,
2015
, “
An Investigation of Grid Convergence for Spray Simulations Using an LES Turbulence Model
,”
SAE
Paper No. 2015-01-0768, Vol.
2
.10.4271/2013-01-1083
42.
Parrish
,
S.
, “
Gasoline Spray (Spray G) Drop Size Measurements
,” ECN 3 Workshop Proceedings, Engine Combustion Network, accessed Dec. 16, 2022, https://ecn.sandia.gov/ecn-workshop/ecn3-proceedings/
43.
Paredi
,
D.
,
Lucchini
,
T.
,
D'Errico
,
G.
,
Onorati
,
A.
,
Pickett
,
L.
, and
Lacey
,
J.
,
2020
, “
Validation of a Comprehensive Computational Fluid Dynamics Methodology to Predict the Direct Injection Process of Gasoline Sprays Using Spray G Experimental Data
,”
Int. J. Engine Res.
,
21
(
1
), pp.
199
216
.10.1177/1468087419868020