Large Eddy Simulations of the fuel injection and mixing process of the ECN GDi Injector Spray G

Understanding unsteady effects of the Gasoline Direct injection (GDi) process, due to short injection duration, plays a major role in the analysis of the mixture formation and so combustion efﬁciency in the spark ignition engines. Focusing on the phase when the needle is at its max-imum lift, there still are some uncertainties. For instance, the time evolution of the upstream pressure, together with the detailed geometry of the needle and the seat, shape of the rate of injection could affect the uniformity of the ﬂow between oriﬁces. Experimentally addressing these issues nowadays remains challenging, if not impossible, so Computational Fluid Dynamics (CFD) are used. Homogeneous Relaxation Model (HRM) is employed to consider the mass exchange between liquid and vapor phases of the fuel inside the nozzle. Different Large Eddy Simulations (LES) sub-grid models are employed in order to account for the unsteady effects of turbulence and analyze the predictive capabilities in resolving the spatial-temporal scales. Results are validated against experimental data. Vortices are generated inside the counter-bore and enhance the liquid dispersion near the injector. Spray parameters, after a proper time-averaging, also accurately match the results reported in the literature. Root mean square velocity (Urms) ﬂuctuation has higher levels of spatially located ﬂuctuations in areas with signiﬁcant levels of small-scale turbulence and noticeable liquid-vapor mixture. tion (LES) models. This work is based on the comparison of two different sub-grid models (both one-equation). On the one hand the Dynamic Structure, which achieves good approximations for engine applications, and on the other hand, the Viscous One-Equation model, which allows coarser meshes. These models add one additional transport equation for sub-grid kinetic energy ( k sgs ). For both approaches, the turbulent boundary layer is modeled using the Werner and Wengle wall function.


Introduction
The evolution of the automotive industry has been part of the causes of the climate change and therefore, stringent emission regulations that have emerged. Conventional engine operation begins when the fuel is introduced into the combustion chamber. The way in which it is delivered directly affects the combustion process and the pollutants generated [1]. Within the concept of cleaner and more efficient transport, the importance of liquid atomization and mixture formation in reducing emissions and achieving high efficiency in spark ignited (SI) engines is highlighted. For this purpose, the technology of Gasoline Direct injection (GDi), which stands out in terms of reducing harmful emissions, is gaining interest and is the main focus of the article. Given the background, the Engine Combustion Network (ECN), an international network of researchers working in the field of sprays and combustion in Internal Combustion Engines (ICEs), made available a multi-orifice injector known as Spray G. The aim of this group is to share the findings of their studies and facilitate the understanding of the phenomena that take place in the injection process and that significantly affect the engine performance. A great deal of effort has been dedicated to studying the geometry of this injector in detail, in particular the work by Matusik et al. [2] who obtained the main characteristic parameters of the injector such as hole diameter, hole inlet corner radius and counter-bore diameter using x-rays and linear regression techniques. The distribution of the orifices and the arrangement of 5 bumps along the nozzle favor the hole-to-hole variation and the interaction between sprays [3]. It is also worth mentioning the efforts made to estimate certain macroscopic and microscopic parameters such as spray shape, jet penetration, axial velocity at the centre of the discharge volume and particle diameter, among others [4].
Researchers have always worked in parallel by experimentally studying injection behavior and computationally analyzing the phenomena observed. In the case of the near-field spray, close to the injector tip, the experiments have limited analytical capabilities due to the complexity of the geometry, commonly designed as multi-hole valve covered orifice (VCO), and the small size of the orifices (in the order of a few hundred microns). For this reason, Computational Fluid Dynamic (CFD) tools have acquired importance over the years and have become a powerful resource for the study of the in-and near-nozzle regions. However, CFD analysis can go even further in terms of accuracy, which is why the use of LES models is becoming more and more widespread [5]. During the injection process the spray is influenced by numerous phenomena such as flow oscillations, vortex structures, shear stresses or cavitation resulting in a faster atomization of the spray. Due to the importance of this behaviour, the transient parts of the injection, opening and closing phases, are currently of scientific interest [6]. Despite the relevance of the transient phases, observation of the steady-state phenomenon is still useful to understand the underlying physics, what is typically done by using a long pulse duration and thus focusing the analysis on the quasi-stable main injection phase. Therefore, the motivation of this article is to study the characteristics of the flow during the steady state of the injection process in the ECN Spray G injector at nominal operating conditions. For this purpose, LES approximations will be used to model the turbulence. Two different sub-grid models will be analyzed and compared. The first part of the study will involve the quality of the mesh which is directly related to the accuracy of the results. The behavior of the injector will be evaluated in terms of macroscopic variables such as mass flow rate and momentum flux and compared with experimental results from previous studies of the group available in the literature. The second part of this study will analyze the deviation in the plume direction in comparison with the geometric angle (drill angle) and the spray angle using a new processing methodology.

Numerical Methods
For the present study, the code used is CONVERGE v2.4 from the Convergent Science group. The analysis of the multi-phase fluid inside and near the nozzle has been carried out within an Eulerian framework using a single-fluid approach governed by the classical conservation equations for mass, momentum and energy [7]. The overall set of equations are solved through the finite volume method. A Pressure Implicit with Splitting of Operators (PISO) is selected to iteratively solve the transport equations. In order to prevent the decoupling of the pressure-velocity, the Rhie-Chow algorithm is activated. First-order upwind was selected as discretization scheme for computing the convection flux in density, energy, species and passives transport equations, while second-order central difference discretization scheme is used for the momentum. In an attempt to promote better stability, the successive over-relaxation (SOR) algorithm has been used. The time-step is controlled by the Courant-Friedrichs-Lewy (CFL) and results in the range between 5 · 10 −7 and 10 −10 s. Values for velocity-based Courant-Friedrichs-Lewy (CFL) are below 0.5 and for viscosity-based CFD below 2.5. This research encompasses the presence of three different species: the liquid fuel, the vapor fuel and the non-condensable ambient gas (N 2 ). When solving for these species, first the species mass fraction transport equation is employed, and then the void fraction is calculated, which means that it is not directly transported. The two phase (liquid and gas) flow in the domain is simulated through an interface-capturing Volume-of-Fluid method which is explained in detail in the work from Battistoni et at. [8]. A non-equilibrium Homogeneous Relaxation Model (HRM) is selected for the mass exchange between the liquid and vapor phases of the same species. This approach accurately predict not only cavitation but also flash-boiling [9]. The model is setup following the recommendations of Saha et al. presented on the basis of the sensitivity analysis developed [10]. Regarding the turbulence, the Navier-Stokes equation contains an unclosed non-linear term which is modeled decomposing the fields into resolved and sub-grid using Large Eddy Simula-tion (LES) models. This work is based on the comparison of two different sub-grid models (both one-equation). On the one hand the Dynamic Structure, which achieves good approximations for engine applications, and on the other hand, the Viscous One-Equation model, which allows coarser meshes. These models add one additional transport equation for sub-grid kinetic energy (k sgs ). For both approaches, the turbulent boundary layer is modeled using the Werner and Wengle wall function.

Geometry Details, Computational Domain and Operating Conditions
As mentioned before, this research is focused on the study of the behavior of GDi injectors. In order to participate and contribute to the creation of a broad benchmark in direct injection systems, the analyzed geometry is the one belonging to the ECN and known as Spray G [11]. The Spray G injector is a solenoid type designed with a VCO nozzle and characterized by 8 counterbore orifices and 5 needle guides on the nozzle wall as Figure 1 illustrates. Notwithstanding the aforementioned, the injector is considered symmetric in the plane A-A that includes holes 1 and 5, a feature used in this case to reduce the computational cost. The analyzed geometry has been obtained from x-rays together with their corresponding post-processing. It stands out for being an ideal geometry in which manufacturing defects have been eliminated and for being the first geometry model developed (Generation 1) [3]. The outlet plenum is defined as semispherical with a diameter of 6 mm, which is large enough to avoid the influence of the boundary conditions on the solution. The main features of this injector are its nominal diameters of 165 µm and 388 µm respectively, the ratio between the length and diameter of the counter-bored nozzle holes of 1 to 1.2 and the inclination of 37 • of the geometrical axis of the holes. Simulations carried out in the analysis are based on the operating condition proposed by ECN and named after the injector, Spray G (explained in detail in [12]). Concerning the initialization and boundary conditions, same procedure presented in previous works has been applied in this study [12]. The only difference relates to the initialization of the turbulence where the use of LES models, one equation in this particular case, allows to initialize only the turbulent kinetic energy defined through the turbulent intensity with a value of 0.01.  The objective of this study is the steady state of the injection process whose variation in lift and wobble profile is 3 microns. Nevertheless, the needle movement and wobble during this stage of the injector operation have been introduced from the complete lift profile obtained by Argonne National Laboratory [3]. The modeling of turbulence from LES approximations is highly dependent on cell size. For this reason, the present study uses a fixed mesh whose sizes have been based on the mesh quality criteria explained in the next section. The mesh selected for this approach, as illustrate in Figure 2, is a cartesian grid with hexahedral elements. The use of fixed embedding allows a minimum size of 4.22 µm to be defined in the seat from a base of 135 µm.

LES quality assesment
By the use of LES models for turbulence modeling, the numerical method employed to model the small scales and the resolution of the mesh influence the turbulent resolution. In particular, the grid resolution not only affects the contribution of the sub-grid scale but also greatly influences the numerical discretization error. Therefore, these models require a quality assessment to ensure adequate resolution of the turbulent flow energy and thus accuracy of the results. From this necessity arises the defined quality index in terms of both numerical and model accuracy such as the proposed by Celik et al. [13], one of the most widespread on the basis of viscosity. The criteria evaluates the contribution relative to the laminar viscosity (υ), the subgrid viscosity (υ sgs ) and the numerical viscosity (υ num ) according to Equation 1. It is considered for High-Reynolds-number flows that the resolution of 75% of the turbulent kinetic energy is enough to ensure a good quality of the LES analysis.
The constants α υ = 0.05 and n = 0.53 are grounded in DNS outcomes [13]. The estimation of such quality index requires turbulent statistics, therefore time-averaged variables will be considered. The simulated time ranges from 0.15 to 0.51 ms (360 µm). For this study, 19 residence times with an acquisition frequency of 1e6 Hz have been used.   Figure 3 presents the results obtained from the mesh quality analysis for the two sub-grid models analyzed. It can be observed how the Dynamic Structure model has a much lower mesh quality, being well below the ranges established as optimal (75% of the energy resolved). On the contrary and mentioning its capacity of favoring the analysis of coarser meshes, the sub-grid Viscosity model meets the quality criteria stipulated by Celik et al. Apart from the quality index, the difference between models in the response with respect to the computational cost was analyzed. Although both models resolve an additional equation to compute the subgrid kinetic energy, the Viscous One-Equation model, that uses the turbulent viscosity to model the sub-grid tensor, has higher CPU cost than the Dynamic Structure approach, 48141 CPU-h compared to 39917 CPU-h.

Results
This section summarizes the results obtained from the study of the near field of the proposed injector and compares them with experimental data from previous studies carried out at CMT-Motores Térmicos. It is important to note that the computationally obtained mass flow rate profiles have been generated by measuring the value in perpendicular planes to the counterbored orifices. In order to follow the same experimental methodology, the momentum flux has been measured at a distance downstream, taking as a reference the injector tip.
Respecting to the mass flow results, in the stationary part of the injection, the changes in the values are minor and mainly due to the presence of vortices inside the counter-bore. Although  2.10 ± 0.062 2.14 ± 0.057 the quality of the mesh of both approached was different and even the Dynamic Structure model was far below the selected criterion, Figure 4 displays similar trends. If zoomed in as shown in the image, a slight difference between models is visible at certain instants of the stationary, leading to a 2% variance (Table 1). Quality compliance makes the Vicous One-Equation model more reliable. Nevertheless, these results tend to overpredict the experimental results by around 9%. In order to explain this effect, the time-averaged mass flow rate of each orifice has been calculated and presented in Table 1. The experimental stationary ROI value is 14.71 g/s, assuming 1.84 g/s for each hole. Comparing this result with that presented in the table, it is observed minimal differences between experimental and computational in orifices 2, 3 and 4, with the greatest discrepancy in orifices 1 and 5. These are the holes where symmetry has been applied which means that the flow condition is not correctly capturing the behavior of the fluid, resulting in a higher mass flow rate and therefore an overestimation of the results. With reference to momentum flux, displayed in Table 1, an overestimation of about 6% is also found. Due to the computational discharge volume is 6 mm in diameter, the momentum measurements have been made in a plane at 1 mm while the experimental data are collected at 3 mm. Although momentum rate is supposed to be conserved, previous work by this group [14] has shown how the data acquisition distance affected the results. This fact coupled with the symmetry condition may be the main reasons for the disparity in results. The root-mean-square (RMS) velocity is one of the leading indicators of the resolved energy and is expressed by the use of the resolved instantaneous velocity vectorũ and its timeaveraged value noted <ũ >:ũ rms = <ũũ > − <ũ > 2 . Figure 6 represents the Urms comparison between sub-grid models. The central part of the spray near the nozzle does not indicate the presence of Urms and refers to the area where the liquid remains intact. As the flux moves away from the injector outlet, the liquid mixes with the ambient gas producing turbulent phenomena which in turn increase the mixing and atomization of the spray. In the Dynamic Structure approach, the spray begins to interact with the environment a few millimeters after the orifice exit whereas, for the Viscous One-Equation model, this intact part extends downstream to near the outlet. The Urms has higher levels of spatially located velocity fluctuations in areas which are characterized by higher levels of small-scale turbulence and where the liquid-vapor mixture is noticeable. The Viscous One-Equation has larger oscillations compared to the Dynamic Structure generating more mixing and resolving a larger amount of energy in the inner part of the spray. The opposite effect is observed on the external side. The cone angle can be defined as the jet opening measured from nozzle tip to some distance downstream of the flow while the plume direction is the angle between the vertical injector axis and the spray axis ( Figure 5). These are significant parameters and determinants of the results accuracy when studying the external flow from, for example, Droplet Discrete Model (DDM) approximations. To date, and due to the lack of experimental data because of the complexity of the calculation, the external flow studies have been carried out starting from an initial ad hoc value for the cone angle and calibrate it to match experimental data, thus leading to incorrect approximations. This calculation has already been carried out in previous works of this group [12], although the used methodology was different (explained below). The new technique avoid dependence on image quality and take into account the influence of geometry on jet development [9]. Therefore, this methodology has been applied to each of the orifices separately. The starting point for the angle calculation is to choose a variable (velocity, Liquid Volume Fraction, etc.) that will define the spray itself. For the present analysis, the liquid mass fraction has been used. When the liquid comes into contact with the environment it starts to mix, so that the boundaries of the spray are a mixture of liquid and gas. For this reason, it is necessary to choose a threshold that determines the limits of the spray being in this case a mass fraction greater than 0.25. An upper limit (z = 0.6 mm) and a lower limit (z = 1.6 mm) have been chosen to define the sample space, avoiding the influence of the outlet, if any ( Figure 5). Once the boundaries are determined, the planes that cut the holes through the center of each hole are defined. Then, the domain is vertically discretized with 20 microns samples to obtain the limits of the spray. For each of these limits, the slope is calculated taking as a reference point the center of the inner hole as defined in the ECN angle specification. Once all the slopes have been calculated, they are averaged to obtain the spray angle as final value. The spray direction is determined in the same way by making, at each point, the bisector of the spray angle and averaging them at the end.   Figure 7 provides the results obtained from the calculation of the cone angle and plume direction for both sub-grid models. Both approaches, Viscous One-Equation and Dynamic Structure, show a slight variability of the values for hole 3 while the other two, holes 2 and 4, have more differences. On the one hand, the orientation of the mesh, since hole 3 is oriented with the cartesian mesh, could influence the results as seen in other works [15]. On the other hand, the studied geometry has 5 bumps and 8 orifices, so the flow distribution is not the same in all of them. Therefore, the turbulence generated upstream of the injector may affect some directions more than others. The difference between orifices can also be seen in the direction of the spray. Thus, orifice 3 is deflected inwards, deviating from the geometrical axis of the orifice, while the other orifices show a minimal deviation. This deflection effect is displayed in Figure 5, where the counter-bore of orifice 2 (on the left) is practically full of liquid while orifice 3 (on the right) is largely occupied by non-condensable gas. This recirculation of air is the main reason of the deflection of the spray, causing it to deviate from the geometric axis of the orifice. The time-averaged and standard deviation cone angle and plume direction for both sub-grid models studied are plotted in Table 2. The results in both models are quite similar. A narrower spray is seen in hole 3 while holes 2 and 4 are around 20 • of cone. With respect to the plume direction, a small deviation of approximately 1-1.5 • is observed in holes 2 and 4 while hole 3, in both approaches, has more significant deviation. Predicted trends from both approaches estimate the liquid plume direction and cone angle close to literature values and contribute to the GDi spray behavior understanding.

Conclusions
The current research carried out a numerical study of the internal and near-nozzle flow for the multi-hole GDi injector from ECN, using two different LES approaches for the turbulence modeling. Macroscopic variables such as mass flow rate and momentum flux have been calculated and compared with data available in the literature. The spray angle and plume direction, essential parameters for understanding spray behavior, have also been obtained. The main findings of this report are summarized below, • Similar results were observed for ROI and ROM variables in both approaches, although the Dynamic Structure tends to provide slightly higher results than the Viscous One-Equation, namely 2%. Comparing these results with the experimental data an overprediction of 9% in ROI results and 6% in ROM values have been found. Most of this error comes from the fact that the applied symmetry condition does not correctly capture the fluid behavior. The distance of data acquisition can also affect the momentum flux results.
• Referring to the spray cone angle and plume direction, results have been shown as a function of time, something that, to the author's best knowledge, has not been studied to date. Mesh orientation or the asymmetric effect of upstream turbulence due to the distribution of bumps and holes in the geometry affects the spray wide as well as the direction.
Comparing the sub-grid models analyzed, both offer similar results and behavior.
• Applying the presented methodology to the analysis of the transient phases of the injection (opening and closing phases) could give interesting results. On the other hand, the calculation of the angle as a function of time could be used for coupling Eulerian-Lagrangian approximations, giving more accurate results of the spray behavior.