Stochastic sub-grid scale model for LES of secondary atomization – assessment and evaluation for ECN Spray-A conditions

At high Reynolds and Weber numbers typical of diesel-like sprays, the finest turbulent structures are highly intermittent in nature, resulting in intense fluctuations in gas phase velocities and can significantly contribute to the atomization process. In order to account for their influence on breakup of spray droplets, we introduce a new stochastic breakup model to be used in conjunction with large eddy simulation (LES) for the gaseous flow. The model is based on a stochastic parent-to-child relaxation of droplets, whose parameters are linked to the viscous dissipation rate on residual scales “seen” by a droplet along its trajectory. In order to introduce the intermittency effects on the droplet breakup, this dissipation rate is simulated stochastically, in the framework of log-normal process. The non-reacting “Spray-A” experiment from Engine Combustion Network (ECN) is used to assess the performance of the new stochastic breakup model in comparison to the standard hybrid KH-RT breakup model in terms of evolution of liquid penetration length and droplet size statistics. The results clearly show that in comparison to the hybrid KH-RT model the stochastic breakup model gives a more accurate prediction of different parameters with relatively less sensitivity to the grid resolution.


Introduction
In order to ensure efficient fuel-air mixing on short time scales in diesel engines, the liquid fuel jet is atomized into very small droplets by injecting at relatively high velocities of the order of 300-500 m/s. The numerical modelling of spray atomization is largely based on Lagrangian description of the spray used in conjunction with the Eulerian modelling of the turbulent gas phase [1][2]. The spray formation is modelled by injecting computational droplets or parcels typically of the size of the nozzle diameter, which then undergo a series of breakup events by a presumed breakup mechanism. The most commonly used hybrid KH-RT model [3] is based on breakup from fastest growing instability on the droplet surface which then results in formation of a single-scale of child droplets. Since in Diesel-like conditions, the high speed injection generates strong turbulence, with highest frequencies of order of characteristic times in KH-RT instabilities, the conditions for droplets breakup become stochastic, and thereby the atomization can be viewed as a stochastic random process. To this end, in [4] the child droplet size is generated from the numerical solution of the master equation for the size distribution function. Having as initial condition, the delta pick around the size of parent drop, the evolution of this master equation has been constructed in a way to assure the Rosin-Rammler type distribution on long times. Also, introducing the inertial response of droplet to turbulent velocity fluctuations, a new formulation for the critical radius was proposed in [4], as an alternative to the Kolmogorov's critical size. In [5][6], the approach is different. Namely, it was shown that at large times under scaling symmetry at the constant fragmentation frequency, the fragmentation equation reduces exactly to the Fokker-Plank equation with two model parameters i.e., the first and second logarithmic moments of the fragmentation intensity spectrum. Based on these principles [7][8] formulated a stochastic breakup model for LES modelling of sprays with the closures for the logarithmic moments expressed in terms of local flow field properties. In order to introduce the intermittency effects of gas phase turbulence on droplet breakup [9] proposed a new stochastic breakup rate expression for relaxation of parent droplet size towards a maximum stable child droplet size. In this model the breakup frequency and the relaxation radius of child droplets are defined as a function of the local viscous dissipation rate "seen" by the droplet along its trajectory. Within the framework of RANS modelling of the gas phase turbulence, the fluctuations in the dissipation rate were randomly sampled from a log-normal distribution with model parameters expressed in terms of local flow Reynolds number. However, the random sampling procedure does not account for the spatiotemporal correlations of dissipation rate field along droplet trajectory. Therefore, in this paper, in the context of LES modelling of the gas phase, the stochastic log-normal process [10] is used to model the evolution of dissipation rate field along each droplet trajectory.
In the next section a brief description of the governing equations for the gas phase and droplet motion are presented. This is followed by an elaborate presentation of the formulations for both the hybrid KH-RT model and the new stochastic breakup model is given. Finally, a detailed assessment of the performance of new breakup model in comparison to the nonevaporating ECN Spray-A experiment and the KH-RT model is presented.

Governing Equations
In this study the carrier gas phase is computed using large eddy simulation and solving the individual droplet dynamics using the Lagrangian approach. The interaction between the two phases due to mass and momentum transfer are modelled by the source terms in the governing equations of the gas phase. In LES, the instantaneous velocity ( ) of ambient gas phase is decomposed into filtered (�) and SGS components ( ′ ).
= � + ′ (1) where, � is solved by the continuity and momentum equation obtained by filtering the Navier Stokes equations of the instantaneous velocity field.
In Eq. (2-3), is the viscous stress tensor, is the source term accounting for the interaction between the gas and liquid phases. On the other hand, Γ is the sub-grid shear stress tensor which is modelled using sub-grid scale (SGS) eddy viscosity ( ) : where is the mean strain rate. In this study, the one-equation eddy viscosity model [11] is used, where SGS viscosity is calculated from the subgrid kinetic energy ( ) and the filter width (Δ) i.e., = Δ � . Further is solved using a transport equation accounting for its production, dissipation and convection as shown in Eq. (5).
where ϵ = ϵ 3 2 Δ � , is the rate of dissipation at the resolved scales. On the other hand, the dynamics of the individual Lagrangian droplets is modelled using the equation of motion accounting for the drag force acting on it due to the relative velocity between the droplet and local gaseous flow. The filtered droplet drag force equation is given by Eq. (6), where τ is the droplet relaxation time scale expressed further in terms of Eq. (7). and is the drag coefficient which is a function of the droplet Reynolds number ( ). The source term given by Eq. (8) accounts for the contribution of the droplet drag forces on the gas-phase momentum.
and , are the mass, acceleration and the number of actual droplets in a given parcel , is volume of the computational cell. Since the objective of this study is to assess the relative influence of breakup models on spray structure, the effects of droplet dispersion and collision are accounted in the computations.

Spray breakup models Hybrid KH-RT breakup model
In [3] the breakup rate and the size of resulting child droplets after each breakup event are assumed to be proportional to the frequency (Ω) and wavelength (Λ) of the fastest growing instabilities on droplet surface. In the near-nozzle region, the growth of surface instabilities is characterized by Kelvin-Helmholtz (KH) mechanism. A parent droplet with radius r upon breakup after a characteristic time scale results in stripping off child droplets from of size , which are then added as new parcels to the computation. The correlation between the model parameters i.e., , and the parameters of fastest growing KH-instability i.e., Ω , Λ is given by Eq. (10).
The model parameters 0 and 1 are used to vary the breakup rate to fit the experimental data. Further breakup of droplets in the downstream region is modelled using the Rayleigh-Taylor (RT) instability mechanism. In this case, the breakup time is found to be equal to the inverse of the frequency of the fastest growing wave i.e., = Ω −1 . At a time = , the parent droplet completely breaks down into small droplets whose size is proportional to the wavelength this disturbance i.e., = Λ . In either case the analytical solutions [12][13] for the maximum growth rate (Ω / ) and the corresponding wavelength (Λ / ) are expressed in terms of parameters, as gaseous Weber number, Ohnesorge number, ratio of the densities of liquid/gas phases and the droplet acceleration.

Stochastic breakup rate model
In [9], the droplet breakup rate is modelled as a random stochastic process by assuming that both the model parameters i.e., the breakup frequency * and the relaxation radius * in the breakup rate expression given by Eq. (11) to be random functions of a stochastic variable represented by the local instantaneous viscous dissipation, .
= − * * (11) As outlined in [4] the expression for the relaxation radius * , is obtained by accounting for the droplet inertia in force balance between the turbulent shear forces acting on the droplet and the capillary forces: where ν is the viscosity of the gas phase, and is the surface tension coefficient. An analytical expression for the breakup frequency is obtained by dimensional analysis of the three important parameters influencing droplet breakup i.e., droplet density , surface tension and dissipation rate .
In [9] the dissipation rate "seen" by the droplet is randomly sampled from a log-normal distribution once over a breakup time. This random sampling approach does not account for the temporal correlations of the instantaneous dissipation rate field ( ) along the particle trajectory. So instead of the random sampling of , in this paper we propose to use the lognormal stochastic process of [10] as given by Eq. (14-15). In difference to [10], the model parameters expressed in terms of the local flow field parameters, and the local filtered dissipation rate ε � is used instead of its mean value 〈ε〉 : where ( ), is the increment of a standard Brownian process, 〈 〉=0, 〈 2 〉=dt. Here Δ is the filter width proportional to computational grid size, η is the Kolmogorov length scale.

Experiments and numerics
The standardized spray experiments [15][16][17] from the Engine Combustion network (ECN) group are used in this study to evaluate the spray breakup models. The n-dodecane fuel jet is issued from an injector with a nozzle diameter of 90 µm at an injection pressure of 150MPa into a constant volume spray chamber filled with an ambient gas at a density of 22.8 kg/m 3 and a temperature of 440K. This experimental configuration is usually referred to as "ECN Spray-A". A cylindrical domain of length 100mm and a diameter 50mm is used to computationally represent the spray chamber. The base coarse mesh referred to "C-Grid" consists of a uniform cell size of 250µm both in axial and radial directions. Another finer mesh referred to as "F-Grid" with a cell size linearly varying from 125µm to 250µm both in axial and radial directions is used to study the effects of mesh resolution. For solving the Eulerian gaseous flow, a weakly-compressible flow solver with an implicit pressure treatment based on the PISO-algorithm is used. While for discretization of spatial gradients second order numerical schemes are used, an implicit first-order Euler scheme is used for the time integration. For the Lagrangian modeling of spray droplets, we used the classical "blob" approach of [2], wherein the initial size of all computational droplets or parcels is assumed to be the same as the nozzle diameter. Initial blob velocity is calculated from the mass flow rate profile obtained from the experimental data. The orientation of the initial velocity is defined randomly within a user-specified "spray cone angle" which is assumed to be 12 degrees. The number of parcels injected is determined by the computational time step (Δt) and the total number of parcels to be injected per second (PPS) which is pre-defined manually. For both the spray breakup models, the number of actual droplets per parcel after breakup is calculated based on principle of mass conservation. The time step is controlled by the user defined where Δx is the grid size. In this study a fixed PPS of 2 × 10 7 and a of 0.3 are used. All the calculations are performed using the opensource software OpenFOAM version 6.0.

Results and Discussion
The liquid spray tip "penetration length" is the most commonly used global parameter to characterize temporal evolution of spray. It is defined as the distance where the accumulated liquid droplet mass reaches 95% of the total liquid mass injected at any given instance of time. First the influence of breakup model parameters and the grid size on the overall spray evolution predicted by different breakup models is presented. In case of hybrid KH-RT model the key parameter controlling the breakup rate and the overall spray characteristics is the coefficient 1 and a wide range of values between 5 to 100 have been reported [18,19]. So, two different values of 10 and 40 are used to characterize the influence of this parameter on spray evolution in case of KH-RT model. From the comparison of the penetration length evolution shown on left hand side of Figure 1, it can be seen that the model constant has no influence on the spray evolution. Moreover, it can also be seen that the hybrid KH-RT model tends to largely over-predict the liquid penetration length even on finer grid resolution of Δ=125µm. Next the temporal evolution of the spray penetration length predicted by the stochastic breakup rate model on two grid resolutions is shown on the right-hand side of Figure 1. The results show that even on coarser grid resolution, the new breakup model gives a much better prediction of spray evolution compared to the other two models. Moreover, the spray evolution predicted by the new stochastic breakup model is relatively less sensitive to the grid resolution. The evolution of jet penetration length is directly correlated to the size distribution statistics of parcels resulting from spray breakup through the droplet relaxation time scale τ given by Eq. (7). So, next a quantitative evaluation of the droplet size statistics predicted by different models in comparison to the experimental data is presented. Figure 2 shows the spatial evolution of sauter mean diameter (SMD), which is defined as the ratio of total volume to total surface area of all the droplets located at any given axial position. While the rate of decrease of SMD along the spray centerline quantifies the breakup rate of droplets, the quasi-steady values of SMD farther downstream reflect the resulting droplet sizes after spray breakup. In case of hybrid KH-RT (on left hand side of Figure 2) it can be seen that while the grid resolution has little effect on evolution of SMD, decreasing the modeling constant 1 from 40 to 10 results in faster breakup and smaller droplet sizes. However even with 1 = 10, the SMD of droplets farther downstream of injection are still much larger (~4µm) compared to the experimental values (~1µm). In case of Stochastic breakup rate model, it can be seen that while the breakup rate is faster than hybrid KH-RT model (on right hand side of Figure 2), the resulting droplets sizes are much smaller than the KH-RT model. However, the resulting droplet sizes from the new model (~1-2 µm) are still slightly larger than the experimental values of (~1µm). Usually, the sauter mean diameter statistics are weighted averaged by the number of droplets per parcel and hence do not always give a full picture of correlation between spray breakup and its temporal evolution. In order to identify the underlying differences in spray evolution, the instantaneous spray structure predicted by the two breakup models at the time t=1.5ms are shown next in Figure 3 and 4. Here the droplets are scaled and coloured by their relative sizes i.e., dark regions indicate presence of large droplets. As shown in Figure 3 in case of KH-RT model, we notice a very slow cascade of droplet sizes moving away from the nozzle exit. This results in the jet penetration being characterized by relatively mono-dispersed small parcels present in the leading edge of the spray. On the other hand, as shown in Figure  4, the stochastic model gives a more realistic statistical representation of the spray with wide spectrum of parcel sizes randomly distributed all over the spray domain and therefore result in completely different spray dynamics. However, the presence of these large parcels is filtered out in the droplet statistics, as they have very few droplets in them.

Conclusions
A new stochastic sub-grid scale mode for spray atomization within the framework of LES for the gas phase has been described and assessed in comparison with the KH-RT spray breakup model and the non-reacting ECN Spray-A experiment. The model involves a stochastic breakup rate expression with parameters expressed in terms of viscous dissipation rate field, whose evolution along the droplet trajectory is modelled using the log-normal process with parameters correlated to local flow Reynolds number. Thus, the model aims to introduce the residual scale turbulence effects on spray breakup. Two important parameters quantifying fuel-air mixing i.e., temporal evolution of the liquid spray penetration length and the spatial evolution of the droplet size are used as a measure to evaluate these breakup models. It was shown that the stochastic breakup rate model, being linked to viscous dissipation through the stochastic process on residual scales, and not through the resolved velocity gradients, has demonstrated the ability to give fairly good representation of the spray parameters even on relatively coarse grid and the results are less sensitive to the grid resolution compared to other two models.