© S. Kaundinya Oruganti et al., published by IFP Energies nouvelles, 2019

1 Introduction

With increasing restrictions on emissions, there has been a continuous evolution towards high pressure fuel injection systems, resulting in highly turbulent in-cylinder flows.

When the spray combustion takes place in highly turbulent conditions and when the chemical reactions are fast, it was shown in [1, 2] that the dispersion of vaporizing droplets is a strong factor controlling the combustion rate. When modeling this dispersion, the effect of turbulent structures of very different length scales on the droplet dynamics have to be considered. Direct Numerical Simulation (DNS), where the entire spectrum of turbulent scales are resolved, have been used to study the spray-turbulence interactions in simple flow configurations. Elghobashi et al. [3–5] have studied the mechanisms of two-way coupling in particle laden isotropic turbulence. They have characterized the modification of the turbulence by inertial particles depending on the magnitude of particle response time (τ _p) in comparison to Kolmogorov time scale (τ _η). Squires and Eaton [6] and Rouson and Eaton [7] have studied the preferential concentration of different particles i.e., accumulation of dense particles because of their inertia in regions of low vorticity and high strain. Toschi et al. [8] and Cencini et al. [9] have characterized the effects of intermittency at small scales on the acceleration statistics for inertial particles. They showed the presence of highly non-Gaussian acceleration Probability Density Functions (PDF) with broad tails indicating high acceleration events. But DNS can be computationally very expensive and not feasible for real engine applications. Over the last decade, there has been a significant increase in application of Large Eddy Simulations (LES) in Internal Combustion Engines (ICE) because of its ability to directly resolve the large flow structures using spatially filtered equations. Usually, the small unresolved scales, often referred to as Sub-Grid Scales (SGS) are modeled using SGS turbulence models [10–13]. As stated by Rutland [14], it is a common practice in the context of Lagrangian modeling in LES to either use laminar correlations or extend sub-models from RANS for modeling the interaction of SGS with the droplets. These approaches usually require very high grid resolutions. Hence Rutland [14], has emphasized on the need to develop SGS models specific to the context of LES modeling of liquid sprays.

To this end, Bharadwaj et al. [15, 16], have developed a model to account for spray induced turbulence, by reconstructing the SGS gas velocity at droplet positions using by deconvolution of the resolved velocity field. And the effect of unresolved scales on droplet dispersion is modeled decomposing the instantaneous gas velocity into the resolved velocity and the velocity at the SGS. The later is calculated from the subgrid kinetic energy (k _sgs) and assumed to be piece-wise continuous in time [17]. Tsang et al. [18] further extended the model by decomposing the SGS velocity into deterministic and stochastic components, where the former is evaluated from deconvolution [15] and the later from the dispersion model [17]. The performance of the SGS model under diesel spray conditions has been evaluated. It was concluded that the dispersion only influences the spatial distribution of liquid droplets and does not have any effect on statistics of the carrier gas phase. Pozorski and Apte [19] have modeled the SGS velocity seen by the droplets using a Langevin stochastic equation to study preferential concentration in particle laden flows. This model has been successful in modeling the dynamics of only large-inertia particles but did not work well with small inertial particles which are concentrated in the small scales usually unresolved by LES. Bini and Jones [20, 21] have modeled the particle velocity increments using the Langevin equation with the variable intensity of the Wiener process. The model has been shown to capture the broad tails in the velocity increment distribution.

All the above proposed SGS models are based on reconstructing the fluctuating flow velocity along the particle trajectories from the subgrid kinetic energy (k _sgs). On the other hand, it was shown in the book of Kuznetzov and Sabel’nikov [22] that the statistics of the droplet velocity relative to the flow velocity is also controlled significantly by the mean viscous dissipation. The latter characterizes the small-scale dynamics in turbulence. This was also confirmed in DNS by Zamansky et al. [23, 24]. Consequently, in the droplet acceleration models, the instantaneous dissipation rate seen by the droplet along its trajectory is to be considered as an important parameter. However in LES this parameter is usually under-resolved, thereby the droplet acceleration is under-resolved as well. This motivated Gorokhovski and Zamansky [25, 26] to propose a new approach, in which the inertial particle acceleration on SGS is simulated by a stochastic process for viscous dissipation rate seen by the droplet along its trajectory. The parameters of those processes are defined from the locally resolved velocity field. Two SGS models are proposed, one for a particle above the Kolmogorov scale (referred to as “stochastic drag” model), and another one for the particle below the Kolmogorov scale as described in [25, 26]. This approach is assessed by DNS [27] for heavy particles immersed in the homogeneously sheared turbulence, and is referred to as LES-STRIP (Stochastic Response of Inertial Particle). Another parameter of this model is the stochastic direction of the droplet acceleration. The model for orientation vector of droplet acceleration is based on random walk over a unit sphere. While in [28], the stochastic equation for the orientation vector on the unit diffusion sphere is based on Ornstein-Uhlenbeck (OU) process, this approach is retained in our current study.

In the previous works of Gorokhovski et al. [25–28] particle laden flows were studied with one-way coupling i.e., dynamics of particles do not influence the turbulence. But in direct injection engines, the momentum transfer from spray to the ambient gas flow contributes significantly to turbulence production. This motivated us to assess the LES-STRIP approach in diesel-like conditions by applying two-way coupling between the spray and the surrounding gas flow. This paper is organized as follows: In Section 2, the computational framework i.e., the governing equations of both Eulerian-Lagrangian phases and the standard dispersion model are described. In Section 3, the LES-STRIP model is recalled in the modified formulation. In Section 4, the experimental test cases and related numerical modeling issues are listed. In Section 5, validation of the dispersion model using the temporal evolution of spray structure is presented. Additionally, an analysis of the statistics of droplet and gas flow fields is provided to illustrate the capabilities of the new model to account for the two-way coupling of spray-turbulence interaction.

2 Governing equations

In this study the two-phase problem is addressed by computing the carrier gas phase using LES and solving the individual droplet dynamics using the Lagrangian approach. The interaction between the two phases due to mass and momentum transfer are modeled by the source terms in the governing equations of the gas phase.

2.1 Gas phase equations – LES

In LES, the instantaneous velocity (u) of ambient gas phase is decomposed into filtered ($\overline{u}$) and SGS components (u′):

$$u=\overline{u}+\mathrm{u\text{'}},$$(1)where, $\overline{u}$ is solved by the continuity and momentum equation obtained by filtering the Navier Stokes equations of the instantaneous velocity field:

$$\frac{\mathrm{ }\mathrm{\partial }{\rho }_f}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\rho }_f{\overline{u}}_i}{\mathrm{\partial }{x}_j}=0,$$(2)

$$\frac{\mathrm{\partial }\left({\rho }_f{\overline{u}}_i\right)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left[{\rho }_f{\overline{u}}_i{\overline{u}}_i\right]}{\mathrm{\partial }{x}_j}=-\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }{\rho }_f{\tau }_{\mathrm{ij}}}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\mathrm{\Gamma }}_{\mathrm{ij}}}{\mathrm{\partial }{x}_j}-S_{\mathrm{liq}}.$$(3)

In equation (3), τ _ij is the viscous stress tensor, S _liq is the source term accounting for the interaction between the gas and liquid phases. On the other hand Γ_ij, is the sub-grid shear stress tensor which is modeled using sub-grid scale eddy viscosity (ν _sgs):

$${\mathrm{\Gamma }}_{\mathrm{ij}}=2{\rho }_f{\nu }_{\mathrm{sgs}}{S}_{\mathrm{ij}}-\frac{2}{3}{\rho }_{f}k{\delta }_{\mathrm{ij}},$$(4)where S _ij is the mean strain rate. In this study, the one-equation eddy viscosity model [9] is used, where SGS viscosity (ν _sgs) is calculated from the sub-grid kinetic energy (k _sgs):

$${\nu }_{\mathrm{sgs}}=C_k\Delta \sqrt{k_{\mathrm{sgs}}},$$(5)where Δ is the filter width and k _sgs is the sub-grid kinetic energy. The k _sgs is solved using a transport equation accounting for its production, dissipation and convection:

$$\frac{\mathrm{\partial }\left({\rho }_f{k}_{\mathrm{sgs}}\right)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left({\rho }_f{\bar{u}}_j{k}_{\mathrm{sgs}}\right)}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left({\rho }_f{\nu }_{\mathrm{sgs}}\frac{\mathrm{\partial }{k}_{\mathrm{sgs}}}{\mathrm{\partial }{x}_j}\right)+{\rho }_f{\mathrm{\Gamma }}_{\mathrm{ij}}{S}_{\mathrm{ij}}+{\epsilon }_{\mathrm{sgs}},$$(6)where ${\epsilon }_{\mathrm{sgs}}=C_{\epsilon }\frac{k^{\frac{3}{2}}}{\Delta }$ is the rate of dissipation at the resolved scales.

2.2 Liquid phase equations – Lagrangian approach

Usually the liquid phase is modeled as discrete parcels, where each parcel represents a ensemble of droplets with same properties. The dynamics of these parcels are obtained by individually solving the momentum equation of each parcel:

$$a_p=\frac{\mathrm{\partial }{u}_p}{\mathrm{\partial }t}=\frac{{\rho }_f}{{\rho }_p{d}_p}\frac{3}{4}{C}_d\left|u-u_p\right|\left(u-u_p\right),$$(7)

$$a_p=\frac{\mathrm{\partial }{u}_p}{\mathrm{\partial }t}=\frac{u-u_p}{{\tau }_p},$$(8)

$${\tau }_p=\frac{4}{3}\frac{1}{C_d|u-u_p|}\frac{{\rho }_p{d}_p}{{\rho }_f},$$(9)where a _p is the droplet acceleration, u _p is the droplet velocity, τ _p is the droplet response time to fluid solicitations and C _d is the drag coefficient. The drag coefficient C _d is a function of droplet Reynolds number ($R{e}_p=\rho |u-u_p|d_p/\nu$). For (Re _p < 1000), it is assumed that the boundary layer around the droplet experiences a transition from laminar to turbulent and the drag coefficient is strongly dependent on Re:

$$C_d=\frac{24}{R{e}_p}\left(1+0.15R{e}_p^{0.687}\right).$$(10)

For much higher velocities (Re _p > 1000), the boundary layer is highly chaotic with vortices of many different scales being shed from the droplet. In this regime the drag coefficient is roughly constant i.e., C _d = 0.42.

2.3 Source terms – two-way coupling

In the case of high pressure diesel sprays, the relative velocities of the two phases is significantly high (usually in the order of 100–400 m/s). So the impact of the momentum transfer from the liquid droplets on the dynamics of the carrier gas phase has to be considered. To this end, a source term S _liq is added to the momentum equation of the gas. In case of non-vaporizing sprays, it accounts only for the contribution of the drag force (F _p = m _pa _p) between the two phases. For a given computational cell, the net momentum transfer is the volume average of the contributions of all the droplets in the cell:

$$S_{\mathrm{liq}}=\frac{1}{V_{\mathrm{cell}}}{\sum }_{m=1}^{\mathrm{Np}}\mathrm{}{n}_{\mathrm{dm}}{F}_{\mathrm{dm}},$$(11)where n _dm is number of droplets in a given parcel m, V _cell is volume of the computational cell and F _dm is the drag force acting on a single droplet in a given parcel m.

2.4 Standard SGS model for dispersion

Most commonly used SGS dispersion model for diesel spray studies [15, 16] is extended from RANS [17], where each component of SGS gas velocity is randomly sampled from a Gaussian distribution:

$$G\left({\mathrm{u\text{'}}}_p\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left[-\frac{{\mathrm{u\text{'}}}_p^2}{2{\sigma }^2}\right],$$(12)where $\sigma =\frac{2}{3}{k}_{\mathrm{sgs}}$, k _sgs is interpreted as the subgrid kinetic energy. And the fluctuating SGS velocity u′_p is sampled once over a turbulent time scale, t _turb. The sampling time t _turb usually represents the time taken by the droplet to traverse through an eddy.

$$t_{\mathrm{turb}}=\min \left[\frac{k_{\mathrm{sgs}}}{{\epsilon }_{\mathrm{sgs}}},0.162\frac{k_{\mathrm{sgs}}^{1.5}}{{\epsilon }_{\mathrm{sgs}}}\frac{1}{|\overline{u}-u_p|}\right].$$(13)

3 LES-STRIP approach for droplet acceleration

The main proposal in LES-STRIP approach [25, 26] is as follows. The acceleration of a finite size particle with size larger than the Kolmogorov length scale but smaller than the filter width is given by:

$$\begin{array}{cc}{a}_p=\frac{u-u_p}{{\tau }_{p,t}}& \mathrm{for} ∆>d_p>\eta \end{array},$$(14)where τ _p,ti.e., the droplet response time is a fluctuating variable due to the turbulence “seen” at the length scale of its diameter.

$${\tau }_{p,t}=\frac{{\rho }_p}{{\rho }_f}\frac{d_p^2}{{\nu }_{p,t}}.$$(15)

Using the Kolmogorov scaling for turbulent viscosity, $\left({\nu }_{p,t}\sim \mathrm{u\text{'}}{d}_p\sim {\epsilon }^{\frac{1}{3}}{d}_p^{\frac{4}{3}}\right)$, the droplet acceleration in equation (14) is characterized by a “fluctuating drag”. The statistics of the latter is conditioned by random viscous dissipation (ϵ) seen by a droplet along its trajectory:

$$a_p=\frac{{\rho }_p}{{\rho }_f}\frac{{\epsilon }^{\frac{1}{3}}}{d_p^{\frac{2}{3}}}\left(u-u_p\right).$$(16)

Introducing the unit vector of the droplet acceleration direction,

$${\vec{e}}_p=\frac{u-u_p}{|u-u_p|},$$(17)and from the Kolmogorov scaling, $\left|u-u_p\right|\sim (\epsilon d_p{)}^{\frac{1}{3}}$, the droplet acceleration in equation (16) can be expressed by the product of two independent random values i.e., the dissipation rate (ϵ) and the acceleration direction vector (${\vec{e}}_p$). Each of these random variables is represented by a stochastic process in the LES-STRIP approach:

$$a_p=\frac{{\rho }_p}{{\rho }_f}\frac{{\epsilon }^{\frac{2}{3}}}{d_p^{\frac{1}{3}}}{\vec{e}}_p.$$(18)

It is worthwhile to note that the experimental studies on the droplet acceleration statistics [29, 30] have shown that the variance of acceleration scales with $\left(d^{\frac{1}{3}}\right)$, as it is the case in equation (18).

For small droplets with size below the Kolmogorov length scale, the LES-STRIP approach [25, 26] is based on the decomposition of the instantaneous acceleration of droplet into two parts i.e., response to resolved flow and the contribution on SGS:

$$a_p=\frac{\overline{u}-u_p}{{\tau }_p}+a_p^{\mathrm{sgs}}\hspace{1em}\mathrm{for} d_p<\eta ,$$(19)where the residual part is given by:

$$a_p^{\mathrm{sgs}}=\sqrt{\frac{\epsilon }{{\tau }_p}.}$$(20)

Equations (18)–(20) require stochastic equations for unit vector of the droplet acceleration direction (${\vec{e}}_p$) and for powers of viscous dissipation i.e., ${\epsilon }^{\frac{1}{2}}$ and ${\epsilon }^{\frac{2}{3}}$. The Ornstein-Uhlenbeck process for ${\vec{e}}_p$ and the method of its integration is taken from [28]. While the stochastic equations for ${\epsilon }^{\frac{1}{2}}$ and ${\epsilon }^{\frac{2}{3}}$ are derived by Ito transformation of stochastic log-normal process [31] as described in [32].

$$\mathrm{d}\epsilon =-\epsilon \left(\ln \frac{\epsilon }{{\epsilon }_{\mathrm{sgs}}}-\frac{{\sigma }_{\chi }^2}{2}\right)T_{\chi }^{-1}\mathrm{d}t+\epsilon \sqrt{2{\sigma }_{\chi }^2{T}_{\chi }^{-1}}\mathrm{d}W\left(t\right),$$(21)where dW(t) is the increment of standard Brownian process i.e., 〈dW(t)〉 = 0 and 〈dW(t)²〉 = dt. In equation (21), ϵ _sgs represents the locally resolved dissipation rate, and the parameters ${\sigma }_{\chi }^2$ and T _χ are chosen as functions of the mean eddy-viscosity (ν _t), the Kolmogorov length scale η and the spatial filter width (Δ):

$$\begin{array}{cc}{\sigma }_{\chi }^2=\ln \frac{\Delta }{\eta },& T_{\chi }=\frac{{\nu }_t}{{\Delta }^2}.\end{array}$$(22)

3.1 LES-STRIP model for diesel sprays

Because of the high injection velocities in diesel sprays, the flow in the near-nozzle region is strongly sheared. This results in very high levels of viscous dissipation along with very small magnitudes of local Kolmogorov’s dissipative scales. At the same time, the liquid is not yet well atomized, with liquid elements typically larger than the local Kolmogorov length scale. Equation (18) assumes the drag to be independent of the laminar viscosity. It is coherent with the classical observation for the sphere drag if Re > 1000. Therefore the condition for equation (18) is Re > 1000. Finally, the reformulated LES-STRIP model for diesel spray applications is given by:

For Re _p > 1000 and Δ > d _p > η

$$a_p=\frac{3}{4}{C}_d\left(\frac{{\rho }_f}{{\rho }_p}\right)\left(\frac{{\epsilon }^{\frac{2}{3}}}{d_p^{\frac{1}{3}}}\right){\vec{e}}_p.$$(23)

Otherwise for Re _p < 1000

$$a_p=\frac{\overline{u}-u_p}{{\tau }_p}+\sqrt{\frac{\epsilon }{{\tau }_p}}{\vec{e}}_p.$$(24)

4 Experimental and numerical setup

4.1 Experimental data

The non-vaporizing spray data from Engine Combustion Network, ECN [33, 34] is used . Three different test cases with do-decane (C₁₂H₂₆) as surrogate of the diesel fuel are used in this study. The parameters for different test cases are listed in Table 1. Apart from the comparison of the spray penetration length and spreading angle, the schlieren images of the spray evolution for test case-3, as shown in Figure 1 are used to validate the LES-STRIP model.

4.2 Numerical models

The simulations are performed in OpenFOAM, using a two-way coupled Eulerian-Lagrangian approach with exchange of mass and momentum between the phases. In Lagrangian spray simulations, the dynamics of the droplets representing the liquid fuel are tracked in a fixed Eulerian framework. A constant volume, cylindrical vessel with a radius of 25 mm and length of 100 mm is used in the simulations as shown in Figure 2. The vessel has a rectangular cross section (region 1) of 10 mm close to the injector to discretize the mesh with a fine grid to better resolve the spray dynamics. Initially a uniform grid size of 0.25 mm (referred to as Grid A) is used in both the regions to discretize the geometry. An additional fine mesh with grid size of 0.125 mm in region-1 and a non-uniform grid with cell size geometrically expanding from 0.125 mm to 0.25 mm in region-2 (referred to as Grid B) is simulated for studying the mesh sensitivity. The parcels are injected using a solid-cone injection method, where the parcels are sampled randomly from the user-defined initial size distribution. In this study Rosin-Rammler distribution model is used to represent the initial size distribution of injected parcels. The parcel injection velocity is calculated from the mass flow rate profile and its initial direction is randomly sampled to lie within a user-specified cone angle. Because of the high injection velocities, the secondary breakup occurs close to the injector and thereby the choice of an appropriate breakup model significantly effects the predictions of the spray evolution. In several studies [35–37], KHRT breakup model [38] has indicated its superior performance in predicting the spray structure under high pressure diesel conditions. The KHRT model includes two modes of breakup: KH breakup, accounts for the growth of unstable waves on the liquid jet due to differences in velocity between the ambient gas and the liquid; and RT breakup, accounts for growth of waves on the droplet’s surface due to aerodynamic drag at the droplet-gas interface. An Euler scheme for temporal discretion and a second order scheme for spatial discretization are used in the simulation. The injection and breakup model constants have been selected based on [39, 40], where influence of these parameters on the spray results have been reported in detail.

5 Results

5.1 Global spray characteristics

In order to evaluate the dispersion models, a quantitative comparison of the liquid penetration lengths, spreading angles for the different experimental conditions is presented in this section. In this study, the liquid penetration length 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. A comparison of liquid penetration length predicted for different test cases using Grid A is shown in Figures 3–5. The experimental penetration lengths presented in this study have an error of about 1–2%.

Initial penetration is strongly controlled by high momentum and secondary breakup of liquid droplets. So the initial under-prediction of the penetration lengths by both models in Figures 3–5 are overlooked in this study. But further downstream (i.e., after an axial distance of 20 mm away from injector), the LES-STRIP model gives more accurate prediction of the evolution of penetration length. Also it is clearly evident from Figures 3–5, that as the injection pressure increases the deviation in the penetration length predicted by the standard SGS model from the experiment increases significantly. Figures 6 and 7 present the grid density analysis of the SGS models for test case-3 using two different filter sizes i.e., Grid A and B.

Figure 6 shows that the standard SGS model requires much finer grid resolution to capture the evolution of penetration length. As grid resolution is increased, more and more scales are directly resolved by LES through $\overline{u}$. Therefore, the effect of SGS velocity component ($u\mathrm{\prime }\sim \sqrt{k_{\mathrm{sgs}}}\sim 0)$ is negligible, which clearly implies that the effect of unresolved scales is not well accounted by the model. On the other hand, from Figure 7 it can be seen that the STRIP model captures the evolution of penetration length even on coarse grid size. This is because, the formulation of STRIP model is based on modeling the droplet acceleration in terms of the ϵ which is a characteristic of the inertial scales. The model parameters of the stochastic equation for ϵ in equation (22) account for the effects of grid size on the evolution of ϵ, thereby making the dispersion model less sensitive to the grid size. The variation in the penetration length for different grid sizes as shown in Figure 7, is attributed primarily to the grid sensitivity of the KHRT breakup model. In addition to the penetration lengths, the spreading angle is a most commonly used parameter to quantify spray dispersion. The spreading angle is calculated from the isosceles triangle built considering 75% of the spray penetration, similar to [41]. For a much simplified analysis, an average spreading angle is calculated by averaging only the angles obtained when the penetration was greater than 35 mm. For all the simulations, the droplet parcels are injected within a cone angle of 20. From Figure 8 it can be seen that the LES-STRIP model consistently predicts much larger spreading angle (~4–5) than the standard SGS model for all the experimental test conditions.

Finally a direct comparison of the spray structure with the experimental schlieren images of test case-3 is presented in Figures 9–11.

Figure 9 shows the schlieren images of spray structure for test case-3 at different time instances. Figure 10 shows the spray structure for the corresponding time instances obtained using LES-STRIP model on Grid A. Even though dispersion is not evident close to the injector (upto 20 mm), overall spray structure is in good agreement with the experiment. On the other hand, Figure 11 represents the spray structure predicted by standard SGS model on Grid A. It can be noticed that the standard SGS model predicts a highly penetrative spray core surrounded by large number of parcels dispersed very close to the injector.

5.2 Droplet acceleration statistics

In order to explain the underlying physical differences in the results predicted by the two approaches, a more detailed analysis of the droplet acceleration statistics are presented in this section. All the results presented in this section are for the test case-3 using Grid A. In recent experiments in homogeneous statistically stationary turbulence [42–47], it was shown that the norm of acceleration is correlated on large times comparable with the integral time, while the direction is correlated on short times of order of the Kolmogorov time (τ _η). In this way the intermittency is manifested i.e., the vortical filaments of small scales are as much energetic as the large-scale turbulent structures. So a comparison of the autocorrelation function (${\rho }_{a_p}\left(\tau \right))$ of the droplet acceleration (a _p) for the two dispersion models is presented in Figures 12–13. The autocorrelation function of the droplet acceleration is calculated using equation (25):

$${\rho }_{a_p}\left(\tau \right)=\frac{\langle a_{p,k}\left(t+\tau \right)a_{p,k}\left(t\right)\rangle }{\langle a_{p,k}\left(t\right)a_{p,k}\left(t\right)\rangle },$$(25)where k = 1, 2 represents the axial and the radial components of acceleration respectively. The correlation time (τ) is normalized by the the Kolmogorov time scale (τ _η).

In the case of standard SGS model (Fig. 12), the axial component of droplet acceleration vector has much longer correlation (~8τ _η) than its radial component and norm. This shows that the droplets retain their axial acceleration component for a longer duration implying less radial dispersion. This is the consequence of the model employed in equation (12), where the norm of the acceleration obtained from sampling of the turbulent velocity (u′) is non-correlated in time. At the same time the direction of the droplet acceleration is retained on times defined by large turbulent structures. Therefore, a large number of droplets in the spray core only penetrate axially without any radial dispersion. Moreover, only the secondary droplets from the KHRT breakup which acquire an additional radial velocity component are dispersed radially as shown in Figure 11. On the other hand, in the case of LES-STRIP model (Fig. 13) both the axial and radial components of droplet acceleration correlate with the Kolmogorov time scales (τ _η), whereas the norm has much longer time correlation (5τ _η). The shorter time correlations of the components is because of the stochastic model, where the frequency of the fluctuations in orientation vector of droplet acceleration are scaled in the order of Kolmogorov time scales. Therefore, the LES-STRIP model is more effective in modeling the effects of intermittency at small scales on droplet dispersion.

A comparison of the PDF of the droplet acceleration at two different time instances i.e., t = 1.0 and 1.5 ms is shown in Figures 14 and 15. With increasing injection pressure, the probability of the high acceleration events increases. But as seen from Figures 14 and 15, the standard SGS model usually under-resolves the droplet acceleration (a _p). On the other hand, the acceleration PDF of LES-STRIP model has much broader tails indicating better resolution of the large acceleration events resulting from interaction of droplet with inertial scales.

5.3 Impact on the gas phase dynamics

Finally, the impact of the SGS dispersion models on resolved gas phase dynamics is investigated. A qualitative comparison is made by plotting iso-surfaces of the Q-criterion to visualize the vortex structures. A positive value of Q representing the local rotational motion of the fluid was chosen for plotting the iso-surfaces to visualize coherent vortices. Figure 16 shows the iso-surfaces of the Q criterion predicted by the different SGS models. Similar to the notable difference of the liquid spray structures as shown in Figures 7 and 8 the structure of the gas jets predicted by the SGS models are clearly different.

Next a comparison of radial profiles of the gas phase velocity at two different cross-sections i.e., 20 mm and 30 mm downstream are presented in Figures 17 and 18. The results represent an ensemble average of five different realizations of each model. The average velocity profiles from simulations without any dispersion model and standard SGS model are of the same order of magnitude. On the other hand the average velocity profiles from the LES-STRIP simulations differ by an order of magnitude. This qualitatively implies that in the case of LES-STRIP approach the spray dynamics and gas-flow field are strongly coupled. Since, LES-STRIP model directly models the droplet drag force, the effect of spray dispersion on the ambient gas dynamics is implicitly modeled through the source term (S _liq) represented by equation (11).

6 Conclusion

In this paper, the recently proposed stochastic equations (i.e., LES-STRIP) for the droplet motion in turbulent flows are used to model spray dispersion in diesel like conditions. The model uses the concept of stochastic drag for droplets larger than Kolmogorov scale and large droplet Reynolds numbers. On the other hand, for droplets smaller than Kolmogorov length scales the concept of stochastic forcing on droplet motion is employed. A detailed analysis of the performance of LES-STRIP in comparison with the most commonly used SGS dispersion model for three different non-vaporizing spray conditions is performed. In comparison with the standard dispersion model, LES-STRIP model gave good liquid spray penetration predictions for varying parameters such as injection pressure, and nozzle size. The grid-sensitivity analysis has shown that the LES-STRIP model is less sensitive to the grid size and can resolve the spray dynamics even on relatively coarse grid. In order to quantify spray dispersion, a comparison of the time averaged spreading angles for different test cases is presented. In comparison with the standard dispersion model, LES-STRIP has consistently predicted larger spreading angles. The higher dispersion in case of LES-STRIP is explained using the temporal correlation of the droplet acceleration. The shorter correlation of droplet acceleration components represent large fluctuations in orientation of droplet acceleration due to intermittency effects of small scales. Accounting for the intermittency effects by LES-STRIP results in larger spray dispersion. On the other hand, the higher spray penetration with standard dispersion model is attributed to longer correlation of the axial component of the droplet acceleration. Finally, to study the impact of the dispersion on two-way coupling of spray-turbulence interaction, a comparison of the gas flow velocities is presented. It was found that, while the results predicted without any dispersion model and with standard model are of same order of magnitude, the results of LES-STRIP differ by an order of magnitude. This qualitatively shows that the LES-STRIP model is more effective in modeling the two-way coupling of spray-turbulence interaction.

Acknowledgments

The authors would like to thank the Association Nationale de la Recherche et de la Technologie (ANRT) and the Volvo group for sponsoring the project. All the simulations in this study are performed on internal clusters of Volvo group and PMCS2I cluster at Ecole Centrale de Lyon. The authors would like to express their gratitude to Engine combustion network and Professor Raul Payri from CMT – Motores Térmicos, Universitat Politècnica de València for allowing access to some of their experimental data base.

References

All Tables

All Figures