Impact of the Primary Break-Up Strategy on the Morphology of GDI Sprays in 3D-CFD Simulations of Multi-Hole Injectors

: The scientiﬁc literature focusing on the numerical simulation of fuel sprays is rich in atomization and secondary break-up models. However, it is well known that the predictive capability of even the most di ﬀ used models is a ﬀ ected by the combination of injection parameters and operating conditions, especially backpressure. In this paper, an alternative atomization strategy is proposed for the 3D-Computational Fluid Dynamics (CFD) simulation of Gasoline Direct Injection (GDI) sprays, aiming at extending simulation predictive capabilities over a wider range of operating conditions. In particular, attention is focused on the e ﬀ ects of back pressure, which has a remarkable impact on both the morphology and the sizing of GDI sprays. 3D-CFD Lagrangian simulations of two di ﬀ erent multi-hole injectors are presented. The ﬁrst injector is a 5-hole GDI prototype unit operated at ambient conditions. The second one is the well-known Spray G, characterized by a higher back pressure (up to 0.6 MPa). Numerical results are compared against experiments in terms of liquid penetration and Phase Doppler Anemometry (PDA) data of droplet sizing / velocity and imaging. CFD results are demonstrated to be highly sensitive to spray vessel pressure, mainly because of the atomization strategy. The proposed alternative approach proves to strongly reduce such dependency. Moreover, in order to further validate the alternative primary break-up strategy adopted for the initialization of the droplets, an internal nozzle ﬂow simulation is carried out on the Spray G injector, able to provide information on the characteristic diameter of the liquid column exiting from the nozzle.


Introduction
Spray modeling capabilities represent a crucial aspect in 3D-Computational Fluid Dynamics (CFD) simulations of internal combustion engines. For example, focusing on Gasoline Direct Injection (GDI) engines, spray directly affects air-fuel mixing [1][2][3] and, thus, combustion [4,5], knock [6,7], and emissions [8]. Moreover, knock suppressor techniques such as water (or water/methanol) injection are more and more diffused [9][10][11][12]. In this case, a detailed numerical representation of the water spray is mandatory to properly predict water evaporation or even phenomena such as liquid film formation on the intake ports, which impact on the effectiveness of water injection itself.
For 3D-CFD simulations of fuel injection, the Arbitrary Lagrangian-Eulerian (ALE) approach is the most widespread, mainly due to its computational efficiency. In the scientific literature, a very wide variety of models and methods is available to mimic the numerous physical processes which characterize the spatial and temporal evolution of the injected liquid. Among such processes, atomization of the exiting liquid column and primary break-up into droplets play a fundamental role. Detailed reviews of some of the most popular models are available in [13][14][15][16][17].
In the ALE framework, computational parcels are introduced directly at the nozzle exit, hence the choice of the initial conditions (i.e., droplet diameter and velocity) is extremely important to achieve fully representative simulations of the spray temporal and spatial evolution [18]. Despite coupling of internal nozzle flow simulations to Lagrangian ones has been recently introduced and applied to directly infer the characteristics of exiting droplets based on nozzle flow and turbulence [19][20][21][22], a vast majority of the simulations of fuel spray relies on primary breakup phenomenological models. As for the initial droplet diameter, an alternative approach to both internal flow simulations and phenomenological models consists in the adoption of droplet diameter distribution functions inferred by experimental Phase Doppler Anemometry (PDA) data, which is a very common practice in the industrial field. In fact, in the industrial framework, a rapid spray calibration is often preferred at the expense of the predictive capabilities of the numerical model. However, for the sake of experimental data validity, PDA measurements are usually carried out at least 15 ÷ 20 mm far from the injector tip, where droplets have already undergone massive secondary break-up. Therefore, PDA-based droplet initial diameters to be adopted in Lagrangian simulations are often characterized by values even one order of magnitude smaller than nozzle hole diameters. In the present paper, such an approach (even if questionable) provides relatively consistent results (i.e., in line with measurements) if the backpressure is close to the ambient one (typical condition for spray calibration).
To this aim, Lagrangian simulations at different injection pressures are carried out on a 5-hole GDI prototype injector [23]. Simulations are characterized by initial droplet diameters lower than 10 µm, inherited from "far field" experimental PDA data. The same atomization strategy is then applied to the well-known Spray G. The analyzed operating condition is the most investigated one in literature and it is characterized by a backpressure equal to 0.6 MPa. In this case, the adoption of small diameters (nearly 10 µm) coming from PDA reveals unacceptable misalignment with experimental outcomes. In fact, unlike experimental evidences, numerical spray tends to collapse similarly to the behavior exhibited under flash boiling conditions [24] (which are considerably far from the investigated ones). Therefore, larger droplet diameters (closer to the hole size of 165 µm) are tested, able to provide a representation of the spray in line with the experiments. In order to further confirm the importance of droplet diameters comparable with hole dimension, an internal nozzle flow simulation on the Spray G injector is carried out. This is able to provide information on droplet initial conditions (to be applied to Lagrangian simulations) and, hence, on droplet diameters. These last are found to be slightly greater than 130 µm, thus similar to nozzle hole dimension and one order of magnitude larger than values provided by PDA measurements at 15 mm from the injector tip. This leads us to conclude that, apart from particular conditions such as flashing ones (usually characterized by small droplets exiting from the nozzle), a simpler blob model with droplet dimension equal to hole diameter may perform much better in terms of 3D-CFD numerical results than diameter distribution functions.
As stated earlier, for the Lagrangian simulations of both 5-hole prototype and Spray G, numerical results are validated against experimental data in terms of liquid length, particle diameter, velocity provided by PDA and, finally, in terms of imaging. As for the internal nozzle flow simulation on the Spray G injector, simulation outcomes are compared to experiments in terms of mass flow rate, hydraulic coefficients and spray pattern.
At the end of the present introduction, investigated operating conditions and numerical setup are presented. Then results on both 5-hole injector and Spray G are discussed. After that, both numerical setup and results of the internal nozzle flow simulation are shown. Finally, conclusions on the activity are drawn.

Spray G Injector
Liquid length, PDA data and hydraulic coefficients are available on the Engine Combustion Network (ECN) website [25]. As for the coefficients, they are obtained from mass flux and momentum flux measurements carried out at Càtedra Motores Térmicos (CMT) in Valencia [26]. PDA measurements are available at different locations on a plane 15 mm far from the injector tip, as depicted in Figure 1. For the sake of validity of the experimental data, only three locations are considered for the comparison with numerical outcomes.

Spray G Injector
Liquid length, PDA data and hydraulic coefficients are available on the Engine Combustion Network (ECN) website [25]. As for the coefficients, they are obtained from mass flux and momentum flux measurements carried out at Càtedra Motores Térmicos (CMT) in Valencia [26]. PDA measurements are available at different locations on a plane 15 mm far from the injector tip, as depicted in Figure 1. For the sake of validity of the experimental data, only three locations are considered for the comparison with numerical outcomes.

5-Hole Injector
Experimental tests over a 5-hole injector are carried out at the SprayLAB of the University of Perugia. It is a side-mounted, counter bore, GDI prototype, whose hole diameters are equal to 125 µm. The injector is experimentally investigated in a wide range of injection pressures, from 5 up to 60 MPa [27]. For the sake of brevity, only three pressure levels (20,40, and 60 MPa) are considered for the comparison with numerical simulations.
A complete hydraulic characterization of the prototype is carried out. In order to obtain statistically significant results, Injection Rate (IR) profiles of 300 consecutive shots are measured. The mean IR profiles, for an energizing time (ET) of 1.5 ms and for the examined injection pressures, are reported in Figure 2. For each one, also the coefficient of variation of the injected mass is reported, to provide an estimation of the shot-to-shot variability. Global evolution of the spray plumes is investigated by means of high-speed imaging. The resulting average tip penetration curves for the most advanced plume are reported in Figure 3. Droplet sizing and velocity are investigated at different locations. The ones used for the numerical-experimental comparisons are reported in Figure  4. As an example, PDA raw data (along with their average), measured at x = −1.5 mm, y = −1 mm and z = 50 mm, are reported in Figure 5.

5-Hole Injector
Experimental tests over a 5-hole injector are carried out at the SprayLAB of the University of Perugia. It is a side-mounted, counter bore, GDI prototype, whose hole diameters are equal to 125 µm. The injector is experimentally investigated in a wide range of injection pressures, from 5 up to 60 MPa [27]. For the sake of brevity, only three pressure levels (20,40, and 60 MPa) are considered for the comparison with numerical simulations.
A complete hydraulic characterization of the prototype is carried out. In order to obtain statistically significant results, Injection Rate (IR) profiles of 300 consecutive shots are measured. The mean IR profiles, for an energizing time (ET) of 1.5 ms and for the examined injection pressures, are reported in Figure 2. For each one, also the coefficient of variation of the injected mass is reported, to provide an estimation of the shot-to-shot variability. Global evolution of the spray plumes is investigated by means of high-speed imaging. The resulting average tip penetration curves for the most advanced plume are reported in Figure 3. Droplet sizing and velocity are investigated at different locations. The ones used for the numerical-experimental comparisons are reported in Figure 4. As an example, PDA raw data (along with their average), measured at x = −1.5 mm, y = −1 mm and z = 50 mm, are reported in Figure 5. provide an estimation of the shot-to-shot variability. Global evolution of the spray plumes is investigated by means of high-speed imaging. The resulting average tip penetration curves for the most advanced plume are reported in Figure 3. Droplet sizing and velocity are investigated at different locations. The ones used for the numerical-experimental comparisons are reported in Figure  4. As an example, PDA raw data (along with their average), measured at x = −1.5 mm, y = −1 mm and z = 50 mm, are reported in Figure 5.
A single operating condition is simulated for the Spray G, whose main characteristics are reported in Table 1. This is the most investigated operation in the ECN community and a lot of experimental data are available, provided by different institutions [30][31][32]. As for the 5-hole prototype injector, operating conditions selected for the numerical investigations are reported in Table 2. Table 2. Simulated operating conditions for the 5-hole injector.
A single operating condition is simulated for the Spray G, whose main characteristics are reported in Table 1. This is the most investigated operation in the ECN community and a lot of experimental data are available, provided by different institutions [30][31][32]. As for the 5-hole prototype injector, operating conditions selected for the numerical investigations are reported in Table 2. The computational domains consist in block-shaped vessels whose characteristic dimensions are 40 × 40 × 80 mm for the Spray G and 130 × 140 × 180 mm for the 5-hole injector, respectively. Numerical grids depicted in Figures 6 and 7 consist of hexahedral cells and they are characterized by the presence of cone-shaped refinements in the spray core region, whose minimum cell size is equal to 0.4 mm for the Spray G and 0.7 mm for the 5-hole injector [23]. As known, main spray characteristics such as penetration are deeply affected by the numerical grid size. This is the reason why a grid sensitivity analysis is carried out for both the investigated injectors and minimum cell sizes reported above represent the best trade-off between reasonable computational cost and reduced grid-dependency [33].   On the left, a numerical grid adopted for the 5-hole injector simulations; on the right, a section is reported along with a detail.
A Reynolds-averaged Navier-Stokes (RANS) approach to turbulence is adopted for all the simulations. The widely diffused k-ε Renormalization Group (RNG) two-equation turbulence model [34] is adopted. A combined Eulerian-Lagrangian approach allows to properly account for both the vessel gaseous ambient and the dispersed liquid phase [35]. The second order Monotone Advection and Reconstruction Scheme (MARS) numerical scheme is adopted for momentum, temperature, and turbulent quantities transport equations. For both the injectors, 15 parcels are injected from each single nozzle at every time-step; this last is fixed to nearly 1e-6 s to keep the maximum Courant number well below unity. The computational domain is initialized with experimental pressure and temperature values; since vessel is supposed to be quiescent, k and ε initial values are set equal to 0. Apart from the top of the domain, which is modeled as a non-slip adiabatic wall, all other boundaries are set as pressure outlets. Single-component Lagrangian parcels are injected, whose properties are inherited from the National Institute of Standards and Technology (NIST) database both for the liquid   On the left, a numerical grid adopted for the 5-hole injector simulations; on the right, a section is reported along with a detail.
A Reynolds-averaged Navier-Stokes (RANS) approach to turbulence is adopted for all the simulations. The widely diffused k-ε Renormalization Group (RNG) two-equation turbulence model [34] is adopted. A combined Eulerian-Lagrangian approach allows to properly account for both the vessel gaseous ambient and the dispersed liquid phase [35]. The second order Monotone Advection and Reconstruction Scheme (MARS) numerical scheme is adopted for momentum, temperature, and turbulent quantities transport equations. For both the injectors, 15 parcels are injected from each single nozzle at every time-step; this last is fixed to nearly 1e-6 s to keep the maximum Courant number well below unity. The computational domain is initialized with experimental pressure and temperature values; since vessel is supposed to be quiescent, k and ε initial values are set equal to 0. Apart from the top of the domain, which is modeled as a non-slip adiabatic wall, all other boundaries are set as pressure outlets. Single-component Lagrangian parcels are injected, whose properties are inherited from the National Institute of Standards and Technology (NIST) database both for the liquid and the vapor phase [36]. For both injectors, experimental injection rate profiles are adopted as mass A Reynolds-averaged Navier-Stokes (RANS) approach to turbulence is adopted for all the simulations. The widely diffused k-ε Renormalization Group (RNG) two-equation turbulence model [34] is adopted. A combined Eulerian-Lagrangian approach allows to properly account for both the vessel gaseous ambient and the dispersed liquid phase [35]. The second order Monotone Advection and Reconstruction Scheme (MARS) numerical scheme is adopted for momentum, temperature, and turbulent quantities transport equations. For both the injectors, 15 parcels are injected from each single nozzle at every time-step; this last is fixed to nearly 1e-6 s to keep the maximum Courant number well below unity. The computational domain is initialized with experimental pressure and temperature values; since vessel is supposed to be quiescent, k and ε initial values are set equal to 0. Apart from the top of the domain, which is modeled as a non-slip adiabatic wall, all other boundaries are set as pressure outlets. Single-component Lagrangian parcels are injected, whose properties are inherited from the National Institute of Standards and Technology (NIST) database both for the liquid and the vapor phase [36]. For both injectors, experimental injection rate profiles are adopted as mass flow rates.
Moreover, the syringe-like effect described in [23] is taken into account to improve numerical outcomes during the first stage of injection. As a consequence of this effect, both initial mass flow rates and velocities are characterized by non-zero values. Primary break-up is replaced by a simplified Blob model. As for the secondary break-up, only Reitz's [37] model is used for the 5-hole injector simulations. For the Spray G, both Reitz's and Kelvin-Helmholtz Rayleigh-Taylor (KHRT) [38] models are used. For all the simulations proposed in the present work, no secondary break-up model constants are modified compared to the reference papers.

5-Hole Injector Numerical Results
As mentioned in the introduction, the approach adopted for the spray modeling is mainly based on experimental outcomes. In particular, droplet initialization relies on PDA measurements at 50 mm from the tip. At this distance, geometrical average droplet diameters are nearly 7, 6, and 5 µm for injection pressures of 20, 40, and 60 MPa, respectively. Therefore, similar values are used to initialize droplets in the Lagrangian simulations. In particular, 8, 7, and 6 µm are droplet initial diameters for injection pressures of 20, 40, and 60 MPa, respectively. As for the initial droplet velocities at start of injection, values are chosen to match penetration curves while, for the static initial velocities, they are estimated from experimental static mass flow rates; for each injection pressure, single hole steady state actual mass flow rate . m r,i (which is the actual injector mass flow rate divided by hole number) is divided by both liquid density ρ l (at injection temperature) and geometrical area of the nozzle hole A, as follows: In this case, an approximation is made since the geometrical area is used and no vena contraction is accounted for. For all the investigated conditions, the cone angle of each single plume is set equal to 30 • , as indicated by the experiments.
The first comparison between numerical and experimental data deals with liquid penetrations. droplets in the Lagrangian simulations. In particular, 8, 7, and 6 µm are droplet initial diameters for injection pressures of 20, 40, and 60 MPa, respectively. As for the initial droplet velocities at start of injection, values are chosen to match penetration curves while, for the static initial velocities, they are estimated from experimental static mass flow rates; for each injection pressure, single hole steady state actual mass flow rate ̇ , (which is the actual injector mass flow rate divided by hole number) is divided by both liquid density (at injection temperature) and geometrical area of the nozzle hole , as follows: In this case, an approximation is made since the geometrical area is used and no vena contraction is accounted for. For all the investigated conditions, the cone angle of each single plume is set equal to 30°, as indicated by the experiments.
The first comparison between numerical and experimental data deals with liquid penetrations. Comparisons at 20, 40 and 60 MPa are reported in Figures 8,9, and 10 respectively.     It is useful to point out that numerical spray penetrations are computed as the distance from the injector nozzle outlet section at which 95% of the plume mass is met, consistently with a widely diffused practice [39]. Moreover, for all the conditions, despite liquid penetration along the injector axis considers the whole spray, it is always determined by the central plume whose axis is almost coincident with the injector one.
That said, numerical penetrations closely reproduce experimental outcomes for all of the investigated conditions.
As for the comparison in terms of PDA data, experimental values represent an average over a specific time interval. For a proper comparison, the same time average is considered also for the numerical results.
PDA comparison is carried out at 20, 30, 40 and 50 mm downstream the injector tip, as reported in Figure 4. From the tip up to a 20 mm distance, experimental data spherical validation is too low to ensure reliable outcomes. This is the reason why, at 20 mm, differences between experimental and numerical outcomes are more evident, as visible in Figures 11 and 12.
At 30, 40 and 50 mm, experimental mean droplet diameters (D10) and velocities are reasonably matched by numerical results for each analyzed injection pressure. In fact, both experiments and simulations show a decrease of the geometric diameter for increasing injection pressures. It is worthwhile to point out that initial droplet diameters and values measured at 20, 30, 40, and 50 mm away from the tip are very similar, which proves that secondary break-up poorly affects results for the chosen numerical setup in terms of initial diameters.
As for the average velocities, it is interesting to note that, despite an increasing initial velocity with injection pressure, an inversion of the trend can be noticed for injection pressures higher than It is useful to point out that numerical spray penetrations are computed as the distance from the injector nozzle outlet section at which 95% of the plume mass is met, consistently with a widely diffused practice [39]. Moreover, for all the conditions, despite liquid penetration along the injector axis considers the whole spray, it is always determined by the central plume whose axis is almost coincident with the injector one.
That said, numerical penetrations closely reproduce experimental outcomes for all of the investigated conditions.
As for the comparison in terms of PDA data, experimental values represent an average over a specific time interval. For a proper comparison, the same time average is considered also for the numerical results.
PDA comparison is carried out at 20, 30, 40 and 50 mm downstream the injector tip, as reported in Figure 4. From the tip up to a 20 mm distance, experimental data spherical validation is too low to ensure reliable outcomes. This is the reason why, at 20 mm, differences between experimental and numerical outcomes are more evident, as visible in Figures 11 and 12.
worthwhile to point out that initial droplet diameters and values measured at 20, 30, 40, and 50 mm away from the tip are very similar, which proves that secondary break-up poorly affects results for the chosen numerical setup in terms of initial diameters.
As for the average velocities, it is interesting to note that, despite an increasing initial velocity with injection pressure, an inversion of the trend can be noticed for injection pressures higher than 40 MPa. This behavior is well captured by both simulations and experiments.     As for the average velocities, it is interesting to note that, despite an increasing initial velocity with injection pressure, an inversion of the trend can be noticed for injection pressures higher than 40 MPa. This behavior is well captured by both simulations and experiments.
The last comparison between experimental and numerical data is carried out in terms of imaging. For the sake of brevity, only a few snapshots are chosen and they are reported in Figures 13-15.       The main spray characteristics such as the overall cone angle and relative penetrations between the different plumes are well captured by simulations.    The main spray characteristics such as the overall cone angle and relative penetrations between the different plumes are well captured by simulations. The main spray characteristics such as the overall cone angle and relative penetrations between the different plumes are well captured by simulations.

Spray G Injector Lagrangian Simulation Results
In the light of the promising results obtained with the 5-hole injector, the same numerical approach is adopted for the Spray G. In particular, the initial droplet diameter is set equal to nearly 10 µm, which roughly corresponds to the experimental Sauter Mean Diameter (SMD) measured on a plane 15 mm away from the tip along the injector axis. As for the static initial droplet velocity, compared to the previous case, a more reliable estimation can be obtained thanks to the availability of C d , C a and C v [25], which are discharge, contraction and velocity coefficients, respectively, and whose definitions are reported hereafter.
Initial droplet velocity can be obtained as follows.
A e f f is the effective area computed as C a A. Alternately, velocity can be obtained as C v v th , v th being the theoretical velocity provided by the Bernoulli equation. Initial droplet velocity at the start of injection is determined, similarly to the previous injector, in order to match penetration curves at the first stage of injection. Finally, the cone angle of each single plume is set equal to 30 • , as reported by experiments.
Before moving to the discussion of the results, it is useful to anticipate at this stage that the adoption of a small initial droplet diameter does not allow to obtain an acceptable numerical representation of the Spray G, if compared to the experimental evidence. For this reason, larger initial diameters are also investigated; their values are reported in Table 3. A first numerical-experimental comparison is proposed in Figure 16, which shows liquid penetrations for the different investigated droplet initial diameters. Initial droplet diameter is able to affect liquid penetration both at the beginning and at the end of the injection. In particular, if initial diameter decreases, initial liquid penetration decreases as well, due to a lower droplet momentum. Even if it is not so evident because of the reduced injection duration, moving towards the end of the injection, a more complicated trend can be noticed: starting from 110 µm and reducing droplet diameter, penetration reduces as well. However, moving from 50 µm to 20 µm liquid length increases. Finally, adopting a diameter of 10 µm penetration remarkably increases.
Energies 2019, 12, x FOR PEER REVIEW 11 of 23 condition (0.1 MPa). In fact, as visible in Figure 18, with low backpressure (5-hole injector), the pressure difference between the core of the spray and the outer region is almost negligible. Conversely, at higher backpressure (Spray G), the difference is remarkable. Because of the lower pressure in the spray core region, plumes are attracted by each other. Only larger diameters ensure a droplet momentum able to overcome attraction between plumes.  Two observations are necessary: the first one regards the slight over-penetration obtained with the initial diameters of 80 and 110 µm. This may be solved with a dedicated tuning of the secondary break-up. In fact, it is useful to remind that simulations are carried out on equal secondary break-up model, namely the Reitz's one and no specific tuning is carried out as it is not the main focus of the present activity. The second remark deals with the poor results obtained with 10 µm as initial droplet diameter, as anticipated at the beginning of the paragraph. This is inconsistent with expectations as, in the light of previous results on the 5-hole injector, such droplet sizing was considered reasonable.
In order to understand the reasons behind such unexpected behavior, a comparison between numerical snapshots and experimental imaging 0.5 ms after the start of injection is shown in Figure 17.
While larger initial droplet diameters lead to a consistent numerical representation of the Spray G, smaller values produce large deviations from the experimental outcomes. Moving from 110 to 10 µm, a collapse of the spray plumes is increasingly evident. Such phenomenon, clearly visible in the case with droplet initial diameter of 10 µm, is the main responsible for the over-penetration in Figure 16 and it is mainly due to the higher back pressure (0.6 MPa) compared to the 5-hole injector operating condition (0.1 MPa). In fact, as visible in Figure 18, with low backpressure (5-hole injector), the pressure difference between the core of the spray and the outer region is almost negligible. Conversely, at higher backpressure (Spray G), the difference is remarkable. Because of the lower pressure in the spray core region, plumes are attracted by each other. Only larger diameters ensure a droplet momentum able to overcome attraction between plumes.   For droplet initial diameters higher than 50 µm, numerical spray closely resembles the actual one. Small differences can be noticed considering the overall cone angle. However, this is mainly due to the fact that numerical images account for all the droplets, even the smallest ones. On the contrary, the experimental acquisition system is characterized by a mass threshold of 1%, so the smallest droplets are excluded. It is useful to note that, not by chance, the collapsing tendency of the spray with small initial diameters closely resembles the behavior of a spray under flash boiling conditions. Under these conditions, which are more and more investigated by experiments, it is widely demonstrated that droplets exiting from the nozzle are characterized by diameters much smaller than those under non-flashing conditions [40].
Even numerical-experimental comparison in terms of droplet sizing (SMD) confirms the importance of larger diameters. For the sake of brevity and validity of the experimental data, comparison is carried out only at three locations. Figures 19-21 show comparisons on a plane 15 mm far from injector tip at different radial locations (10, 11 and 12 mm). It is worthwhile to specify that For droplet initial diameters higher than 50 µm, numerical spray closely resembles the actual one. Small differences can be noticed considering the overall cone angle. However, this is mainly due to the fact that numerical images account for all the droplets, even the smallest ones. On the contrary, the experimental acquisition system is characterized by a mass threshold of 1%, so the smallest droplets are excluded. It is useful to note that, not by chance, the collapsing tendency of the spray with small initial diameters closely resembles the behavior of a spray under flash boiling conditions. Under these conditions, which are more and more investigated by experiments, it is widely demonstrated that droplets exiting from the nozzle are characterized by diameters much smaller than those under non-flashing conditions [40].
Even numerical-experimental comparison in terms of droplet sizing (SMD) confirms the importance of larger diameters. For the sake of brevity and validity of the experimental data, comparison is carried out only at three locations. Figures 19-21 show comparisons on a plane 15 mm far from injector tip at different radial locations (10, 11 and 12 mm). It is worthwhile to specify that numerical SMDs are calculated using the following expression: D i is the diameter of the i-th droplet passing through the measurement point during the whole injection process.
to the fact that numerical images account for all the droplets, even the smallest ones. On the contrary, the experimental acquisition system is characterized by a mass threshold of 1%, so the smallest droplets are excluded. It is useful to note that, not by chance, the collapsing tendency of the spray with small initial diameters closely resembles the behavior of a spray under flash boiling conditions. Under these conditions, which are more and more investigated by experiments, it is widely demonstrated that droplets exiting from the nozzle are characterized by diameters much smaller than those under non-flashing conditions [40]. Even numerical-experimental comparison in terms of droplet sizing (SMD) confirms the importance of larger diameters. For the sake of brevity and validity of the experimental data, comparison is carried out only at three locations. Figures 19-21 show comparisons on a plane 15 mm far from injector tip at different radial locations (10, 11 and 12 mm). It is worthwhile to specify that numerical SMDs are calculated using the following expression: is the diameter of the i-th droplet passing through the measurement point during the whole injection process.
It is clearly visible that only larger initial droplet diameters allow to match experimental values at measurement stations. On the contrary, small sizing leads to a not acceptable under-estimation of the experimental data.   Before moving to the inner-nozzle flow simulation, a few considerations are drawn on the importance of secondary break-up. In fact, the same collapsing effect obtained with reduced initial It is clearly visible that only larger initial droplet diameters allow to match experimental values at measurement stations. On the contrary, small sizing leads to a not acceptable under-estimation of the experimental data.
Before moving to the inner-nozzle flow simulation, a few considerations are drawn on the importance of secondary break-up. In fact, the same collapsing effect obtained with reduced initial droplet diameters can be achieved by means of larger diameters along with an enhanced secondary break-up. As an example, initial droplet diameter equal to 80 µm is considered. Secondary break-up model is switched from Reitz's one to KHRT. This latter provides a faster and more effective secondary break-up. Therefore, as shown in Figure 22, plume collapsing and a rapid increase of the liquid length at the end of the injection are noticed, similarly to the previous case with Reitz's model for the secondary break-up and 10 µm as initial droplet diameter.   Before moving to the inner-nozzle flow simulation, a few considerations are drawn on the importance of secondary break-up. In fact, the same collapsing effect obtained with reduced initial droplet diameters can be achieved by means of larger diameters along with an enhanced secondary break-up. As an example, initial droplet diameter equal to 80 µm is considered. Secondary break-up model is switched from Reitz's one to KHRT. This latter provides a faster and more effective secondary break-up. Therefore, as shown in Figure 22, plume collapsing and a rapid increase of the liquid length at the end of the injection are noticed, similarly to the previous case with Reitz's model for the secondary break-up and 10 µm as initial droplet diameter. The same can be stated for SMD; with higher initial diameters but a promoted break-up, numerical droplet sizing remarkably under-estimate experimental data, as visible in Figure 23.  Before moving to the inner-nozzle flow simulation, a few considerations are drawn on the importance of secondary break-up. In fact, the same collapsing effect obtained with reduced initial droplet diameters can be achieved by means of larger diameters along with an enhanced secondary break-up. As an example, initial droplet diameter equal to 80 µm is considered. Secondary break-up model is switched from Reitz's one to KHRT. This latter provides a faster and more effective secondary break-up. Therefore, as shown in Figure 22, plume collapsing and a rapid increase of the liquid length at the end of the injection are noticed, similarly to the previous case with Reitz's model for the secondary break-up and 10 µm as initial droplet diameter. The same can be stated for SMD; with higher initial diameters but a promoted break-up, numerical droplet sizing remarkably under-estimate experimental data, as visible in Figure 23. The same can be stated for SMD; with higher initial diameters but a promoted break-up, numerical droplet sizing remarkably under-estimate experimental data, as visible in Figure 23. Therefore, besides a proper atomization of the liquid column (which is the main focus of the present activity), a reliable secondary break-up model is also mandatory to match experiments.
In light of the results shown above, the diffused approach based on the adoption of small diameters to initialize droplets seems to be acceptable (even if wrong) only with reduced backpressures (i.e., nearly equal to ambient pressure), while at higher backpressures such an approach provides numerical results extremely different from experimental outcomes. As a further confirmation of such evidence, an inner-nozzle flow simulation on the Spray G is proposed hereafter, able to provide information about the liquid column exiting from the injector nozzle. It is mandatory to point out that the adoption of droplet diameters comparable to injector holes cannot be considered as a universal approach. Even if correct for most of the industrial applications (such as gasoline and Diesel injectors under standard internal combustion engine operations), there are some situations in Therefore, besides a proper atomization of the liquid column (which is the main focus of the present activity), a reliable secondary break-up model is also mandatory to match experiments.
In light of the results shown above, the diffused approach based on the adoption of small diameters to initialize droplets seems to be acceptable (even if wrong) only with reduced backpressures (i.e., nearly equal to ambient pressure), while at higher backpressures such an approach provides numerical results extremely different from experimental outcomes. As a further confirmation of such evidence, an inner-nozzle flow simulation on the Spray G is proposed hereafter, able to provide information about the liquid column exiting from the injector nozzle. It is mandatory to point out that the adoption of droplet diameters comparable to injector holes cannot be considered as a universal approach. Even if correct for most of the industrial applications (such as gasoline and Diesel injectors under standard internal combustion engine operations), there are some situations in which initial droplet diameters have to be carefully investigated (such as in the case of gasoline injectors under flash boiling conditions).

Inner-Nozzle Flow Simulation Setup
Inner-nozzle flow simulations represent a crucial aspect of spray modeling, since they are able to provide essential information for the droplet initialization, such as diameter and velocity of the liquid column. In this activity, an inner-nozzle simulation is carried out on the Spray G injector at the operating condition reported in Table 1. This is done to infer information about the droplet initial diameter for the validation of the approach proposed in the previous paragraphs. Among the different CAD files available on the ECN website, the so called "Generation 1" is considered, which represents the nominal geometry of the Spray G injector, with nozzle diameters equal to 165 µm. Static condition of the injector (namely maximum needle lift equal to 45 µm) is simulated [25].
As for the computational mesh, polyhedral cells are adopted in the core region, while prism layers are preferred at the walls, to improve near-wall modeling. Mesh refinements are locally introduced, such as in the needle seat region. Minimum and maximum cell sizes are 10 µm and 50 µm, respectively [41]. The total number of fluid cells is nearly seven million. The computational domain and related numerical grid are reported in Figure 24. In order to prevent numerical instabilities, properties (such as density and viscosity) are constant and calculated at the injection temperature. Boundary conditions are specified in Figure 25. As for the turbulence modeling, the two-equation K-Omega SST model [44] is adopted to close the set of Reynolds-Averaged Navier-Stokes equations. The so called "All-" wall treatment is preferred for the near-wall modeling [45][46][47][48]. This model relies on a hybrid approach able to work as a Low-Reynolds wall treatment for values belonging to the viscous sub-layer and as a High- Concerning the CFD setup, an unsteady Eulerian multiphase approach is adopted, along with the VOF (Volume of Fluid) model [42], which is able to predict the distribution of the interface between immiscible phases. Phase distributions and interface position are described by phase volume fractions αi, defined as: where V i is the volume of i-th phase in the cell and V is the cell volume. The Schnerr-Sauer's cavitation model is used to predict vapor formation [43]. In order to prevent numerical instabilities, properties (such as density and viscosity) are constant and calculated at the injection temperature. Boundary conditions are specified in Figure 25. In order to prevent numerical instabilities, properties (such as density and viscosity) are constant and calculated at the injection temperature. Boundary conditions are specified in Figure 25. As for the turbulence modeling, the two-equation K-Omega SST model [44] is adopted to close the set of Reynolds-Averaged Navier-Stokes equations. The so called "All-" wall treatment is preferred for the near-wall modeling [45][46][47][48]. This model relies on a hybrid approach able to work as a Low-Reynolds wall treatment for values belonging to the viscous sub-layer and as a High-Reynolds one for values pertaining to the fully turbulent region. Experimental available data from ECN [26] include hydraulic coefficients of the injector (Ca, Cv and Cd), total mass flow rate and plume centroid locations on a plane located 50 mm far from the tip and normal to the injector axis.
In order to compare numerical and experimental outcomes, sections are purposely created at hole exits (before the hole steps), as reported in Figure 26. Such sections are used to measure quantities such as mass flow rate and hydraulic coefficients.  As for the turbulence modeling, the two-equation K-Omega SST model [44] is adopted to close the set of Reynolds-Averaged Navier-Stokes equations. The so called "All-y + " wall treatment is preferred for the near-wall modeling [45][46][47][48]. This model relies on a hybrid approach able to work as a Low-Reynolds wall treatment for y + values belonging to the viscous sub-layer and as a High-Reynolds one for y + values pertaining to the fully turbulent region.
Experimental available data from ECN [26] include hydraulic coefficients of the injector (C a , C v and C d ), total mass flow rate and plume centroid locations on a plane located 50 mm far from the tip and normal to the injector axis.
In order to compare numerical and experimental outcomes, sections are purposely created at hole exits (before the hole steps), as reported in Figure 26. Such sections are used to measure quantities such as mass flow rate and hydraulic coefficients. In order to prevent numerical instabilities, properties (such as density and viscosity) are constant and calculated at the injection temperature. Boundary conditions are specified in Figure 25. As for the turbulence modeling, the two-equation K-Omega SST model [44] is adopted to close the set of Reynolds-Averaged Navier-Stokes equations. The so called "All-" wall treatment is preferred for the near-wall modeling [45][46][47][48]. This model relies on a hybrid approach able to work as a Low-Reynolds wall treatment for values belonging to the viscous sub-layer and as a High-Reynolds one for values pertaining to the fully turbulent region. Experimental available data from ECN [26] include hydraulic coefficients of the injector (Ca, Cv and Cd), total mass flow rate and plume centroid locations on a plane located 50 mm far from the tip and normal to the injector axis.
In order to compare numerical and experimental outcomes, sections are purposely created at hole exits (before the hole steps), as reported in Figure 26. Such sections are used to measure quantities such as mass flow rate and hydraulic coefficients.  A dedicated methodology is developed for the evaluation of the hydraulic coefficients. For each single nozzle, C d is evaluated considering both actual and theoretical mass flow rates, the first being measured at the nozzle exit section. Ca is calculated at the same section as surface average of the liquid volume fraction. Finally, C v is derived as ratio between C d and C a , so that the effective injection velocity can be obtained. Moreover, plume directions are obtained by means of the resulting velocity vectors at the nozzle exits. In particular, given the (red colored) coordinate system reported in Figure 27, components (v x , v y and v z ) of such vectors are calculated as surface averages on the portions of the nozzle exit sections where velocity is directed outside. Figure 27 shows that, once components are known, angles α x and α y can be calculated for each nozzle as in Equations (6) and (7), providing the plume directions in the global reference system.
liquid volume fraction. Finally, Cv is derived as ratio between Cd and Ca, so that the effective injection velocity can be obtained. Moreover, plume directions are obtained by means of the resulting velocity vectors at the nozzle exits. In particular, given the (red colored) coordinate system reported in Figure 27, components ( , and ) of such vectors are calculated as surface averages on the portions of the nozzle exit sections where velocity is directed outside. Figure 27 shows that, once components are known, angles and can be calculated for each nozzle as in Equations (6) and (7), providing the plume directions in the global reference system. = atan = atan ( , ) ( , ) = atan = atan ( , ) ( , ) In Equations (6) and (7), , is the x-component of the velocity in the global reference system of the i-th cell of the nozzle exiting section.
It is useful to point out that numerical results are to be intended as time average quantities. Calculation of the velocity vector direction is depicted in Figure 27.

Inner-Nozzle Flow Simulation Results
Comparisons between experimental and numerical outcomes for Spray G injector start with the total injected mass flow rate, as reported in Figure 28.  In Equations (6) and (7), v i,x is the x-component of the velocity in the global reference system of the i-th cell of the nozzle exiting section.
It is useful to point out that numerical results are to be intended as time average quantities. Calculation of the velocity vector direction is depicted in Figure 27.

Inner-Nozzle Flow Simulation Results
Comparisons between experimental and numerical outcomes for Spray G injector start with the total injected mass flow rate, as reported in Figure 28.
A dedicated methodology is developed for the evaluation of the hydraulic coefficients. For each single nozzle, Cd is evaluated considering both actual and theoretical mass flow rates, the first being measured at the nozzle exit section. Ca is calculated at the same section as surface average of the liquid volume fraction. Finally, Cv is derived as ratio between Cd and Ca, so that the effective injection velocity can be obtained. Moreover, plume directions are obtained by means of the resulting velocity vectors at the nozzle exits. In particular, given the (red colored) coordinate system reported in Figure 27, components ( , and ) of such vectors are calculated as surface averages on the portions of the nozzle exit sections where velocity is directed outside. Figure 27 shows that, once components are known, angles and can be calculated for each nozzle as in Equations (6) and (7), providing the plume directions in the global reference system.
= atan = atan ( , ) ( , ) In Equations (6) and (7), , is the x-component of the velocity in the global reference system of the i-th cell of the nozzle exiting section.
It is useful to point out that numerical results are to be intended as time average quantities. Calculation of the velocity vector direction is depicted in Figure 27.

Inner-Nozzle Flow Simulation Results
Comparisons between experimental and numerical outcomes for Spray G injector start with the total injected mass flow rate, as reported in Figure 28.  Numerical mass flow rate closely reproduces the static value detected in the experiments, i.e., approximately 14 g/s. Moving to the hydraulic coefficients, Figure 29 shows the numerical discharge coefficients of each individual injector hole, their average and its corresponding experimental value. A similar comparison is proposed also for C a and C v in Figures 30 and 31, respectively. Numerical mass flow rate closely reproduces the static value detected in the experiments, i.e., approximately 14 g/s. Moving to the hydraulic coefficients, Figure 29 shows the numerical discharge coefficients of each individual injector hole, their average and its corresponding experimental value. A similar comparison is proposed also for Ca and Cv in Figures 30 and 31, respectively.   Thanks to the agreement in terms of global mass flow rate shown Figure 28, numerical mean Cd is significantly close to its experimental counterpart. Similar satisfactory results are obtained also for Cv and Ca; differences between numerical and experimental data are lower than 3%. A Ca equal to nearly 0.7 means that part of the section area is not exploited by the flow. Figure 32 confirm this statement since only a portion of the exiting sections is occupied by liquid. Similarly, only a portion of the area is characterized by positive axial velocity (i.e., exiting flow), as it is visible in Figure 33. Numerical mass flow rate closely reproduces the static value detected in the experiments, i.e., approximately 14 g/s. Moving to the hydraulic coefficients, Figure 29 shows the numerical discharge coefficients of each individual injector hole, their average and its corresponding experimental value. A similar comparison is proposed also for Ca and Cv in Figures 30 and 31, respectively.   Thanks to the agreement in terms of global mass flow rate shown Figure 28, numerical mean Cd is significantly close to its experimental counterpart. Similar satisfactory results are obtained also for Cv and Ca; differences between numerical and experimental data are lower than 3%. A Ca equal to nearly 0.7 means that part of the section area is not exploited by the flow. Figure 32 confirm this statement since only a portion of the exiting sections is occupied by liquid. Similarly, only a portion of the area is characterized by positive axial velocity (i.e., exiting flow), as it is visible in Figure 33. Numerical mass flow rate closely reproduces the static value detected in the experiments, i.e., approximately 14 g/s. Moving to the hydraulic coefficients, Figure 29 shows the numerical discharge coefficients of each individual injector hole, their average and its corresponding experimental value. A similar comparison is proposed also for Ca and Cv in Figures 30 and 31, respectively.   Thanks to the agreement in terms of global mass flow rate shown Figure 28, numerical mean Cd is significantly close to its experimental counterpart. Similar satisfactory results are obtained also for Cv and Ca; differences between numerical and experimental data are lower than 3%. A Ca equal to nearly 0.7 means that part of the section area is not exploited by the flow. Figure 32 confirm this statement since only a portion of the exiting sections is occupied by liquid. Similarly, only a portion of the area is characterized by positive axial velocity (i.e., exiting flow), as it is visible in Figure 33. Thanks to the agreement in terms of global mass flow rate shown Figure 28, numerical mean C d is significantly close to its experimental counterpart. Similar satisfactory results are obtained also for C v and C a ; differences between numerical and experimental data are lower than 3%. A Ca equal to nearly 0.7 means that part of the section area is not exploited by the flow. Figure 32 confirm this statement since only a portion of the exiting sections is occupied by liquid. Similarly, only a portion of the area is characterized by positive axial velocity (i.e., exiting flow), as it is visible in Figure 33.           Finally, thanks to both C a and C v , it is possible to obtain essential information for the droplet initialization (i.e., velocity and diameter).
Such values are reported in Table 4. The most important numerical outcome for the present analysis is the estimation of the initial droplet diameter. This last is higher than 130 µm and, in the light of a nominal hole diameter of 165 µm, the approach proposed in the previous paragraphs, i.e., the choice of an initial droplet diameter comparable to the injector hole, seems to be solid. Therefore, trying to generalize the results of the present activity, compared to poorly predictive atomization models or droplet distribution functions inferred from far-field experiments, a simpler blob model based on the hole diameter may be able to provide 3D-CFD numerical results in much better agreement with the experimental counterpart. However, it is necessary to point out that such conclusion cannot be considered of universal validity. For example, as mentioned before, flash boiling conditions are characterized by diameters much smaller than ones under non-flashing conditions [50][51][52].
It is useful to point out that diameters reported in Table 4 are slightly larger than the largest one tested in the previous paragraphs. The adoption of such a large diameter (nearly 130 µm) or even to use of a blob model based on the hole size (165 µm) may require a small tuning of the Reitz's model to emphasize secondary break-up and, in turn, to match experimental data.

Conclusions
In the present paper, the impact of the primary break-up modeling strategy in GDI sprays is addressed. Atomization modeling in common industrial practice often relies on the adoption of droplet diameter distribution functions inherited from experimental PDA data. However, for the sake of validity of the latter, measurements are usually carried out at least 15 ÷ 20 mm away from the injector tip, where droplets have already undergone the secondary break-up. Therefore, droplet initial diameters adopted in Lagrangian simulations are characterized by values even one order of magnitude smaller than nozzle hole diameters. Even if questionable, such an approach is able to provide numerical results in line with experimental ones if the backpressure is close to the ambient one. For this purpose, Lagrangian simulations at ambient conditions and different injection pressures are carried out on a 5-hole GDI prototype injector. Such simulations are characterized by initial droplet diameters lower than 10 µm and inherited from experimental PDA data. The same atomization strategy is also applied to the well-known Spray G injector. The analyzed operating condition is the most investigated one in literature and it is characterized by a backpressure equal to 0.6 MPa. In this case, the adoption of small diameters (nearly 10 µm) coming from tests is responsible for relevant misalignments with reference to the experimental outcomes. In fact, unlike experimental evidences, the numerical spray tends to collapse similarly to flash boiling conditions. Therefore, larger droplet diameters (closer to the hole size of 165 µm) are tested, able to provide a representation of the spray in line with the experiments.
In order to confirm the importance of choosing droplet diameters comparable with the hole size, an internal nozzle flow simulation of the Spray G is carried out. Such a simulation provides relevant information on droplet initial conditions to be applied in the Lagrangian simulations, particularly in terms of droplet diameters. These last are found to be slightly larger than 130 µm, thus similar to nozzle hole dimension and one order of magnitude larger than values provided by PDA measurements 15 mm away from the injector tip. Therefore, compared to droplet distribution functions inferred from experiments (or even poorly predictive atomization models), a simple blob model with droplet size equal to hole diameter may better perform in terms of 3D-CFD numerical results. Acknowledgments: Authors would like to thank both SprayLAB of Perugia University and ECN for experimental data on 5-hole prototypal injector and Spray G injector respectively. Moreover, special thanks are reserved also to SIEMENS PLM for granting licenses.

Conflicts of Interest:
The authors declare no conflict of interest.