Effect of Particle Diameter on Primary Breakup of High-Pressure Diesel Spray Atomization: A Study Based on Numerical Simulations Using the Eulerian–Lagrangian Model

Abstract

The coupling of Eulerian and Lagrangian methods in the Eulerian–Lagrangian Spray Atomization (ELSA) approach is critical. This study proposes an equation for the primary breakup particle diameter D of a diesel fuel spray and adopts it as a key transition criterion for coupling. A three-dimensional diesel spray is modeled by the large-eddy simulation (LES) approach. This improved ELSA simulation was conducted using various transition criteria for particle diameter D_cr. The results show that fuel spray experiences two stages: stage I, when a liquid column appears without a dispersed phase, and stage II, when primary breakup occurs with many discrete particles. Although D_cr has little influence on the macro-spray characteristics, such as top penetration distance S and spray cone angle θ, it has significant effects on discrete particles, such as their number, average diameter, distribution and location, and spray cone area. D_cr should be determined on the basis of actual operating conditions.

1. Introduction

Diesel vehicles, such as heavy-duty trucks, commercial vehicles, and engineering machinery, have been widely applied in all professions and trades. The diesel engine is one of the most important power devices because of its unique advantages, such as high thermal efficiency, powerful output, and good reliability. However, an improvement in the fuel injection performance of diesel engines is necessary to meet more stringent emission regulations. Fuel injection and spray features have a direct impact on air–fuel mixing, subsequent combustion process, and, thus, emissions (Yang et al. and Rakopoulos et al.) [1,2].

Many researchers have conducted studies focusing on fuel spray characteristics using optical experiments and numerical simulations to obtain a better understanding of fuel injection penetration and spray processes. Li et al. [3] used CONVERG software to simulate and analyse the fuel injection mixing and ignition process in a cylinder of a high-pressure direct injection diesel micro-ignition liquefied natural gas (LNG) engine and found that the timing of natural gas injection is critical for natural gas ignition, and the shorter the duration of natural gas injection, the more beneficial it is to the combustion of diesel or natural gas. Sorusbay et al. [4] used the KIVA-3v algorithm to investigate diesel spray with a Kelvin–Helmholtz/Rayleigh–Taylor spray break-up model. The experimental results verified the model, which estimated fuel tip penetration with good accuracy. Pei et al. [5] used a phase Doppler anemometry (PDA) system to investigate the free and impingement spray characteristics of fuel spray of various fuels, such as methanol, ethanol, isooctane toluene reference fuels, and gasoline, and reported that many droplets appeared during the spray process, and the Reynolds number with droplet diameter as the characteristic length reached 10,000. Gad et al. [6] measured diesel fuel spray characteristics using a digital camera with particle image velocimetry (PIV), and the images revealed that liquid fuel bulks appeared in the spray, as the spray was not stabilized or formed. Si et al. [7] optically observed a near-field diesel spray and found that it contained liquid, large droplets, and numerous tiny droplets. The macroscopic fuel spray characteristics revealed that liquid diesel fuel contains a continuous liquid phase and liquid droplets. In diesel fuel injection, the liquid column first jets out of the nozzle and then transforms into discrete droplets. As the fuel spray develops further, the liquid column breaks up into numerous droplets. Therefore, both a continuous liquid jet and dispersed droplets existed in the fuel spray. Moreover, most of the liquid columns and films are present in the near field, but the liquids in the near field are quite different from those in the far field.

The microcosmic characteristics of fuel sprays have been investigated to determine the details of the fuel spray process, especially the primary breakup of the fuel near the nozzle. Crua et al. [8] successfully captured breakup images, including those of the first emerging liquid fuel, which had the shape of an undisrupted oblate spheroid and dispersed tiny bubbles and droplets in the far field. Lai et al. [9] conducted an optical test of diesel fuel injection and reported that, at the very beginning of a diesel spray, the first jets out of the nozzle are a liquid column, which is a continuous liquid jet. Numerous droplets appeared in the spray owing to breakup; however, all the details of the spray could not be ascertained through optical methods alone. Some of the processes involved in phenomena, such as primary atomization and nozzle cavitation, are still not fully understood, despite the fact that the full set of equations describing the physical situation is known. Thus, spray computational fluid dynamics (CFD) has become one of the most important tools for understanding and improving diesel spray development in modern diesel engines. M. Ghiji [10] et al. used LES (large-eddy simulation) with VOF model and successfully simulated the process of fragmentation development of liquid column in near field. Yin [11] studied the cavitation phenomenon inside the nozzle and found that the cavitation flow can increase the effective velocity at the nozzle exit and that the selection of the cavitation model is important for the influence of the near-field liquid column fragmentation process in high-pressure jetting.

Recently, several studies have been devoted to the accurate modeling of the spray atomization phenomenon. One useful numerical model for fuel sprays is the Eulerian–Lagrangian Spray Atomization (ELSA) model, which can be used for modeling the discrete phase in high-pressure injection sprays and can realistically describe the dense zone of the spray and its atomisation. Vallet et al. [12] proposed an Eulerian model for the atomization of a liquid jet using direct numerical simulation (DNS) models based on fixed modeling parameters and constants. They reported that the ELSA model could capture the influence of changes in gas temperature. The ELSA method is an improved method based on the Eulerian–Lagrangian method, and the coupling of the Eulerian and Lagrangian methods is critical for the accuracy of the ELSA model. Many studies have been conducted on this coupling in the ELSA model by setting some specific parameters that switch from the Euler to the Euler–Lagrange approach. Blokkeel et al. [13] built the ELSA method for fuel spray and set the ratio between the average distance between droplets and their average diameter to complete the coupling between the Eulerian and Lagrangian phases. Beau et al. [14] realized the coupling by using the liquid volume fraction as a threshold value (1%), and the switch between the Eulerian and Lagrangian phases occurs if the liquid volume fraction is below 1%. However, the transition zone of this method may be disturbed. Anez et al. [15] conducted a simulation of an LES model on the basis of the ELSA framework and used the volume formulation with different meshes by considering a subgrid turbulent liquid flux that depends on the local flow condition. Yan et al. [16] proposed an automatic identification method for high-pressure spray droplet characteristics based on the SD-ELSA model, using sphericity (S) and particle size (D) as critical criteria and studied the initial atomization of a high-pressure common rail injector. The results show that this method has good prediction accuracy and realizes a coupling between Eulerian and Lagrangian approaches.

The above studies show that the ELSA model can capture important spray characteristics. Comparisons of the numerical results from the ELSA model with experimental measurements show encouraging agreement, and the ELSA model is suitable in the near-injector region. Currently, in most studies on the ELSA model, the Eulerian and Lagrangian models are calculated on the basis of the same grids but different space partition meshes. This calculation method requires every calculation condition to be revised to identify the Eulerian continuous phase interface (He et al.) [17]. The coupling of the Eulerian and Lagrangian methods is critical for the accuracy of the ELSA model. Many studies have been conducted on the coupling of the ELSA model by setting up certain specific parameters that switch from the Euler to the Euler–Lagrange approach. Ning et al. [18] developed a high-pressure diesel spray and atomization model on the basis of the ELSA method and reported that the droplet size, droplet number, and droplet distributions were determined using the local liquid mass fraction and local liquid surface density. Thus, diesel fuel droplet size is a key parameter for fuel spray and atomization. Furthermore, Dasappa et al. [19] measured the soot particle mobility size in premixed hydrocarbon flames and reported that the first inception of particles and persistence of nucleation-sized particles with time depend on the structure of the parent fuel. Therefore, the molecular structure and size of the fuel have important effects on fuel spray, combustion flame, and emissions.

The aim of this study was to investigate the effect of fuel particle diameter on a fuel spray model based on the ELSA method. Both experimental and numerical studies were carried out to better understand the primary break-up of diesel sprays. An optical test rig consisting of a high-pressure common-rail injector was designed and tested for fuel spray characteristics. Furthermore, on the basis of the ELSA method, a numerical method was developed to simulate high-pressure diesel sprays, with a greater focus on the primary break-up. The proposed numerical method uses two key transition criteria: sphericity (Sp) and particle diameter (D). An improved ELSA model, validated on the basis of optical test results, was adopted to numerically simulate different transition criteria and particle diameter conditions. The primary break-up of the high-pressure diesel spray was analyzed on the basis of the simulation results of the improved ELSA model.

2. Primary Breakup Analysis

2.1. ELSA Method

In the Eulerian–Lagrangian method, discrete droplets are tracked through the domain. The particle tracker uses the physical properties of the individual droplets to account for the exchange of mass, momentum, heat, species, etc., with the continuous phase. However, ignoring the gas volume displacement of droplets in spray-dense regions can affect the accuracy of the simulation.

In the volume-of-fluid (VOF) approach based on the Eulerian model, the volume fraction of the liquid is stored in each cell. The gas–liquid interface can be tracked using an explicit discretization scheme. When using the VOF model, the mesh of the phase interface around each droplet must be smaller than that of the smallest droplet, resulting in a higher mesh resolution and smaller time steps. However, this approach is expensive and requires significant computational resources.

The VOF-to-DPM model transition mechanism combines both approaches in a reasonable way. The initial jet and primary break-up are predicted using the VOF model on a sufficiently fine grid, and the bulk dilution region of the dispersed spray is simulated using DPM. The VOF-to-DPM transition is evaluated by specifying the particle diameter. If the compact liquid clumps separated from the liquid core meet the specified conversion criteria, they are converted into particle parcels in the Lagrangian formulation. The equations of motion of the particles are given below.

ANSYS Fluent predicts the trajectory of a discrete phase particle (or droplet or bubble) by integrating the net force on the particle, which is expressed in a Lagrangian reference frame. The net force is obtained by equating particle inertia with the forces acting on the particle and can be written as

$$m_p\frac{d{\overrightarrow{u}}_p}{d_t}=m_p\frac{\overrightarrow{u}-{\overrightarrow{u}}_p}{{\tau }_r}+m_p\frac{\overrightarrow{g}\left({\rho }_{p-\rho }\right)}{{\rho }_p}4+\overrightarrow{F}$$

(1)

where $m_p$ is the particle mass; $\overrightarrow{u}$ is the fluid phase velocity; $\overrightarrow{u_p}$ is the particle velocity; $\rho$ is the fluid density; ${\rho }_p$ is the density of the particle; $\overrightarrow{F}$ is an additional force; $m_p\frac{\overrightarrow{u}-{\overrightarrow{u}}_p}{{\tau }_r}$ is the drag force; and ${\tau }_r$ is the droplet or particle relaxation time calculated by:

$${\tau }_r=\frac{{\rho }_p{d}^2{}_p}{18\mu }\frac{24}{C_{d}Re}$$

(2)

Here, $\mu$ is the molecular viscosity of the fluid; $d_p$ is the particle diameter; and Re is the relative Reynolds number, which is defined as:

$$R_e\equiv \frac{\rho d_p|\overrightarrow{u_p}-\overrightarrow{u}|}{\mu }$$

(3)

The trajectory equations and any auxiliary equations describing the heat or mass transfer to or from the particle can be solved by stepwise integration over discrete time steps. The integration of time in Equation (1) yields the velocity of the particle at each point along the trajectory, with the trajectory itself predicted by

$$\frac{d_x}{d_t}{=u}_p$$

(4)

Equations (1) and (4) are a set of coupled ordinary differential equations, and Equation (1) can be cast into the following general form:

$$\frac{{\mathit{du}}_p}{d_t}{=u}_p\left(u-u_p\right)+a$$

(5)

where the term a includes accelerations due to all forces except drag force.

This set can be solved for constants u, a, and ${\tau }_p$ by analytical integration. For the particle velocity at the new location $u_p{}^{n+1}$ we obtain:

$$u_p{}^{n+1}{=u}^n+e^{-\frac{\Delta t}{{\tau }_p}}\left(u_p{}^n-u^n\right)-a{\tau }_p\left(e^{-\frac{\Delta t}{{\tau }_p}}-1\right)$$

(6)

2.2. LES Model

In high-pressure diesel injection, fuel injection is a transient flow with an unstable turbulent phase interface, which results in strong anisotropy throughout the injection break-up process. In this study, the large-eddy simulation (LES) method was used as a turbulence model. The Eulerian gas phase is described by the compressible Navier–Stokes equations. For an incompressible fluid, the basic equations of LES are obtained by spatially filtering the continuity, momentum, species, and energy equations as follows:

$$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial {\overline{u}}_i{\overline{u}}_j}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}+\nu \frac{{\partial }^2{\overline{u}}_i}{\partial x_j\partial x_j}+\frac{\partial {\overline{\tau }}_i{}_j}{\partial x_j}$$

(7)

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(8)

Here, $u_i$ and $u_j$ are the velocities in i and j directions, respectively; p is the pressure; ρ is the density; ν is the kinetic viscosity; and ${\overline{u}}_i$ is the filtered velocity. The superscript indicates filtering.

${\tau }_{ij}=\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j$ is defined as the subgrid Reynolds stress (SGS). To realize the LES model simulation, unclosed terms can be constructed according to the SGS. With the definition of $u_i={\overline{u}}_i+u_i{}^{\prime }$, SGS τ_ij is derived on the basis of the eddy viscosity turbulence model, as shown in Equation (9):

$${\tau }_{ij}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}=2{\nu }_T{\overline{S}}_{ij}$$

(9)

where ${\nu }_T$ is the SGS eddy viscosity coefficient; ${\delta }_{ij}$ is the Kronecker symbol; ${\delta }_{ij}=\{\begin{array}{c}0 \left(i\ne j\right)\\ 1 \left(i=j\right)\end{array}$; and ${\overline{S}}_{ij}$ is the rate of the deformation tensor, ${\overline{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$.

The model is closed as shown in Equation (10):

$$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial }{\partial x_j}\left({\overline{u}}_i{\overline{u}}_j\right)=-\frac{1}{\rho }\frac{\partial \left(\overline{p}+\rho \frac{{\delta }_{ij}}{3}\overline{u_k^{\prime }{u}_k^{\prime }}\right)}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\nu +{\nu }_T\right)\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)\right]$$

(10)

For single phase, the volume fraction is adopted in the continuity equation:

$$\nabla \left({\alpha }_l{\rho }_l\overline{u}\right)=-\frac{\partial }{\partial t}\left({\alpha }_l{\rho }_l\right)+S_M$$

(11)

where ${\rho }_l$ is the density of the gas phase; ${\alpha }_l$ is the vapor volume fraction; and S_M is the source term of the mass.

The energy equation is shown as follows:

$$\frac{\partial }{\partial t}\left(\rho e\right)+\nabla \left(\rho e\overline{u}\right)=\nabla \left(k\nabla T\right)+\nabla \left(p\overline{u}\right)+S_E$$

(12)

with

$$e=\frac{{\alpha }_l{\rho }_l{e}_l+{\alpha }_v{\rho }_v{e}_v}{{\alpha }_l{\rho }_l+{\alpha }_v{\rho }_v}$$

(13)

where k is the heat conductivity coefficient; e is the internal energy; S_E is the source item of energy; ${\alpha }_v$ is the volume fraction of the liquid phase; and ${\rho }_v$ is the density of the liquid phase.

2.3. Spray Model

For the spray model, the ELSA model takes advantage of the Eulerian description of near-nozzle flow. For each phase of the two-phase flow near the nozzle zone, the Eulerian method is employed to solve the mass, momentum, and energy equations. The governing equations are provided in the literature (Yousefifard et al.) [20]. For dispersed particles, the Lagrangian method is used to solve the mass, momentum, and energy equations. The breakup model used was the KH–Rayleigh–Taylor (RT) model.

2.4. Transition of the Phase Interface

During the primary breakup process, two zones can be simulated using the ELSA model: the liquid column zone and droplet zone. In the context of the ELSA model for two-phase flows, an Eulerian approach was used to describe the gas phase, whereas a Lagrangian approach was used for the liquid phase. During the total fuel spray and atomization process, both Eulerian and Lagrangian processes exist, with no interaction between the Lagrangian liquid phase and Eulerian vapor phase. For the ELSA model simulation, grid partitioning was first performed to define the liquid phase domain, whereas the disperse phase domain was used for the Lagrangian model. To realize effective convergence between Eulerian and Lagrangian approaches, an improved ELSA model is developed in this study, which employs two key parameters based on the transient spray process: the sphericity (Sp) of the breakup droplet and average stochastic particle diameter (D). These two parameters were adopted as the transition criteria between the Eulerian and Lagrangian domains. The sphericity (Sp) was used to determine whether the phase interface was a sphere. The model has greater accuracy for the particle transition to the dispersed phase at the spherical phase interface. The average particle diameter (D) was used to determine the size of the dispersed particles, which were ready for the secondary breakup model. Furthermore, the sphericity of the particles is related to the particle diameter, which directly determines sphericity. Thus, this study focused mainly on the particle diameter, and the sphericity was set to a constant.

Average Particle Diameter D

Currently, all the equations for the average particle diameter (D) are deduced on the basis of the disturbance breakup theory of the liquid volume surface wave. Classic liquid droplet breakup models include Taylor analogy breakup (TAB), WAVE, and KH-RT models. In the TAB model, a liquid droplet deformation equation is built by modeling the vibration of the droplet on the basis of the vibration of a spring, and the breakup is determined according to the deformation. In both the WAVE and KH-RT models, the droplet is considered to peel off the spray surface owing to the enhancement of the unsteady surface wave.

In contrast to the classical model of droplet rupture, an analysis of the forces acting on the droplet during the rupture of a nonstationary liquid film is presented here. The liquid first breaks through the ruptured membrane. At this time scale, when the air resistance is greater than the surface tension, the first large droplet formation is caused by the instability of the liquid film rupture, as shown in the following equation:

$$p_l-p_{gj}-2{\mu }_l\frac{\partial v_l}{\partial y}\pm {\sigma }_l\frac{{\partial }^2{\xi }_j}{\partial x^2}<0$$

(14)

where p_l and p_gj are the liquid pressure and gas pressure owing to the disturbance, respectively; v_l is the tangential velocity along the liquid jet direction; and x is the displacement of the liquid film along the jet direction.

The amplitude of the surface wave increases with time and film expansion and is described as follows:

$$\xi ={\xi }_{0}exp\left(\varpi t\right)$$

(15)

Droplet diameter D is derived from a static liquid film under the action of gravity, as described in Equation (16):

$$D=3.3\sqrt{\frac{{\sigma }_l}{{\rho }_{l}g}}$$

(16)

where σ_l is the liquid surface tension; ρ_l is the liquid density; and g is the acceleration due to gravity.

Only gravity is considered in Equation (16). However, for small droplets, air resistance plays a significant role in liquid fuel breakup, as shown in Figure 1. Considering the air resistance on the spherical droplet along with gravity, the droplet diameter in Equation (16) should be changed to

$$D=3.3\sqrt{\frac{{\sigma }_l}{{\rho }_l{a}_r}}$$

(17)

where a_r is the acceleration induced by the resultant force, F_r.

The resultant force, F_r, acting on the liquid droplet surface is:

$$F_r=F_f-mg=m{a}_r$$

(18)

where F_f is the air resistance force, which can be deduced using Equation (19)

$$F_f=\frac{1}{2}C{\rho }_{g}A{v}^2$$

(19)

where C denotes the air resistance force coefficient; ρ_g denotes the air density; and v denotes the relative fluid velocity. A is the windward area. Here, the particle is assumed to be a spherical droplet, and its windward area is equal to the projected area of the droplet surface, that is, A = πR². R is the radius of the spherical droplet.

The acceleration induced by the air resistance force is calculated using Equation (20)

$$a_f=\frac{F_f}{m_l}=\frac{F_f}{{\rho }_l{V}_p}$$

(20)

where m_l is the mass of the liquid. The particle volume is defined as follows:

$$V_p=\frac{4}{3}\pi R^3=\frac{4}{3}\pi {\left(\frac{A}{\pi }\right)}^{\frac{3}{2}}=\frac{4}{3}\pi {\left(\frac{A_p}{4\pi }\right)}^{\frac{3}{2}}$$

(21)

Substituting Equations (19) and (21) into Equation (20) yields the following:

$$a_f=\frac{3C{\rho }_g{v}^2}{8{\rho }_l\sqrt{\frac{A}{\pi }}}$$

(22)

When both the air resistance force and gravity are considered, Equation (17) is converted into Equation (23).

$$D=3.3\sqrt{\frac{{\sigma }_l}{{\rho }_l{a}_r}}=3.3\sqrt{\frac{{\sigma }_l}{{\rho }_l\left(a_f-g\right)}}=3.3\sqrt{\frac{{\sigma }_l}{{\rho }_l\left(\frac{3C{\rho }_g{v}^2}{8{\rho }_l\sqrt{\frac{A}{\pi }}}-g\right)}}$$

(23)

Thus, the average diameter of the liquid droplet can be estimated using Equation (23) under the action of air resistance force and gravity.

2.5. Model Calculation Process

The total field was calculated using the Eulerian model until the liquid phase entered the Lagrangian domain, where the dispersed phase was calculated using the Lagrangian model until the end of the simulation. For the ELSA model, the predefined grid partitions (defining the Eulerian and Lagrangian domains) were constants for each operational boundary condition. However, every time the operational boundary changed (e.g., when the inlet pressure or back pressure changed), the grid had to be rebuilt to create a new grid partition, rendering this method inconvenient.

In this study, an improved approach based on the ELSA model was adopted, and two key parameters, namely the sphericity S_pcr (here set as a constant, S_pcr = 0.5) and random particle mean diameter D_cr, were used as transition criteria between the Eulerian and Lagrangian domains. In the improved ELSA model, droplets in the primary and secondary crushing are transferred from the Eulerian liquid phase block to the dispersed phase particles, which achieves dynamic convergence between the Eulerian and Lagrangian equations. Figure 2 illustrates the simulation process of the improved ELSA model. The improved ELSA model allows for flexible identification of the phase interface using these two criteria. Similar to the ELSA simulation procedure, the model calculations were initially based on the Eulerian model. Once the dispersed particles satisfy these two criteria, the improved ELSA model uses the Lagrangian method to calculate the dispersed-phase domain. During fuel injection, the Lagrangian phase is transformed from the Eulerian phase. Furthermore, the Eulerian phase was actively transformed into the Lagrangian phase. During this active transformation, the position and timing of initial Lagrangian particles may change. Therefore, the transition from the Eulerian to the Lagrangian phase is transient rather than constant, as in the ELSA model. The improved ELSA method allows for dynamic identification of the phase interface, even if the operational boundary conditions change.

3. Numerical Model Building and Validation

3.1. Geometry and Mesh Modeling

The high-pressure diesel spray was tested using a constant volume pump (CVP) test rig under varying injection and back pressures, as shown in Figure 3a. The CVP consists of a cuboid with dimensions of 330 mm × 300 mm × 300 mm and a cylindrical core chamber. A high-pressure diesel multi-hole injector was located at the top center of the CVP. To realize single-hole injection, a matching nozzle cover was designed, which had a vertical hole with a diameter of 0.3 mm. A 3D model mesh based on the diesel injection prototype CVP system is shown in Figure 3b. In this CVP model, an orifice with a diameter of 0.3 mm and length-to-diameter ratio of 5 mm is present in the middle of the model. The zone around the cylindrical orifice entrance and outlet corners should be given the most attention owing to the complex flow caused by the variable cross-section and steep gradients. Therefore, the meshes in these zones were refined.

3.2. Physical Setup

3.2.1. Model Settings

We used the CFD software ICEM for the creation of the CVP model and meshing and solved the model using FLUENT2020R1 for the calculation. For the LES model, the WMLES S-Omega submodel was adopted. Although widely used in the academic community, LES has had a very limited impact on industrial simulations. The main reasons for this are: (1) The LES model requires a high-resolution wall boundary layer. (2) Nevertheless, near the wall, the largest scales in the turbulent spectrum are geometrically very small and require a very fine grid and small-time step. (3) This can be achieved only for flows with low Reynolds numbers.

The wall modeling LES (WMLES) S-Omega model compensates for the high requirements for wall boundary layers in the LES model as well as the need for very dense computational grids and small time steps. This model also does not require any significant increase in mesh resolution at high Reynolds numbers. Thus, it compensates for the shortcomings of the LES model.

The VOF model with an explicit setting was used in the near-nozzle zone. The DDM model was used in the far zone, and the KH-RT model was adopted for the breakup. We investigated the effect of the particle diameter D on the simulation results, which are presented here for only one injection pressure condition with varied transition criteria D_cr.

Table 1. Model physical parameters.

Parameter	Specification
Injection pressure (MPa)	80
Back pressure (MPa)	1
Nozzle diameter (mm)	0.3
Length/diameter ratio	5
Wall temperature (K)	298.15
Fluid	Diesel
Fluid temperature (K)	298.15
Fluid density (kg/m3)	850
Surface tension coefficient (N/m)	0.031
Air	Compressible gas
Air temperature (K)	298.15
Sphericity, Spcr	0.5
Average particle diameter, Dcr (mm)	0.05, 0.07, 0.1, 0.2

presents details of the boundaries and model parameters.

3.2.2. Choice of Solver

Two main solvers are included in FLUENT: the pressure-based solver and density-based solver. The pressure-based solver, which uses a pressure correction algorithm, solves the control equations in scalar form and can efficiently solve for incompressible flows but can also solve for compressible flows. This solver is primarily used for most conventional flow processes, where the effect on the fluid density is small. The density-based solver, which solves the control equations in vector form, is available in the main discrete formats, Roe and AUSM+. The solver adds a solution equation for density, which makes solving high-speed compressible flow processes possible.

For simulation of atomization breaking processes for high-speed incompressible diesel and compressible background gas, the pressure-based solver can be used to meet the computational requirements.

3.2.3. Setting of the Time Step

The time step and Courant number are mutually determined. The Courant number measures the stability of numerical calculations and is also known as the Courant–Friedrichs–Lewy (CFL) number or CFL criterion. The Courant number is the ratio of the distance of fluid motion in a time step to the unit length of the rectangular grid. The two-dimensional Courant number was solved using the following Equation:

$$\mathrm{Courant}=\frac{u\mathsf{\Delta }t}{\mathsf{\Delta }x}$$

(24)

where u is the fluid velocity; Δt is the time step; and $\mathsf{\Delta }x$ is the grid size.

The formula for the time step can be easily derived from Equation (24)

$$\mathsf{\Delta }t=\frac{\mathrm{Courant}\ast \mathsf{\Delta }x}{u}$$

(25)

In CFD, the length of a 3D mesh is characterized by the ratio of the mesh volume to the total surface area of the upper mesh:

$$\mathsf{\Delta }x=\frac{{\mathrm{Cell} \mathrm{Volume}}^{}}{\begin{array}{c}\mathrm{Total} \mathrm{Surface} \mathrm{Area}\\ \end{array}}=\frac{V}{A}$$

(26)

where V is the volume and A is the total surface area of the mesh.

In a three-dimensional mesh, the direction of velocity is not necessarily perpendicular to the face. As the velocity is at an angle to the face, we must consider the velocity component ($U_f·{\stackrel{\wedge }{n}}_f$) in the direction normal to the face. In this direction, the face points to the center of mass of the mesh. The velocity component is then multiplied by the time step Δt to obtain the distance traversed by the fluid. In summary, the three-dimensional grid Courant number was calculated as follows:

$$\mathrm{Courant}=\frac{1}{2}\ast \Delta t\ast \frac{{\displaystyle {\sum }_f\left|(U_f·{\stackrel{\wedge }{n}}_f)\right|}{{\displaystyle A}}_f}{{{\displaystyle V}}_p}$$

(27)

The formula for the time step can be derived from Equation (27)

$$\Delta t=\frac{{{\displaystyle 2V}}_p\mathrm{Courant}}{{\displaystyle {\sum }_f\left|(U_f·{\stackrel{\wedge }{n}}_f)\right|}{{\displaystyle A}}_f}$$

(28)

Because this solver is based on a pressure solver, the Courant number cannot be set directly. The flow velocity u and the minimum grid size can be roughly estimated, and the Courant number is usually set between 1 and 10. The Courant number in this study was set to 5 and the time step to 1 × 10⁻⁹.

3.2.4. Spatially Discrete Formats and Transient Formulation

The LES method is, to some extent, a direct simulation (DNS), and for time integration schemes, a Crankl–Nicoson semi-implicit scheme with at least second-order accuracy should be chosen. For spatially discrete formats based on the finite-volume method, a discrete format with at least second-order accuracy should be chosen to overcome pseudo-diffusion. Considering the limitations of computational time and resources, in this study, the bounded second-order implicit scheme was chosen for the time-integration scheme, and the second-order upwind format was chosen for the spatial discrete format.

3.3. Grid Independence Analysis

The mesh density is a very important parameter in simulation calculations. A grid with too small a density will not give accurate simulation results. A grid with very high density will yield a more accurate solution but will increase computational effort and time. The ELSA model was simulated with different mesh sizes. Figure 4 shows the results of the spray liquid column length for injection pressure P_injection = 80 MPa and back pressure p_b = 1 MPa. The results show that the liquid column length initially increased and then stabilized. During the stable stage, the length decreased slightly as the grid number increased, and this liquid column length curve became constant for grid numbers greater than 30 million, demonstrating that the model solution was relatively independent of the grid density. Hence, for better simulation accuracy, a model with 30 million grid points was acceptable in this simulation.

3.4. Model Validation

An optical diesel spray test rig was used to validate the numerical CVP model. The details of the test rig, test devices, and experimental procedures are presented in our previously published paper (Lei et al., 2019) [21]. In the CVP test rig, a high-pressure diesel injector was located on the center top head of the CVP system, and the fuel spray process of the injector was examined by Schlieren with a high-speed camera. The optical experimental results of the fuel spray are the macro parameters, such as the penetration distance S, spray cone angle θ, and fuel spray cover area A, as shown in Figure 5. To derive the effective spray macro parameters, we developed an image processing program on the basis of MATLAB, which has also been used in the literature (Lei et al., 2019) [22]. Both the optical test images and model simulation pictures were processed using this program, and the edges of all these images and pictures were extracted, as shown in Figure 5b. The top penetration distance S is defined as the vertical distance between the nozzle outlet center and top edge. The spray cone angle was the intersection angle of the penetration cone at 90% of S. The spray-side cover area A was the area of the spray contour edge.

Figure 6 shows a comparison of the model results with experimental data. The top penetration distance S is shown in Figure 6a for different transition standards (D_cr). All the S curves increased with time. The results show that all ELSA results are in good agreement with the optical experimental data for different transition standards (D_cr). The deviations in the models are shown in Figure 6b. The results show that, although the models have a relatively high bias in the early injection phase (time < 0.15 ms), their error is relatively low at injection times larger than 0.15 ms, not exceeding 10%. Furthermore, the error of the model is particularly small, less than 5%, in the range of 0.15 ms to 0.25 ms. Furthermore, there is little difference between the model results for different transition criteria (D_cr).

4. Results and Discussion

In this ELSA model of fuel injection, the transition standard particle diameter D_cr is an important parameter for the simulations. Four different transition standard particle diameters (D_cr) are presented; one value of D_cr is derived from Equation (11), i.e., D_cr = 0.1 mm, different values within a range of 0.1 mm, i.e., D_cr = 0.05 mm, 0.07 mm, and larger than the range of 0.1 mm, i.e., 0.20 mm.

Figure 7 shows the Reynolds number Re of the fuel spray flow at the nozzle outlet, which rapidly increases to its maximum value and then becomes stable. For high-pressure diesel fuel injection, Re is always high, with values of approximately 27,000, which shows that the fuel spray flow at the nozzle outlet is highly turbulent.

Figure 8 shows the Weber number We of the diesel fuel spray. The Weber number curves have similar trends to those of the Reynolds number curves. These results show that fuel injection has a rapidly increasing tendency at the beginning of the injection, which soon becomes stable. For this ELSA model, the stable Weber number with an injection pressure of 80 MPa is approximately 7.8 × 10⁵. Such high We will have a significant influence on diesel fuel spray. A high We means the shear tension is large, causing more deflection and deformation of the fluid flow, which is beneficial for breakup.

Furthermore, both Reynolds number and Weber number are almost constant but large, even when the transition criteria particle diameter, D_cr, varies. Such high Reynolds and Weber numbers illustrate that this high-pressure diesel fuel injection is a turbulent flow process with severe deformation of the liquid droplet.

The simulation results of the CVP model illustrate the diesel fuel spray process, as shown in Figure 9. Here, we present four curves with varied transition criteria for the particle diameter D_cr. All these top penetration distance S curves exhibit an increasing trend. Figure 9a clearly shows that diesel fuel spray penetration distance can be divided into two stages. Stage I: The early injection stage ranges from 0 ms to 0.1 ms. During stage I, all the curves coincide with each other. This is because during the early injection stage, the diesel fuel maintains the liquid, a liquid column appears, and the disperse phase does not occur. Stage II is the period after 0.1 ms. During stage II, primary breakup occurs, and many droplets appear. At 0.1 ms, the front of the liquid column breaks up, and the liquid column splits to form liquid sheets and large droplets surrounding the center liquid column. After that (>0.12 ms), more vapors appear, the spray becomes stable, and the diesel fuel continues to spray in a stable cone with the vapor surrounding the liquid center. During stage II, the differences among these S curves with varied transition criteria for particle diameter D_cr increase. This indicates that the transition criterion particle diameter D_cr may affect the simulation results. The differences, however, are minor, as shown in Figure 9b. Hence, the spray penetration distance is substantially stable for the varied transition criteria of particle diameter D_cr.

Figure 10 shows a comparison between the fuel jet shapes in the test image and model simulation results. During early stage I (e.g., at t = 0.08 ms), the model simulation results were somewhat different from the test images, although the penetration distances were very close to the test data. However, the fuel jet shape in the simulated results was similar to that in the images. For these simulation results, the center of the spray is a liquid column (red area, diesel liquid volume fraction ≈ 1), whereas the dispersed phase is diesel vapor (blue area, diesel liquid volume fraction ≈ 0). The optical image is only the outermost layer of the entire spray, whereas the surface of the fuel spray contains mainly fuel vapor. Therefore, the results of the ELSA model agree well with experimental data, not only in terms of penetration distance but also in terms of spray shape. The ELSA model of diesel spray can be used to analyze the macroscopic and microscopic properties of fuel sprays.

Figure 11 shows the details of the primary break-up of the fuel spray for different transition standard particle diameters, D_cr. At the beginning of the fuel spray process, the diesel spray was almost always a liquid column comprising a continuous phase. In addition, a thin film forms on the liquid column head in the shape of a mushroom. The liquid column then grows in length. At time = 0.08 ms, the liquid column film becomes unstable and starts to break up into large primary droplets at the edge of the film, while the liquid column shrinks. For all D_cr, droplets appeared at the edges. The transition model is then further broken up with continued jetting (from 0.1 ms), resulting in the transformation of the primary large droplets into dispersed particles. At the same time, the results show that the atomized particles not only penetrate the vertical spray direction but also expand in the horizontal direction, implying that the phase interface changes as the spray develops. The shape of the spray becomes long in the vertical direction and wide in the horizontal direction, resulting in a semi-conical profile. When the film front encounters the liquid column, the column breaks up, creating transition particles around the column. These results demonstrate that the improved ELSA model can effectively capture the dynamic phase interface. Furthermore, the entire development of the liquid column is similar for different transition standard particle diameters, D_cr, with differences in the droplet position and distribution.

Figure 12 shows the details of the spray cone angle, θ. For both the test data and model results, the spray cone angle rises to a maximum early on and then decreases. For different transition standard particle diameters, D_cr, the spray cone angles become different. Furthermore, the value of the transition standard particle diameter D_cr = 0.10 mm is closest to the test data.

Several planes perpendicular to the spray direction were defined to compare the model results in more detail, as shown in Figure 13. These four planes, numbered 1 to 4, were arranged along the vertical direction at distances of 15, 20, 25, and 30 mm from the nozzle outlet. These planes intersect the spray cone to form the intersection plane, which is projected onto the vertical plane to obtain the segments. The length L_s of the spray projection segments is defined as the maximum line distance.

The results for spray coverage area A are shown in Figure 14. This shows that the spray coverage area for all model results is slightly lower than the test data for A. Furthermore, for smaller transition standard particle diameters (D_cr = 0.05 mm), the difference in spray area A between the model results and test data becomes significant. As the transition standard particle diameter increases (D_cr ≥ 0.07 mm), the results of the model for spray area A match the test data. The smaller the transition standard particle diameter, the more the liquid column tends to shift towards the discrete phase, causing the spray cone to cover a smaller area. Therefore, the transition standard particle diameter is important for the simulation of the ELSA model, particularly for primary rupture.

Figure 15 shows the details of the discrete particles for the transition criteria particle diameter D_cr. Until time = 0.1 ms, no discrete particles appear. This means at the very beginning of the fuel spray (stage I, time < 0.1 ms), the spray remains liquid, and the primary breakup occurs at time = 0.1 ms, with an average particle diameter of approximately 0.045 mm (Figure 15b). Thenceforth (stage II, > 0.1ms), the discrete particle number increases sharply for all transition criteria particle diameter D_cr (Figure 15a); additionally, the average particle diameter rapidly decreases and eventually tends to a stable value (approximately 0.01mm). For too small a transition criterion particle diameter (i.e., for D_cr = 0.05 mm), both the discrete particle number and average particle diameter are smaller than those for other values of D_cr. The discrete particle number, in particular, is considerably smaller. This further explains why the spray cover area for D_cr = 0.05 mm is relatively smaller (Figure 14). For larger D_cr (i.e., for D_cr > 0.07 mm), the discrete particle number and diameter remain almost constant for different D_cr values.

Table 2. Average error.

	Average Error (%)
	Penetration Distance S		Spray Cone Angle θ		Spray Area A		Ls
	Total	Stage II	Total	Stage II	Total	Stage II
Dcr = 0.05 mm	10.91	7.07	24.79	19.17	38.61	33.33	16.84
Dcr = 0.07 mm	11.36	7.69	15.81	7.19	26.29	19.17	9.09
Dcr = 0.10 mm	8.60	4.33	9.56	1.63	28.91	22.57	4.21
Dcr = 0.20 mm	9.13	4.94	12.26	5.23	23.73	16.59	12.17

shows the average error in macro-spray parameters, including top penetration distance S, spray cone angle, spray area A, and L_s, for different transition criteria particle diameters. The results reveal that the model has the largest average error for the spray area compared with the test data, while it has the smallest error for the penetration distance S. In addition, the ELSA model has a relatively smaller average error during stage II compared with the total fuel spray process. This result indicates that the ELSA model is suitable for the simulation of the primary breakup of the early fuel injection process. Furthermore, D_cr = 0.10 mm had the lowest average error of all D_cr values.

Therefore, for the ELSA fuel injection model, the transition standard particle diameter D_cr is an important parameter. Although D_cr has little influence on macroscopic spray characteristics, such as the top penetration distance S and spray cone angle θ, it has a significant influence on the discrete particles in a single break-up. D_cr affects the number of discrete particles, average particle diameter, and distribution of particles, which ultimately determines the position of the particles and the area they cover (spray cone area). Thus, setting the correct D_cr is crucial for spray models and is related to spray parameters, such as the injector nozzle size, fuel properties, and fuel injection operating conditions (pressure and temperature).

5. Conclusions

In this study, an improved modeling method for diesel injection was proposed on the basis of the classical ELSA model. A three-dimensional model of high-pressure diesel injection was developed and validated on the basis of the experimental results obtained by optical tests in the CVP. The main findings of this study are summarized below.

The fuel spray had a two-stage feature during the early injection process. Stage I: At the very beginning of fuel injection, the diesel fuel maintains the liquid, and a liquid column appears, but the disperse phase does not occur. Stage II: The primary breakup occurs, many discrete particles appear, and the spray becomes a stable cone with vapor surrounding the liquid center.

Particle diameter is critical for the numerical model. The transition criterion particle diameter D_cr is an important parameter for the ELSA model. Although D_cr has little influence on the macro-spray characteristics, such as the top penetration distance S and spray cone angle θ, it significantly affects the discrete particles of the primary breakup. D_cr affects the discrete particle number, average particle diameter, and particle distribution, which determine the particle location and cover area (spray cone area).

Thus, setting a suitable D_cr that is related to injection parameters, such as the injector nozzle size, fuel properties, and fuel injection operation conditions (pressure and temperature), in the spray model is critical. D_cr should be set on the basis of the practical operating conditions of the fuel injection.