Development of a Eulerian Multi-Fluid Solver for Dense Spray Applications in OpenFOAM

: The new generation of internal combustion engines is facing various research challenges which often include modern fuels and different operating modes. A robust modeling framework is essential for predicting the dynamic behavior of such complex phenomena. In this article, the implementation, veriﬁcation, and validation of a Eulerian multi-ﬂuid model for spray applications within the OpenFOAM toolbox are presented. Due to its open-source nature and broad-spectrum of available libraries and solvers, OpenFOAM is an ideal platform for academic research. The proposed work utilizes advanced interfacial momentum transfer models to capture the behavior of deforming droplets at a high phase fraction. Furthermore, the WAVE breakup model is employed for the transfer of mass from larger to smaller droplet classes. The work gives detailed instructions regarding the numerical implementation, with a dedicated section dealing with the implementation of the breakup model within the Eulerian multi-ﬂuid formulation. During the veriﬁcation analysis, the model proved to give stable and consistent results in terms of the selected number of droplet classes and the selected spatial and temporal resolution. In the validation section, the capability of the developed model to predict the dynamic behavior of non-evaporating sprays is presented. It was conﬁrmed that the developed framework could be used as a stable foundation for future fuel spray modeling.


Introduction
To increase the efficiency of internal combustion (IC) engines, which is tightly coupled with the increase of the compression ratio, modern engines are being designed to operate in compression ignition mode [1]. Due to the high combustion temperatures (resulting from non-premixed combustion mode), conventional diesel engines suffer from high nitrogen oxides' emissions. New strategies should lower the pollutant emissions while keeping the efficiencies as high as possible. One of the possible approaches is the partially premixed combustion [2], where the fuel spray is injected directly into the cylinder, but the timing and duration of the start of injection are varied to reach optimal combustion efficiencies at all working conditions. Optimization of fuel-air mixing significantly depends on the nozzle design, which controls the spray penetration length and droplet sizes. The characteristics of the in-nozzle flow affect the spray by causing velocity fluctuations, which enhance the mechanical breakup of the liquid jet, and consequently the formation and collapsing of cavitating bubbles [3,4]. Furthermore, the increase of injection pressure (up to 3000 bars) promotes effective breakup and atomization of liquid fuels [5,6]. Therefore, understanding these complex physical phenomena of spray dynamics at high pressures is crucial for improving the efficiency of IC technology.

Formulation of the Mathematical Model
This section presents a Eulerian multi-fluid model specialized for spray applications, which is generalized for an optional number of incompressible fluids. The first fluid is the continuous gas phase, and the remaining phases describe the liquid fuel phase. The liquid phase is sub-divided into an arbitrary number of fluids (using the classes method). The model equations are conditionally averaged using the procedures described by [23,[35][36][37][38]. The proposed model is an upgrade of the work given in [39], which was developed for monodisperse bubbly flows. The model is enhanced with advance interfacial momentum transfer models specialized for deforming droplet flows at high phase fractions, which can reproduce thick spray effects in near nozzle regions. Furthermore, the proposed model now includes breakup functionality for high Weber number flows (We > 100) using the WAVE breakup model [32][33][34]. Within this work, the proposed model is tested for non-evaporating spray conditions. Therefore, the presented model does not contain an evaporation model.
Linking the population balance equation (PBE) [40] with the standard continuity and momentum equations enables the model to predict polydisperse flows. Moreover, the multi-fluid formulation (in comparison with the standard two-fluid formulation) allows the model to capture velocity and spatial variance because the interfacial momentum transfer models are strongly dependent on the droplet size. In this work, the PBE is discretized using the classes method, which means that droplets are divided into a finite number of droplet classes. This approach is similar to the Multiple Size Group (MUSIG) [41] or Inhomogeneous MUSIG [42,43] model, but it offers a higher resolution and precision (each droplet class has its momentum and phase continuity equation, i.e., there are no velocity groups).
The droplet diameters are discretized using the equal diameter distribution, i.e., the i-th droplet diameter d i is calculated from: where d max and d min are the maximum and minimal droplet diameter and n droplets is the total number of droplet classes. The solver employs the RANS approach using the single-phase k -turbulence model for the continuous gas phase [44]. The dispersed phase turbulence is evaluated using the turbulence response coefficient.

Phase-Intensive Momentum Equation
The conditionally averaged phase-intensive momentum equation for phase ϕ is given by the following expression: where U ϕ gives the averaged phase velocity, α ϕ is the phase fraction, R eff ϕ gives the joined viscous and turbulent stress, p denotes the mixture pressure, ρ ϕ is the phase density, g is the gravitational acceleration, M ϕ is the averaged interfacial momentum transfer term, and S Mϕ is the net momentum source term due to breakup, which is caused by the transfer of mass between the droplet classes.
The interfacial momentum transfer term M ϕ reads: where n phases is the total number of fluids, and M ϕ,i = −M i,ϕ gives the momentum transferred between phases ϕ and i. The presented work limits the model to droplet flows, where the gas is described by only one continuous phase, and the fuel phase is divided into an arbitrary number of classes depending on the droplet diameter. Therefore, n phases = n droplets + 1, where n droplets is the total number of droplet classes. The momentum between the droplet phases and the continuous gas phase is exchanged via the turbulent dispersion force and drag (other forces such as virtual mass and lift can be neglected): where the subscript d indicates the dispersed phase, and the subscript c signifies the continuous phase. The relative velocity term is calculated as U r,i = U c − U d,i . C d,i and C td,i are the drag and turbulent dispersion force coefficient of the i-th droplet phase. The k c term represents the turbulence kinetic energy of the continuous phase. The diameter of the i-th dispersed phase is given with d i .
In Equation (5), the turbulent dispersion force is implemented following the approach presented by Reeks [45] and Bertodano [46]. The coefficient C td,i can be treated as a constant value, but it can also be linked to the time scales associated with droplets, using the following expression: where τ c,i is the time constant of the particle, and τ d,i is the effective time constant of the fluctuating force acting on the particle. τ d,i is calculated as: and τ c,i is given by: In Equation (7), ν c indicates the kinematic viscosity of the gas (continuous) phase, and Re d,i is the Reynolds number for the i-th droplet class given by: and, in Equation (8), c is the dissipation of turbulence energy of the continuous phase. In Equation (5), the drag coefficient of i-th droplet class is implemented following the procedure described by Liu et al. [47]. Due to the large deformations of droplets in engine-like conditions (which also lead to droplet breakup), Liu et al. [47] suggested blending the drag coefficient between an ideal sphere and a disc (which is approx. 3.6 times greater): where C d,sphere,i is the drag coefficient of ideally spherical particle with the diameter d i , and y i is the normalized distortion parameter (of the i-th droplet class) calculated with the Taylor-Analogy (TAB) model [48]. The TAB model assumes that the droplet distortion can be described as a one-dimensional mass spring system, where the droplet viscosity ν d,i is the damping force and the surface tension σ is the restoring force, which leads to the following expression (when defined using the droplet diameter): Integration of Equation (11) gives the time-dependent normalized distortion equation, which is used for the evaluation of the additional drag term.
The drag coefficient of an ideally spherical particle C d,sphere,i in Equation (5) can be calculated with the following relations: which include the influence of the local phase fraction on the droplet drag presented by O'Rourke and Bracco [49]. The momentum transfer term for the gas phase is calculated as:

Phase Continuity Equation
For incompressible flows, the phase continuity equation (for phase ϕ) can be written in the following form: where S ϕ denotes the net source term due to breakup mass transfer between droplet classes. In this work, the phase continuity equation is implemented following the formulation given by Weller [36], which contributes to the conservativeness and boundedness of the solution. The generalisation for the multi-fluid formulation is described in [50,51]. Consequently, the modified phase continuity equation for polydisperse flows can be written in the following form: where U denotes the mixture velocity, which is defined as: The net source term S i is evaluated using the WAVE breakup model, which is presented in Section 2.3.

WAVE Breakup Model
The aerodynamic interaction between the high-speed droplets and the gas phase introduces the development and growth of disturbances on the droplet surface. The generated deformations of the droplets are practically the dominant cause of droplet breakup, especially in regions further away from the injector nozzle. Reitz and co-workers made a great effort in deriving [32][33][34] a continuous and unified breakup model, often referred to as the WAVE or the Kelvin-Helmholtz model, which was used for modeling of high-speed diesel jets [52,53].
The derived model assumes that a cylindrical liquid jet penetrates a stationary incompressible gas through a round opening. The surface of the liquid jet is subject to initial perturbations which are further increased by the liquid-gas interaction. It is also assumed that only the fastest growing disturbances (denoted with the growth rate Ω, which matches the wavelength Λ) will cause the breakup. Furthermore, Reitz [53] simplified the problem by fitting the numerical results to analytical expressions which give the maximum growth rate Λ i (for i-th droplet class): and its wavelength Ω i : In Equations (17) and (18), Z i gives the Ohnesorge number defined as: T i is the Taylor number: We d,i is the liquid Weber number: We c,i is the gas Weber number: and Re d,i (in Equation (19)) defines the liquid phase Reynolds number (defined using the droplet radius, and not the diameter as in previous models): The size of droplets (stable radius r s,i ) which are formed by the breakup process is usually linearly coupled to the most unstable surface disturbance, i.e., to the wavelength Λ i : where the proportionality coefficient B 0 is of order unity, and, in this work, the standard value of 0.61 is employed. Due to the breakup process and the generation of new smaller droplets, the parent droplets lose mass, i.e., the radius of parent droplets is reduced with the following expression: where the breakup time τ i is calculated as: In Equation (26), B 1 denotes a constant which describes the effects of the inner nozzle flow on the breakup time because those effects cannot be resolved directly with the model [52].
When using the Eulerian multi-fluid approach, the reduction of the parent droplet diameter needs to be converted into a phase sink term (in the parent phase continuity equation) and a corresponding source term (in the child phase continuity equation). Therefore, the net source term S i in Equation (15) is divided in the following manner: where B B,d,i is the droplet birth rate due to breakup from larger droplets (into class i), and D B,d,i is the droplet death rate due to breakup (from class i) into smaller droplets.
Following the procedure described in [54], the rate of change of parent class radius, given in Equation (25), can be reformulated in mass loss per unit volume of phase i, i.e., it can be converted into D B,d,i : More details about the numerical implementation of the model and details regarding the calculation of the droplet birth rate B B,d,i will be discussed in Section 3.1.

Numerical Model
This section gives an overview of the numerical procedures utilized for the implementation and solving of the previously described mathematical model. The collocated Finite Volume Method (FVM) is used for the solution of the previously given equations [55,56]. The proposed solution procedure uses the PISO algorithm [57] and the implemented procedure per each time step is given in Algorithm 1.

Algorithm 1 Solution algorithm per each time step.
Calculate the source/sink terms due to breakup.
Construct and solve the phase continuity equations.
Calculate the interfacial momentum transfer terms.
Construct the phase momentum equations and predict fluxes.
Construct and solve the mixture pressure equation.
Correct fluxes and reconstruct phase velocities.
Construct and solve the turbulence model equations.

Implementation of the WAVE Breakup Model and Phase Continuity Equations
This section gives details about the implementation of the phase continuity equations and the WAVE breakup model. As previously described, the presented model is limited to droplet flows, where the gas phase is represented with only one continuous phase, and the droplets of various sizes are described with an arbitrary number of droplet phases n droplets . The continuous gas phase does not undergo breakup, and since the model does not account for evaporation, the net source term in the phase continuity equation is equal to zero, i.e., S c = 0. The phase continuity equation for the continuous phase is implemented as: The dispersed phase continuity equations are implemented in the following form: In this work, the droplet diameters are discretized using the equal diameter distribution using Equations (1) and (2). Consequently, the temporal change in the parent droplet radius is implemented as: where the temporal change of the droplet radius and the droplet death rate are greater than zero only if the stable radius is smaller than the lower bound of the i-th droplet class. Therefore, the smallest droplet class does not undergo breakup. The droplet death rate D B,d,i defined by Equation (28) can introduce negative solutions of the droplet phase continuity equation, especially when larger time steps are enforced. This work suggests a limiter, which keeps the solution bounded. The proposed limiter compares the local droplet death rate predicted by the model to the maximal allowed value: where ∆t gives the time step value. However, the limiter requires implicit treatment (in terms of α d,i ) of the advection terms in the droplet phase continuity equation Equation (30). The corresponding droplet birth rate of phase j (from phase i) is implemented as: where the mass is transferred from i-th to j-th class only if the stable radius of phase i is within the bounds of the droplet class j. Considering that the mass transfer due to breakup always goes from larger to smaller droplets, the total droplet birth rate of phase j is given by: It is required that the implementation satisfies the conservation criterion, i.e., the total source needs to be zero when summed over all droplet classes:

Numerical Procedure
In this work, all presented calculations given in Section 4 used identical linear solver and discretization settings. Any differences in the case set-up are explicitly mentioned.
The turbulence model equations and the phase continuity equations were solved with a Bi-Conjugate Gradient Method preconditioned by DILU [58], and the pressure equation used the selection algebraic multigrid algorithm [59] with the Gauss-Seidel smoother [60]. All equations used the same absolute tolerance for the normalized residual value of 10 −10 . For a matrix system Ax = b, the normalized residual r is evaluated as [61]: and the normalization factor n is calculated as: where x is the current solution vector, and x denotes the average value of x.
The phase fractions variables were advected using the linear upwind-biased approximation with a limiter for stronger bounding. The upwind scheme is used for the advection of the turbulence model variables. The momentum variables employed the Gamma scheme [62], which is a member of the normalized variable diagram family. All (first) time derivative terms were evaluated using the Crank-Nicholson scheme. Gradients, Laplacians, and cell-to-face interpolations were assessed using linear interpolation.

Results
In this section, a detailed verification study is presented, where the spatial and temporal resolution were systematically varied. Furthermore, the implemented model was tested with different numbers of droplet classes to examine the sensitivity of the model to the droplet class resolution. The last sub-section deals with the validation, where the results are compared with the available experimental measurements such as spray penetration, spray angle, and droplet size distribution.
The testing of the presented and implemented model is done for non-evaporating conditions, where the liquid fuel is injected into a pressurized (2.1 MPa) constant volume vessel filled with carbon dioxide. The diesel fuel is injected through a Mini-Sac nozzle with a diameter of 140 µm and bore length of 0.8 mm. The experimental measurements are available in [63][64][65], where the data were used for testing of various numerical approaches for predicting spray behavior. The physical properties of the gas and liquid phase, which were employed in the following simulations are given in Table 1. The given results utilize the blob injection model [53], i.e., through the duration of injection, large blobs (droplets which are the same size as the nozzle hole) are being added at the inlet boundary, and the inlet velocity is calculated from the corresponding fuel flow-rate. Immediately after the blobs enter the computational domain, the WAVE breakup model shears off smaller child droplets from the surface of the blobs.
The selected fuel injection flow-rate is shown in Figure 1, which was obtained by fitting the curve to the available experimental measurements available at [65].

Verification
The performed verification analysis follows the guidelines for unsteady flows given by [66]. The analysis was carried out by systematically varying spatial and temporal resolution, i.e., five structured grids with uniformly varied refinement levels and four different time step sizes were employed in the study. The grid density was increased towards the nozzle, both in the radial and axial direction. The initial coarse grid was constructed to have two cells per nozzle diameter, and for finer grids, the cell density was uniformly increased. The selected three-dimensional cylindrical computational domain is shown in Figure 2. The outer dimensions of the domain were quite large (the cylinder is 50 mm in radius and 80 mm in length) in comparison with the nozzle diameter, to minimize the influence of the boundary conditions on the solution.
To reduce the computational load of the verification study, and, due to the too high Courant number, when using larger time steps for finer grids, only the smallest time step is used for all meshes. A visual representation of the employed computational grids (with a detailed view of the refinement area near the nozzle) and the corresponding number of cells are given in Figure 3 (sub-figures 3a-e).
The uncertainty and the achieved accuracy in space and time were estimated using the freely available ReFRESCO application [67]. The study is conducted for the spray tip penetration length after 0.5 ms (after the start of the fuel injection) with 14 droplet classes. The input values for the ReFRESCO application are given in Table 2. Table 2. Test matrix for the verification analysis. The values denote the spray tip penetration length (in mm) after 0.5 ms for various spatial and temporal resolutions.  The results of the estimator are given in Table 3, where φ 0 is the extrapolated exact solution, φ 1 is the finest level solution, U φ is the uncertainty estimate, p is the achieved accuracy in space, and q is the achieved temporal accuracy. Table 3. Results of the uncertainty estimation.

Number of Cells
Spray tip penetration 22.4 22.5 0.9% 2.00 2.00 The achieved second order accuracy (both in space and time) was expected, considering the employed numerical methods described in Section 3.

Sensitivity to the Selected Number of Droplet Classes
The sensitivity of the implemented model to the employed number of droplet classes is tested for the previously described flow conditions using the second finest computational grid with 21,465 cells. The sensitivity of the model is tested for the droplet size distribution and the spray tip penetration.
The droplets size distributions were calculated by integrating the fluxes of the individual droplet phases through time, i.e., counting the number of droplets, passing through the predefined circular sampling surface (1 mm in radius and located 62 mm downstream in the axial direction from the nozzle exit). The comparison of the droplet size distributions for 7, 14, and 28 classes is given in Figure 4. In Figure 4, the left sub-plots give the predefined droplet population at the inlet boundary, i.e., they present the employed blob population. The right sub-plots show the droplet population obtained with the previously described sampling surface. Figure 4 shows that the model behavior is consistent in terms of the selected number of droplet classes. However, as expected, the increase in the selected number of classes improves the resolution of the solution. The increased resolution predicts the distribution peak around 7.5 µm, and not in the smallest droplet class, which is not visible from the lower resolution results. In Section 4.3, the presented droplet size distribution (obtained with 28 classes) is compared to available experimental measurements.

Validation
The validation of the model is performed using the results obtained with the second finest grid (21,465 cells) and using 28 droplet classes. The numerical results are compared to available experimental measurements [63][64][65]. The validation is performed in terms of spray angle, spray tip penetration, and droplet size distribution. Figure 6 shows a comparison of the experimentally measured spray penetration curve (denoted by the dashed line) and the one obtained with the previously presented numerical model (indicated by the solid line). The numerical results capture the spray behavior quite well, but there is a significant lag in the penetration between 0.1 and 0.25 ms. The slowdown is a consequence of drag overprediction in the near-nozzle region. In the remaining time interval, the two curves are practically parallel, which suggest that the spray dynamics is captured adequately. The presented spray penetration curves indicate that the blob injection model should be replaced by a more advanced modeling approach, e.g., primary atomization modeling or initialization using the high-fidelity atomization simulations. In future work, the presented Eulerian multi-fluid model is planned to be initialized using the DNS results. The comparison of the droplet size distributions is given in Figure 7. The left sub-plot presents the employed blob population in the numerical simulation, and the right one compares the droplet size distribution obtained with the experimental measurements (dashed line) and the numerical model (denoted by solid bars). The numerical simulation overpredicts the generation of smaller droplets in comparison with the measurements. The smallest droplets are generated immediately after the blobs enter the computational domain. In the near-nozzle region, the relative velocity between the blob droplet class and the gas phase is quite large, which results in tiny values of the stable radius. The numerical model correctly predicts the range of occurring droplets, but the distribution peak is shifted towards smaller droplets, due to the previously described issue. Therefore, the accuracy of the predicted droplet size distributions would benefit from a more advanced modeling approach in the nozzle exit region. The spray angle is calculated as the droplet spreading angle at 70% of the spray penetration tip at the end of the fuel injection, i.e., after 2 ms [65]. The comparison of the experimentally measured spray angle and the numerical prediction is shown in Table 4. The implemented model successfully predicts the cone shape of the spray for the given flow conditions.
Overall, the developed numerical model is capable of describing the dynamic spray behavior, but the employed initialization procedure coupled with the WAVE breakup model is not ideal for capturing all near-nozzle phenomena.

Conclusions
The objective of this work is to provide a new modeling approach to the publicly available simulation framework for fuel spray applications, which is vital for research and advancement in the internal combustion technology. This work presents a detailed description of the developed and implemented simulation framework for predicting the dynamic behavior of dense sprays using the Eulerian multi-fluid model. Special attention was given to the numerical implementation of the breakup model using the Eulerian approach. The presented work employs the WAVE breakup model and the blob injection approach. The implemented model was thoroughly tested to determine the achieved accuracy in space and time. Additional tests presented the sensitivity of the model to the chosen number of droplet classes. Furthermore, the numerical results were also compared to the available experimental data. The tests showed that the developed solver is giving stable and consistent results in terms of the selected number of droplet classes and employed computational grids. The validation section showed that the implemented model is capable of predicting the dynamic behavior of non-evaporating sprays in terms of the spray shape, penetration, and droplet size distribution, but with some limitations. The blob injection approach, coupled with the presented model, introduces some issues in the near-nozzle region. The numerical results overpredicted drag and generation of small droplets (due to breakup of blobs) in the vicinity of the nozzle exit. The calculation procedure is not taking into account the complex nozzle flow (e.g., cavitation and local turbulence) when injecting the blobs into the computational domain. Consequently, the droplet breakup is directly influenced only by aerodynamic forces, which undoubtedly reduces the accuracy of the solution. In future work, the Euler-Euler simulations will be initialized using the DNS results of the spray atomization [19] to test the impact on the solution accuracy. Furthermore, the model is planned to be expanded with evaporation functionality. The implementations presented within this work were done within foam-extend (a community-driven fork of OpenFOAM).