Direct Numerical Simulation of Atomization Characteristics of ECN Spray-G Injector: In-Nozzle Fluid Flow and Breakup Processes

Abstract

In order to sustain applications dealing with various liquid fuels in internal combustion engine (ICE), it is essential to make available prediction methodologies that allow an early evaluation of the potential usefulness of such fuels in terms of favourable mixture preparation process already in realistic configurations. Since the air-mixture formation and subsequent processes are predominantly governed by the fuel injection, a DNS based numerical investigation coupled with VOF as an interface tracking method is carried out in this paper to gain better insight on the fuel injection from an industrial injector "Spray G" configuration. Chosen from Engine Combustion Network (ECN), this is a gasoline direct injector (GDI) featuring 8-holes orifices and operating with high injection pressure (200 bar). Under consideration of the required computational cost associated with DNS, only the 1/8 of the nozzle geometry including one orifice is used. The numerical simulation is accomplished for the quasi-steady injection condition with nozzle needle fully opened. The obtained results are first validated with available experimental data for nozzle mass flow rate and spray spread angle showing a good agreement. Then, a detailed numerical analysis is provided for the in/near nozzle flow evolution especially for flow turbulence, primary and secondary atomization. Furthermore, droplet statistics in terms of droplet shape, and droplet size-velocity distribution together with a breakup regime map are reported. Finally, a 2-D data curation technique is proposed to extract the droplet statistics along selected planes and evaluated by direct comparison with three-dimensional droplet data, which may allow handling of the DNS data in more feasible and economical way especially for time series data with higher frequency. The comprehensive DNS data generated by this DNS-VOF approach enable not only to carry out detailed numerical analysis of in- and near-nozzle physical phenomena for which experimental data are still scarce, but also to provide a hint of more reliable injector boundary conditions useful for lower order spray injection method based on Lagrangian particle tracking.

1 Introduction

Notwithstanding the pressing interest in CO$_2$-free energy sources, liquid fuels will remain of great importance for transportation sectors as well as for power generation in stationary applications in the foreseeable future. In order to meet the increasingly stringent emission legislation, modern automobile internal combustion engines (ICE) that are powered by conventional hydrocarbon based fuel have to overcome the major challenge, namely, operating in a manner of best fuel efficiency with the lowest possible carbon footprint. Thereby, two of the most important strategies are, on one side, to replace the fossil fuels with renewable fuels, like bio-ethanol, bio-butanol, etc, on the other side, to improve the fuel efficiency with regard to develop in-cylinder technologies that allow a realization of engine with higher efficiency. These technologies essentially deal with the conducive charge-intake, preparation of optimal air-fuel mixture, proper combustion and mitigation/control of pollutant formation. Because the air-mixture formation and subsequent processes in ICE are predominantly governed by the fuel injection, it is essential to make an early evaluation of the potential usefulness of such fuels in terms of favourable mixture preparation process already in realistic configurations. For that purpose, the design of numerical methodologies that shall allow comprehensive analysis of the relevant processes close to the nozzle exit is highly required. In fact, for ICE operated with Gasoline Direct Injection (GDI) system, the fuel is directly injected in the cylinder chamber during the engine intake stroke, thus helping in spreading and formation of fuel spray in the combustion chamber in order to achieve an air-fuel mixture ideal for an optimum engine performance under varied loading conditions (Sadiki et al., 2016). During the injection and spray formation, the fuel spray absorbs the heat from the surrounding gas resulting evaporation and mixing of the fuel with air. Thus reducing the chamber temperature, which allows more air intake and thus increases the volumetric efficiency and reduces the chances of knocking. Based on engine loading conditions, the fuel is often injected for a shorter duration in the order of magnitude of milliseconds, and with fewer crank angle degree rotations. Consequently, the time scale for the entire process of injecting the liquid fuel and mixing it with the air to form a properly combustible mixture is very small. Accordingly, the issue of achieving a suitable air-fuel mixture becomes one of the most critical bottleneck for an efficient and clean combustion in such engines.

As stated above the quality of the air-fuel mixture largely depends upon the liquid fuel injection strategy (timing and mode). During a full load condition, the fuel is injected at the start of the suction stroke, and in part load conditions, it is injected during the start of the compression stroke, when the intake valve is closed. Thereby, the fuel injection has the direct influence on the engine operation and performance via the fuel disintegration (primary and secondary breakup) processes. This results in a dispersed spray which further undergoes complex liquid–gas interaction (flow drag, turbulence dispersion and modulation, evaporation along with air-mixture formation) and subsequent turbulent combustion process, radiation, and pollutant formation (soot, NOx, CO) (Sadiki et al. 2016; Neroorkar 2011; Serras-Pereira et al. 2010; Devassy et al. 2020). It is also reported that the fuel injection and the resulting spray dynamics are strongly correlated with the processes taking place upstream in- and near-nozzle flow (Serras-Pereira et al. 2010).

As Duronio pointed out in their study (Duronio et al. 2020), the process of spray breakup plays a crucial role in fuel injection into the combustion chamber. This phenomenon has been extensively investigated in the literature using both numerical models and experimental imaging techniques. Over the last few decades, numerical simulations have gained increasing relevance in providing a comprehensive understanding of the ongoing processes. However, the focus of the most works using DNS or highly resolved numerical simulations has been on the development of numerical methods that made it possible to monitor the jet breakup during the atomization process, e.g. (Demoulin et al. 2013; Desjardins et al. 2013, 2008; Gorokhovski and Herrmann 2008; Herrmann 2008, 2011; Jiao et al. 2017; Ménard et al. 2007; Olsson and Kreiss 2005; Sussman and Puckett 2000; Sussman 1994). Thereby, the interface capturing methods Volume of Fluid (VOF)(Hirt and Nichols 1981) and Level-Set (LS) (Sussman 1994) and their combinations (Ménard et al. 2007; Gorokhovski and Herrmann 2008) have emerged and are mainly applied due to their attractive properties in terms of accuracy and reliability. Both procedures introduce an indicator function, which does not directly indicate the phase boundary, but marks rather the fluids involved. This is in contrast to interface tracking methods, in which particles mark the phase boundary and the calculation of strong topology changes require particularly high effort (see Gorokhovski and Herrmann 2008). All these methods have been developed and applied to generic configurations. The first applications under technical conditions (higher pressure) of this combination, known as CLSVOF (Coupled LS/VOF) methodology, can be found in Arienti et al. (2010); Ménard et al. (2007). Rather, a more comprehensive DNS of the two-phase mixing layer between parallel gas and liquid streams is carried out by Ling et al. (2017, 2019) in quasi-3D domain by using mesh up to 4.23 billions control volumes with smallest grid size of 3.125 μm. In this study, mixing layer is generated by continuous and parallel streams of liquid and gas at different velocities. This way, it allows to carry out simulation at sufficient long physical time scale for detailed turbulence and statistical investigation of two-phase flow under study. Thereby correlating the turbulence quantities to the formation of surface instabilities, the breakup of liquid jet into sheets, then ligaments, and finally into droplets. The reported mesh sensitivity study highlighted the importance of mesh size in resolving correctly not only the smaller droplets but also the larger droplets. Moreover, the most complex DNS so far of primary breakup of the direct injection of a liquid jet into a still environment was performed by Shinjo and Umemura (2011, 2010). The authors coupled the LS with a modified VOF methodology and used 6 billion control volumes to address the impact of perturbations on the atomization of the liquid jet tip. No turbulence initialized by the injector on the liquid jet has been taken into account. Detailed analyzes of the physical mechanisms, liquid structures and droplet distributions were presented. Unlike Shinjo and Umemura, Bode et al. (2017) provided turbulent inlet data from a compressible LES of the in–nozzle flow to carry out DNS of the primary breakup in a direct injection of a 3-hole injector. They considered 400 million control volumes. The applied CLSVOF solvers are from Le Chenadec and Pitsch and carefully documented in Chenadec and Pitsch (2013a, 2013b). In order to specifically determine the influence of the flow within the injector on the atomization, Arienti and Sussman (2014) extended their CLSVOF method by an additional LS function for handling the moving walls. The additional LS-equation makes it possible to account for the injector geometry including the opening and closing of the valves, and thus to consider a direct coupling between the injector inner flow and the atomization (Arienti and Sussman 2017). Both publications reported a use of mesh counts of 350–400 million cells.

Furthermore, due to the piston cyclic motion, the in-cylinder processes are characterized by cyclic fluctuations. Apart from influencing several other transient phenomena such as engine characteristic flow (tumble, swirl), turbulence, ignition, flame propagation, combustion etc., the cyclic variability of the internal engine flow can affect the morphology of the injected fuel jet and spray evolution. While the fuel injection can also induce or influence the engine cyclic variability, such interactions have not been fully investigated and not yet well understood. To this purpose, only few works dealt with the effects of the injection spray on the flow and combustion characteristics in the cylinder. However, they rather used the VOF-LES method coupled directly (Herrmann 2008, 2010) or indirectly (Edelbauer 2017; Grosshans et al. 2014) to the Lagrangian droplet tracking technique. In highly resolved numerical simulation or DNS, the realization of the two–way interactions can be addressed comprehensively at least in three steps. First, investigations shall be carried out for quasi-steady nozzle flow condition in which the needle is not moved. This way, the effects of injection pressure and in–nozzle turbulent flow on the liquid and relevant spray dynamics can quantified adequately. This will also allow to generate reliable reference data for validation and injector boundary conditions for lower order methods (e.g. LES, LPT, etc.). A second step shall then examine the additional transient effect of motions of the injector needle valve on the fluctuations in the fuel density and early atomization as it emerges from the nozzle. Finally, the mutual interaction and influence of fuel injection and engine CCV shall be included to DNS study which will provide an extensive data set for both model development and validations. In fact, this is one of the main objectives of the research activities within the Deutsche Forschungsgemeinschaft (DFG, German Research Council) project FOR2687.

However, the present work focuses on the issue representing the first step mentioned above, in which only the effects of in–nozzle turbulent flow on the liquid injection and disintegration are considered, in order to show how DNS is able to provide early information of the characteristics of the fuel disintegration in complex configuration. With this in mind, an industrial gasoline injector (GDI) Spray G3 from the Engine Combustion Network (ECN) (ECN Workshop 2022) is chosen to investigate numerically using pseudo-DNS/VOF and to generate comprehensive reference data in the near nozzle region. In particular, droplet statistics and spray characteristics along with joint droplet size distribution–velocity at various times and positions will be provided including initial droplet size distributions as typically required in frequently used numerical approaches based on the Eulerian-Lagrangian (EL) method. It should be noted that in industrial GDI the fuel is injected for fraction of milliseconds ($\approx$0.8ms) with highly transient fuel mass injection rate. Therefore, it is not possible in this study to carry out simulation for sufficient longer physical time to achieve the so called statistically developed flow that can allow detailed turbulence analysis and its impact on overall breakup processes.

This paper is organized as follows. Section 2 presents the pseudo-DNS/VoF approach used for the description of the turbulent multiphase flows. Section 3 introduces the numerical setup along with the investigated ECN spray G3 test case. The achieved results are discussed in Sect. 4. Section 5 is devoted to summary and concluding remarks.

2 Coupled Direct Numerical Simulation/ Volume-of-Fluid Approach for Multiphase Flows with Cavitation

In this work, a direct numerical simulation (DNS) is proposed to carried out detailed numerical investigation of gasoline injection and resultant spray dynamics for a Spray G configuration listed in Engine Combustion Network (ECN) (ECN Workshop 2022). Thereby, the flow is considered as incompressible single continuum, which is calculated by the DNS approach, while the phases corresponding to the liquid and the gases are treated as one fluid with material properties that change abruptly across the interface. This interface is simultaneously resolved by using a surface capturing technique by applying the Volume of Fluid (VOF) approach. The variant of VOF approach descirbed in Ahmed et al. (2020); Yu et al. (2017) is implemented into OpenFOAM (https://www.openfoam.com/documentation/user-guide) that allows to track the liquid and gaseous phases with cavitation. Therefore, the influence of phase changes due to cavitation on the atomization could be captured in this context (Li et al. 2020)

In the context of coupled DNS/VoF method, it is assumed that the liquid and gaseous flows are mixed homogeneously as a single continuum in the solution domain. Thus, the balance equations of mass and momentum read

$$\begin{aligned} \begin{aligned} \frac{\partial \rho }{\partial t} +\nabla \cdot ({\rho }{{\textbf{U}}})&=0,\\ \frac{\partial \rho {\textbf{U}}}{\partial t} +\nabla \cdot \left( \rho \textbf{UU}\right)&= \nabla p + \nu \nabla \cdot \left( (\nabla {\textbf{U}})^T + \nabla {\textbf{U}} \right) +{\textbf{f}}, \end{aligned} \end{aligned}$$

(1)

where $\textbf{U}$ the velocity, p the pressure and $\nu$ the kinematic viscosity. While, ${\textbf{f}}$ represents the source term due to the surface tension between liquid phase and the gas phase. This force is only active at the location of the phase interface and calculated by the method of Villiers et al. (2004).

As a mixture property the mass density $\rho$ is computed based on the property of individual phases as

$$\begin{aligned} \rho = \alpha _l\rho _l + \alpha _v\rho _v + \alpha _g\rho _g. \end{aligned}$$

(2)

Hereby, $\rho _l$, $\rho _v$ and $\rho _g$ represent respectively the densities of liquid fuel, vaporous fuel and non-condensable gas phase, which does not take part in the phase change process. Whereas, $\alpha _l,\, \alpha _v$ and $\alpha _g$ are the volume fraction of liquid phase, vapor phase and non-condensable gas phase within the control volume. The equations describing the volume fractions are given as:

$$\begin{aligned} \begin{aligned} \frac{\partial \alpha _l}{\partial t} +\nabla \cdot (\alpha _l {\textbf{U}})&+ \nabla \cdot \left( \alpha _l (1-\alpha _g){\textbf{U}}^r_{lg} \right) + \nabla \cdot \left( \alpha _l (1-\alpha _v){\textbf{U}}_{lv}^{r} \right) \\&= {\dot{m}}\left( \alpha _l \left( \frac{1}{\rho _v}-\frac{1}{\rho _l} \right) +\frac{1}{\rho _l} \right) + \alpha _l \nabla \cdot {\textbf{U}}, \\ \frac{\partial \alpha _v}{\partial t} +\nabla \cdot (\alpha _v {\textbf{U}})&+ \nabla \cdot \left( \alpha _v (1-\alpha _l){\textbf{U}}^r_{lv} \right) \\&= -{\dot{m}}\left( \alpha _v \left( \frac{1}{\rho _l}-\frac{1}{\rho _v} \right) +\frac{1}{\rho _v} \right) + \alpha _v \nabla \cdot {\textbf{U}},\\ \alpha _l + \alpha _v + \alpha _g&= 1. \end{aligned} \end{aligned}$$

(3)

Thereby, ${\textbf{U}}^r_{lg}$ and ${\textbf{U}}^r_{lv}$ are the relative velocity between liquid-phase and gas-phases, namely, non-condensable gas phase and vapor phase, respectively, and the involved “artificial compression” terms in the transport equations (see Weller 2008) are added to preserve the sharpness of the liquid–gas interface. Notice that, in the present work, the surface tension of gas–gas phase is assumed to be negligible, thus, the sharping terms for the gas–gas interface is not presented in the above phase transport equations. In the present work, the mass transfer term ${\dot{m}}={\dot{m}}_{c}+{\dot{m}}_{v}$ due to cavitation is modeled by the Schnerr–Sauer model (Schnerr and Sauer 2001), where terms ${\dot{m}}_{c}$ and ${\dot{m}}_{v}$ represent the rate of condensation and evaporation liquid phase on phase interface, respectively, exprressed as follows:

$$\begin{aligned} \begin{aligned} {\dot{m}}_{c}=\frac{3 \rho _v \rho _l}{\rho }\alpha _{l}\alpha _{v}R_b\sqrt{{\frac{2}{3\rho _l \left( \left| p - p_v\right| +0.001 p_v\right) }}}max\left[ p-p_v, 0\right] \\ {\dot{m}}_{v}=\frac{3 \rho _v \rho _l}{\rho }\alpha _{l}\alpha _{v}R_b\sqrt{{\frac{2}{3\rho _l \left( \left| p - p_v\right| +0.001 p_v\right) }}}min\left[ p-p_v, 0\right] \end{aligned} \end{aligned}$$

(4)

the term $R_b$ represents a cavitation nuclei radius, related with $\alpha _l$and and bubble density n by

$$\begin{aligned} R_b = \left( \frac{4\pi n \alpha _l}{3\left( \alpha _v-\alpha _{Nuc}\right) }\right) ^{\frac{1}{3}} \quad \text{ with }\quad \alpha _{Nuc}=\frac{n\pi d_{Nuc}^3}{6+n\pi d_{Nuc}^3} \end{aligned}$$

(5)

where, $\alpha _{Nuc}$ is the nucleation site volume fraction corresponding to bubble diameter $d_{Nuc}$. Notice that, though the equation and source terms due to cavitation are included into the numerical framework and simulation, the detail discussion on flow cavitation and their impact on spray atomization has been reported by author (Li et al. 2020), but it is not the scope of this study.

3 Engine Combustion Network “Spray G” Configuration and Numerical Setup

To support detailed investigations of liquid jet injection using numerical simulations and relevant model development, the workshop series are regularly organized jointly by researcher from experiment and numerical simulation communities. These workshops are primarily dedicated to select the common operating parameters and physical phenomena to be investigated in the context of fuel injection and spray combustion. Most notably the workshops on Turbulent Combustion of Spray (TCS) (http://wwwaggutheil.iwr.uni-heidelberg.de/tcs/tcs8.html), mostly focusing on atmospheric and low–pressure spray burners and the ECN workshops, addressing high-pressure and dense sprays which operate under representative engine conditions. This way, it provides opportunity and a common platform to researcher from both experimental and numerical modeling groups to exchange the relevant data for collaborative comparisons of measured and/or modeled result. Thereby, identify priorities and key processes for further experimental and computational research, thus promoting an intense mutual interaction among themselves. In the present work, one of the gasoline injector configuration known as “Spray G” injector from the ECN is chosen, which was originally manufactured by Delphi Technologies.

3.1 Experimental Configuration

The Spray G injector represents an 8-hole conventional pressure assisted GDI (see Fig. 1 (left)), for which various operating points with well defined initial and boundary conditions are available as summarized in Table 1. In the present study, the case “Spray G” of iso-octane injection with a injection pressure P_inj= 200bar and temperature T_inj = 363K is selected, while the pressure and temperature in the chamber are kept at P_amb= 1bar and T_amb = 333K, respectively. The non-dimensional flow conditions of the Spray G3 case are listed in Table 2 with Reynolds number $Re_{nozzle} = U_{bulk}L/\nu _l$, Weber number $We_{nozzle}=\rho _l\left( U_{bulk}^{l} -U_{bulk}^{g}\right) ^2\,L/\sigma$ and Ohnesorge number $Oh_{nozzle}= \rho _l \nu _l/\sqrt{\rho _l\sigma L}$. Hereby, the quantities are calculated based the local value at the nozzle exit, e.g. the characteristic length L is the orifice diameter at nozzle exit.

Table 1 Spray G: various operating conditions

In order to make the available experimental data for numerical modeling and validation, various spray and injection characteristic data are uploaded at the ECN database, such as mass flow rate, spread angle, droplet size distribution etc. However, the fluid flow processes inside the nozzle and the near–nozzle atomization processes are highly complex and very difficult to carry out any measurements. The difficulty arises partly due to the inaccessibility of region and partly due to the flow is very dense. Thereby, not many experimental studies are reported and measurement data are still scarce for these regions. In this context, a high quality numerical investigation can be of good help for in depth analysis and understanding of in-nozzle and near-nozzle flow phenomena. The data obtained by such study can also be very useful as a reference data for other numerical investigations.

Table 2 Non-dimensional flow conditions of Spray G3 case

3.2 Numerical Configuration

As previously pointed out, in Spray G nozzle the gasoline is injected through 8 holes orifice with relatively high pressure resulting very high injection velocity. A high quality numerical investigation of such flow system essentially requires to resolve properly all the important flow structures (e.g. time and length scales) in the context of spray breakup processes. Therefore, a very fine spatial resolution together with ultra low time steps are indispensable, and hence the associated extremely high computational expenses are unavoidable. In order to ensure that all the relevant scales are properly resolved with feasible and economical computational costs, a 1/8 domain is considered in this study (see Fig. 1 (right)) with the non-dimensional flow conditions listed in Table2. Regarding the small computational domain and short injection duration as well as the focus of this study, to investigate the processes very close to nozzle, the heat transfer and associated mass transfer is not considered and the simulation of the fluid flow is treated as iso-thermal, whereas the physical properties of individual species are specified according to their initial temperature, namely liquid or vapor fuel at 363 K and nitrogen at 333K, respectively. Further, to investigate the fully evolved fuel spray and to circumvent the numerical complexity associated with needle movement, a injection scenario is considered with the nozzle needle fully opened.

Fig. 1

8-hole spray G configuration (left) and computational domain with operating conditions (right)

The resulting computational domain with relevant boundary conditions is depicted in Fig. 2 (left), in which the numerical domain for nozzle geometry is taken in accordance with ECN workshop. A cyclic boundary condition is applied for the circumferential direction, the no-slip condition is utilized at the walls, and the constant total pressure boundary condition is specified for both inlet and outlet with P_inlet=200 bar and P_outlet=1 bar. In order to resolve properly the dynamics of gas–liquid interface using VoF method, the mesh resolution is accomplished to resolve the sub-micron size atomized individual liquid mass (or droplets). Figure 2 (right) shows the mesh size distribution on the mid-section of the computational domain. The mesh density is very fine ($\approx$ 0.5–1.0 μm) in the injector orifice and entertainment region and with mesh size of 1.0–2.0 μm in the atomization region with gradual increment up to $\approx$ 50 μm in the outer region. This way, the numerical domain consists of a total cell number of 98 M. Throughout the simulation, maximal Courant-Friedrichs-Lewy number is set to CFL = 0.6, resulting in a time-step size of $\Delta t \approx 5\times 10^{-10}s$.

The simulation is carried out using a modified OpenFOAM (https://www.openfoam.com/documentation/user-guide) solver for multi–phase flow with cavitation. A detailed implementation can be found in Ahmed et al. (2020); Yu et al. (2017) and elsewhere. Thereby, second-order numerical schemes are employed for transient, convection and diffusion terms. All the computations are performed on the Hawk at High-Performance Computing Center Stuttgart (HLRS) requiring a total simulation time of $\sim$ 8mil. CPUHs. That constitutes an usage of 3072 processors running continuously for about four months.

Fig. 2

Numerical setup of the Spray G3 configuration, nozzle geometry with boundary conditions (left) and grid size distribution along the mid-sectional plane (right)

4 Results and Discussions

As already highlighted that the DNS of such a high velocity nozzle flow inherently involves very high computational expenses due to very fine mesh together with ultra low time step–size. Therefore, it would be highly computational intensive and very costly to obtain DNS results that are averaged for multiple injection events. In this work, DNS results from a single injection event is rather used especially for the comparison with experimental data and analysis of droplet statistics (e.g. droplet size distribution, droplet velocity etc.) as addressed in the following subsections.

4.1 Comparison with Experimental Data

In order to first verify the adopted pressure inlet and outlet boundary conditions in achieving the correct transient nozzle discharge rate, Fig. 3 shows the transient evolution of nozzle flow rate for both DNS and experiment data. It should be noted that the flow rate representing a quasi-steady needle lift position. The compared result shows very good agreement of DNS with experiment suggesting the reliability of the applied fixed inlet and outlet pressure boundary conditions. Noting that, the numerical results are presented with applying the axial symmetrical condition.

Fig. 3

Comparison of the transient evolution of nozzle mass flow rate during the injector needle fully opened phase

Fig. 4

Comparison of the transient evolution of spray angle during the injector needle fully opened phase

In fuel injection and spray combustion processes, the spray angle is of particular importance. It determine the spatial distribution of droplets and air-fuel mixture inside the combustion chamber. Apart from the injection parameters, the spray angle is greatly influenced by the in-cylinder flow properties (tumble, swirl, turbulent kinetic energy etc.). In this regard, Fig. 4 shows the spread angle of DNS results, which are evaluated at 1.2 mm downstream to the nozzle exit, that agrees well with experimental data obtained in further downstream. Thereby, the successful validation of DNS results with experimental measurement confirms that the adopted numerical approach is reliable and allows further DNS based numerical analysis of the spray G3 configuration under study.

4.2 Turbulent Nozzle Flow

Inside injector nozzle, a high speed fuel jet is generally made to pass through a complex flow passages with varying cross-section and flow directions. To substantiate this, an instantaneous velocity profile and fuel distribution with $\alpha _f \ge 0.5$ along the mid-sectional plane are shown in Fig. 5. The entrainment of the non-condensible gas from the combustion chamber into the nozzle hole, i.e. sudden change in flow direction, results in flow separation near the orifice base (see Fig. 5a marked in red). Additionally, the high speed fuel jet encounters sudden and adverse changes in flow direction at the orifice base (see Fig. 5a marked with green circles). These two combined phenomena are reckoned to expedite and enhance the overall primary atomization process.

Fig. 5

Flow and spray evolution along the mid-sectional plane a Instantaneous velocity, b fuel distribution (liquid attached to liquid core in red and detached liquid in blue), c closer view of vortex structure around the liquid interface in the region marked with red color d Closer view of vortex structure and liquid jet disintegration in the region marked with green color

In particular, the liquid jet is accelerated locally up to 300 m/s due to reduction of the effective cross-sectional area leading to local pressure drop in this region. This essentially leads to a flow cavitation which is not the scope of this study rather left for the future research. Additionally, the gas entrainment into the orifice (see the region marked in red) intensifies the surface instability by imparting additional energy to the liquid jet. Thus helping in intense and early breakup of liquid jet and during this processes the detached liquid mass can be seen with higher velocity. The corresponding profile for the liquid fraction is depicted in Fig. 5b,c. Thereby, the region of flow separation can be clearly identified (see the red marking) as there is almost no presence of liquid phase in this region. Additionally, due to the gas entrainment, the early breakup of the liquid jet can also be observed as the distinct pocket of liquid mass can be seen inside the nozzle orifice itself. Subsequently downstream, these liquid pockets further disintegrate into smaller one due to aerodynamic effects (see in Fig. 5d). Similar fuel distribution and phenomena are reported in other numerical studies, e.g. Chen and Oevermann (2018); Baldwin et al. (2016); Saha et al. (2016); Yue et al. (2020).

In order to analyze the turbulence evolution along the spray and the fuel jet interface, the vortex visualization technique based on the Q-criteria is adopted in this study, in particular the Q is calculated for the gas/carrier phase only. An global view of the vortex structures in the nozzle hole and near downstream in the chamber along the liquid jet is provided by Fig. 6, in which vortices are visualized by iso-surfaces of Q with a value of $Q=2\times 10 ^{14} s^{-2}$ and colored by the velocity magnitude. The spray G features a step-hole counter-bored orifices to allow the gas entertainment which is confirmed by the sparse distribution of vortices in the hole. In the near-downstream to the nozzle as also pointed out by Blume et al. (2019), the stream-wise directed vortex structures are observed around the jet core, while at the jet-gas interface, vortices with circumferential direction orientation are more dominant (see Section A in Fig. 6). Further downstream inside the chamber, the vortex structures are relatively smaller and denser. This can be attributed to the flow evolution during the intense spray breakup of ligaments in to droplets due to aerodynamic instabilities. Moreover, the gas flow also appears to be accelerated after the breakup along flow direction as shown in Fig. 6 (see also Fig. 5, top).

Fig. 6

Iso-surface of Q-criterion with $Q = 2\times 10^{14}\,s^{-2}$ in the gas phase colored by velocity magnitude

4.3 Primary and Secondary Breakup

As already pointed out based on earlier observations in Fig. 5 that the on-set of early instabilities are the combined influence of sudden change in flow direction, gas entrainment and subsequent increase in local jet velocity. The resulting three-dimensional morphology of the jet core together with droplets profile on selected planes are shown for liquid structures with volume fraction $\alpha _l \ge 0.5$ in Fig. 7. The selected planes are placed along the center line of the nozzle with a distance of 0.2 mm, where the first plane (P1) is located next to the interface between inner hole and counter-bore to the chamber along the centre line of the nozzle. Notice that, the blue marked region represents the intact structures that are still attached on the liquid core, while the separated and isolated liquid volumes are shown in red. Due to the special design of this spray G injector and sudden change in flow direction in inner hole, the early liquid jet tends to be asymmetric. Subsequently, the gas entrainment augments the surface instability resulting in the expansion of the liquid jet with the formation and breakup of ligaments along the counter-bore region (see also Fig. 5). Thereby, the primary breakup takes place in the counter-bore, and most liquid are still attached on the liquid core even at location P4 in the chamber, while the detached liquid structures are observed mostly further downstream along plane P5 and beyond. As also explained in Fig. 5 (bottom), the presence of few isolated liquid mass or droplets in the upstream region (see planes P1, P2) are due the early primary breakup associated with surface stripping induced by gas etrainment.

Fig. 7

Snapshot of droplet dynamics along the various planes in the jet direction: liquid core in blue, detached liquid structures in red

Fig. 8

Evolution of primary breakup processes and typical secondary breakup evolution of detached liquid structures at six different times with $\Delta t= 5 \mu s$

Additionally, a large amount of liquid can be seen impinging on the wall of counter-bore region and stay attached on the wall until the nozzle exit, which results in a dense distribution of smaller detached liquid structures around the jet-gas interface (see plane P5) while few structures are observed in the core region. Whereas, the liquid structures (individual isolated liquid mass) is getting further smaller from P5 to P8 suggesting the on-set of secondary atomization processes, the distribution of the isolated liquid mass is becoming more homogeneous, as well.

Fig. 9

Quantitative evolution of a selected typical detached liquid structure located at the red marked region in Fig. 8 at each time from (2) to (6) (top to bottom). M for Multi-mode breakup, B for Bag breakup and O/D for Oscillation/ Deformation regime

Figure 8 depicts typical primary breakup processes and a subsequent bag breakup evolution (color marked regions) at six different times with $\Delta t= 5 \mu s$. The present liquid structure in Fig. 8 (1) is part of the liquid core at the position close to plane P4, in which bag formation can be found (red marked). Subsequently, this liquid structure disintegrates into three large fragments and several droplets as shown in Fig. 8 (2). Focusing on the red marked region, from (2) to (3), the red marked large fragment burst in line with bag breakup mechanism, resulting in a dense distribution of ligament formation (shown in (3)). The breakup of these ligament produces a large number of droplets and a few of relative larger fragments (shown in (4)). Thereafter, since the relative velocity in gas and liquid phases is getting small, these relative larger fragments either decomposes into smaller droplets in line with oscillation breakup mechanism or deform only their shapes without further breakup (in (5) and (6)). Note that, the center of all droplets at (6) is in the near upstream of P8. Similar breakup evolution occurs in the green and blue marked region with a delay of $~5\mu s$ comparing with the red marked region, e.g. bag formation (green) and multi-bag formation (blue) are presented at (2) while bag breakup and the ligament breakup occur at (3) and (4), respectively. These qualitative evolution and breakup mechanisms of the detached liquid structures is confirmed by the quantitative relation of We- and Oh-numbers in Fig. 9, which is evaluated for a typical selected liquid structure in the red marked region at each time from (2) to (6) (from top to bottom). The regime map is obtained as functions of Weber and Ohnesorge numbers with multimode breakup, bag breakup and oscillation/deformation regimes clearly demarcated. Here, the detached structures are assumed to be spherical and the their Weber numbers and Ohnesorge numbers are calculated by

$$\begin{aligned} We = \rho _g\left( U^{l} -U^{g}\right) ^2d/\sigma , ~~ Oh = \mu _l/\sqrt{\rho _l\sigma d}, \end{aligned}$$

(6)

where d represents the hydraulic diameter of the liquid structures, $U^{l}$ the volume averaged liquid velocity, $U^{g}$ the spatial averaged gas velocity around the selected liquid structures.

Fig. 10

Map of the droplet breakup regimes for detached droplet structures on the selected planes (P3-P8) at four physical time intervals of $\Delta t = 50\mu s$ (in different colors). M for Multi-mode breakup, B for Bag breakup and O/D for Oscillation/Deformation regime

In order to have a global view of the secondary break up in the ECN Spray G configuration, the map of the droplet breakup regimes for detached droplet structures on the selected planes (P3-P8) are depicted quantitatively in Fig. 10 for four different time steps with $\Delta t = 50\mu s$, with which the aliasing effects are negligible. It is worth mentioning that in this case the Ohnesorge number is calculated according to Eq. (6), while the $U^g$ of a selected plane is averaged spatially over this plane for the gas phase. It can be clearly seen in the scattered plot that in planes P3-P5, the liquid structures experience both bag and multi-mode breakup in addition to oscillation/deformation breakup, as observed in Fig. 8. Further downstream (P6-P8), the most structures undergo oscillation/deformation breakup occurring at low Weber number. This suggests the completion of secondary breakup and generation of more stable droplets ($\textrm{We}\le$1).

Fig. 11

Evolution of spatial averaged sphericity $\varphi$ for the selected planes in the chamber at four physical time intervals of $\Delta t = 50\mu s$ (in different colors). $\Delta l$ represents the distance to the nozzle exit

The breakup of ligaments may not always lead to formation of spherical primary droplets, while aerodynamic induced surface instabilities can also deform the droplets considerably. Moreover, the droplet shape (sphericity) has great impact on the local interfacial exchanges of mass, momentum, and energy. In this study, the droplet shapes are characterized by an averaged sphericity $\varphi$, which represents the max-longitudinal-tangential ratio of liquid structures, are shown for these planes (P4–P8) in Fig. 11. Additionally to evaluate the transient evolution of droplet shapes at these locations, the data ($\varphi$) for four physical times with $\mathrm {\Delta t} = 50 \mathrm {\mu s}$ are also plotted in different colors. The comparison confirms a time independent evolution of droplet shapes. Notice that, $\Delta l=0$ is located at the nozzle exit. With an increasing distance from the nozzle exit, the dominant ligaments break up into much smaller and stable droplets, which are getting more spherical and therefore more stable as observed in Fig. 11 with $\varphi$ values progressively converging towards unity.

4.4 Droplet Statistics

The probability density function (PDF) of the droplet distribution in this configuration is evaluated in terms of the representative droplet diameter and volume averaged velocity, which are rarely available in reported literature by experimental studies or other numerical studies for a spray G configuration under study.

Fig. 12

Number based droplet size distribution for four physical time intervals of $\Delta t = 50\mu s$ (in circle) with fitting function by log-normal distribution (in black solid line) and Weibull distribution (in black dashed line)) along the selected planes (P1-P8) (see Fig. 7)

In Fig. 12 the number based droplet size PDF obtained by the DNS simulation is plotted for the selected planes for four physical times with an interval of $\Delta t = 50 \mu s$ in circle with different colors. Notice that, the stationary droplet size PDFs can be represented by a log-normal distribution $\ln (l_d) \sim N(\epsilon , \psi ^2)$ (black solid line in Fig. 12), which reads:

$$\begin{aligned} f(l_d) = \frac{1}{l_d\psi \sqrt{2\pi }}\exp \left( -\frac{\left( \ln (l_d) -\epsilon \right) ^2 }{2\psi ^2}\right) . \end{aligned}$$

(7)

Here, the mean related coefficient $\epsilon$ and variance associated coefficient $\psi$ are given in Table 3. It is worth mentioning that the Weibull distribution (black dashed line in Fig. 12) failed to represent the droplet size PDF in this ECN spray G configuration.

Table 3 The fitting coefficients of log-normal distribution for number based droplet size distributions

The global max value of $\psi$ and $\epsilon$ is found on P2, and $\epsilon$ increases in the chamber (P4 to P8) that might be caused by the change of grid size see fine case in Table 5, while the decrease of $\psi$ in chamber suggests a more homogeneous droplet distribution along flow direction. It is worthwhile to mention that the peek position of the number based PDF is significantly smaller than the Sauter Mean Diameter (SMD) of the droplets passing through the plane, e.g. on P8 the peak position is located at 4.8 μm, while the SMD of these droplets is 11.4 μm that is comparable with the measured SMD in further downstream under the operating condition of spray G1, reported in Weiss et al. (2020).

Figure 13 shows the droplet-velocity correlations for above mentioned four physical times (in circle) with best-fitting function by a modified exponential distribution (in black solid line) for the selected planes as given by:

$$\begin{aligned} f(l_d) =\lambda - \lambda \exp {(-\beta l_d)}, \end{aligned}$$

(8)

Fig. 13

Droplet size-velocity correlations for four physical time intervals of $\Delta t = 50$ μs (in circle) with best-fitting function by a modified exponential distribution (in black solid line) for the selected planes (see Fig. 7)

where the coefficients $\beta$ and $\lambda$ are listed in Table 4. Similar to the stationary droplet size distribution the droplet size-velocity correlation depends weakly on the time. Whereas the liquid losses its momentum along the flow direction and all droplets on P8 have similar order of velocity magnitude. Moreover, the larger liquids have a higher speed than the smaller droplets, especially inside the nozzle (P1–P3). The droplet size-velocity correlation is also provided as joint PDFs for planes P4 in Fig. 14 and P8 in Fig. 15. It is also evident here that (see also Fig. 13) the collectively droplets at plane P4 has higher velocity as compared to plane P8 (see Bauckhage and Floegel 1984).

Table 4 Fitting coefficient of the modified exponential distribution for droplet size-velocity correlation

Fig. 14

The droplet size-velocity correlation represented by joint PDFs of droplet size with respective droplet velocity along plane P4

Fig. 15

The droplet size-velocity correlation represented by joint PDFs of droplet size with respective droplet velocity along plane P8

4.5 2D Data Curation Technique for Droplet Statistics from VoF/DNS

In a DNS study, apart from it being very computational intensive, the processing and managements of the obtained numerical results are very demanding task due to sheer number of data storage and RAM capacity required. It becomes even more challenging in the need of frequent data storage. Thereby, it is always not feasible to have the droplet size distribution based on the series of three-dimensional results, even if a three-dimensional simulation is carried out. In this regard, a data curation technique based on the obtained two-dimensional results can be more practical and economical to store and post-process the time series data with higher frequency.

Fig. 16

A schematic representation of 2D data curation technique to obtain the droplet statistics

Fig. 17

Number based droplet size distribution obatined using 2D data curation technique for four physical time intervals of with $\Delta t = 50\mu s$ (in circle) with best-fitting function by log-normal distribution (in black solid line) along the selected planes (P1-P8) (see Fig. 7 for post-processing planes)

Fig. 18

Droplet size-velocity correlations obtained using 2D data curation techniques for four physical time intervals of $\Delta t = 50$ μs (in circle) with best-fitting function by a modified exponential distribution (in black solid line) for the selected planes P1-P8 (see Fig. 7 for post-processing planes)

Thereby, two-dimensional droplet size distributions on required planes are extracted in 5 steps (see Fig. 16 for schematic representation of these individual steps), as follows:

• 1. a threshold of $\alpha _f \ge 0.5$ is used to characterize liquid regions for each slice P1–P8;

• 2. to separate connected liquid fragments, an erosion operation with an operator size of $\Delta _F = 2\Delta _G$ ($\Delta _G$ grid size) is performed;

• 3. dilation is performed with the same operator size $\Delta _F$ in step 2;

• 4. droplets are counted, sorted by the diameter and scaled in order to be mass conservative;

• 5. Finally, number based PDFs are calculated.

Fig. 19

Comparison of log-normal distribution coefficients

Fig. 20

Comparison of exponential distribution coefficients

The droplet size distribution obtained using the 2D data extraction method is shown in Fig. 17 along the same planes for which 3D droplets statistic were collected (see also Fig. 12). The comparisons among various time intervals also suggest a time independent droplet size distribution. Additionally, the droplet velocity distributions are depicted in Fig. 18. In order to carry out direct comparison with the three-dimensional data, the droplet size distributions extracted from two-dimensional data as represented by the log-normal distribution and the fitting coefficients are tabulated in Table 3, while the two-dimensional velocity-droplet size distributions are represented by exponential distribution function with coefficient listed in Table 4. Notice that, in two-dimensional data frame, the detached droplet structure can not be exclusively separated from the attached structures and, similarly, a large ligament can be counted as several independent droplets. Therefore, in the primary breakup region (P1–P5) the two-dimensional droplet statistics differs significantly from those provided by three-dimensional data. This can be ascertained further by the comparison of the fitting coefficients, as shown in Fig. 19 and Fig. 20, where the coefficients of two-dimensional procedure are plotted in red cross, while the reference three-dimensional coefficients are shown in circle. In the primary breakup region the liquid core can not be separated from detached liquid by using 2D data extraction method and many large ligaments are present in this region, as well. Therefore, the coefficients provided by 2D data extraction method are in poor agreement with the reference 3D data. Whereas, in the secondary breakup region with less large ligaments (P5–P8), the fitting coefficients provided by two-dimensional procedure represent well the reference three-dimensional data for both mean (or max value) and variation relevant coefficients. Moreover, from P5 to P8 less large structures relevant to multi-mode breakup or bag-mode breakup regimes are present. In the opposite, more fine droplets that can be found in the oscillation or deformation regime appear. Thereby, a trend of better representation of the reference 3D data by the 2D data extraction method can be observed. It can therefore be concluded that in the dilute region all the coefficients represented by the 2D data curation technique should match the 3D data well.

5 Conclusion

A DNS based numerical investigation coupled with VOF as an interface tracking method is carried out for detailed understanding of processes evolving during the fuel injection in one of the Spray G configuration from Engine Combustion Network (ECN). Which is a gasoline direct injector (GDI) featuring 8-holes orifices and operating with high injection pressure (P_inj= 200 bar). Thereby, by taking into account the computational cost associated with DNS, only 1/8 of the nozzle geometry including one orifice is considered.

The numerical simulation is accomplished for the quasi-steady injection condition with nozzle needle fully opened. The DNS results is first validated for experimentally measured nozzle mass flow rate and spray spread angle showing good agreement. Subsequently, a detailed numerical analysis is carried out to gain the deep understanding of transient flow, turbulence and spray evolution inside the nozzle-orifice, counter-bore and chamber regions. In particular, the jet and/or spray morphology, the breakup regime are discussed in detail together with droplet statistics in terms of droplet shapes, size and velocity. Finally, a 2-D data curation technique is proposed to extract the droplet statistics along selected planes and evaluated by direct comparison with three-dimensional droplet data.

Based on the comprehensive numerical analysis, and comparison against the experimental data, the following conclusions can be drawn from this work:

• The sudden change in flow direction and the counter-bore design of nozzle exit induce a significant gas entrainment inside the nozzle orifices, that is reckoned to expedite and enhance the overall primary atomization process.

• The flow separation at orifice base and gas entrainment are responsible for an asymmetric jet-core combined with surface instability results in early and enhanced primary atomization.

• The larger droplets experience further instability during primary atomization combined with aerodynamic instability resulting in more spherical and stable smaller droplets.

• The obtained droplet statistics appears to be time independent especially in terms of droplet size distributions and size-velocity correlations.

• The DNS data generated make it easy to provide more reliable injector boundary conditions useful for lower order spray injection method based on Lagrangian particle tracking as already demonstrated by Lien et al. (2024).

• The new 2-D data curation technique proposed and evaluated against the 3-D data allows handling of the DNS data in more feasible and economical way especially for time series data with higher frequency.

In overall, the DNS demonstrates its great potential to enable detailed numerical analysis of in- and near nozzle complex liquid and flow phenomena for which the experimental data are still scarce.

Appendix A: Effect of Grid Size on the Droplet Statistics

To analyze the influence of grid size on the droplet statistics, an additional simulation with a coarser mesh is carried out with the grid size globally enlarged with a factor 1.5. The variation of grid size on these selected planes are listed in Table5 for the fine and coarse mesh, respectively. The mass flow rate and the spread angle barely depends on the grid size. The droplet statistics on the selected planes using coarse mesh are evaluated using the 2D data curation technique as described in Sect. 4.5, and compared with the data of the fine case. The influence of the grid size on the droplet statistics is shown in Fig. 21 for the number based PDF exemplary on plane P8. The PDF of both fine and coarse cases are represented by the log-normal distribution, whereas the peak position differs slightly with the peak value showing the grid size dependency. Ling et al. (2017) highlighted the importance of mesh size in resolving correctly not only the smaller droplets but also the larger droplets. Additionally, they pointed out that the droplets size smaller than $2\Delta$ ($\Delta$ is the grid size) are most likely not to be well captured. This is confirmed by the observation in the present study, as the droplet size distribution shift toward the larger diameter with coarser mesh size (see, Fig. 21).

Table 5 Averaged cell sizes (in $\mu m$) on the selected planes P1 to P8 for the fine and the coarse case

Fig. 21

Number based droplet size distribution obtained using 2D data curation technique (see, Sect. 4.5) for plane P8