Evaluation and validation of large-eddy simulation sub-grid spray dispersion models using high-fidelity volume-of-fluid simulation data and engine combustion network experimental data

Abstract

A sub-grid model accounting for the interaction of spray and sub-grid turbulence was developed and tested. The model predicts the sub-grid scale dispersion velocity used for calculating the slip velocity in Lagrangian–Eulerian Large-eddy simulation spray models. The dispersion velocity is assumed to be decomposed into a deterministic and a stochastic part, and it is updated in every turbulence correlation time for each computational parcel. The model was validated against two datasets: volume-of-fluid simulations and Engine Combustion Network experiments. The volume-of-fluid data showed that dispersion velocities at the centerline are anisotropic. This qualitative feature is well captured by the current model. For the Engine Combustion Network Spray A cases, it was found that sub-grid scale dispersion has profound impact on the prediction of the spatial distribution of liquid mass. Neglecting the sub-grid scale dispersion model results in underprediction of the width of the lateral projected liquid mass density profiles. Also, the prediction of the projected liquid mass density is sensitive to the two model constants determining the sub-grid scale dispersion velocity magnitude and turbulence time scale. However, the predictions of resolved gas-phase statistics are relatively insensitive to different sub-grid scale dispersion model setups. The primary reason for this was investigated. It was found that the motion of high-momentum liquid blobs in the near-nozzle region leading to air entrainment and subsequent gas jet development is minimally influenced by sub-grid scale dispersion. The importance of sub-grid scale dispersion inversely correlates with drag force magnitude: the larger the drag force, the less critical the sub-grid scale dispersion. Moving further downstream, quasi-equilibrium between the two phases is established, resulting in relatively small slip velocity and drag force.

Introduction

In internal combustion engines, the fuel injection process determines the mixing of fuel vapor and air, which greatly impact engine efficiency and emissions. Turbulent dispersion of fuel droplets is one of the many important physical mechanisms influencing the fuel–air mixing process. The existence of droplets may change the characteristics of turbulent flows, for example, change of the turbulent kinetic energy spectrum. In turn, a range of length and time scales of turbulent flow may have important effects on the dynamics of droplets. Fundamental studies on spray–turbulence interaction have been carried out by direct numerical simulations (DNS) and experiments in simplified flow configurations. Balachandar and Eaton ¹ reviewed some important aspects, such as preferential concentration of droplets and turbulence modulation from fundamental DNS and experimental studies.

Large-eddy simulation (LES) has become a popular approach in practical engine spray simulations in recent years. As pointed out by Rutland, ² there are several expected characteristics of an LES of engine flows. The primary expectation is that more vortices, eddies, and flow structures can be resolved. That is, LES provides the “large-eddy” flow field. More predicted flow structures might imply that LES is more suitable to be used to study more complex phenomena due to turbulence, such as cycle-to-cycle variability and engine knock. To computationally study these phenomena, accurate prediction of spatially and temporally evolving turbulent flows is required. Unlike DNS, which resolves flow length scales down to the Kolmogorov scales, LES only resolves filtered flow length scales. This means that the effect of flow length scales below the filter size on spray dynamics needs to be modeled. This model is commonly called a sub-grid scale (SGS) dispersion model.

A number of SGS dispersion models have been proposed and tested in non-engine flows. They can be categorized into stochastic and deterministic models. For the stochastic models, Bini and Jones ³ modeled the velocity increment of Lagrangian particles as the Wiener vector process pre-multiplied by a diffusion coefficient matrix, which is a function of the particle–turbulence interaction time and the sub-grid kinetic energy. This model was validated in dilute particle-laden mixing layer ⁴ and evaporating acetone sprays. ⁵ Pozorski and Apte ⁶ modeled the increment of the SGS dispersion velocity itself by the Wiener process. The coefficients in the stochastic partial differential equation (PDE) of the SGS dispersion velocity are related to the sub-grid kinetic energy. Preferential concentration of Lagrangian particles was quantified by a radial distribution function (RDF). This model was tested in particle-laden forced isotropic turbulence. It was found that the RDF predicted by this model is in good agreement with DNS data in a large Stokes number case, but it did not work well in a small Stokes number case. Also, the prediction is sensitive to the tunable constant determining the time scale of residual motion.

One of the deterministic models was proposed by Okong’o and Bellan. ⁷ The SGS dispersion velocity has a magnitude of the filtered standard deviation which is determined from the SGS stress model and the filtered velocity, and its sign is opposite to the Laplacian of the filtered velocity. Assumptions needed in this model were verified by the DNS data of an evaporating mixing layer laden with evaporating particles in a priori test. Later, they found that the performance of this model is much better than the stochastic Gaussian distribution method by comparing the DNS mixing layer data of velocities and droplet source terms in a priori test. ⁸ Another deterministic model is to determine the SGS dispersion velocity from the approximate deconvolution of a filtered velocity. The main idea is an approximation of the unfiltered field by truncated series expansion of the inverse filter operator.^9,10 This model was later used by Bharadwaj et al. ¹¹ to calculate the spray source term in the sub-grid kinetic energy transport equation.

The most commonly used SGS dispersion model for engine spray simulations is inherited from the Reynolds-averaged Navier–Stokes equations (RANS) dispersion model used in the KIVA engine simulation codes.^12,13 The SGS dispersion velocity is sampled once every turbulence correlation time from a Gaussian distribution function with variance proportional to the sub-grid kinetic energy. This type of model is widely used in LES of engine sprays.^14,15 However, studies on the impacts of the SGS dispersion model on engine spray simulations are few. One of the few studies was done by Jangi et al. ¹⁶ investigating the effects of the stochastic Gaussian dispersion model on the LES prediction of Diesel spray. They found that turning off the stochastic dispersion model results in overpredicted liquid penetration and the rate of air entrainment is slower. They also showed that including the stochastic dispersion model tends to reduce the sub-grid kinetic energy near the nozzle. The explanation is that excessive sub-grid kinetic energy produced by the spray source term is transformed into the liquid phase.

The objectives of this study are to evaluate and validate different SGS dispersion modeling strategies and to understand their effects on the prediction of spray–air mixing under Diesel engine conditions. Models are validated against not only available experimental data from Engine Combustion Network ¹⁷ but also data generated from high-fidelity volume-of-fluid (VoF) simulations. Recently, Eulerian–Eulerian high-fidelity spray simulations have emerged.^18–20 These types of studies not only provide insight into detailed atomization and breakup mechanisms in the near-nozzle region but also can provide useful data for model development. In this study, the modeled dispersion velocity is compared to the dispersion velocity directly computed from the solutions of the VoF simulations. In the context of LES Lagrangian–Eulerian spray simulations, it is known that several factors play a significant role in predicting the gas–liquid momentum exchange process. These include physics-based models such as atomization and breakup model^21–23 and sub-grid stress model ²⁴ and numerical schemes and mesh resolution.^15,25–27 However, as mentioned above, the effect of SGS dispersion modeling on the gas–liquid momentum exchange process is unclear and related studies are few. Thus, this is investigated in this study.

This article is organized as follows. The LES governing equations for liquid and gas phases and sub-grid stress models are first presented, followed by the description of the SGS dispersion model used in this study. The VoF method and equations are then described, and the test cases for validation and associated computational setups of LES and VoF are explained. The last two sections are “Results and discussion” and “Summary and conclusion.”

LES governing equations

In this section, LES equations governing spray dynamics are presented. There are several unclosed terms resulting from filtering operation that need models to close. These terms represent different physics. One of the terms is the sub-grid dispersion velocity, which is the main focus of this study. Development of the model to close this term is described in a separate section (section “SGS dispersion model”). Other important sub-grid models used in this study are introduced in this section.

This study uses the discrete parcel Lagrangian–Eulerian approach to simulate two-phase flow problems. Liquid droplets having the same properties are grouped into a parcel, and properties (position, velocity, mass, etc.) of each parcel are tracked. Gas-phase quantities are obtained from solving continuous Navier–Stokes equations. Based on the assumption that the parcel does not occupy any volume with respect to the gas phase, the parcel can be treated as point sources of mass, momentum, and energy for the gas phase.

In this section, gas and liquid momentum equations and sub-grid stress models are presented. Gas and liquid mass, species, and energy equations and heat transfer and vaporization models are detailed in Tsang and colleagues.^24,28

Gas-phase governing equations and sub-grid models

In LES, velocity $\mathit{u}(\mathit{x},t)$ is decomposed as

$$\mathit{u}(\mathit{x},t)=\overline{\mathit{u}}(\mathit{x},t)+\mathit{u}\prime (\mathit{x},t)$$

(1)

where

$$\overline{\mathit{u}}(\mathit{x},t)=\int G(\mathit{r},\mathit{x})u(\mathit{x}-\mathit{r},t)d\mathit{r}$$

(2)

$G$ is the low-pass filter to filter out high-frequency components, and $\mathit{u}\prime$ is the sub-grid velocity. For the compressible Navier–Stokes equations, the Favre averaging is frequently used. It is defined as $\tilde{\mathit{u}}=\overline{\rho \mathit{u}}/\overline{\rho }$ . Applying the filtering operation to the Navier–Stokes equations under the Lagrangian–Eulerian framework leads to the following momentum equation ¹¹

$$\frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial {\overline{\tau }}_{\mathrm{ij}}}{\partial x_j}-\frac{\partial \overline{\rho }{\mathrm{\Gamma }}_{\mathrm{ij}}}{\partial x_j}+{\overline{\stackrel{·}{S}}}_{\mathrm{momentum},i}$$

(3)

The subscripts $i$ and $j$ are index notation indicating Cartesian coordinates. Repeated indices of $i$ and $j$ imply summation. The symbol $\rho$ is the fluid density, $u$ the velocity, and p the thermodynamic pressure. The filtered viscous stress tensor ${\overline{\tau }}_{\mathrm{ij}}$ in the filtered momentum equation is written as

$${\overline{\tau }}_{ij}=\overline{\rho \nu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)}\approx \overline{\rho }\nu \left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}-\frac{2}{3}\frac{\partial {\tilde{u}}_k}{\partial x_k}{\delta }_{ij}\right)$$

(4)

where $\nu$ is the kinematic viscosity and ${\delta }_{\mathrm{ij}}$ is the Kronecker delta.

The momentum coupling between liquid and gas is represented by the source/sink term, ${\overline{\stackrel{·}{S}}}_{\mathrm{momentum}}$ , in equation (3). There are two contributions to this term, one from droplet vaporization and the other from drag force between the two phases. The contributions in a computational cell is summed over all Lagrangian parcels in that cell. This is given as

$${\overline{\stackrel{·}{S}}}_{\mathrm{momentum},i}=\frac{1}{V_f}\sum _{m=1}^{\mathrm{np}}{n}_{d,m}{\stackrel{·}{m}}_{d,m}{u}_{d,i,m}-n_{d,m}{F}_{d,i,m}$$

(5)

where $n_p$ is the number of parcels in the filter volume, $V_f$ ; $n_d$ is the number of droplets in a parcel; ${\stackrel{·}{m}}_d$ is the total vaporization rate of a droplet; $u_{d,i}$ is the droplet velocity; and $F_{d,i}$ is the drag force exerted on the droplet. The calculation of the droplet vaporization rate is detailed in Tsang and colleagues.^24,28 The drag force calculation is described in the next section.

The ${\mathrm{\Gamma }}_{\mathrm{ij}}$ term in the gas-phase momentum equation, equation (3), is the unclosed term resulting from the filtering operation. It is termed as the sub-grid stress tensor, representing interactions between resolved/filtered and unresolved/SGSs. It is of the form

$${\mathrm{\Gamma }}_{\mathrm{ij}}=\widetilde{u_i{u}_j}-{\tilde{u}}_i{\tilde{u}}_j$$

(6)

This study employs the dynamic structure sub-grid stress tensor model ²⁹ with the addition of a near-nozzle viscosity model. This model has been well validated under conditions of high-speed gas jet and spray jet. ²⁴ Prediction of fuel vapor penetration was greatly improved. It is given as

$${\mathrm{\Gamma }}_{\mathrm{ij}}=c_{\mathrm{ij}}{k}_{\mathrm{sgs}}-2{\nu }_{\mathrm{noz}}({\tilde{S}}_{\mathrm{ij}}-\frac{1}{3}{\tilde{S}}_{\mathrm{kk}}{\delta }_{\mathrm{ij}})$$

(7)

where ${\tilde{S}}_{\mathrm{ij}}$ is the filtered strain rate tensor

$${\tilde{S}}_{\mathrm{ij}}=\frac{1}{2}(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i})$$

(8)

The first term on the right-hand side in equation (7) was developed by Pomraning and Rutland. ²⁹ The form of the dynamic tensor coefficient $c_{\mathrm{ij}}$ can be derived following Germano’s ³⁰ identity, and the resultant form is given as

$$c_{\mathrm{ij}}=\frac{2L_{\mathrm{ij}}}{L_{\mathrm{kk}}}$$

(9)

where $L_{\mathrm{ij}}$ is the Leonard stress tensor, which is defined as

$$L_{\mathrm{ij}}=\widetilde{{\stackrel{~}{u}}_i{\stackrel{~}{u}}_j}-{\widetilde{\stackrel{~}{u}}}_i{\widetilde{\stackrel{~}{u}}}_j$$

(10)

More features of this model can be referred to Rutland and colleagues.^2,24,29,31 The sub-grid kinetic energy, $k_{\mathrm{sgs}}$ in equation (7), is obtained by solving its modeled transport equation

$$\frac{\partial \overline{\rho }{k}_{sgs}}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_j{k}_{sgs}}{\partial x_j}=-\overline{\rho }{\Gamma }_{ij}{\tilde{S}}_{ij}-C_e\overline{\rho }\frac{k_{sgs}^{3/2}}{\Delta }+\frac{\partial }{\partial x_j}\left(\overline{\rho }{\nu }_{sgs}\frac{\partial k_{sgs}}{\partial x_j}\right)+{\stackrel{.}{W}}_{s,sgs}$$

(11)

where $\Delta$ is the filter size chosen as cube root of the computational cell volume. The value of the SGS dissipation rate constant, $C_e$ , is taken to be 0.5. The sub-grid viscosity, ${\nu }_{\mathrm{sgs}}$ , in the transport equation is modeled as

$${\nu }_{\mathrm{sgs}}=C_k\Delta k_{\mathrm{sgs}}^{1/2}$$

(12)

where the value of the constant, $C_k$ , is taken to be 0.1. The last term on the right-hand side is the source/sink term due to the interaction of spray and sub-grid motions

$${\stackrel{·}{W}}_{s,\mathrm{sgs}}=\frac{1}{V_f}\sum _{m=1}^{\mathrm{np}}{n}_{d,m}{F}_{d,i}u{\prime }_{i,m}$$

(13)

Modeling of the SGS dispersion velocity, $u{\prime }_i$ , the main focus of this study, is described in section “SGS dispersion model.” Bharadwaj et al. ¹¹ found that the effect of this term cannot be negligible. Neglecting the ${\stackrel{·}{W}}_{s,\mathrm{sgs}}$ term resulted in underpredicted spray dispersion and overpredicted liquid penetration.

The second term in equation (7) is an additional viscosity term added to the original dynamic structure model. ²⁴ Unlike conventional viscosity-based models such as the Smagorinsky ³² model, the viscosity term here is only activated in the region where the strain rate is high. For sprays, the viscosity is only added in the near-nozzle region. The purpose of adding this term is to enhance sub-grid mixing in the near-nozzle region where a portion of energy-containing motions may not be adequately resolved under moderate grid resolution. This term is written as

$${\nu }_{\mathrm{noz}}=\{\begin{array}{c}\hfill C_{\mathrm{noz},1}\Delta k_{\mathrm{sgs}}^{1/2},\frac{\left|\tilde{S}\right|{\mathrm{\Delta }}^2}{\mathcal{U}\mathcal{L}}⩾C_{\mathrm{noz},2}\hfill \\ \multicolumn{1}{c}{0.0,\mathrm{otherwise}}\end{array}$$

(14)

where $\mathcal{U}$ and $\mathcal{L}$ are the characteristic velocity and length, respectively. For sprays, $\mathcal{U}$ is the injection velocity and $\mathcal{L}$ is the nozzle diameter. The filtered strain rate magnitude is written as

$$\left|\tilde{S}\right|={\left(2{\tilde{S}}_{\mathrm{ij}}{\tilde{S}}_{\mathrm{ij}}\right)}^{1/2}$$

(15)

The model constant $C_{\mathrm{noz},1}=0.25$ controls the magnitude of the nozzle viscosity, and the constant $C_{\mathrm{noz},2}=0.3$ controls the region where the viscosity is added.

Liquid-phase governing equations

The liquid phase is described by discrete Lagrangian parcels, with each parcel containing droplets having the same properties. Velocity of each parcel is updated by solving momentum equation for the drops

$$\frac{d{\mathit{u}}_d}{\mathrm{dt}}=\frac{{\mathit{F}}_d}{m_d}$$

(16)

and the parcel position ${\mathit{x}}_d$ is updated by

$$\frac{d{\mathit{x}}_d}{\mathrm{dt}}={\mathit{u}}_d$$

(17)

Mass and energy equations are described in Tsang and colleagues.^24,28 The drag force on a droplet, ${\mathit{F}}_d$ , is a function of relative velocity between gas and liquid, drag coefficient, and the frontal area of a droplet. It is written as

$${\mathit{F}}_d=\frac{1}{2}{\rho }_c{U}_{\mathrm{rel}}\pi r_d^2{C}_D({\mathit{u}}_c+\mathit{u}\prime -{\mathit{u}}_d)$$

(18)

where ${\rho }_c$ is the carrier phase density, $r_d$ the droplet radius, and ${\mathit{u}}_c$ the carrier-phase velocity evaluated at the droplet position obtained from an interpolation using adjacent cell values. $U_{\mathrm{rel}}$ is the magnitude of the relative velocity between the droplet and the gas

$$U_{\mathrm{rel}}=|{\mathit{u}}_c+\mathit{u}\prime -{\mathit{u}}_d|$$

(19)

Note that the drag force calculation depends on the SGS dispersion velocity, $\mathit{u}\prime$ , and the drag force, ${\mathit{F}}_d$ , appears in both gas and liquid momentum equations (equations (3) and (16)). This implies that the SGS dispersion is a two-way coupled dispersion. The drag coefficient in equation (18) is determined from

$$C_D=C_{D,\mathrm{sphere}}(1+2.632y_{\mathrm{distortion}})$$

(20)

where the drag coefficient is

$$C_{D.\mathrm{sphere}}=\{\begin{array}{c}\hfill \frac{24}{R{e}_d}(1+\frac{1}{6}{\mathrm{Re}}_d^{2/3}),R{e}_d<1000\hfill \\ \multicolumn{1}{c}{0.424,R{e}_d>1000}\end{array}$$

(21)

and the Reynolds number, $R{e}_d={\rho }_c{U}_{\mathrm{rel}}2r_d/{\mu }_c$ , where ${\mu }_c$ is the carrier phase viscosity and $y_{\mathrm{distortion}}$ is the distortion of a drop modeled by Taylor’s analogy between a drop and a spring–mass system. ³³

SGS dispersion model

The effect of sub-grid turbulent motions on spray dynamics are taken into account by the modeling of the SGS dispersion velocity, $\mathit{u}\prime$ . The SGS dispersion velocity is used extensively in the spray models. As seen from the previous section, it is used to calculate the slip velocity, ${\mathit{u}}_c+\mathit{u}\prime -{\mathit{u}}_d$ , which is used to determine droplet drag and trajectory, and it is used in the calculation of the sub-grid kinetic energy spray source/sink term (equation (13)). It also appears in the interphase transport term, equation (5), which means that the effect of the SGS dispersion has two-way coupling between the liquid and gas phases. In addition, non-dimensional numbers characterizing breakup and vaporization, such as the Reynolds number and the Weber number, are based on the slip velocity. In this section, detailed formulation of the SGS dispersion model used in this study is presented.

It is assumed that the dispersion velocity can be decomposed into the deterministic and the stochastic parts

$$\mathit{u}\prime ={\mathit{u}}_{\mathrm{sgs}}+{\mathit{u}}_{\mathrm{sto}}$$

(22)

The deterministic part, ${\mathit{u}}_{\mathrm{sgs}}$ , is modeled by the approximate deconvolution method^9–11

$${\mathit{u}}_{\mathrm{sgs}}=C_{\mathrm{sgs}}(2\tilde{\mathit{u}}-3\widetilde{\stackrel{~}{\mathit{u}}}+\stackrel{~}{\widetilde{\stackrel{~}{\mathit{u}}}})$$

(23)

where $C_{\mathrm{sgs}}$ is a model constant. The explicit filter is numerically implemented as the cell face area-weighted averaging. Theoretically, the approximate deconvolution method can only recover the SGSs on the order of the cutoff wave number. Information from smaller scales cannot be recovered, as discussed in previous studies.^6,7 In this study, the small-scale motion is represented by the stochastic part of the dispersion velocity. It is assumed to be isotropic and Gaussian distributed with zero mean

$$f\left(u_{\mathrm{sto},i}\right)=\frac{1}{\sqrt{2{\sigma }_{\mathrm{sgs}}^2\pi }}\exp (-\frac{u_{\mathrm{sto},i}^2}{2{\sigma }_{\mathrm{sgs}}^2})$$

(24)

where the subscript i is denoted as the ith component of a vector, and the variance ${\sigma }_{\mathrm{sgs}}^2$ is proportional to the sub-grid kinetic energy, that is

$${\sigma }_{\mathrm{sgs}}^2=C_{\mathrm{sig}}\frac{2}{3}{k}_{\mathrm{sgs}}$$

(25)

where $C_{\mathrm{sig}}$ is a model constant. From the equation above, equations (7) and (11) suggests that in the near-nozzle region where sub-grid mixing/kinetic energy is enhanced by the addition of the artificial viscosity, the sub-grid turbulence has larger effect on the fluctuation of the stochastic dispersion velocity in this region.

The SGS dispersion velocity of a parcel is sampled once every turbulence correlation time, $t_{\mathrm{turb}}$ . This concept was also used in the RANS dispersion model. ¹² In RANS, the turbulence correlation time is the minimum of the eddy breakup time and the time for a droplet to traverse a turbulent eddy. The length and the time scales in the RANS model are determined from the turbulent kinetic energy and its dissipation rate. The model can be extended to be used in the LES framework by simply replacing the turbulent kinetic energy by the sub-grid kinetic energy and the turbulent kinetic energy dissipation rate by the $k_{\mathrm{sgs}}$ dissipation rate, ${\epsilon }_{\mathrm{sgs}}$ , which is modeled as the second term on the right-hand side of equation (11) divided by the filtered density. Thus, $t_{\mathrm{turb}}$ can be written as

$$t_{\mathrm{turb}}=\min (\frac{k_{\mathrm{sgs}}}{{\epsilon }_{\mathrm{sgs}}},C_{\mathrm{ps}}\frac{k_{\mathrm{sgs}}^{3/2}}{{\epsilon }_{\mathrm{sgs}}}\frac{1}{|{\mathit{u}}_c+\mathit{u}\prime -{\mathit{u}}_d|})$$

(26)

where $C_{\mathrm{ps}}=0.16432$ is an empirical constant derived from the RANS $k-\epsilon$ model. Note that this constant could be re-derived for LES. However, the integral length scale concept in RANS must be replaced by a length scale in LES such as the filtered length scale, which has different physical significance. It is more common to not modify this model and model constant in the LES framework.^14–16

Another possible formulation of $t_{\mathrm{turb}}$ can be related to the deterministic part of the SGS dispersion velocity

$$t_{\mathrm{turb}}=C_{\mathrm{turb}}\frac{2\Delta }{|{\mathit{u}}_{\mathrm{sgs}}-{\mathit{u}}_d|}$$

(27)

where $C_{\mathrm{turb}}$ is a model constant, and $2\Delta$ represents the smallest eddy LES that can resolve. The physical meaning of equation (27) is the time for a droplet traveling through the largest unresolved eddy moving with the velocity of ${\mathit{u}}_{\mathrm{sgs}}$ . The underlying assumption of this formulation is that droplets are smaller than the largest unresolved eddy. Typically, the grid size is on the order of (usually larger than) the nozzle hole diameter in LES Lagrangian–Eulerian simulations, so the assumption is reasonable.

It is of interest to compare individual droplet trajectories predicted by equations (26) and (27). Figure 1 shows the trajectories of one liquid blob and its child droplets due to the KH breakup under the Engine Combustion Network (ECN) “Spray-A” conditions, as described in section “ECN Spray-A.” The trajectories predicted by equation (27) are apparently much more reasonable than those predicted by equation (26). We can see that in Figure 1(a), all the child droplet trajectories are straight lines. By examining the magnitude of $t_{\mathrm{turb}}$ of child droplets predicted by equation (26), we found that it is in the order of 10⁻³ s, which is the same order of magnitude of the injection duration. Thus, after a child droplet is born due to the Kelvin–Helmholtz (KH) breakup with a randomly chosen velocity vector orthogonal to the slip velocity vector of its parent blob, its dispersion velocity almost never gets updated or only gets updated once, resulting in the straight trajectories. The child droplets move radially outward and eventually leave the computational domain without any random walk behaviors. However, the turbulence correlation time predicted by equation (27) is on the order of 10⁻⁶ s, so the dispersion velocity of child droplets does get updated multiple times, resulting in the droplet trajectories more similar to the expected stochastic diffusion process. This observation does not restrict to a specific mesh resolution. Two more different grid resolutions (0.5 and 0.125 mm grid) were tested (results not shown here), and similar droplet trajectories as shown in Figure 1 were observed.

Figure 1.

Trajectories of one of the liquid blobs and its child droplets due to the Kelvin–Helmholtz breakup predicted by (a) equation (26) and (b) equation (27) with $C_{\mathrm{turb}}=0.5$ . Parent droplets/blobs are shown in solid trajectories, and child droplets are shown in dashed lines.

This comparison suggests that the simple modification from RANS-type spray models to LES may not work well without consideration of the different modeling methodologies between RANS and LES. For example, simple substitutions of RANS turbulent kinetic energy with LES sub-grid kinetic energy may not work well in all models. In the subsequent results, equation (27) is used to calculate the turbulence correlation time.

The effect of the variance of the Gaussian distribution and the turbulence correlation time on particle trajectories is studied in a simplified case before running full simulations. Consider a particle located at $(z,y)=(0,0)$ initially, moving with a constant mean speed $U$ in the z-direction, and the dispersion velocities in z- and y-directions are $\zeta$ and $\xi$ , respectively. The dispersion velocities are sampled once every constant time interval $h$ by the normal distribution with zero mean and a constant variance of ${\eta }^2$ . Assume that there is no drag, no collision, and no breakup. Then, after time T, the z- and y-positions of the particle are simply

$$z\left(T\right)=\mathrm{UT}+\sum _{i=1}^{T/h}{\zeta }_{i}h$$

(28)

and

$$y\left(T\right)=\sum _{i=1}^{T/h}{\xi }_{i}h$$

(29)

The y-position, $y(T)$ , is also Gaussian distributed with mean zero and variance of

$$\mathrm{var}\left(y\left(T\right)\right)=\frac{T}{h}{h}^2{\eta }^2=\mathrm{Th}{\eta }^2$$

(30)

Hence, larger variance, ${\eta }^2$ , or longer duration between two sampling events, $h$ , leads to larger variance of y-positions and thus larger spreading angles. This is verified by plotting sample paths of particles whose positions are determined by equations (28) and (29), as shown in Figure 2.

Figure 2.

Particle sample paths calculated from equations (28) and (29) with different variances, ${\eta }^2$ , and time intervals between two sampling events, $h$ : (a) $\eta =100\mathrm{m}/\mathrm{s},h=0.05\mathrm{ms}$ , (b) $\eta =100\mathrm{m}/\mathrm{s},h=0.0125\mathrm{ms}$ , and (c) $\eta =25\mathrm{m}/\mathrm{s},h=0.05\mathrm{ms}$ . The mean speed, U, is equal to 500 m/s.

To investigate how different modeling methodologies of SGS dispersion impact the prediction of spray dynamics, simulations with different model setups as listed in Table 1 were run. Case 1 does not use any dispersion model. Namely, the dispersion velocity $\mathit{u}\prime$ is zero. Case 2 only uses the deterministic model in which $\mathit{u}\prime$ is equal to the sub-grid velocity computed from equation (23) with $C_{\mathrm{sgs}}=1.0$ . Both Cases 3 and 4 use the deterministic and the stochastic models, but Case 4 has larger model constant, $C_{\mathrm{sgs}}$ , for the deterministic part. Cases 5 through 8 are used to study the effects of the magnitude of the variance of the Gaussian distribution and the turbulent correlation time. These values of the model constants were chosen to do the sensitivity study as they provide reasonable results compared to the ECN data, as presented in the subsequent sections. These values are not larger than one order of magnitude, which indeed suggests that the formulation of the dispersion model does have physical and mathematical significance.

Table 1.

Different SGS dispersion model setups to be tested.

Case no.	Model	$C_{\mathrm{sig}}$	$C_{\mathrm{turb}}$	$C_{\mathrm{sgs}}$
1	No model $\mathit{u}\prime =0$	N/A	N/A	N/A
2	Deterministic $\mathit{u}\prime ={\mathit{u}}_{\mathrm{sgs}}$	N/A	0.5	1.0
3	Deterministic + stochastic $\mathit{u}\prime ={\mathit{u}}_{\mathrm{sgs}}+{\mathit{u}}_{\mathrm{sto}}$	0.5	0.5	1.0
4		0.5	0.5	5.0
5		0.75	0.5	1.0
6		0.25	0.5	1.0
7		0.5	1.0	1.0
8		0.5	0.25	1.0

VoF method

The VoF simulations reported in this article are performed with an algebraic solver, interFoam, which forms a part of a larger open-source distribution of computational mechanics solvers and C++ libraries of OpenFOAM^®. ³⁴ The solver is based on a finite volume discretization on collocated grids for the solution of two-phase incompressible flow. A thorough evaluation of solver performance with respect to a broad range of two-phase flows is reported in our previous publication. ³⁵ The evaluation was based on the performance with respect to kinematics of advection, dynamics in inertia-dominated regime, and dynamics in the surface tension–dominated regime. An abbreviated description is provided here; a more detailed explanation can be found in Deshpande et al. ³⁵

The first part of the solution consists of advancing the liquid fraction field, $\alpha$ , by solving the following conservation equation

$$\frac{\partial \alpha }{\partial t}+\nabla ·\left(\mathit{u}\alpha \right)=0$$

(31)

The liquid fraction represents the volume fraction of liquid occupying a given computational cell, ${\mathrm{V}}_i,i\in [1,N_{\mathrm{cells}}]$ . The discrete version of this equation is

$$\frac{{\alpha }^{n+1}-{\alpha }^n}{\Delta t}+\frac{1}{\left|V_i\right|}\sum _{f\in \partial {\mathrm{V}}_i}{F}_u+{\lambda }_M{F}_c=0$$

(32)

where the fluxes are defined as

$$F_u={\varphi }_f^n{\alpha }_{f,\mathrm{upwind}}^n\mathrm{and}{F}_c={\varphi }_f^n{\alpha }_f^n+{\varphi }_{rf}^n{\alpha }_{rf}^n(1-{\alpha }_{rf}^n)-F_u$$

(33)

here, $n$ denotes time level, subscript $f$ refers to a cell-face quantity, ${\varphi }_f^n={\mathit{u}}_f^n·{\mathit{S}}_f$ , and ${\mathit{S}}_f$ is the outward normal vector corresponding to a given cell (not normalized). Since velocity (as well as $\alpha$ ) is cell-centered quantity, ${\mathit{u}}_f^n$ is obtained by weighted-averaging from cells sharing the given face. In the relevant region of the domain, where the spray exists, the grid size is uniform and isotropic; hence, the weight is 0.5. In the flux term, $F_u$ , the upwind value for the liquid fraction is denoted by ${\alpha }_{f,\mathrm{upwind}}^n$ . With respect to $F_c$ , ${\alpha }_f^n$ is determined from the second-order van Leer scheme. ³⁶ The remaining quantities represent the compressive flux, that is, ${\varphi }_{\mathrm{rf}}^n{\alpha }_{\mathrm{rf}}^n(1-{\alpha }_{\mathrm{rf}}^n)$ , where

$$\begin{array}{c}\multicolumn{1}{c}{{\varphi }_{\mathrm{rf}}^n=\underset{f\prime \in \partial {\mathrm{V}}_i}{min}(\frac{\left|{\varphi }_{f\prime }^n\right|}{\left|{\mathit{S}}_{f\prime }\right|},U_{\mathrm{rf},\max })({\mathit{n}}_f·{\mathit{S}}_f)\mathrm{and}{U}_{\mathrm{rf},\max }}\\ \multicolumn{1}{c}{=\underset{f\in \mathrm{V}}{\max }\left[\frac{\left|{\varphi }_f^n\right|}{\left|{\mathit{S}}_f\right|}\right]}\end{array}$$

(34)

This compressive flux is used to mitigate the effects of numerical diffusion that would occur as a result of the sharp gradients in $\alpha$ in the interfacial region. Finally, the variable ${\alpha }_{\mathrm{rf}}^n$ is obtained using the interface compression scheme native to OpenFOAM.^34,35 In numerical tests concerning the advection of a discontinuous profile, such as the $\alpha$ field, the treatment given above performs noticeably better than total variation diminishing (TVD) schemes with regard to the preservation of the sharpness of the discontinuity.

With respect to momentum, the following equation is solved

$$\frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·(\rho \mathit{u}\otimes \mathit{u})=-\nabla p_d+[\nabla ·(\mu \nabla \mathit{u})+\nabla \mathit{u}·\nabla \mu ]-\mathit{g}·\mathit{x}\nabla \rho +{\displaystyle \underset{\Gamma (t)}{\int }\sigma }\kappa \delta (\mathit{x}-{\mathit{x}}_s)nd\Gamma (x_s)$$

(35)

where the surface tension coefficient is given by $\sigma$ , local curvature by $\kappa$ , the gas–liquid interface by $\mathrm{\Gamma }(t)$ , the 3D Dirac Delta function by $\delta (\mathit{x}-{\mathit{x}}_s)$ , and ${\mathit{x}}_s$ is the integration variable over $\mathrm{\Gamma }(t)$ . The Continuum Surface Tension Model of Brackbill et al. ³⁷ is employed, giving

$$\underset{\mathrm{\Gamma }{\cap }^{V_i}}{\int }\sigma \kappa \delta (\mathit{x}-{\mathit{x}}_s)\mathit{n}\mathrm{d}\mathrm{\Gamma }\left({\mathit{x}}_s\right)=\underset{V_i}{\int }\sigma \kappa \nabla \alpha \mathrm{dV}$$

(36)

In the predictor step, the density and viscosity fields are regularized according to

$$\rho ={\rho }_l\alpha +{\rho }_g(1-\alpha )\mathrm{and}\mu ={\mu }_l\alpha +{\mu }_g(1-\alpha )$$

(37)

The solution of the momentum equation is obtained via a Pressure-Implicit with Splitting of Operators (PISO) ³⁸ iteration procedure. A predictor velocity is first constructed and then corrected to ensure momentum balance and mass continuity. Explicit formulation of the predictor velocity is a two-step process, where first the viscous, advective, and temporal terms in the momentum equation are used to generate a cell-centered vector field, which is then projected to cell faces using a second-order scheme. Contributions from surface tension and gravity terms are then added, concluding the predictor formulation. This procedure enforces a consistent discretization of surface tension and pressure gradient.^35,39

Within the correction procedure, the pressure distribution is added to the flux of predictor velocity, and mass conservation is invoked to yield a Poisson equation for pressure. The linear system is then solved using a preconditioned conjugate gradient method, with diagonal incomplete Cholesky as the preconditioner. In this work, we have used three PISO steps to arrive at predictions for $({\mathit{u}}^{n+1},p_d^{n+1})$ .

Cases for validation

Two cases were selected to validate the LES models and to study the effects of the SGS dispersion. The first case is the VoF case with lower injection velocity. Under this condition, accuracy of the VoF calculation has been validated and grid convergence achieved. LES-predicted SGS dispersion velocities and VoF-derived SGS dispersion velocities were compared in this case, as discussed in section “Low injection velocity case.” Operating conditions of the second case, Engine Combustion Network Spray-A, are representative of typical Diesel engines. This case is used to validate key global quantities such as spray tip penetrations and fuel vapor spatial distribution. Detailed operating conditions are described in the following sections.

Low injection velocity case

Operating parameters and fluid properties of the VoF case are summarized in Table 2. The fuel properties correspond to Iso-octane at 303 K and 4.5 MPa, and the air properties correspond to the same temperature and pressure. Vaporization is not considered in this case.

Table 2.

Operating parameters and fluid properties of the LIV case.

Nozzle diameter (µm)	90
Injection velocity (m/s)	50
Fuel density (kg/m3)	688.03
Fuel surface tension (N/m)	0.02
Fuel viscosity (Pa s)	4.78e−04
Air density (kg/m3)	50
Air viscosity (Pa s)	1.88e−05

ECN Spray-A

The ECN database, ¹⁷ maintained by Sandia National Laboratories, is an online library collecting a large amount of well-documented spray experiments from different institutions for computational model validations. In this study, experiments done by Sandia National Laboratories and Argonne National Laboratory were used to validate the models. In the Sandia experiments, a single-hole injector is mounted in a constant-volume cube. In-chamber air temperature, pressure, and density are well controlled and are typical of Diesel engine operating conditions. The root-mean-square (RMS) velocity is approximately 0.7 m/s due to the fan inside the chamber, which is much smaller than the spray injection velocity (∼500 m/s).

Several important quantities are measured by different experimental techniques.^40–42 These include temporally resolved liquid and vapor penetrations and ensemble-averaged spatial fuel vapor distribution. The Advanced Photon Source facility at Argonne National Laboratory provides data in the near-nozzle region ⁴³ and inside the nozzle. ⁴⁴ Ensemble-averaged liquid projected mass density (PMD) is one of the valuable data used for validating the models in the near-nozzle region. The PMD data were obtained from X-ray radiography path-integrated measurements. The data provide a two-dimensional projection of liquid fuel mass onto a plane whose normal vector is orthogonal to the spray axis. In this study, the “Spray-A” experimental data were used to validate the models. The Spray-A operating conditions are listed in Table 3.

Table 3.

ECN Spray-A operating conditions.

Fuel	100 % n-dodecane
Oxygen concentration	0%
Ambient gas temperature (K)	900 (vaporizing) and 300(non-vaporizing)
Ambient gas density (kg/m3)	22.8
Nozzle hole diameter(µm single hole)	90
Nozzle discharge coefficient	0.89
Fuel injection pressure (MPa)	150
Injection duration (ms)	6 (vaporizing) and 1.5(non-vaporizing)
Fuel temperature (K)	373

LES and VoF computational setup

Physical models and numerical setup (mesh and numerical methods/schemes) are equally essential to have a successful spray simulation. The computational setup and comparison details between the LES and VoF methodologies are thus described in this section.

LES computational setup

A two-dimensional cut-plane of the three-dimensional domain and grid sizes are shown in Figure 3. Note that all cells are isotropic, namely, $\Delta x=\Delta y=\Delta z$ . The grid size is refined to 250 µm in the main spray development region. Note that the grid size is larger than the nozzle diameter (90 µm). The authors’ previous work ²⁴ compared performance of different sub-grid stress models with moderate grid resolution and found that different sub-grid stress models predicted significantly different gas-phase quantities. This suggests the importance of sub-grid models in this type of simulations.

Figure 3.

2D cut-plane of the LES computational mesh. This mesh was used in the both cases described in section “Cases for validation.”

In this study, the computational fluid dynamics (CFD) toolbox, OpenFOAM, ³⁴ was used to solve the governing equations with the SGS and spray models implemented. The PISO pressure–velocity coupling algorithm ³⁸ is used for solving the LES Navier–Stokes equations. OpenFOAM uses collocated grid system with a remedy similar to the Rhie–Chow scheme ⁴⁵ for pressure/velocity oscillations. Euler implicit method is the time integration scheme for the gas-phase governing equations. The central cubic interpolation scheme is used for the convection term in the gas momentum equation. This combination was shown to be suitable for LES of free shear flows and sprays. ²⁶ The computational time step size was fixed at 2.5e−07 s, resulting in maximum Courant number of ∼0.3 to ensure numerical stability.

The time step size, $\Delta t_l$ , for solving the Lagrangian governing equations is restricted by the Lagrangian Courant number, which is the maximum fraction of the current cell that can be traversed by a parcel in a Lagrangian time step. The Lagrangian Courant number is taken to be 0.2. Additional restriction is applied to the time step size for the droplet momentum equation (equation (16)). The time step size is a minimum of $\Delta t_l$ and ${\tau }_d/10$ where ${\tau }_d$ is the droplet response time scale

$${\tau }_d=\frac{8}{3}\frac{1}{C_d}\frac{{\rho }_d}{{\rho }_c}{r}_d{U}_{\mathrm{rel}}^{-1}$$

(38)

The reason for this additional restriction is that the solution may become unstable due to the dramatic change of the dispersion velocity, $\mathit{u}\prime$ , in one time step resulting from the random sampling.

The stochastic Kelvin–Helmholtz/Rayleigh–Taylor (KH-RT) atomization and breakup model ²⁸ is used. Model constants such as the KH time constant and the RT size constant are determined stochastically and dynamically according to local conditions. Compared to the original KH-RT model, ⁴⁶ the stochastic model gave satisfactory simulation results in a wider range of operating conditions without tuning model constants. ²⁸ Initial velocity of Lagrangian parcels is determined by the synthetic eddy method.^28,47 This method simulates turbulence at nozzle outlet and generates time- and space-correlated injection velocity data without simulating flow inside the nozzle.

VoF computational setup

The VoF method is performed by the incompressible two-phase flow solver interFoam, as discussed in section “VoF method.” The computational domain along with the grid structure, which consists of hanging nodes, is shown in Figure 4. The domain extent is $49D\times 49D\times 360D$ along the respective $(x,y,z)$ directions, where $D$ is the jet diameter (90 µm). The mesh resolution in the spray region is $D/\Delta x=24$ , which is 67 times finer than the LES mesh. Away from this region, the mesh resolution systematically drops to $D/\Delta x=3$ , as illustrated in Figure 4. The grid dependency has been estimated through the calculation of intact liquid core length. What we have found is that for grid sizes below $\Delta x=3.75\mathrm{\mu }\mathrm{m}$ (the current resolution), the intact core length no longer changes appreciably.

Figure 4.

2D cut-plane of the VoF computational mesh; zoom-in view from top to bottom.

The identical physical setting used in LES is implemented in the VoF simulation. The fuel is iso-octane ( ${\nu }_f=6.95e-07m^2/s,{\rho }_f=688.03\mathrm{kg}/m^3,\sigma =0.02N/m$ ), and it is injected into a high-pressure air environment at 303 K and 44.85 bar $({\nu }_a=3.75e-07m^2/s,{\rho }_a=50\mathrm{kg}/m^3)$ through a 90-μm circular jet nozzle. The injection velocity $u_z$ is 50 m/s with uniform velocity distribution.

Method to compare LES and VoF data

In the low injection velocity (LIV) case, mean and RMS values of the SGS dispersion velocity modeled by LES and computed directly from VoF were compared. This section introduces how to obtain SGS dispersion velocity from VoF solutions and how to perform averaging in VoF and LES to obtain mean and RMS values.

Data processing in VoF

Computational volumes (shown in Figure 5) are employed to volume average the VoF data for use in LES. These volumes are of the same size as the LES mesh $(250\mathrm{\mu }\mathrm{m}\times 250\mathrm{\mu }\mathrm{m}\times 250\mathrm{\mu }\mathrm{m})$ , and their centroids reside along the spray axis (Y/D = 0) and off the spray axis (Y/D = 1).

Figure 5.

Computational volumes ( ${\mathrm{\Omega }}_1$ to ${\mathrm{\Omega }}_4$ ) with size $\Delta x_{\mathrm{LES}}^3$ , that is, $(250\mathrm{\mu }\mathrm{m}{)}^3$ , used in the volume-averaging of VoF simulation data for use in LES. The coordinates of the centroids are ${\mathrm{\Omega }}_1:(Y/D,Z/D)=(0,44)$ ; ${\mathrm{\Omega }}_2:(Y/D,Z/D)=(0,50)$ ; ${\mathrm{\Omega }}_3:(Y/D,Z/D)=(0,56)$ ; and ${\mathrm{\Omega }}_4:(Y/D,Z/D)=(0,61)$ . Note that additional four volumes (not shown here) are located at Y/D = 1 with the same axial locations of ${\mathrm{\Omega }}_1$ to ${\mathrm{\Omega }}_4$ .

The gas velocity, ${\mathit{u}}_g({\mathit{x}}_{\mathit{i}},t)$ , is first volume-averaged by performing the following operation

$${\mathit{u}}_g({\mathit{x}}_i,t)=\frac{\underset{{\mathrm{\Omega }}_i}{\int }\left(\rho I_g{\mathit{u}}_z\right)\mathrm{dV}}{\underset{{\mathrm{\Omega }}_i}{\int }\left(\rho I_g\right)\mathrm{dV}}$$

(39)

where ${\mathit{x}}_{\mathit{i}}$ is the coordinate of centroid for each computational volume ${\mathrm{\Omega }}_i$ , $t$ is the time, and $I_g$ is the gas indicator field, which is one if locally the computational cell is filled with gas and zero otherwise. Based on these values for velocity, the following first- and second-order statistical quantities are calculated, namely

$$〈{\mathit{u}}_g\left({\mathit{x}}_i\right)〉=\frac{1}{T}\underset{t_0}{\overset{t_0+T}{\int }}{\mathit{u}}_g({\mathit{x}}_i,t)\mathrm{dt}$$

(40)

$${\mathit{u}\prime }_g({\mathit{x}}_i,t)={\mathit{u}}_g({\mathit{x}}_i,t)-〈{\mathit{u}}_g\left({\mathit{x}}_i\right)〉$$

(41)

$$〈{\mathit{u}\prime }_g〉=\frac{1}{T}\underset{t_0}{\overset{t_0+T}{\int }}{\mathit{u}\prime }_g\mathrm{dt}$$

(42)

$${\mathit{u}\prime }_{g,\mathrm{var}}=\frac{1}{T}\underset{t_0}{\overset{t_0+T}{\int }}{({\mathit{u}\prime }_g-〈{\mathit{u}\prime }_g〉)}^2\mathrm{dt}$$

(43)

The beginning of the time integration is set $t_0=961\mathrm{\mu }\mathrm{s}$ , which occurs well beyond the initial liquid injection transient phase. The time integration period, $T=24.7\mathrm{\mu }\mathrm{s}$ , is equal to $99\Delta t_{\mathrm{LES}}$ , where $\Delta t_{\mathrm{LES}}$ is the LES time step size.

Data processing in LES

The parcel-based dispersion velocity, $\mathit{u}\prime$ , is converted into CFD cell-based quantities, ${\mathit{u}\prime }_p$

$${\mathit{u}\prime }_p({\mathit{x}}_p,t)=\frac{\sum _{m=1}^{\mathrm{np}}{\mathit{u}\prime }_m}{n_p}$$

(44)

where ${\mathit{x}}_p$ is the cell position and $n_p$ the number of parcels in a cell. Time averaging of ${\mathit{u}\prime }_p$ is performed to obtain mean and variance

$$〈{\mathit{u}\prime }_p〉=\frac{\underset{t_0}{\overset{t_0+T}{\int }}{\mathit{u}\prime }_p\mathrm{dt}}{T}$$

(45)

and

$${\mathit{u}}_{p,\mathrm{var}}=\frac{\underset{t_0}{\overset{t_0+T}{\int }}{({\mathit{u}\prime }_p-〈{\mathit{u}\prime }_p〉)}^2\mathrm{dt}}{T}$$

(46)

where $t_0$ and $T$ are the same as those in the VoF calculations. The RMS value is the square root of ${\mathit{u}}_{p,\mathrm{var}}$ . To evaluate the model performance, modeled mean (equation (45)) and RMS (equation (46)) values are compared against the VoF mean (equation (42)) and RMS (equation (43)) values.

Results and discussion

This section presents results from the two test cases described in section “Cases for Validation.” For the LIV case, SGS dispersion velocities modeled by LES and calculated by VoF are compared. For the ECN spray cases, LES-predicted liquid PMD and liquid and vapor penetrations are compared to the available experimental data. Also, the effects of the SGS dispersion on liquid and gas motions are discussed.

LIV case

Figures 6 and 7 show the comparison of the three components of the mean and RMS dispersion velocities at Y/D = 0 calculated by VoF (equations (42) and (43)) and predicted by LES (equations (45) and (46)) Case 3. The comparison of the LES Case 4 results and the VoF results at Y/D = 0 and Y/D = 1 are shown in Figures 8–11. The VoF data have been shifted to the left by 20 nozzle diameters for comparison. This is to compare LES and VoF data in the same breakup regime where breakup of liquid core and dramatic increase in gas–liquid interface area occur. In the VoF calculations,^20,48 this occurs in the region between the upstream initial instability development and downstream dispersed droplets. In LES, this region is modeled by the concept of the breakup length in the stochastic KH-RT model. Within the breakup length, only the primary KH breakup can occur; beyond the breakup length, both KH and the secondary RT breakup which tends to produce smaller drops can occur. It was found that the VoF calculation predicted longer liquid core, and thus, the breakup of liquid core and dramatic increase in gas–liquid interface area occurred further downstream than LES by 20 nozzle diameters. This may be due to the different boundary conditions between VoF and LES. The VoF simulation used a top-hat uniform velocity profile, while the LES simulation employed the blob injection with the synthetic eddy perturbation.

Figure 6.

Mean values of the SGS dispersion velocities atY/D = 0. The LES results are from Case 3 $(C_{\mathrm{sgs}}=1.0)$ listed inTable 1.

Figure 7.

RMS values of the SGS dispersion velocities atY/D = 0. The LES results are from Case 3 $(C_{\mathrm{sgs}}=1.0)$ listed in Table 1.

Figure 8.

Mean values of the SGS dispersion velocities atY/D = 0. The LES results are from Case 4 $(C_{\mathrm{sgs}}=5.0)$ listed in Table 1.

Figure 9.

RMS values of the SGS dispersion velocities at Y/D = 0. The LES results are from Case 4 $(C_{\mathrm{sgs}}=5.0)$ listed in Table 1.

Figure 10.

Mean values of the SGS dispersion velocities atY/D = 1. The LES results are from Case 4 $(C_{\mathrm{sgs}}=5.0)$ listed in Table 1.

Figure 11.

RMS values of the SGS dispersion velocities atY/D = 1. The LES results are from Case 4 $(C_{\mathrm{sgs}}=5.0)$ listed in Table 1.

Several observations from the VoF data are discussed as follows. The mean values in the x- and y-directions are close to zero at the spray axis due to the axisymmetric jet, which is expected, but the axial component is greater than zero. At the off-axis locations (Y/D = 1), the x and the y components of the mean dispersion velocities are also close to zero predicted by VoF and LES as shown in Figure 10. The Y/D = 1 location may be still close to the centerline at axial locations greater than Z/D = 20. Statistics at further off-axis locations were not computed due to the restricted region of grid refinement in the VoF calculation. The RMS values in the two lateral directions are not zero at the spray axis and off the spray axis (Y/D = 1) and are smaller than the RMS values in the axial direction. Note that the x and y components of the RMS values computed in VoF were not exactly matched. The difference stem from the fact that in the near field, the momentum interaction between phases is quite a bit more complicated than an interaction described by a populations of spherical elements as captured in the LES Lagrangian–Eulerian simulations. This complexity is the result of the dynamic shape of the interface in this region of the domain, which unavoidably affects the resulting turbulence and is reflected in different components for the SGS dispersion velocities. The VoF results revealed that the SGS dispersion velocity is anisotropic at and near the spray axis. The current model is able to capture the anisotropy. Note that the anisotropy must result from the deterministic part of the model (equation (23)) since the stochastic part (equation (24)) is assumed to be isotropic. However, Case 3 $(C_{\mathrm{sgs}}=1)$ underestimated the z-component mean and RMS values. By increasing $C_{\mathrm{sgs}}$ from 1 to 5 (Case 4), the z-component dispersion velocities can be effectively increased, as shown in Figures 8–11.

Note that the VoF results are dependent on the choice of the volume size and shape for the averaging operation (equation (39)). In this study, the volume is assumed to be equal to the LES cell volume. In the present implicit LES approach, discretization of derivatives acts as a low-pass filter, but the size and the shape of the filter volume are not explicitly defined. Thus, computational cell volume in LES being equal to the LES filter volume is indeed an assumption. This assumption is used in calculating model terms in which filter length scales are required, such as the SGS viscosity (equation (12)). Since there may be a reasonable range of volume sizes for the averaging operation in VoF, this study is more focused on qualitative comparison rather than tuning the dispersion model constants to perfectly match the current VoF data. For Case 4, qualitative features of the SGS dispersion can be well captured. These include anisotropy and x and y components of the RMS values being smaller than the z-component of the RMS values. Also, the LES and VoF results are in the same order of magnitude. In the next section, the role of the SGS dispersion model in practical Diesel spray simulations is investigated.

ECN Spray-A cases

Non-vaporizing and vaporizing ECN spray results are discussed in this section.

Non-vaporizing sprays

Simulated spatial distributions of liquid mass are compared against the X-ray ECN data as shown in Figures 12 and 13. The experimental results were time- and ensemble-averaged. Simulation results were time-averaged and spatial-averaged in the azimuthal direction which is the statistically homogeneous direction. The duration of the time averaging, 0.4–1.2 ms, is the same in the simulations and in the experiments. As shown in Figures 12 and 13, the liquid spray width is ∼0.8 mm at z = 5 mm and ∼2 mm at z = 10 mm. As shown in Figure 1(a), however, the child droplets can reach as far as 10 mm in the lateral direction. This is an overprediction. This comparison shows that the RANS turbulence correlation time model (equation (26)) predicts the existence of spray droplets at the outer region of the spray where no droplets were detected by experimental measurements. Again, this confirms the unreasonable droplet trajectories predicted by the RANS model.

Figure 12.

Lateral PMD profiles at z = 5 mm. Case 1: $u\prime =0$ ; Case 2: $u\prime =u_{\mathrm{sgs}}(C_{\mathrm{sgs}}=1)$ ; Case 3: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}}(C_{\mathrm{sgs}}=1)$ ; and Case 4: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}}(C_{\mathrm{sgs}}=5)$ .

Figure 13.

Lateral PMD profiles at z = 10 mm. Case 1: $u\prime =0$ ; Case 2: $u\prime =u_{\mathrm{sgs}}(C_{\mathrm{sgs}}=1)$ ; Case 3: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}}(C_{\mathrm{sgs}}=1)$ ; and Case 4: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}}(C_{\mathrm{sgs}}=5)$ .

The SGS dispersion model has impact on the prediction of liquid spray width as shown in Figures 12 and 13. Without using the stochastic dispersion model (Cases 1 and 2), the PMDs at y = 0 are significantly overpredicted, and the spray width is underpredicted in the downstream region (z = 10 mm). The stochastic model (Cases 3 and 4), however, predicts wider spray and improves the prediction at y = 0. These differences can be explained by Figure 14, showing the mass-weighted and spatially averaged mean and RMS values of the dispersion velocity at different axial locations at 0.6 ms. Mass-weighted and spatially averaged mean and RMS values of any properties of Lagrangian parcels, $\varphi$ , are computed as

$${\varphi }_{\mathrm{mean}}\left(Z\right)=\frac{\sum _{m=1}^N{m}_{d,m}{n}_{d,m}{\varphi }_m\chi (z_{d,m},Z)}{\sum _{m=1}^N{m}_{d,m}{n}_{d,m}\chi (z_{d,m},Z)}$$

(47)

and

$${\varphi }_{\mathrm{rms}}\left(Z\right)=\sqrt{\frac{\sum _{m=1}^N{m}_{d,m}{n}_{d,m}{({\varphi }_m-{\varphi }_{\mathrm{mean}}\left(Z\right))}^2\chi (z_{d,m},Z)}{\sum _{m=1}^N{m}_{d,m}{n}_{d,m}\chi (z_{d,m},Z)}}$$

(48)

where

$$\chi (z_d,Z)=\{\begin{array}{c}\hfill 1,Z-\Delta Z⩽z_d⩽Z+\Delta Z\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}}\end{array}$$

(49)

where $N$ is the total number of parcels in the computational domain, $z_d$ the z-coordinate of a droplet position, and $\Delta Z$ is 0.5 mm. Note that in Figure 14, intervals between data points are 1 mm, which means that all parcels in an interval were used for sampling since $\Delta Z$ is 0.5 mm. The RMS values of the dispersion velocity, shown by the error bars in Figure 14, in the two lateral directions (x and y) predicted by the deterministic model (Case 2) are much smaller than those predicted by the deterministic and the stochastic model (Case 3). Hence, it is less probable for a droplet to move further away from the centerline, resulting in narrower liquid sprays for the deterministic model. The stochastic part of the dispersion model enhances the diffusion of droplets, resulting in the improved PMD profiles. Next, the effect of the model constants, $C_{\mathrm{sig}}$ and $C_{\mathrm{turb}}$ , in the stochastic part of the model is investigated.

Figure 14.

Mass-weighted and spatially averaged dispersion velocities predicted by (a) Case 2 $(u\prime =u_{\mathrm{sgs}})$ and (b) Case 3 $(u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}})$ . The symbols are mean values, and the error bars are RMS values.

Figure 15 shows the qualitative comparison of the instantaneous PMD contours predicted by the different values of $C_{\mathrm{sig}}$ and $C_{\mathrm{turb}}$ . The display range of the PMD is from 0 to 1 µ/mm². The value of 1 is chosen as suggested by Pickett et al. ⁴¹ suggesting that the spray edge observed from optical diffused backlighting measurements correspond to the PMD of 1 µg/mm² measured from the X-ray measurements. Also, the PMD contours are shown up to z = 20 mm to compare the predictions in the near-nozzle region, although the penetration is ∼35 mm at 0.55 ms. As shown in Figure 15, larger values of $C_{\mathrm{sig}}$ or $C_{\mathrm{turb}}$ predict wider liquid sprays. This is consistent with Figure 2 in which the simplified case is considered. Note that these five cases in Figure 15 predict different spray widths and hence spray mass in the near-nozzle region (<20 mm), suggesting that the downstream liquid penetrations and spray mass (>20 mm) are also different since the total mass injected are the same for each case.

Figure 15.

Instantaneous PMD contours at 0.55 ms predicted by the different values of the model constants, $C_{\mathrm{sig}}$ and $C_{\mathrm{turb}}$ , in the stochastic part of the SGS dispersion model. Case 3: $C_{\mathrm{sig}}=0.5,C_{\mathrm{turb}}=0.5$ ; Case 5: $C_{\mathrm{sig}}=0.75,C_{\mathrm{turb}}=0.5$ ; Case 6: $C_{\mathrm{sig}}=0.25,C_{\mathrm{turb}}=0.5$ ; Case 7: $C_{\mathrm{sig}}=0.5,C_{\mathrm{turb}}=1.0$ ; and Case 8: $C_{\mathrm{sig}}=0.5,C_{\mathrm{turb}}=0.25$ .

Figures 16 and 17 show the quantitative comparison of the PMD profiles predicted by the different values of $C_{\mathrm{sig}}$ or $C_{\mathrm{turb}}$ . Different values of the two constants predict similar PMD values at y = 0. Comparing the lateral profiles at the two different axial distances, the prediction at z = 10 mm is more sensitive to the two constants than at z = 5 mm. Larger values of $C_{\mathrm{sig}}$ or $C_{\mathrm{turb}}$ predict more diffusion of droplets and thus wider profiles. This suggests that the SGS dispersion model has larger impact on the spatial distribution of liquid droplets having small inertia in the downstream region.

Figure 16.

PMD profiles at (a) z = 5 mm and (b) z = 10 mm predicted by the different values of the SGS dispersion model constant, $C_{\mathrm{sig}}$ . The constant determines the variance of the Gaussian distribution (equation (23)). Case 6: $C_{\mathrm{sig}}=0.25$ ; Case 3: $C_{\mathrm{sig}}=0.5$ ; and Case 5: $C_{\mathrm{sig}}=0.75$ .

Figure 17.

PMD profiles at (a) z = 5 mm and (b) z = 10 mm predicted by the different values of the SGS dispersion model constant, $C_{\mathrm{turb}}$ . The constant determines the turbulence correlation time (equation (25)). Case 8: $C_{\mathrm{turb}}=0.25$ , Case 3: $C_{\mathrm{turb}}=0.5$ , and Case 7: $C_{\mathrm{turb}}=1.0$ .

The results above show that SGS dispersion has a profound effect on the spatial distribution of liquid mass. Here, the effect of the SGS dispersion model on gas-phase solutions are examined. Qualitative comparison is made first by visualizing vortex structures by plotting iso-surfaces of the Q-criterion. The Q-criterion ⁴⁹ is defined as $1/2({\mathrm{\Omega }}^2-S^2)$ where $\mathrm{\Omega }$ is the rotational rate tensor and $S$ is the strain rate tensor. A positive value of Q, which means that the rotational motion dominates locally, was chosen for plotting iso-surfaces to visualize coherent vortices. Figure 18 shows the iso-surfaces of the Q-criterion predicted by the different model setups listed in Table 1. Unlike the notable difference of the liquid spray angles shown in Figure 15, the spreading angles of the gas jets predicted by the different dispersion models and the model constants are similar to each other. Unlike previous studies showing that different sub-grid stress models (model for the ${\mathrm{\Gamma }}_{\mathrm{ij}}$ term in equation (3)) lead to significant different predictions of flow structures and statistics,^24,50 the SGS spray dispersion model has relatively small effect on the prediction of gas-phase motion.

Figure 18.

Q-criterion iso-surfaces at 0.4 ms, $Q=1.5\times {10}^8{\mathrm{s}}^{-2}$ , colored by the vorticity magnitude predicted by the different model setups: Case 1: $u\prime =0$ ; Case 2: $u\prime =u_{\mathrm{sgs}},C_{\mathrm{sgs}}=1.0$ ; Case 3: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0$ ; Case 4: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0$ ; Case 5: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0$ ; and Case 6: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0$ .

Quantitative comparisons are made by plotting mean and RMS streamwise gas velocity radial profiles at three different axial locations, as shown in Figures 19 and 20. Time averaging from 0.4 to 1.2 ms and spatial averaging in the azimuthal direction were performed to obtain the mean and RMS values. There are some noticeable differences of the mean and RMS values predicted by the different model setups at the centerline and at z = 5 and 15 mm. At z = 5 mm, the RMS value at centerline predicted by Case 4 is lower than those predicted by the other cases. This might correspond to less small-scale structures predicted by Case 4 as shown in Figure 18. Nevertheless, the overall shapes of the profiles predicted by the different model setups are close to each other. At 25 mm, these profiles are almost on top of each other. Again, compared to the sub-grid stress model which has profound impacts on the prediction of gas entrainment as shown in the previous studies, these results suggest that the SGS dispersion model plays a less critical role in predicting air entrainment under Diesel spray conditions.

Figure 19.

Mean streamwise gas velocity profiles predicted by the different SGS dispersion model setups at three axial locations: (a) 5 mm, (b) 15 mm, and (c) 25 mm. Case 1: $u\prime =0$ ; Case 2: $u\prime =u_{\mathrm{sgs}}$ ; Case 3: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0$ ; Case 4: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0$ ; Case 5: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0$ ; and Case 6: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0$ .

Figure 20.

RMS streamwise gas velocity profiles predicted by the different SGS dispersion model setups at three axial locations:(a) 5 mm, (b) 15 mm, and (c) 25 mm. Case 1: $u\prime =0$ ; Case 2: $u\prime =u_{\mathrm{sgs}}$ ; Case 3: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0$ ; Case 4: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0$ ; Case 5: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0$ ; and Case 6: $u\prime =u_{\mathrm{sgs}}+u_{\mathrm{sto}},C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0$ .

The fact of the SGS dispersion model being less important in predicting gas-phase quantities is further examined. Figure 21 shows the time series of the convection, sub-grid stress, pressure gradient, and spray momentum source terms in the gas-phase z-momentum equation, predicted by Case 6 at different positions. The drag force term predicted by Cases 4, 5, and 6 are quite close to each other, meaning that the momentum budget predictions should also be similar. Hence, Case 6 should be sufficiently representative for the momentum budget analysis. In the upstream region, 5 mm from the injector, the magnitude of the spray momentum source term is comparable to the convection term, meaning a significant amount of momentum transfer from liquid to gas. Figure 22 also shows a sharp decrease in the drag force within 5 mm. Moving further downstream to 10 and 15 mm from the injector, the magnitude of the spray momentum source term is much smaller than the convection, sub-grid stress, and pressure gradient terms, and it is negligible at 15 mm. This corresponds to the small drag force downstream as shown in Figure 22. Thus, air entrainment and subsequent gas jet development due to spray injection is mainly initiated in the near-nozzle region (within ∼5 mm). In this region, liquid injected from the nozzle carry significant amount of momentum and experience little deceleration as shown in Figure 24. Figure 23 shows that within 5 mm, the dispersion ratio, $|\mathit{u}\prime |/\left||\tilde{\mathit{u}}-{\mathit{u}}_d|\right|$ , is smaller than that in the downstream region. This suggests that the motion of the liquid parcels upstream leading to air entrainment and subsequent gas jet development is less influenced by the SGS dispersion. This may be one of the reasons for the similar gas velocity profiles predicted by the different SGS dispersion model setups.

Figure 21.

Time series of the z-momentum budgets at (a) $(x,y,z)=(0,0,5)$ , (b) $(x,y,z)=(0,0,10)$ , (c) $(x,y,z)=(0,0,15)$ , and(d) $(x,y,z)=(0,2,15)$ predicted by Case 6. The injection position is $(x,y,z)=(0,0,0)$ . The labels “conv_z,”“sgs_z,”“gradP_z,” and “Smom_z” represent the convection term, the sub-grid stress tensor term, the pressure gradient term, and the spray momentum source term in the gas momentum equation (equation (1)), respectively.

Figure 22.

Mass-weighted and spatially averaged drag force magnitude versus distance from the injector predicted by Case 1 $(u\prime =0)$ , Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ . Case 1 without the stochastic forcing results in smaller drag force magnitude downstream.

Figure 23.

Mass-weighted and spatially averaged dispersion ratio versus distance from the injector predicted by Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ . At the upstream region, the dispersion ratio increases with $C_{\mathrm{sgs}}$ or $C_{\mathrm{sig}}$ , indicating that the resolved slip velocity does not change proportionally with the dispersion velocity.

In the downstream region (z > 10 mm), the dispersion velocity magnitude is comparable to the resolved slip velocity as shown in Figure 23, indicating that the dispersion model is more important in determining the droplet motion. Figure 24 shows that the larger the dispersion velocity magnitude, the larger the acceleration magnitude. This is also consistent with Figures 13 and 16(b) indicating that larger dispersion velocity results in more diffusion of droplets. Another note is that in the 15 mm downstream region, a quasi-equilibrium has been established between the two phases. This is evidenced by large acceleration of droplets (Figure 24) equivalent to short response time to ambient flows and negligible drag force term in the gas momentum equation (Figure 21(c) and (d)) indicating that slip velocity approaches zero.

Figure 24.

Mass-weighted and spatially averaged droplet accelerations versus distance from the injector predicted by Case 1 $(u\prime =0)$ , Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ .

The main point from Figures 21–24 is that the SGS dispersion model does not have significant effect on the air entrainment process, but has significant impact on the droplet motion and droplet spatial distribution downstream of the spray where the momentum transfer from liquid to gas is small. In fact, a correlation between the dispersion ratio, $|\mathit{u}\prime |/\left||\tilde{\mathit{u}}-{\mathit{u}}_d|\right|$ , and the drag force as shown in Figure 25 can be clearly seen: the more important the SGS dispersion, the smaller the drag force.

Figure 25.

Drag-dispersion ratio scatter plot predicted by Case 6 $(C_{\mathrm{sig}}=0.25)$ and Case 5 $(C_{\mathrm{sig}}=0.75)$ .

Additional note is on the singularity at a dispersion ratio of 1 in Figure 25. This occurs where droplets are at periphery of the spray where the resolved gas velocity, $\tilde{\mathit{u}}$ , is close to zero and is much smaller than the stochastic component of the SGS dispersion velocity ${\mathit{u}}^{\prime }$ . That is, the dispersion ratio, $|\mathit{u}\prime |/\left||\tilde{\mathit{u}}-{\mathit{u}}_d|\right|$ , can be approximated as $|\mathit{u}\prime |/\left|\left|{\mathit{u}}_d\right|\right|$ . The ratio of 1 indicates that droplet velocities at this region equilibrate with the sub-grid gas motion predicted by the model and thus small drag force (i.e. $\mathit{u}\prime -{\mathit{u}}_d\approx 0$ ).

Vaporizing sprays

Simulation results of the vaporizing Spray-A are compared against the available experimental data. As shown in Figure 26, larger $C_{\mathrm{sig}}$ or $C_{\mathrm{sgs}}$ results in shorter liquid penetration. This is consistent with the PMD results showing that larger $C_{\mathrm{sig}}$ or $C_{\mathrm{sgs}}$ results in larger liquid spray angle. Nevertheless, the difference in the liquid penetrations predicted by the different model setups is within 15%. This observation is consistent with the work done by Irannejad and Jaberi ⁵¹ also showing that the SGS dispersion model has small effect on the prediction of the liquid penetration. The prediction of Diesel spray liquid penetration is more sensitive to the atomization and breakup model, as shown in a number of studies.^21–23,52 Vapor penetration data as shown in Figure 27 do not show correlation with $C_{\mathrm{sig}}$ or $C_{\mathrm{sgs}}$ , and the different cases gave very similar results. Figure 28 compares mean fuel vapor mass fraction profiles. Except for the centerline fuel vapor mass fraction within z = 15 mm, the spatial distributions of fuel vapor mass predicted by the different models and the model constants are close to each other, which is consistent with the similar velocity profiles. From these vaporizing spray results, the SGS dispersion modeling being less important in predicting gas entrainment is again observed.

Figure 26.

Liquid penetrations predicted by Case 1 $(u\prime =0)$ , Case 3 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0)$ , Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ and measured by the ECN experiment.

Figure 27.

Vapor penetrations predicted by Case 1 $(u\prime =0)$ , Case 3 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0)$ , Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ and measured by the ECN experiment.

Figure 28.

Mean fuel vapor mass fraction profiles at (a) centerline, (b) 25 mm, and (c) 40 mm axial distances predicted by Case 1 $(u\prime =0)$ , Case 3 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=1.0)$ , Case 4 $(C_{\mathrm{sig}}=0.5,C_{\mathrm{sgs}}=5.0)$ , Case 5 $(C_{\mathrm{sig}}=0.75,C_{\mathrm{sgs}}=1.0)$ , and Case 6 $(C_{\mathrm{sig}}=0.25,C_{\mathrm{sgs}}=1.0)$ and measured by the ECN experiment. To compare the simulation results to the ensemble-averaged experimental data, time averaging from 1.5 to 2.5 ms and spatial averaging in the azimuthal direction were performed.

Summary and conclusion

This work is focused on the development and evaluation of the SGS dispersion model used in LES of engine sprays. The model computes the dispersion velocity term needed in the calculation of the slip velocity used in a number of spray models. It is assumed that the dispersion velocity is decomposed into a deterministic part and a stochastic part. The deterministic part is calculated by the approximate deconvolution method, which recovers the largest unresolved scales from the solution of the filtered velocity. Small-scale sub-grid motions are taken account by the stochastic part of the model. The stochastic part of the dispersion velocity is assumed to be isotropic and Gaussian distributed. The dispersion velocity is sampled once every turbulence correlation time, the time needed for a droplet to traverse an eddy in SGSs. A new form of the turbulence correlation time for LES was proposed to replace the RANS turbulence correlation time, which was shown to predict unreasonable droplet trajectories.

Mean and variance values of the SGS dispersion velocity were directly computed from the solutions of the high-fidelity VoF simulations. The injection velocity for the VoF case is approximately 1/10 of typical Diesel injection velocity. The data were compared to LES results. The VoF data revealed that both mean and variance of the SGS dispersion velocity are not isotropic, with magnitudes in the streamwise direction larger than those in the other two lateral directions. The anisotropic feature is well captured by the SGS dispersion model. The VoF and LES Case 4 results showed a good match.

Effects of SGS dispersion on Diesel spray conditions were studied using ECN data. It was found that neglecting the SGS dispersion resulted in the poor predictions of the spatial distribution of liquid mass. The stochastic part of the SGS dispersion model enhanced the diffusion of droplets and improved the prediction. In particular, larger values of the constants $C_{\mathrm{sig}}$ and $C_{\mathrm{turb}}$ associating with larger variance of the SGS dispersion velocity and longer turbulence correlation time led to wider lateral PMD profiles (larger liquid spray angle). However, different SGS dispersion model parameters predicted similar gas jet spreading angle, length scales of vortices, and mean and RMS values of gas velocity. The reason for this was further investigated. It was found that the SGS dispersion model is more important in determining the motions of droplets having small inertia and thus contributing less to momentum transfer from liquid to gas. However, for the liquid parcels having large inertia and initiating gas entrainment and subsequent gas jet development, the dispersion velocity magnitude is approximately 40% smaller than the resolved slip velocity. This ratio is large enough to influence centerline gas velocity prediction upstream, but overall the gas jet spreading and the downstream gas velocity profiles were little influenced by the SGS dispersion. Similarly, in vaporizing sprays, the predictions of vapor penetration and fuel vapor profiles are little influenced by the SGS dispersion model. These results suggest that SGS dispersion may have even smaller effect on subsequent ignition and combustion event.

It is suggested that one cannot neglect the SGS dispersion model in LES of Diesel engine sprays since it provides more reasonable prediction of the spatial distribution of liquid mass, although the model plays a less important role in predicting gas entrainment due to high-speed spray injection. One direction of future research is to examine whether the same conclusion is still applied to gasoline sprays.