Numerical investigation on transient internal cavitating flow and spray characteristics in a single-hole diesel injector nozzle: A 3D method for cavitation-induced primary break-up

Fuel 233 (2018) 778–795

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Numerical investigation on transient internal cavitating flow and spray characteristics in a single-hole diesel injector nozzle: A 3D method for cavitation-induced primary break-up ⁎

T

⁎

Donghua Lia, Shenghua Liua, , Yanju Weia, , Renzhi Lianga, Yonghong Tangb a b

School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China Department of Fluid Machinery and Engineering, Xi’an Jiaotong University, Xi’an 710049, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Cavitation characteristics Standard method Critical stability position Primary break-up

The spray development and atomization are strongly influenced by the internal nozzle flow. In this study, a 3D cavitation-induced primary break-up method is proposed based on cavitation characteristics studies. The reason why to develop this advanced standard method and abandon the total Eulerian approach has been discussed. The transient cavitation characteristics inside a single-hole diesel injector nozzle were numerically investigated. Combined with previous experimental studies, a dynamic mesh method and its related numerical model were utilized to investigate the inner cavitating flow characteristics for cases of different maximum needle lift values, injection pressures and ambient pressures. It is found that the mass flow rate increases significantly with the increase of maximum lift position for low needle lift cases. The cavitating flow for low needle lift cases requires necessary simulations and the traditional standard method with only empirical equations is not enough. It is also found that for high needle lift cases, the internal flow keeps stable for over 90 percent of the total injection process. The 3D method aims to incorporate as much flow information as possible to avoid using the total Eulerian approach. The 3D method allows for nozzle cross section with both two-dimensional and three-dimensional data. The spray experiment in a constant chamber with different injection pressure and the same ambient pressure is employed to validate the 3D method. The 3D method with transient internal flow profiles is implemented into the AVL Fire CFD codes. The 3D method, together with the typical standard method (with or without 1D method) are applied to the experimental data, including spray tip penetration (STP), spray shape and droplet diameter. The results show that the 3D method has good predictive ability on the spray development and shape, and it has better performance on droplet diameter prediction.

1. Introduction Higher requirements for emissions are put forward to current diesel engines. The spray atomization is one of the most important parameters in combustion and emission characteristics for diesel engine. To get better performance and meet stringent regulations on emissions, it is necessary to give through investigations on the spray atomization of the injector. Previous studies [1–5] found that the internal flow inside the nozzle has a close connection with the spray characteristics. But the contribution of internal flow and ambient environment to spray breakup has not been measured or analyzed precisely. Additionally, the cavitation phenomenon in the internal flow make obscure researches conducted on nozzle geometry influence on spray characteristics because the cavitation phenomenon was neglected during the analysis. One

⁎

evidence may be that the spray angle at the instant the fuel comes out from the nozzle exit changes greatly for nozzles of different geometry. Arai [1–3] and Bode [4] found that the spray behavior changes dramatically with the presence of cavitation. They have verified that the spray angle increases significantly and the jet breakup length shortens when cavitation occurs. Zhifang Chen [5] has investigated the cavitation characteristics using a methanol injector and found that as the temperature increases, the cavitation region enlarges and the droplet Sauter mean diameter reduces. This is a sign that the cavitation has similar trends with spray effect. Due to lack of practical measurement equipment on the experiment of internal flow, the computational fluid dynamics (CFD) can be a preferable and useful tool to gain flow characteristics inside the injector nozzle. A combined experimental and computational study was carried out by F. Payri [6] and they have drawn the conclusion that the increase of injection velocity was

Corresponding authors at: School of Energy and Power Engineering, Xi’an Jiaotong University, Xianning West Road 28#, Xi’an 710049, China (Shenghua Liu). E-mail address: [email protected] (S. Liu).

https://doi.org/10.1016/j.fuel.2018.06.103 Received 29 March 2018; Received in revised form 6 June 2018; Accepted 25 June 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.

Fuel 233 (2018) 778–795

D. Li et al.

Nomenclature

RB ρv ρl ρ → Vv Pt P Pv αv R αnuc Fvap Fcond k ṁ pin pout Cc Cd Ca

Cv σ Ugeo Ueff Ageo A Aeff veff αi Ai vi u v w Dini a r ∧ Ω τ · ml i

Bubble radius Vapor phase density Liquid phase density Mixture density Vapor phase velocity Turbulent pressure Local pressure Saturation pressure Vapor volume fraction Phase change rate nucleation volume fraction of the bulk liquid Vaporization coefficient Condensation coefficient Turbulent kinetic energy Mass flow rate Inlet pressure Outlet pressure Contraction coefficient Discharge coefficient Area contraction coefficient

Velocity coefficient Cavitation number Theoretical velocity of a laminar flow Effective velocity Area of the nozzle exit cross section Effective area Effective velocity Diesel liquid fraction for the i-th type Area of the i-th type region Velocity of the i-th type parent drop Axis in the nozzle exit cross section plane Axis in the nozzle exit cross section plane Axis normal to the nozzle exit cross section plane Initial diameter Wave width (ligament width) New drop radius Wave length Wave growth rate Break-up time Liquid mass flow rate Parent type label

methods for spray modeling, namely, the standard method and the total Eulerian approach. The standard method is a Eulerian-Lagrangian approach. First the flow characteristics is obtained and then the Discrete Phase Model is employed for spray modeling based on the flow data at the nozzle exit region. Recent years the total Euler approach has been a hot topic for its main advantage of complete flow information incorporated [18–23]. Similar experimental study has also focus on the dense spray transition process [24,25]. F. Oley has used the total Eulerian approach for the cavitating nozzle flow and primary jet break-up with Large-Eddy simulation [18]. Similar way has also been conducted by Lebas using DNS method [20]. The total Euler approach with DNS method is much time consuming and this restricts its wide applications. Besides, these methods for total Eulerian approach cannot well involve the Sauter mean diameter (SMD) and this is important for spray atomization on diesel engines. Only turbulence-induced air entrainment or mixture distribution result is not enough. Berg [22] has proposed the total Eulerian multi-fluid approach by separating sets of conservation laws for each phase. The author aimed to gain the SMD result by assuming different levels of droplets as different phases that formed the continuous multi-phase fluid. This method can also bring assumption problem, because that the assumed droplets diameter may even approach the size of nozzle diameter. In the simplest approach the droplet size classes are fixed and may not change in space and time. But fixed size classes cause an underestimation of the evaporation and decreases the resolution of the droplet size distribution. When the upper limit level of droplet size is too big, the underestimation of evaporation process may be further compounded and produce bigger errors. Besides, due to the fixed size class, the accuracy of secondary break-up model is decreased and this may get worse when the upper limit level of droplet size is very big. In addition, the total Eulerian approach is not well compatible with the present combustion model yet. Due to these shortcomings, the standard method (Eulerian-Lagrangian approach) is preferable for spray modeling with further improvement. It is obvious the main weakness of the standard method is the lack of flow information or simple handling of internal flow data. The aim of this study is to develop a 3D cavitation-induced primary break-up method based on incorporating utmost flow information acquired by cavitating internal flow simulations. The typical empirical equations used in standard method is replaced with a three-dimensional data distribution in order to retain the original advantage and weaken the advantage of complete flow information advocated by the total Eulerian

induced by cavitation at both cavitating and non-cavitating cases for one-hole diesel injector nozzle. All of this has been attributed to the presence of vapor, which consequently reduces the viscosity. Hadi Taghavifar [7] conducted a numerical simulation to figure out the effect of cavitation on the spray characteristics, in which the RANS [8] method was employed to explain cavitation-turbulence interaction. They found that turbulence induced sprays are more susceptible to be evaporated rapidly. Zhixia He and Wenjun Zhong [9] investigated the effect of nozzle geometrical and dynamic factors on cavitating and turbulent flow inside diesel nozzle. In their study, it is found that the cavitation inception is caused not only by geometrical factors but also by dynamic factors, including injection pressure and injector needle lift. They also insisted that the needle lift has significant effect on the flow characteristics. Indeed, the needle lift influences the flow characteristics through the injection process, not simply the needle lift values for diesel injector. To get a further understanding on the relationship between needle motion and the transient cavitating flow characteristics, how the stability position for flow characteristics changes with different needle movement profiles need to be investigated. It has been reported by many researches that the injection process may influence the cavitation evolution in the internal flow [10–16]. Jose Maria Desantes [10] analyzed the influence of needle lift on the cavitation using Large-Eddy simulation method. The simulations show that the needle lift significantly influences the momentum flux and effective velocity. Similar conclusions have been drawn by Zeidi’s research [11]. These related researches attributed such phenomenon to the consequence of the varied needle lift and conducted the research at fully opened needle position for varied needle lift cases. Many numerical researches assumed the internal flow always keeps stable at varied needle lift position, but the internal flow may never come up to stability throughout the entire injection process for low needle lift cases. From previous visualized experiments [12,13,17], interesting cavitation phenomenon was observed, but these experiments were not close to real applications due to slight structure modifications and lack of highquality photographic equipment. Thus, to get transient nature and the critical stability position of the cavitating flow, it is necessary to take the dynamic needle motions into consideration and give a thorough investigation on the transient cavitation characteristics in an injection process. Coupling the cavitating flow characteristics is the first step for the standard method for spray modeling. Indeed, there are two main 779

Fuel 233 (2018) 778–795

D. Li et al.

approach. First the transient cavitating flow characteristics is investigated with a modified and validated Zwart model. The turbulence effect is included. Then all the flow data at the nozzle exit region is imported with a 3D distribution form into the 3D method. The Ansys Fluent and AVL Fire CFD codes are used for the cavitating flow characteristics and spray formation studies, respectively. 2. Cavitation model

the turbulence effect due to the fact that the turbulence significantly influences cavitating flow characteristics. Thus, the estimation of local turbulent pressure fluctuation can be defined as [27]: Thus, the local pressure P can be modified as: This practice has been found to be much simple, robust and almost as good as the more rigorous practice of Reference [28] and this is strongly recommended by Singhal [29]. Then the phase change rate is modified as:

2.1. Model profile

R=

All the cavitating flow simulations were conducted by the commercial software Ansys Fluent 14.0. To model cavitating flow, the mixture model was utilized in this study. The mixture model solves for the mixture momentum equation and prescribes relative velocities to describe the dispersed phases. The internal flow is of large pressure gradient, it is preferable to use the realizable k-epsilon turbulent model which has shown substantial improvements over the standard k-epsilon turbulent model. The flow features of this model have included strong streamline curvature, vortices and rotation. The near-wall treatment used the non-equilibrium wall functions because the critical dimensions, including initial needle lift position and nozzle hole diameter, were much small. The key elements in the non-equilibrium wall functions are as follow: Launder and Spalding’s log-law for mean velocity is sensitized to pressure-gradient effects; The two-layer-based concept is adopted to compute the budget of turbulence kinetic energy in the wallneighboring cells. It can be seen from the key element that the nonequilibrium wall function is much suitable for high pressure-gradient cases and flows through a small gap. In this study, the nozzle channel is extremely small, which right suits this condition. The gravity effect is neglected, because the high level of turbulence does not allow large bubble growth and the flow is of high speed. To predict this type of cavitating flow, a two-phase cavitation model named Zwart model is utilized. The details of the Zwart model can be found in Reference [26]. The vapor transportation, bubble dynamics consideration and phase change rate are given as follows:

→ ∂ (α v ρv ) + ∇ ·(α v ρv Vv ) = R ∂t

(1)

2 PB−P 3 ρl

(2)

DRB = Dt R=

R=

3Fvap·αnuc·(1−α ) RB 3Fcond ·α ρv RB

ρv

2 Pv−P , if Pv > P 3 ρl

2 −Pv + P , if Pv < P 3 ρl

R=

3Fvap·αnuc·(1−α ) RB

3Fcond ·α ρv RB

ρv

2 Pv−P + 0.195ρk , if Pv > P−0.195ρk ρl 3

2 −Pv + P−0.195ρk , if Pv < P−0.195ρk ρl 3

(7)

(8)

The parameter values in the cavitation model are determined by performing several series of computations for the tunnel and the sharpedged orifice. The assessment criteria include comparisons of mass flow rate and cavitation zone location and extent. Due to the linear relationship between the coefficients and other parameters, we only need to adjust the Fvap and Fcond values. A large number of Fvap and Fcond values were tried for several flow conditions (e.g., for the tunnel, pressure drop = 19, 45, 60, 70, 80 bar). After many permutations and combinations, the most satisfactory values of Fvap and Fcond were found to be 50 and 0.01, respectively. These values were also tested on other problems, including hydrofoils. All these simulations gave satisfactory results with flow pattern and pressure distribution in good agreement. The nucleation formation volume fraction αnuc is fixed at 5 × 10−4. In this model, it is assumed that all the bubbles in a system have the same size, and the total interphase mass transfer rate per unit volume is calculated using the bubble density numbers(n), and the mass change DR rate of a single bubble: R= n× 4π RB 2ρv DtB . The bubble number increases when the vapor volume fraction increases, but the size of bubble never changes. For example, the vapor volume is doubled, all the bubble size keeps at RB , but the bubble number is doubled to conserve the mass balance. Thus RB should be set at a very small value. The bubble radius RB in this study is set at 0.001 mm.

(

)

2.2. Cavitation model validation Several simulations for different ambient pressures were utilized to validate the accuracy of the present computational model. The experiment results used for validation come from research of Winklhofer [30]. The experimental results are often served as validation data for describing basic cavitation flow phenomenon in many numerical researches. The experiment was carried out on a visualized nozzle under stationary flow conditions at inlet pressure of 10 MPa. The flow channel dimensions were selected to provide flow velocities and velocity gradient. The flow pattern is close to those found in fuel injections of diesel

(3)

(4)

In this study, the Zwart model has been slightly modified to include

Fig. 1. Tunnel geometry sketch and computational mesh grid. 780

Fuel 233 (2018) 778–795

D. Li et al.

engine. In the experiment, the cavitation regimes have been performed comprehensively. The sketch of the tunnel is shown in Fig. 1a). The tunnel width is around 0.3 mm and the tunnel inlet rounding radius is of 0.02 mm. Fig. 1b) gives the mesh grid of the tunnel. As shown in Fig. 1b), the mesh is of over 220,000 grids, with over 100,000 grids around the critical middle domain. This approach is proved to have the ability to avoid the grid-dependence of simulations [31]. Further enhance on grid-independence solutions is attached by refining the near-wall grids of the tunnel. It is important that the mesh near the wall is properly sized to ensure accurate simulation of the flow field. The height of the first mesh cell off the wall required to achieve a desired Y+ using flatplate boundary layer theory. The y plus value is around 40. Fig. 2 shows the experimental results and simulations of vapor diesel volume fraction distribution inside the rectangular shape nozzle. The cavitation length is defined as the vapor length with vapor volume fraction over 0.9 along the tunnel axial direction and the cavitation region is defined as the region where the vapor volume fraction is over 10%. As shown in Fig. 2, the simulations of the present computational model have good agreement on cavitation length with the experimental results of all the cavitation regimes, including cavitation inception, developing cavitation and super cavitation. The cavitation onset of the simulation is slightly earlier than the experiment obtained by backlight illumination. This can be explained by the fact that the experimental results are taken from the photographs of a nozzle. The nozzle is captured with backlight, and the intensity of the transmitted light was used to indicated the cavitating region. A high transmitted intensity means that only one phase exists and a weak light means a vapor-liquid mixture exists in this region. Thus, it is impossible to determine the exact vapor volume fraction in the experiment, and only whether the cavitation occurs or not can be revealed. It is weak on experimentally capture the vapor volume fraction with high resolution. The measured cavitation area with a high-vapor volume fraction is much larger than the simulations in the figure. The color label in the figure is only used for simulation. Besides, the real inside surface of the tunnel in the experiment is rough. The roughness of the surface can to some degree influence the cavitation phenomenon inside the tunnel. The large green area represents the area with vapor volume fraction around 0.4–0.6. The numerical scheme used to solve the volume fraction equation is the coupled scheme and the solution method is pressure-velocity coupling. The solution is of second order precision and the currant number is at around 25.

Fig. 3. Comparisons of experiment and simulations for mass flow rate.

The mass flow rate has also been tested for model validation. Fig. 3 illustrates the mass flow rate comparisons between experiment and simulations. The mass flow rates plotted in the Fig. 3 is liquid mass flow rate and the position where the mass flow rates were measured is at the tunnel outlet, as shown in the Fig. 1. The discrepancy is tiny and this verifies that the model gives good prediction on mass flow rate. The velocity profiles in the tunnel entrance (0.053 mm from the orifice inlet), shown in Fig. 4, are also compared with the experiment for the case of 1.5 MPa, 3.3 MPa and 4.5 MPa outlet pressure. The results show that the model has a good prediction on the velocity profile of the experiment. To further validate the turbulence-included Zwart model, the sharpedged orifice experiment by Nurick [32] is employed for comparisons. The orifice profile is depicted in Fig. 5a). Fig. 5b) shows the built mesh for the orifice with the same scale. The mesh is of over 52,000 grids, with dense refinement at the narrow tunnel for over 31,000 grids. Fig. 6 shows the discharge coefficient with varied cavitation number. The ṁ ,where ṁ is the mass discharge coefficient (Cd ) is defined as A 2ρ (pin − pout )

flow rate obtained from the simulation and A refers to the area of nozzle exit cross section. The rest parameters ρ , pin and pout are the fuel density, injection pressure and ambient pressure, respectively. The caP −P vitation number σ is defined as σ = P in− P v , where Pv refers to the in

Fig. 2. Comparisons of experiment and simulations for vapor diesel volume fraction distribution. 781

out

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 6. The discharge coefficient with varied cavitation number.

Fig. 4. Comparisons of experiment and simulations for velocity distribution.

which is also useful for the spray modeling in the next step. According to the characteristics of two-dimensional turbulent boundary layer flow, the spacing of the grid node near the wall is smaller and the distance is gradually increasing that is away from the wall. All the meshes used in this study has been refined in this method. To take the interaction between the internal cavitating flow and the spray jet into consideration, the computational domain was extended to the chamber zone. The computations were performed in the dynamic mesh with the needle moving by the path defined in the varied needle motion profiles. The dynamic mesh was handled with remesh method embedded in the Ansys Fluent. To make the mesh in each step with high quality, the high-quality predesigned meshes were set to replace the working mesh at the corresponding needle position. The predesigned meshes were replaced at needle positions with spacing distance of 0.035 mm. The complete datas distributed on the previous mesh were also transfered to the new mesh. Part of the predesigned meshes are shown in Fig. 8 for explaining purpose. This is the recommended method to keep precision in the CFD codes. All the predesigned meshes have been refined with the similar method mentioned above. The time step set herein is at 0.5 us in the initial 200 us and at 1 us for the rest time up to 700 us.

saturation pressure at the local temperature. In Nurick’s experiment, a phenomenological equation highly consistent with the experiment is derived as follow:

Cd = Cc σ

(9)

where the Cc is the contraction coefficient related to the orifice geometry, usually around 0.62. It can be found that the simulating result is in high agreement with the experiment. The coefficient of determination is over 0.977. From the above comparisons of experiment and simulations, it can be concluded that the computational model can accurately reproduce the behavior of internal flow at cavitating conditions and predict the flow characteristics with high confidence. 3. Cavitation characteristics analysis 3.1. Nozzle geometry and mesh grid The nozzle used in this study is a single-hole diesel nozzle. The physical prototype comes from a high-pressure BOSCH injector used in our previous experimental study [33], with the nozzle diameter of 0.18 mm. The sketch of the nozzle is shown in Fig. 7a). With the employment of dynamic mesh method, the needle is designed to move up and down. The details of the nozzle geometry can be found in Table 1. Fig. 7b) shows the mesh grid for the nozzle. The mesh is of around 200,000 to 350,000 structured grids. The mesh grid number increases with the needle position lifted towards the fully opened position. The gap between the needle seat and the needle is with refinement of over 52,000 grids for calculating precision. The sac and orifice region has also been refined to get more details for the turbulent flow information,

3.2. Grid dependency test To testify the grid-independence of the above mesh refinement method, three meshes with different number of grids were utilized to simulate a steady case for the internal flow. The boundary conditions is set with injection pressure of 90 MPa and ambient pressure of 4 MPa. All the meshes have been refined with the same refinement method mentioned above. The total grid number are 35,612, 94,448 and 377,792, respectively. The case for each mesh is labeled with the grid

Fig. 5. Orifice geometry sketch and computational mesh grid. 782

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 7. Nozzle geometry sketch and computational mesh grid. Table 1 Nozzle geometry. orifices

1

rounding radius (mm) entrance diameter (mm) exit diameter (mm) orifice length (mm)

0.03 0.18 0.18 0.80

number. Fig. 9 shows the initial pressure and mixture distribution of the cases. The chamber environment was initialized with pure nitrogen at 4 MPa pressure. Two valves were set at the sac zone entrance for transient modeling, and it herein kept fully open for the steady case. Fig. 10 illustrates diesel vapor volume fraction of the three cases. It is obvious that the diesel vapor volume fraction distributions are almost the same. This indicates that this mesh refinement method can gain satisfactory grid-independence. The mass flow rates are respectively 7.785, 7.787 and 7.784 mg/ms for the 35,612, 94,448 and 377,792 case, which is in high agreement. This also proves the grid-

Fig. 9. The initial pressure and mixture distribution of the test cases.

Fig. 8. Predesigned meshes at certain needle position. 783

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 10. The diesel vapor volume fraction distribution of the test cases.

indicates that the internal flow may never come to stability for very low needle lift cases. As the internal flow goes into the closing stage, the cavitating flow keeps stable and then goes through great changes. Fig. 12b) gives the details of cavitation development in the closing stage. When the needle approaches the needle seat, the gap between the needle and the seat is narrowed. When the time over 620 micro seconds, the gap is narrowed to a very small value. The flow from the gap to the sac zone is much like a liquid jet. The mass flow rate significantly reduces. The pressure at the sac zone quickly drops and the vapor formed at this zone when the pressure is below the saturation pressure. The cavitation region in the orifice shrinks and the vapor diesel fraction decreases significantly until cavitation occurs in the sac zone. The wake stage is the transient process that the slight-level cavitation that confined in orifice and sac zone shrinks to none. All the characteristics during the opening and closing stage can be further proved by the mass flow rate at the nozzle exit. The normalized mass flow rate is shown in Fig. 13 attached with normalized needle position. The normalized mass flow rate is defined as the real mass flow rate divided by the maximum value of mass flow rate, and the normalized needle position is defined as the real needle position divided by the maximum needle lift value. Thus, the two parameters can be normalized to 0–1 range for convenient viewing. As the main stream of liquid diesel approaches the nozzle exit, the mass flow rate at the nozzle exit increases rapidly. When the needle arrives at certain critical

independence of the mesh refinement method. 3.3. Needle lift profiles and test conditions Fig. 11 shows one of the typical needle lift profiles for the singlehole diesel injector. As shown in the figure, the lift process is divided into three stages: the opening stage, the closing stage and the wake stage. The maximum needle lift value is defined as the maximum distance between the needle and needle seat during the whole needle lift process. All the needle motion profiles in this study share the same average opening and closing speed. The relationship between needle lift and time is almost linear in real engine applications. The slight difference on the needle lift profiles caused by varied injection pressure was neglected in order to exclude the possible effect brought by needle lift speed. The test conditions for numerical simulations are listed in Table 2. The simulations were conducted on varied needle motion profiles with the same opening and closing speed and with varied boundary pressures. The only difference among all the needle motion profiles is the maximum needle lift value. 3.4. Transient cavitating flow evolution A typical needle lift motion profile was used for demonstration of the transient cavitating flow process. The needle motion profile details can be found in Fig. 11, with average opening speed at around 1 m/s and closing speed at around 1.17 m/s. The injection pressure and ambient pressure were set at 60 MPa and 4 MPa, respectively. The maximum lift value was set at 0.35 mm. Fig. 12a–b) illustrates the cavitation phenomenon inside the nozzle when the needle arrives at certain position during one lift process. All the figures were presented in chronological order. In the open stage, the cavitation happens at the corner of the orifice entrance, as shown in Fig. 12a). It is often assumed in many numerical researches that cavitating flow keeps stable and the needle lift is kept at fully opened position throughout the whole injection process. The details during the opening and closing stage were neglected. In many researches it is often assumed that the internal flow is in a relatively steady state during the whole injection process, even for very low lift position cases, but it is away from the reality because the internal flow did not come to stability until the needle moves past certain position. It is found in this simulation that the steady state arrives during the opening stage. This indicates that the cavitating flow comes to steady state as the needle is still moving towards the fully opened position. It also somehow

Fig. 11. Typical needle lift profile. 784

Fuel 233 (2018) 778–795

D. Li et al.

Table 2 Test conditions for numerical simulations. Maximum lift value (mm)

0.021, 0.028, 0.042, 0.140, 0.350

Injection pressure (MPa) Ambient pressure (MPa) Fuel type Fuel temperature (K)

30, 60, 90, 120, 180, 240, 300 0.1, 1, 2, 4, 5, 6 diesel 293.15(@inlet)

position, the mass flow rate stops increasing and keeps stable when the needle continues moving upward. The steady state remains until the needle approaches the needle seat. The two critical positions of starting point and ending point are also shown in the figure, namely SP and EP. To investigate the SP and EP changes can help further understand the transient internal cavitating flow characteristics. The critical stability position refers to SP position in this study when no explanation is attached. The flow coefficients, as shown in Fig. 14, were also analyzed to further investigate the characteristics during the whole injection proA cess. The area contraction coefficient Ca is defined as Al ,where Al refers to the area occupied by liquid diesel in the nozzle exit cross section. The velocity coefficient Cv can be obtained by such formula: Cd = Ca*Cv . Consistent with the mass flow rate trends, the flow coefficients, including Cd,Ca andCv, all come to stability at around the SP position and turn to instability after EP position. The cavitation development has significant influences on the internal flow characteristics. The mass flow rate was also consistent with the cavitation development.

Fig. 13. Normalized mass flow rate at varied needle position.

gives the corresponding mass flow rate for these varied needle motion profiles. For conditions that the needle was lifted over 0.042 mm, they share the same mass flow rate and the only difference is that needle movement profile with short injection time width goes into close stage earlier. This indicates that the position around 0.042 mm is a critical point for nozzle with this type of geometry and motion profiles. To further investigate this critical point, the mass flow rate of 0.021 mm and 0.028 mm cases were also appended. Maximum needle lift values of such cases were designed to be much lower than 0.042 mm in order to make the needle fail reaching the SP position. As shown in 0.021 mm and 0.028 mm cases, the mass flow rate seems to get another critical stability point. It is obvious that comparing with cases with maximum needle lift value over 0.042 mm, the critical stability position for 0.028 mm case moves ahead. The 0.021 mm case almost has not come into stability. This indicates that for low needle lift cases, the critical stability position shifts to an earlier position when the maximum needle lift value decreases. This trend shall last until the stability never exists throughout the injection process. As for higher needle lift cases, the

3.5. Influence of maximum needle lift value To investigate the transient cavitating flow inside the nozzle, varied needle lift motion profiles were utilized. The needle profiles have different maximum lift values. All the profiles share the same opening and closing speed in order to exclude the speed factor. The needle lift profiles are shown in Fig. 15a). The maximum lift values are as follows: 0.021 mm, 0.028 mm, 0.042 mm, 0.140 mm and 0.350 mm. Fig. 15b)

Fig. 12. Cavitation phenomenon of the whole needle lift process. 785

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 14. Discharge coefficient (Cd ), area contraction coefficient (Ca ) and velocity coefficient (Cv ) throughout the needle lift process.

Fig. 16. Critical stability position (SP) and stability time ratio for varied injection and ambient pressures.

critical stability position remains stable and never changes as the maximum needle lift value increases. 3.6. Influence of boundary pressure Fig. 16a) gives the critical stability position (SP) and stability time ratio for varied injection pressures. The stability time ratio is defined as the time width in steady state divided by the time width of the whole needle lift process. The critical stability position SP moves ahead as the injection pressure increases in the 30–90 MPa range. When the injection pressure is over 90 MPa, as shown in the figure, the SP position seems to stabilize down at around 35 mm position. This can be explained by the fact that the internal flow with higher injection pressure shares high velocity and this reduces the time of the flow reaching the exit. And high velocity facilitates the transition from instability to stability. The stability time ratio also increases with the increase of injection pressure and gradually stabilizes down. This indicates that increasing injection pressure can extend the steady flow time width during injection process. Fig. 16b) shows the critical stability position (SP) and stability time ratio for varied ambient pressures at injection pressure of 60 MPa. For practical applications on diesel engine, the ambient pressures were set at 0.1, 1, 2, 4, 5 and 6 MPa, respectively. As shown in the figure, the SP position maintains stable with the increase of ambient pressure. This indicates that the critical stability position is not sensitive to ambient pressure. The stability time ratio also keeps almost constant for varied ambient pressures. It is obvious that the change of ambient pressure hardly influences the steady state of the internal cavitating flow.

Fig. 15. Mass flow rate and corresponding needle lift motion profiles for varied maximum needle lift values.

786

Fuel 233 (2018) 778–795

D. Li et al.

Table 3 Differences of the typical method and 3D method at certain time with time step δt. type

Typical method (AVL Fire)

3D

parcels type initial velocity mass fraction

1 v average

n vi (1ñ)

initial position Nozzle file Injected mass

random or threshold only 3D form M (for δ t)

accumulated mass at fixed positions 2D or 3D form

Ai ·αi·vi ∑in= 1 (Ai ·αi·vi )

t1 + δt

∫ (Ai ·αi·vi) + Mi − loss t1

position (cell), which is obtain by simulation result. The 0 mm position is set at the center of the nozzle exit. Overall, the effective velocity distributions are quite close. With further detailed observation at around −0.09–0.09 mm positions, slight increase can be observed as the ambient pressure decreases. When the ambient pressure decreases to 2 MPa, the discrepancy of effective velocity almost reduced to none. The effective velocity distributions at the nozzle exit for varied injection pressures are illustrated in Fig. 17a). The ambient pressure is set at 4 MPa. Obviously higher injection pressure shares higher effective velocity. The effective velocity distribution shall be used for the spray modeling in the next section.

3.7. Cavitation characteristics discussion From the above analysis, it is easy to find that the transient cavitating flow keeps stable for more than 90% of the total needle lift process. Thus, the flow information at the nozzle exit can be relatively regarded as quasi-steady, which can be employed in the spray modeling. This also facilitates the Lagrangian approach in the next step due to the convenience for initialization. In addition, this can greatly improve the accuracy for Lagrangian approach. It is because that the Lagrangian approach is weak on predicting spray with violent fluctuations. Moreover, the cavitation simulation is necessary to be conducted before spray simulations, because the case with low needle lift can never gain the mass flow rate value of much higher needle lift. The mass flow rate for the low needle lift case is in close relationship with the maximum needle lift value. When the needle lift position is over certain

Fig. 17. Effective velocity distribution for varied injection and ambient pressures.

Similar pattern is also found in the effective velocity distribution at the nozzle exit for varied ambient pressures, as shown in Fig. 17b). The effective velocity is actually the simulated liquid velocity at certain

Fig. 18. Comparisons of typical method and 3D method. 787

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 19. The sketch of the nozzle exit cross section for the single hole nozzle.

importance is further analyzed with comparisons of the total Eulerian approach and the standard method.

4. Spray characteristics analysis 4.1. Previous studies on primary break-up modeling To investigate the spray formation, two typical methods are available among previous studies [34,35], namely, the Eulerian Multi-Fluid approach and the standard method. The Eulerian Multi-Fluid approach is to perform a single calculation covering both the internal nozzle region and the spray domain. This method has been a hot research topic recently, because this method can transfer the complete flow information. This is the main reason to utilize this approach in many numerical research. But this approach is not recommended in this study for these reasons: Firstly, the total Eulerian approach often prefers to use the Direct Numerical Simulation (DNS) method to gain better accuracy. The DNS method is much time consuming due to large amount of calculations. It is not available for many computers with ordinary calculation ability, which largely limits its application. Secondly, the total Eulerian approach may not involve the Sauter Medium Diameter (SMD) of the spray. Only turbulence-induced mixture distribution or turbulence intensity may be not enough. Thirdly, the corresponding cost of gaining the complete flow information is that the simple handling of calculation is allowed. But the simple handling is much inaccurate as the assumptions can bring great errors. For example, the bulk liquid is assumed to consist of gas, liquid and n types of droplets, with each type denotes a droplet diameter level. The droplets diameters of these n types are usually ranging from 5–100 um and the typical upper limit is empirically set at 50 or 100um to cover all the SMD levels in reality. It is too divorced from the reality because in the simplest approach the droplet size classes are fixed and may not change in space and time. However, fixed size classes cause an underestimation of the evaporation and decreases the resolution of the

Fig. 20. The initial condition sketch of the n types of droplets parcels per time step.

Table 4 Initial conditions for i-th type parent drop per time step. type

i

initial diameter velocity magnitude mass fraction

Dini vi

velocity angle total mass at this time step

Ai ·αi·vi ∑in= 1 (Ai ·αi·vi )

arccos

( )) viw vi

t1 + δt

∫ (Ai ·αi·vi) + Mi − loss t1

value (slightly less than 0.1 mm), the mass flow rate shall stabilize around one constant value even with a much higher maximum needle lift value. Luckily the critical needle lift position is often much lower than the typical needle lift values set for real spray and engine applications. But this also reveals the importance for the implementation of cavitating flow simulation in order to exclude the occurrence of cases similar with the low needle lift case, avoiding only using empirical equations for the parcel initialization in the spray modeling. The 788

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 21. The method for flow data changing from 1D to 2D form.

Fig. 22. Facilities for measuring spray development and Sauter mean diameter.

interaction in the spray simulation. As for the collision/coalescence model in the Lagrangian approach, the collision calculation is performed for pairs of particles only if they are in the same computational cell. It is assumed that the particles associated with each parcel are uniformly distributed within the sampling volume (computational cell) in which they are located. In addition, the total Eulerian approach, even for the Semi-Eulerian approach which extends the internal flow domain to a bit part of the spray region, is not well compatible with the present models on combustion engine applications so far. The standard method is to separate the simulation into two steps: First, a two-fluid model is used to get the flow data in the nozzle orifice through Eulerian approach. Then the flow data at the nozzle exit is transferred into the Discrete Droplet Method (DDM) as initial conditions. This method is also known as Discrete Phase Model (DPM). Additionally, a primary break-up model (usually the “blob” injection model) is applied in the DDM model. As for this approach, the two-fluid flow can be considered with or without cavitation characteristics in the internal flow. The present CFD codes, including Fluent or AVL Fire, has already offered the options on consideration of cavitation effect. The cavitation is usually taken into account by simply changing the effective initial diameter. The effective diameter is derived through the 1D method [36]. The details of this 1D method is list below: First, the Nurick’s phenomenological equation Eq. (9) is employed. The effective conditions at the exit are then calculated as:

Fig. 23. comparison of injection rate for measured and simulated result.

droplet size distribution. When the upper limit level of droplet size is too big, the underestimation of evaporation process may be further compounded and produce bigger errors. Besides, due to the fixed size class, the accuracy of secondary break-up model is decreased and this may get worse when the upper limit level of droplet size is very big. Thirdly, the total Eulerian approach may neglect the effect of particle

789

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 24. Mesh for spray modeling. Table 5 Main models used in spray modeling. type

model

primary break-up secondary break-up evaporation coalescence/collision wall interaction

0D, 1D, 3D (blob injection) Wave [38] Dukowicz [39] O’Rourke [40,41] Walljet0 [42]

Pin−Pout

Ugeo = Cd

Ueff =

Ugeo Cc

1 ρ 2 l

−

Aeff = Ageo

Deff =

Pout −Pv ρl ·Ugeo

(10)

(11)

Ugeo Ueff

(12)

4·Aeff π

(13)

The details of the 1D method can be found in the reference [36]. Obviously the 1D method weakens the description of the cavitating flow information. And the key advantage of the total Eulerian approach is the acquisition of the complete flow information, but this approach has also many great physical problems as listed above. Thus, apart from the 1D method, the CFD code has offered a nozzle interface method to utmost utilize the internal flow data. The nozzle file which contains internal flow data shall be completely imported into the spray modeling. The 3D method also uses the similar way to incorporate flow data for primary break-up. The very difference is the 3D method has a very different scheme for introducing droplet parcels. The 3D method allows for nozzle file with 2D form and 3D form distribution. Besides, the 3D method also applies for transient spray modeling with relatively big internal flow changes. The method in CFD code, named as typical method, introduces certain constant number of blobs per time step. The mass fraction shall be equally divided for each orifice cell or one unit

Fig. 25. Comparisons of spray tip penetration for injection pressures of 60, 90 and 120 MPa.

790

Fuel 233 (2018) 778–795

D. Li et al.

measured by the experiment. The parent drops in the 3D method are divided into n types according to the diesel liquid fraction level along the radial direction of the nozzle exit. It is worth noting that the more detailed flow data brought by the cavitating flow simulation, the more types of parent drops can be set in the 3D method. The detailed initial condition of the i (representing 1ñ ) type of parent drop is shown in Table 4. Each of these types of drops represents the flow characteristics of corresponding region, which contains the information of mass fraction, liquid fraction, velocity magnitude distribution and even velocity angle. The total initial number of parent drops (droplets parcel form in CFD codes) should be enough for an even circumferential distribution when the angle profile is attached. Then the information of flow data can be coupled as much as possible, which compensates the loss of flow information by 1D method and avoids the simple and inaccurate handling of calculations by using total Eulerian approach. The injected droplet parcels for each type varies due to the transient internal flow data per time step. If the i-th type is allocated with no sufficient mass, the mass at the fixed position shall accumulate to the next step and no parcel shall occur this time step. The mass per time step is allocated based on mass fraction and the mass is equally divided only when it goes to certain type. Due to the axisymmetric geometry feature of the single hole diesel nozzle, the Eulerian simulation is conducted on a two-dimensional mesh grid. The Eulerian simulation result is corresponding to the onedimensional positions at the nozzle exit. To couple the Eulerian simulation result, the one-dimensional result at the nozzle exit is extended into a two-dimensional form. According to the axisymmetric geometry feature, we have copied the flow characteristics along the circumferential direction, as shown in Fig. 21a). Each color type in Fig. 21a) denotes one type of region with respective velocity and vapor volume fraction characteristics. For each type of region, the initial parcel positions are evenly distributed along the circumferential direction and the initial position number (K) is determined by the mass per time step (Δm ) and the mass fraction for this type of region (βmi ), the calculating equation is as follow:

Table 6 Comparisons of average injection velocity in the initial 20 us. injection pressure

average injection velocity in the initial 20 us (m/s)

60 MPa 90 MPa 120 MPa

0D

1D

3D

271 321 347

275 362 381

325 391 421

blob, and the parent parcel occurs at every cell face or random position at the exit, respectively. But the 3D method introduces different types of parcels per time step to represent the velocity gradient with the initial ligament. The positions of these type of parcels is relatively fixed, but the occurrence of each type depends on the mass fraction for each type and the total mass per step. The main differences of typical method and 3D method are shown in Fig. 18 and Table 3. 4.2. Primary break-up analysis for transient spray modeling In order to overcome the shortcomings of the total Eulerian approach and the standard method (with or without 1D), a compromise scheme is proposed in this study. The 1D method is extended into a 3D form with new modeling scheme to gain as much flow data information as possible that is obtained from the former two-phase flow simulation step. The nozzle exit cross section is the most important region for the spray formation, it is thus seriously considered. Fig. 19 gives the sketch of the nozzle exit cross section. The cross section is divided into n part. Each part labeled as i denotes a region with diesel liquid fraction of αi and the velocity in this region is denoted as vi , and it is vertical to the cross section. All these data come from the cavitation simulation in the above section. Thus, it is easy to obtain the mass flow rate of the diesel liquid in such form: n

ṁ l = ρl ·Aeff ·veff =

∑ (ρl ·Ai ·αi·vi)

(14)

i=1

K= βmi

The n represents n type of region, all of these regions make up the total nozzle exit cross section. The blobs initialized at certain position is obtained by mass component per time step. Suppose the mass component at certain position is mli and the blob diameter is Dini , the blob 1 number and velocity are mli / ρl 6 πDini3 and vi , respectively. In general, the initial blob velocity at the nozzle exit is based on the conservation of mass. Based on this division, the “blob” injection model is incorporated and modified. This is the initial step for the DDM for full cone sprays. A series of large blobs of approximately nozzle diameter, which are close to the nozzle exit, represent the coherent liquid jet. The diameter is subsequently reduced due to the mass detachment calculated from the primary break-up model. In this study, the initial diameter Dini is set at an effective diameter calculated by such equation:

(

Dini =

4·Aeff π

βmi =

ρl ·Ai ·αi·vi n

∑i = 1 (ρl ·Ai ·αi·vi )

(17)

(18)

r = B0 · ∧

4 ∑i = 1 (αi Ai ) π

(16)

Fig. 21b) shows the initial parcel positions along the circumferential direction for one type of the regions. The parent drops that undergo the primary break-up process are also called ligament. When a certain criterion (empirically We< 40 ) is fulfilled, primary break-up is switched off and the ligament is transformed into the assumed regular droplet parcel form, which then goes into secondary break-up [37] process. The liquid break-up is modeled by postulating that new drops are formed (with drop radius, r ) from the n types of parent drops (with radius, a ). The new drop radius is derived with such equation:

)

n

=

Δm 1

ρl · 6 πDini3

where B0 = 0.61. This is a recommended value by Reitz [37]. The change of the radius of a parent drop or “blob” is assumed to follow the rate equation:

(15)

The present application primarily deals with droplets and a spherical shape is assumed as initial step. To couple within the DDM, the blobs release the primary spray droplets with constant initial diameter of the parent drops. The initial condition sketch for these n types of droplets parcels per time step is depicted in Fig. 20. The n types of drops denote the very series of large parent drops. And the u , v and w correspond to the axis in the 3D coordinate system, with w axis perpendicular to the nozzle exit cross section. In the 1D method or CFD nozzle interface (typical method), only one type of drop is initialized with the radial velocity is achievable by specifying the half outer spray angle at 9 degree

da (a−r ) =− (r ⩽ a) dt τ

(19) 3.726·B1·a . ∧·Ω

Ω denotes the where τ is the break-up time, defined as τ = wave growth rate. The B1 is a constant with recommended value of 20 by Reitz. This value is dependent on the injector characteristics and other values. The B1 is set at 26.4 in this study, which is the most satisfactory value after adjustment for the 0D and 1D method. The 3D method is not much sensitive to the B1 value. The basis of this model is the concept that the atomization of the 791

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 26. Spray shape for 60, 90 and 120 MPa injection pressure at the same STP.

792

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 27. Method for extracting droplet diameter data at certain position.

simulated result at 90 MPa injection pressure. From the figure, it can be found that the simulated result is in high similar pattern with the measured data. All the initial data for spray modeling comes from the measured injection rate and transient simulated distribution data for primary break-up. For spray modeling, the three-dimensional structured mesh grid shown in Fig. 24 is utilized. The mesh is of over 30,000 grids and has been refined with hexahedrons about 50 mm around the axis normal to the nozzle exit cross section. The detail of these refinement is also presented in the figure. The mesh is built as a cylinder with bottom diameter of 80 mm and height of 150 mm, which is enough for the spray development. All the above detailed handling can prevent error brought by the mesh grid. For comparison, the standard method without internal flow considered, 1D method and 3D method has been all applied to a series of injection pressure of 60, 90 and 120 MPa. The standard method without internal flow considered is labeled as 0D method for shorthand. For all these three methods, the secondary break-up model is set for the standard wave model [38] and droplet interaction model is set for O’Rourke [40,41]. More details of employed models for spray modeling can be found in Table 5. The time step set herein is at 5 us at the initial 0.5 ms and at 10 us for the rest time up to 2 ms. The time step of initializing parent drops is set at 50 us. The spray tip penetration (STP) result of each method is compared with the experiment for the 60, 90 and 120 MPa injection pressure and 4 MPa ambient pressure. The spray pictures were repeatedly captured

injected liquid and the subsequent break-up of drops are indistinguishable processes within a dense spray. All these n type parent drops shall be initialized with different mass flow data and the process shall be treated as wave growing process in the next step. This is the main difference between the typical standard method and the 3D method. 4.3. Experimental and numerical comparisons To validate the 3D method, the experiment done previously has been employed for comparisons. The experiment is conducted at a constant chamber with constant temperature of 20 degrees Celsius and constant ambient pressure of 4 MPa. The injection pulse is measured at around 2.0 ms and the nozzle orifice diameter is 0.18 mm. The spray development is captured by a high-speed camera (shown in Fig. 22a) and the Sauter mean diameter (SMD) is calculated by PDIA technique with LaVision system (shown in Fig. 22b). The details of the experiment work can be found in Reference [33]. The injection rate has all been measured and imported into the CFD codes for 60, 90 and 120 MPa injection pressure. The measurement facility is a small constant chamber with a Kistler pressure sensor set inside the wall. The chamber is filled with diesel at 4 MPa pressure in order to recreate the same pressure environment. The injection rate is calculated by the measured pressure data based on small level low-pass filtering algorithm and related equations. Fig. 23 gives the typical one of the measured and 793

Fuel 233 (2018) 778–795

D. Li et al.

Fig. 28. The extraction position and the SMD result.

ability on the spray penetration. The average injection velocity predicted by the 0D, 1D and 3D are shown in Table 6. The average injection velocity is calculated by calculating the difference of spray penetration in the initial 20 us. It is obvious that the 3D predicts the highest velocity and 0D predicts the lowest. This also indicates that the cavitation can somehow increase the injection velocity at the nozzle exit. Fig. 26 gives the spray shape for 60, 90 and 120 MPa injection pressure. To make the comparisons reasonable, the figures are extracted at the same STP length for convenient viewing. Each figure in Fig. 26 is made by combining half of the experiment picture and half of the simulating result. The left part refers to the simulation and the right part is the experiment. Each part is cut at the nozzle center axis. Two positions of each injection pressure are employed for spray shape comparisons. The simulation result has been processed with a color inversion and the white dots that represent droplets in the figure are only used to present the spray shape. From the figure, it is obvious that the spray shape of the 3D method is the most similar one among the three methods. The 1D and 0D method seem to overestimate the spray area and form a much fatter spray than the real spray. This is due to the lack of consideration of diesel liquid distribution and velocity distribution for the 1D and 0D method. From the statistic point of view, it is inevitable that with the repeat of the same type parent drop initialized at the nozzle exit, the chance of droplets emerging at certain position shall increase and this in turn forms the local spray shape in this region. This is the common situation for the 1D and 0D method. Besides, the underestimation of STP also can bring errors on the spray shape prediction, because the late arrival of STP at certain position may reduce the velocity to some extent and can bring other further effects caused by velocity reduction. In general, the 3D method overcomes the shortcomings by attaching the distribution information of velocity and diesel liquid fraction. In this way, the mass fraction distribution is also introduced. The 3D method gives good predictive ability on spray shape and this can be much valuable for high efficient combustion and low soot emission engine applications. The droplets density at the center of the spray is very high, which weakens the transmission of the background illuminating light. Therefore, the microscopic measurements were conducted at 40 mm

Table 7 Comparisons of Sauter mean diameter result for spray modeling and experiment. Type

experiment

0D

1D

3D

SMD (um)

28.06

16.37

15.06

24.54

for over three times. The spray tip penetration (STP) was calculated for the picture captured each time and found that the STP at certain time for repeated captured pictures are almost the same, with error length less than 1 mm and error rate less than 0.5%. The injection process of the experiment has a good consistency. The STP of the 3D method is in high degree agreement with the experiment, as shown in Fig. 25. The B1 constant for the 0D and 1D method has been tuned to make sure the present result is the most approaching result for the experimental result, i.e., the B1 constant for each method is set at the most satisfactory value. In the 3D method, all the initial settings are predefined. The 3D method is not sensitive to the B1 constant from our repeated test, but it is much sensitive to the time step set in the predefined profiles for droplet parcels initialization at the nozzle exit. The smaller the time step is, the more-smooth the spray tip penetration curve is. As shown in the figure, the spray tip penetration for 3D method is not much smooth due to the time step set for droplet parcels initialization. But the general accuracy is not affected. When the injection pressure is over 90 MPa, the 1D method shows obvious advantage on spray tip penetration predicting. This is because when the injection pressure is over 90 MPa, super cavitation happens inside the nozzle. The super cavitation denotes that the cavitation develops to the nozzle exit. The flow characteristics changes greatly when the super cavitation happens, which also greatly increases the initial parcel velocity. When the cavitation is considered, the liquid velocity near the axis increases significantly. This also explains why the 0D method underestimates the spray tip penetration and 1D method improves on spray tip penetration predicting over 0D method. But the 0D and 1D method all underestimate the spray development to some extent. The error rate, defined as the standard deviation, is also given in the figure. The 3D method shares the determination coefficient with value over 0.994 for all the injection pressures. This indicates that the 3D method has satisfactory forecasting 794

Fuel 233 (2018) 778–795

D. Li et al.

below the injector nozzle tip and 8 mm from the spray axis at injection pressure of 60, 90 and 120 MPa and ambient pressure of 4 MPa. Fig. 22b) shows the typical image illustrating the local droplets. The data at the corresponding captured time and position in the simulation is extracted. The captured time corresponds to the time of 1.2 ms in the simulation. The SMD result can be directly extracted by picking the result at certain position with output box in AVL Fire, as shown in the Fig. 27. The extracted result corresponds to SMD result of the cells around certain position is shown in Fig. 28. The extraction process is as follow: First, the x, y and z planes are used to find the position of the measurement point; Then, the pick value function is used to extract the SMD result of the demanded position. Table 7 shows the SMD result for the three methods at injection pressure of 90 MPa and ambient pressure of 4 MPa. The microscopic characteristics of the spray is only conducted on injection pressure of 90 MPa. The 3D method is slightly weak on predicting the SMD, but it has more satisfactory result than the 1D and 0D methods. This can be explained by the fact that the initialized velocity of the parent drops in the 3D method is not fixed at a constant and the parent drop with lower initialized velocity can contribute to slightly bigger droplet.

[10] Desantes Jose Maria, Salvador Francisco Javier, et al. Large-eddy simulation analysis of the influence of the needle lift on the cavitation in diesel injector nozzles. J. Automobile Eng. 2015;229(4):407–23. [11] Zeidi SeyedMohammadJavad, Dulikravich George S, Reddy Sohail R. Effects of needle lift and fuel type on cavitation formation and heat transfer inside diesel fuel injector nozzle. Napoli, Italy: ICHMT International Symposium on Advances in Computational Heat Transfer; 2017. [12] He Zhixia, Zhong Wenjun, Wang Qian, et al. An investigation of transient nature of the cavitating flow in injector nozzles. Appl Therm Eng 2013;54:56–64. [13] He Zhixia, Zhang Zhengyang, Guo Genmiao, et al. Visual experiment of transient cavitating flow characteristics in the real-size diesel injector nozzle. Int Commun Heat Mass Transfer 2016;78:13–20. [14] Javier Lopez J, Salvador FJ, de la Garza Oscar A, Arregle Jean. A comprehensive study on the effect of cavitation on injection velocity in diesel nozzles. Energy Convers Manage 2012;64:415–23. [15] Lockett RD, Loverani L, Thaker D, et al. The characterisation of diesel nozzle flow using high speed imaging of elastic light scattering. Fuel 2013;106:605–16. [16] Sun Zuoyu, Li Guoxiu, Yusong Yu, et al. Numerical investigation on transient flow and cavitation characteristic within nozzle during the oil drainage process for a high-pressure common-rail DI diesel engine. Energy Convers Manage 2015;98:507–17. [17] Jiang Guangjun, Zhang Yusheng, Wen Hua, Xiao Gan. Study of the generated density of cavitation inside diesel nozzle using different fuels and nozzles. Energy Convers Manage 2015;103:208–17. [18] Oley F, Trummler T, Hickel S, et al. Large-eddy simulation of cavitating flow and primary break-up. Phys Fluids 2015;27:086101. [19] Berg EV, Edelbauer W, Tatschl R, et al. Validation of a CFD Model for Coupled Simulation of Nozzle Flow, Primary Fuel Jet Break-Up and Spray Formation. Salzburg, Austria: Spring Technical conferences of the ASME Internal Combustion Engine Division; 2003. [20] Lebas R, Menard T, Beau PA, et al. Numerical simulation of primary break-up and atomization: DNS and modeling study. Int J Multiph Flow 2009;35:247–60. [21] Edelbauer W. Numerical simulation of cavitating injector flow and liquid spray break-up by combination of Eulerian-Eulerian and Volume-of-Fluid methods. Comput Fluids 2017;144:19–33. [22] Berg EV, Edelbauer Wilfried, et al. Coupled calculation of cavitating nozzle flow, primary Diesel fuel break-up and spray formation with an Eulerian multi-fluidmodel. J. Eng. Gas Turbines Power 2005;4:127. [23] Pandal A, Payri R, García-Oliver JM, Pastor JM. Optimization of spray break-up CFD simulations by combining Σ -Y Eulerian Atomization Model with a Response Surface Methodology under diesel engine-like conditions (ECN spray A). Comput Fluids 2017;156:9–20. [24] Poursadegh F, Lacey JS, Brear MJ, Gordon RL. On the fuel spray transition to dense fluid mixing at reciprocating engine conditions. Energy Fuels 2017;6:31. [25] Zigan L, Schmitz I, et al. Effect of fuel properties on spray breakup and evaporation studied for a multihole direct injection spark ignition injector. Energy Fuels 2010;24(8):4341–50. [26] Zwart PJ, Gerber AG, Belamri T. A two-phase model for predicting cavitation dynamics. Proceedings of the International Conference on Multiphase Flow. Yokohama, Japan: ICMF-2004; 2004. [27] Hinze JO, Turbulence. 2nd Edition, McGraw Hill, New York. [28] Vaidya Singhal AK., et al. Multi-Dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change. ASME FED Meeting, Paper No. FEDSM’97-3272, Vancouve, Canada, 1997. [29] Singhal Ashok K, Athavale Mahesh M, et al. Mathematical basis and validation of the full cavitation model. J Fluids Eng 2002;124:617–24. [30] Winklhofer E, Kull E, Kelz E, et al. Comprehensive hydraulic and flow field documentation in model throttle experiments under cavitation conditions. In: Proceedings of the ILASS-Europe Conference, Zurich, Schwitzerland; 2001. [31] Zhang XB, Qiu LM, Qi H, et al. Modeling liquid hydrogen cavitating flow with the full cavitation model. Int J Hydrogen Energy 2008;33:7197–206. [32] Nurick WH. Orifice cavitation and its effect on spray mixing. J Fluids Eng 1976;December:681–7. [33] Li Donghua, Gao Yuxin, Liu Shenghua, et al. Effect of polyoxymethylene dimethyl ethers addition on spray and atomization characteristics using a common rail diesel injection system. Fuel 2016;186:235–47. [34] Dukowicz JK. A Particle-Fluid Numerical Model for Liquid Sprays. J. Computational Physics 1980;35:229–53. [35] Durst F, Milojevic D, Schönung B. Eulerian and Lagrangian predictions of particulate two-phase flows: a numerical study. Appl Math Modelling 1984;4:101–15. [36] Kunsberg V, Sarre C, Kong SC, Reitz RD, Modeling the effects of injector nozzle geometry on diesel sprays. SAE, 1999-01-0912. [37] Reitz RD, Diwakar R. Structure of high-pressure spray. SAE Technical Paper Series 1987:870598. [38] Liu, AB, Reitz RD, Modeling the effects of drop drag and break-up on fuel sprays. SAE, 930072. [39] Dukowicz, JK, Quasi-steady droplet change in the presence of convection. Informal report Los Alamos Scientific Laboratory, LA7997-MS. [40] O’Rourke PJ. Statistical properties and numerical implementation of a model for droplet dispersion in turbulent gas. J. Comput Phys 1989:83. [41] O’Rourke PJ, Bracco FV. Modeling of drop interactions in thick sprays and a comparison with experiments. IMECHE 1980. [42] Naber, JD, Reitz RD, Modeling engine spray/wall impingement. SAE, 880107.

5. Conclusions Combined with previous experimental studies, a cavitation-induced 3D primary break-up model is developed. The varied needle motion profiles are used to extend the scope of this study and test the consistency for the internal flow characteristics. For high needle lift cases, the cavitating flow remains stable for more than 90 percent of the time. For low lift cases, the mass flow rate at steady state increases with the increase of needle lift value. All of these testify the implementation necessity of cavitating flow simulations and somehow validate the feasibility of the 3D method. With the cavitating flow data incorporated in three-dimensional distribution with several types of parcels per time step, the 3D method has better predictive ability on the experimental spray penetration and spray shape than the typical standard method. The 3D method is slightly weak on predicting droplet diameter, but it has better performance comparing with the 1D and 0D methods. In general, the 3D method has good performance for spray modeling and constitutes a promising platform for future development. Acknowledgement This work is financially supported by the Key Industry Innovation Chain Project of Shanxi Province and the National Natural Science Foundation of China (No. 51576159). References [1] Hiroyasu H, Arai M, Shimizu M. Break-Up Length of a Liquid Jet and Internal Flow in a Nozzle. Gaithersburg, MD: ICLASS-91; 1991. [2] Arai M, Shimizu M, Hiroyasu H. Similarity Between the Breakup Lengths of a High Speed Liquid Jet in Atmospheric and Pressurized Conditions. Gaithersburg, MD.: ICLASS-91; 1991. [3] Arai Masataka. Break-up Mechanisms of a High Speed Liquid Jet and Control Methods for a Spray Behavior. Hksohima: International Symposium on Advanced Spray Combustion; 1994. [4] Bode, Jiirgen. Max-Planck-Inst fur Stromungsforschung. March, 1991. [5] Chen Zhifang, Yao Anren, Yao Chunde, et al. Effect of fuel temperature on the methanol spray and nozzle internal flow. Appl Therm Eng 2017;114:673–84. [6] Payri F, Bermudez V, Payri R, Salvador FJ. The influence of cavitation on the internal flow and the spray characteristics in diesel injection nozzles. Fuel 2004;83:419–31. [7] Taghavifar Hadi, Khalilarya Shahram, Jafarmadar Samad, Taghavifar Hamid. A RANS simulation toward the effect of turbulence and cavitation on spray propagation and combustion characteristics. Theor Comput Fluid Dyn 2016;30:349–62. [8] Nordin N Complex. Chemistry Modeling of Diesel Spray Combustion. Chalmers University of Technology; 2001. PhD Thesis. [9] He Zhixia, Zhong Wenjun, Wang Qian, Jiang Zhaochen, Shao Zhuang. Effect of nozzle geometrical and dynamic factors on cavitating and turbulent flow in a diesel multi-hole injector nozzle. Int J Therm Sci 2013;70:132–43.

795