<figure><img src="image_1.png" /></figure>

HAL
open science

# Three-component volume of fluid method coupling with interface compression method and Eulerian–Lagrangian spray atomization surface density model for prediction of cavitating sprays

P. Mcginn, G. Tretola, K. Vogiatzaki

► To cite this version:

P. Mcginn, G. Tretola, K. Vogiatzaki. Three-component volume of fluid method coupling with interface compression method and Eulerian–Lagrangian spray atomization surface density model for prediction of cavitating sprays. Physics of Fluids, 2024, 36 (2), pp.022118. 10.1063/5.0187353. hal-04653007

HAL Id: hal-04653007

https://ifp.hal.science/hal-04653007v1

Submitted on 18 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-component volume of fluid method coupling with interface compression method and Eulerian–Lagrangian spray atomization surface density model for prediction of cavitating sprays

P. McGinn ; G. Tretola ; K. Vogiatzaki

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for updates

*Physics of Fluids* **36**, 022118 (2024)

https://doi.org/10.1063/5.0187353

<figure><img src="image_5.png" /></figure>

View
Online

<figure><img src="image_6.png" /></figure>

Export
Citation

<figure><img src="image_7.png" /></figure>

AIP
Publishing

<figure><img src="image_8.png" /></figure>

Physics of Fluids

Special Topic:
Selected Papers from the 2023 Non-Newtonian
Fluid Mechanics Symposium in China

Submit Today

<figure><img src="image_9.png" /></figure>

AIP
Publishing
# Three-component volume of fluid method coupling with interface compression method and Eulerian-Lagrangian spray atomization surface density model for prediction of cavitating sprays F

Cite as: Phys. Fluids **36**, 022118 (2024); doi: 10.1063/5.0187353

Submitted: 13 November 2023 · Accepted: 16 January 2024 ·

Published Online: 23 February 2024

<figure><img src="image_10.png" /></figure>

View Online

<figure><img src="image_11.png" /></figure>

Export Citation

<figure><img src="image_12.png" /></figure>

CiteMark

P. McCinn,<fnref n="1" />,<fnref n="2" />,ᵃ) C. Tretola,<fnref n="3" /> ᵈ and K. Vogiatzaki<fnref n="3" />,⁵ ᵈ

## AFFILIATIONS

## ABSTRACT

In this work, a new three phase cavitation model with Eulerian-Lagrangian Spray Atomization (ELSA) and Interface Capturing Method (ICM) coupling is presented to allow for, in a unified approach, physical insight of the surface evolution of cavitating sprays at the Sub-Grid Scale (SGS). Phase change is accounted for in the framework via mass transfer across a liquid and corresponding vapor phase. The surface density model was validated against direct numerical simulation data of an atomizing jet with a systematic variation of mesh resolutions. The sensitivity of the adjustable parameters such as the critical Weber number was also demonstrated. Then, experimental comparisons were made with a more realistic orthogonal spray geometry within a cavitating and turbulent non-cavitating system. These comparisons include both large scale visualizations and small scale SGS quantities such as the Sauter mean diameter. To our knowledge, this is the first study that presents the performance of a three phase cavitating framework with SGS ELSA-ICM coupling.

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0187353

## I. INTRODUCTION

One of the dominant characteristics of applications relying upon multi-phase flows is the presence of evolving interfaces between the fluids and phases present. A well-known problem of such nature is linked to the challenges of the accurate prediction of spray injection dynamics that phase change occurs through either evaporation or cavitation/flash boiling. For example, in industrial systems, such as high pressure injection systems and pumps, a third fluid-vapor is generated and destroyed from cavitation phenomena taking place locally when high velocity flow, geometric conditions, and foreign gaseous particles present along solid boundaries cause local static pressure to reduce to saturation conditions where vapor is generated. This form of cavitation, manifested commonly as wall-attached cavitation, results in the generation of sharp surfaces between the liquid and the vapor that subsequently grow and evolve due to the immiscibility of the two fluids and local dynamic phenomena such as turbulence. Formed vapor is

immiscible with respect to the main liquid inside the injector but is miscible with the gaseous chamber environment when is transported outside the injector. It has been shown that internal nozzle cavitation has a significant influence on the primary atomization process of the fuel injectors<fnref n="1" />⁻⁵ due to the liquid/gas surface evolution process. When the formed vapor structures collapse violently, the surface energy is conserved and influences the evolution of the surrounding surfaces.<fnref n="3" />,⁶ First, the conserved energy further augments the collapse of other vapor cavities from the shock that is caused within the structure which is due to the bubbly mixture leading to a reduced speed of sound.⁶,⁷ Second, this energy is also transported downstream the nozzle exit where it influences the surface breakup of the atomization process.<fnref n="3" /> Third, the energy that is conserved from the violent collapse also causes material erosion and inevitable failure of the injector.⁸,⁹ It is essential then to be able to track the evolution of the various liquid interfaces at different

<footnote n="1">¹Advanced Engineering Center, University of Brighton, Cockcroft Building, Lewes Road, Brighton, United Kingdom</footnote>

<footnote n="2">²IFP Energies Nouvelles, 1 et 4 Avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote n="3">³University of Oxford, Department of Engineering Science, Oxford, United Kingdom</footnote>

<footnote>ᵃ) Electronic mail: paul.mcginn@ifpen.fr</footnote>

<footnote>ᵇ) Author to whom correspondence should be addressed: konstantina.vogiatzaki@eng.ox.ac.uk</footnote>
scales throughout the injection process to determine more clearly the role of cavitation in the atomization process.

Characteristics of the atomization process and cavitation can be identified at both large and small scales. As fuel atomizes, oscillations due to shear as well as due to cavitation generate pressure waves that wrinkle the interface leading to the detachment of large ligaments that then break into smaller droplets. These droplets have various sizes but as the atomization progresses they mostly exist at the micro-scale. Similarly, with cavitation, large vapor cavities form which typically have small bubble structures contained within. A limitation to modeling the whole process is that most studies isolate the internal nozzle and the spray atomization processes due to resolution restrictions. However, in doing so, a holistic understanding of how cavitation effects the subsequent surface dynamics and what is its wider role in the atomization process is missing and new insight into the interaction of the evolving surfaces at various scales throughout the whole process is needed. Geometric complexities make this challenging to achieve experimentally and compromises are required to do so. High-speed x-ray phase contrast imaging has been conducted on realistic injector nozzle geometries by Sforzo et al.² demonstrating a widening of the spray cone angle due to the presence of cavitation. The visualizations show the regions where cavitation occurs; however, structures cannot be clearly distinguished with no measurements made for any micro-scale structures. Suh and Lee³ attempt to overcome these limitations by obtaining shadowgraph visualizations of the whole flow field of an enlarged orthogonal nozzle of a Diesel injector which presents the large-scale structures. Alongside, droplet Sauter Mean Diameters (SMDs) are obtained by PDPA measurements to analyze the small-scale structures of the external spray. A range of cavitation regimes where explored which showed an increase in flow, and thus cavitation, intensity produces a widening of spray cone angle. Furthermore, smaller SMDs were observed when cavitation was present. A similar nozzle geometry to inject water was adopted by Sou et al.,⁴,¹⁰ which also demonstrated a widening of spray cone angle when the cavitation intensity was increased.

The dynamics of the evolving surfaces between the three fluids/ phases in cavitating spray flows (liquid, vapor, ambient gas), aside from the experimental visualization complexities described above, also present a number of numerical challenges. Numerical modeling of a three phase flow system is required, and the surface evolution must be tracked uniformly throughout the process, from the onset of cavitation internal to the nozzle to the external primary atomization process. A two fluid phase model proposed by Örley et al.¹¹ and further adopted by Mithun et al.¹² considers two transported phases, one as an isothermal atmospheric gas phase and the other as a barotropic mixture of liquid and vapor with a similar phase change mechanism to the Homogeneous Equation Model (HEM).¹³ The pressure wave dynamics from the collapse of vapor cavities can be modeled due to the compressible framework of these models; however, their reliability is limited due to the implicit definition of the cavitation interface, and a sharp, not diffusive, interface between the vapor and gas phases which is nonphysical. The two phase Volume Of Fluid (VOF) method, that is a common model used for spray atomization and phase change flows, can be expanded to consider a third transported phase.⁵ This still takes the one-fluid approach as a single equation for momentum, and pressure is resolved for all phases and is not treated as separate Eulerian phase systems. Such a three fluid model was proposed by Yu et al.¹⁴

which transports the volume fractions of all three phases and models the cavitation with common mass transfer phase change models. A similar solver was developed by Guissani et al.¹⁵ where the model underwent numerical benchmarking and was tested using the experimental geometry provided by Sou et al.¹⁰

The models described above reproduce the generated large-scale turbulent structures adequately; however, any structures that are present at the Sub-Grid Scale (SGS) are lost. In the context of Large Eddy Simulations (LES), a filtering operation is performed over the computational cell that dampens any numerical complexity below the grid scale. The Eulerian Lagrangian Spray Atomization (ELSA) model, developed by Vallet et al.,¹⁶ introduces a quantity that is transported alongside the phases and describes the liquid surface density. By transporting the surface density variable alongside the volume fraction, quantities such as the Sauter Mean Diameter (SMD) and droplet dispersion can be extracted. Comparisons made with Direct Numerical Simulation (DNS) data show that the ELSA model resolved the dense liquid part accurately.¹⁷ This approach was expanded by Navarro-Martinez²⁸ so that the phase volume fraction and surface density transport resolved using an SGS Probability Density Function (PDF) within quasi-iterations of stochastic fields in order to address better the SGS dynamics. This model was further developed and validated by Tretola et al.,¹⁸ which showed that the SGS structures such as the SMDs were better resolved using this approach. Lyras et al.,¹⁹ integrated the ELSA approach with the Homogeneous Relaxation Model (HRM) to model vaporization phase change to simulate flash boiling. All of these models, however, have been adopted within the two-phase VOF approach (only accounting for vapor and liquid) and do not consider the third phase (gaseous ambient phase) which plays an important role for phenomena occurring outside the nozzle.

In the present study, the integration of the ELSA-Interface Capturing Method (ICM) surface density approach into a new multi-fluid, incompressible framework accounting for phase change is presented. This couples the sub-grid generative and destructive surface dynamics with cavitation modeling which allows for the tracking of the evolution of all the surfaces present in atomizing spray systems. The main aspects of the model are

* Taking as starting point the one-fluid VOF method, and work done in our previous publication,²⁰ in the new framework presented here we incorporate the transport of three fluids (liquid, vapor, and non-condensable gas). The Schnerr-Sauer cavitation model is adopted to model the mass transfer across the interface of the liquid and vapor phases. Other mass transport models can be used; thus, the framework can be modified to account for alternative phenomena to cavitation such as flash boiling.

* The ELSA surface density approach is coupled and transported alongside the phase volume fractions. To our knowledge, it is the first implementation of the ELSA model in a three fluid context with phase change and phase miscibility. Most of the parameters of this model have been validated against DNS data though some need to be adjusted on a case-by-case basis such as the critical Weber number; this is explored in the validation.

* The ICM adaptation is adopted here to maintain a sharp liquid interface where necessary which is important when modeling cavitating flows.
## II. NUMERICAL MODELLING

### A. Multiphase modeling

To account for three phases in the context of VOF for liquid, vapor, and non-condensable gas, a transport equation for each phase volume fraction $\alpha$ is independently transported:

$$ 
\begin{aligned}
\frac{\partial \alpha_l}{\partial t} + \nabla \cdot (\alpha_l u) &= +S_{\alpha}, \\
\frac{\partial \alpha_v}{\partial t} + \nabla \cdot (\alpha_v u) &= -S_{\alpha}, \\
\frac{\partial \alpha_{nc}}{\partial t} + \nabla \cdot (\alpha_{nc} u) &= 0,
\end{aligned}
\quad (1) $$

where the subscripts l, v, and nc denote liquid, vapor, and non-condensable gas phases, respectively. The phase change mass transfer source terms are denoted by S which describes the mass transfer across the liquid-vapor interface. The above set of equations in their current form transports explicitly the “large” (i.e., above grid scale) fluid structures, while any sub-grid scale surface evolution is lost. Furthermore, due to discretization limitations, the interface of liquid structures becomes artificially diffusive. In Secs. IIB–IID, a summary of the mass transfer model that is used to predict internal nozzle cavitation is first presented. Then, the fundamentals of the ELSA approach are described, and it is explored how its implementation in the current framework will allow for the SGS surface density evolution to be tracked. Finally, the adoption of the ICM method is presented that allows for the sharpening of liquid structures when appropriate.

### B. Cavitation modeling

Cavitating flows are characterized by highly turbulent and chaotic phenomena that are non-linear in nature so describing the correct physical behavior at a given flow condition is very demanding. A number of different mass transfer cavitation models exist and are typically based on simplified forms of the Rayleigh-Plessent equation. The Rayleigh-Plessent equation estimates the radius growth of a spherical bubble due to the pressure acting on it.21 The Zwart, Singhal, and Schnerr Sauer model are all derived from the Rayleigh-Plessent equation and predict the mass transfer from the evolution of bubble dynamics.22,23 Our framework can accommodate any of the above models. Here, for demonstration purposes, we will focus on the Schnerr-Sauer model.

The onset of cavitation begins with the nucleation of small vapor bubbles. The Rayleigh-Plessent equation relates the radius growth of a perfect spherical bubble to the pressure dynamics.23 The Schnerr-Sauer model assumes that the mass transfer rate is proportional to the bubble radius growth $R_b$

$$ R_b = \sqrt[3]{\frac{1}{n_b} \frac{3}{4\pi} \frac{1-\alpha_l}{\alpha_l}}, \quad (2) $$

where $n_b$ is the number of bubbles $n_b$ in a given volume. The presence of the third non-condensable phase requires modification to the conventional Schnerr-Sauer model to accommodate this as follows:

$$ 
\begin{aligned}
\dot{m}_{-} &= C_{v} \frac{3 \rho_{1} \rho_{2}}{\rho_{m}} \frac{(1-\alpha_{l}-\alpha_{nc})}{R_{b}} \sqrt{\frac{2 p_{sat}-p}{3 \rho_{1}}}, \\
\dot{m}_{+} &= C_{c} \frac{3 \rho_{1} \rho_{2}}{\rho_{m}} \frac{(1-\alpha_{l}-\alpha_{nc})}{R_{b}} \sqrt{\frac{2 p-p_{sat}}{3 \rho_{1}}},
\end{aligned}
\quad (3) $$

where $\dot{m}_{-}$ and $\dot{m}_{+}$ are the volumetric mass transfer rates of vaporization and condensation, respectively. $C_v$ and $C_c$ are the respective tuning coefficients. For cavitating flows, these coefficients can typically be set to a value of 1 for highly cavitating flows.

### C. Eulerian Lagrangian spray atomization (ELSA)

The scales at which cavitation and atomization occur are small [O (μm)–O(nm)] and can be smaller than a typical cell size in LES. Thus, the liquid structures resulting from these processes cannot be accounted for accurately within a simple VOF-LES framework since the SGS liquid interface is averaged over the cell due to the filtering operation. Moreover, droplets break up into increasingly small sizes during atomization, and such structures are difficult to distinguish since VOF methods do not provide a description of the effect the breakup process has on the SGS wrinkling of the liquid interfaces. Distortion of the interface is only captured at the cell scale level. Vallet et al. introduced the idea of the ELSA model that solves, alongside the liquid mass fraction (from VOF), an additional transported quantity for surface density $\Sigma$.16 This quantity defines the amount of liquid surface that is present within the computational cell volume and, thus, gives a description of the interface at the SGS. For instance, if in a cell the liquid volume fraction is low but the surface density is high, this potentially indicates that in this cell there is a large number of small droplets. A high liquid volume fraction with low surface density, on the other hand, indicates that only a few larger liquid structures with a smooth interface are present across the cell. Figure 1 shows a schematic to clarify how to two cells with the same volume of liquid can have different surface quantity.

Although the initial idea of the ELSA model can help increase the accuracy of the modeling of the evolving fluid surface, its numerical implementation has some challenges. A characteristic of wall attached cavitation is the sharp interface present along the surface of the cavity structure. One of the limitations of the conventional ELSA is that introduces a diffusive interface from local fluctuations in the velocity field. To overcome this limitation, a coupled ELSA with Interface Capturing Method (ICM) approach is adopted in this work that introduces criteria to determine locally whether the flow should be resolved using ELSA or conventional VOF with sharp interface compression. When the flow is dilute, this compressive interface definition is switched off and the ELSA modeling is then applied. An Interface Resolved Quality (IRQ) quantity sets a threshold for this switch based on two criteria relevant to the local surface density $\Sigma$ distribution and curvature of the Volume fraction $\alpha$ interface.

A description of the surface density $\Sigma$ along with the volume fraction $\alpha$ allows extraction of relevant parameters. This includes the SMD

<figure><img src="image_13.png" /><figcaption>FIG. 1. Examples of two different surface densities within a computational cell with the same volume fraction: a) shows a low surface density, this case represents $\Sigma_{eq}$ as this is the minimum surface in a cell and b) shows a high surface density that indicates the presence of smaller structures within the cell.</figcaption></figure>
which provides quantitative comparisons with experimental results. Using the surface density $\Sigma$ to determine the SMD is considered in two ways. For bubbly flows with a very small volume fraction, the SMD is defined as $d_{32} = 6(1 - \alpha_l)/\Sigma$. For mono-dispersed spray consisting of spherical droplets, the SMD is defined as $d_{32} = 6\alpha_l/\Sigma$. All these structures can be accounted for using a combined definition for the SMD:

$$d_{32} = 6 \frac{\alpha_l (1 - \alpha_l)}{\Sigma}. \qquad (4)$$

To formulate a transport equation for $\Sigma$, first we need to consider the interface Dirac function or grained surface density

$$\delta_I = -\mathbf{n}_k \nabla \chi_k, \qquad (5)$$

where $\chi_k$ is the binary phase indicator function. $\delta_I$ is a Dirac distribution with units of inverse length, $m^{-1}$, and is sometimes termed local instantaneous interfacial area concentration.²⁴ It has been shown by Lhuillier et al.²⁵ and summarized by Morel²⁴ that the transport equation of $\delta_I$ can be directly derived from this distribution giving

$$\frac{\partial \delta_I}{\partial t} + u_j^I \frac{\partial \delta_I}{\partial x_j} = -\delta_I n_i n_j \frac{\partial u_i}{\partial x_j}. \qquad (6)$$

The assumption that the other volume fields is considered smooth where SGS effects are not considered unless modeled in the surface density transport equation is made. This volume averaging operation is typically applied to a volume which defines the threshold width between the sub-grid and resolved scales, but here, it is considered to be smaller than any length scale. This operation performed on a generic function $\phi$ is defined as

$$V[\phi] = \hat{\phi} = \frac{1}{|\delta\Omega|} \int_{\delta\Omega} \phi dx'. \qquad (7)$$

Applying the volume averaging operation, Eq. (7), to the evolution equation of $\delta_I$, Eq. (6), yields the first form of the interface surface density transport equation

$$\frac{\partial \Sigma}{\partial t} + \frac{\partial u_j^I \Sigma}{\partial x_j} = \frac{1}{|\delta\Omega|} \int_{\delta\Omega} \delta_I n_i n_j \frac{\partial u_i}{\partial x_j} dV. \qquad (8)$$

The volume averaging operation artificially thickens the interface to a thickness that is smaller than all resolved scales. The lengthening of the liquid interface due to local flow gradients and surface effects is described in the term on the RHS of Eq. (8). The interface normal is not a continuous function and discontinuities arise from liquid breakup effects and surface coalescence. If the integral was to be performed over a large volume than this integration requires separation. The interface velocity $u^I$ needs to be closed if the control volume V that the transport equation for $\delta_I$ is larger than the integral volume $\delta\Omega$, $V \gg \delta\Omega$. The problem with the closure of the surface velocity has been summarized by Delhaye.²⁶ A proposed closure formulation from Séro-Guillaume and Rimbert²⁷ is given by

$$u_j^I = u_j - \frac{\sigma}{T} V L \frac{\partial \hat{\Sigma}}{\partial x_j}. \qquad (9)$$

The first term on the RHS is the flow field velocity and the second is a restorative velocity term that diffuses the surface. If the local flow

field is in thermodynamic equilibrium or if the interface velocity is uniform, then this term disappears. Similar simplifications can be considered for very dilute flows where liquid structures are generally much smaller than the volume. The derivation of such simplifications has led to the development of surface density models that consider droplet breakup and bubble coalescence. These events are modeled as source terms, and several have been adopted co-currently in the literature. However, experimental validation for these models is extremely limited. Thus, a common approach has been to limit the number of source terms as much as possible. This can be achieved by having a source term representing surface generation/destruction and a source term for the presence of two co-existing phases. The same general form arises from literature for the surface density transport equation that accounts for these considerations which was adopted by Vallet and Borghi,¹⁶ Navarro-Martinez,²⁸ and Anez et al.²⁹ which is used in this work, is defined as

$$\frac{\partial \Sigma}{\partial t} + \frac{\partial u_j \Sigma}{\partial z_j} = \Sigma_{\text{mix}} + \Sigma_{\text{int}} = \mathbf{S}_{\Sigma}, \qquad (10)$$

where $\Sigma_{int}$ is the production/destruction of the surface caused by turbulent effects, collision events, coalescence, etc. $\Sigma_{mix}$ is the surface generation caused by liquid gas mixing and exits due to the presence of co-existing phases. Previous studies have explored various approaches in which to model this term.¹⁶,¹⁷,²⁹ The common outcome was that the term has a negligible effect on the overall solution of the surface density thus is not considered here. Vallet et al.¹⁶ propose the restoration to equilibrium formulation for the combined source term $\Sigma_{int}$ that has been adopted by a number of other ELSA models that is given by

$$\Sigma_{\text{int}} = C_{\Sigma} \frac{\Sigma}{\tau_{\Sigma}} \left( 1 - \frac{\Sigma}{\Sigma_{\text{eq}}} \right), \qquad (11)$$

where it is assumed that the equilibrium surface density value $\Sigma_{eq}$ is reached within the characteristic timescale $\tau_{\Sigma}$, which is proportional to a flow timescale.³⁰ Assuming the local surface energy is in equilibrium with local kinetic energy $k = \frac{u^2}{2}$ and isothermal flow, the estimation of the equilibrium surface density $\Sigma_{eq}$ can be defined as (neglecting viscous stresses)

$$\Sigma_{\text{eq}} \approx \frac{\rho k \alpha}{\sigma}. \qquad (12)$$

Deviation from the local equilibrium surface density can be characterized by a critical Weber number

$$We^* = \frac{\rho k \alpha}{\sigma \Sigma_{\text{eq}}}. \qquad (13)$$

In dilute flow where a distribution of small droplet sizes is present, the different droplets have to be represented by the source terms for the generation and destruction. The associated phenomenon for these terms include droplet breakup, collision, and coalescence where each has a unique local equilibrium value of $We^*$ and timescale $\tau_{\Sigma}$. These different phenomena are considered within the source terms in conventional ELSA models, where a weighting factor for each is based on the volume fraction to distinguish between dilute and dense regions. Duret et al. proposed that the local equilibrium $We^*$ to be used in combination to consider both dense and dilute flow.
Inspection of Eq. (11) presents no clear definition of the surface at rest, and the source term component $\Sigma_{mix}$ is neglected. Here, there is no kinetic energy present; thus, the equilibrium surface $\Sigma_{eq}$ vanishes and the surface is destroyed instantaneously. To overcome this limitation, Chesnel et al.³¹ proposed a quantity for the minimum surface $\Sigma_{min}$ in the presence of two phases

$$ \Sigma_{min} = C_{min} \alpha_l (1 - \alpha_l), \quad (14) $$

where $C_{min}$ is a constant defined as $C_{min} = 2.4$ based on DNS calculations. The consideration of a minimum surface value $\Sigma_{min}$ allows for the definition of surface density to be decomposed as

$$ \Sigma = \Sigma_{min} + \Sigma', \quad (15) $$

where $\Sigma'$ is the additional surface at the SGS through droplet interaction and surface evolution and is needed to complete the total surface density. Thus, only this term is transported. Combining Eqs. (10) and (15) leads to the transport equation for $\Sigma'$ that is considered in this work

$$ \frac{\partial \Sigma'}{\partial t} + \frac{\partial u_j \Sigma'}{\partial x_j} = S_{\Sigma}. \quad (16) $$

Before a definition is made for coupling the liquid interface at the SGS with the large scale phase volume fraction transport within the LES context, first, we should consider the filtered continuity equations for LES

$$ \frac{\partial \bar{u}_j}{\partial x_j} = 0, \quad (17) $$

$$ \frac{\partial \bar{\rho} \bar{u}_i}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_i \bar{u}_j}{\partial x_j} = - \frac{\partial \bar{p}}{\partial x_i} + \frac{2 \mu \bar{S}_{ij}}{\partial x_j} - \frac{\partial \tau_{ij}^{\text{sgs}}}{\partial x_j} + \sigma \bar{\kappa} \nabla \bar{\alpha} \bar{\delta}_I, \quad (18) $$

where the unknown sub-grid stress tensor, $\tau_{ij}^{\text{sgs}}$, which is a result of splitting the advection term before filtering, $\bar{u}_i \bar{u}_j = \tau_{ij} + \bar{u}_i \bar{u}_j$, is considered proportional to the filtered strain rate tensor and turbulent eddy viscosity, $\tau_{ij}^{\text{sgs}} = -\nu_t \bar{S}_{ij}$. Applying the same LES filtering operation to the phase volume fraction $\alpha$ transport equations and the surface density transport leads to the following set including closer terms

$$ \begin{aligned} & \frac{\partial \bar{\alpha}_l}{\partial t} + \frac{\partial \bar{u}_j \bar{\alpha}_l}{\partial x_j} = \frac{\partial [\bar{\alpha}_l (\bar{U} - \bar{U}_{\alpha_l})]}{\partial x_j} + \mathbf{S}_{\alpha}, \\ & \frac{\partial \bar{\alpha}_v}{\partial t} + \frac{\partial \bar{u}_j \bar{\alpha}_v}{\partial x_j} = \frac{\partial [\bar{\alpha}_v (\bar{U} - \bar{U}_{\alpha_v})]}{\partial x_j} - \mathbf{S}_{\alpha}, \\ & \frac{\partial \bar{\alpha}_{nc}}{\partial t} + \frac{\partial \bar{u}_j \bar{\alpha}_{nc}}{\partial x_j} = \frac{\partial [\bar{\alpha}_{nc} (\bar{U} - \bar{U}_{\alpha_{nc}})]}{\partial x_j}, \\ & \frac{\partial \bar{\Sigma}'}{\partial t} + \frac{\partial \bar{u}_j \bar{\Sigma}'}{\partial x_j} = \frac{\partial [\bar{\Sigma}' (\bar{U} - \bar{U}_{\Sigma'})]}{\partial x_j} + \mathbf{S}_{\Sigma}, \end{aligned} \quad (20) $$

where the superscripts $l$, $v$, and $nc$ denote the liquid, vapor, and non-condensable gas phase, respectively. The first two terms on the RHS of these equations are closure terms required from the filtering operation. These can be thought of the additional fluctuations at the SGS which is unresolved. This term for the volume fraction is known as the turbulent liquid flux where the volume fraction is transported via velocity fluctuations. The closure term for the surface density is considered in the same way where the velocity fluctuations transport the surface density. These turbulent fluxes can be closed using a classical gradient law¹⁶

$$ \bar{\alpha}_n (\bar{U} - \bar{U}_{\alpha_n}) = \mathfrak{R}_{\alpha_n}^{\text{sgs}} = \frac{\partial}{\partial x_j} \frac{\nu_{\text{sgs}}}{S_{\text{sgs}}} \bar{\alpha}_n, \quad (21) $$

$$ \bar{\Sigma}' (\bar{U} - \bar{U}_{\Sigma'}) = \mathfrak{R}_{\Sigma'}^{\text{sgs}} = \frac{\partial}{\partial x_j} \frac{\nu_{\text{sgs}}}{S_{\text{sgs}}} \bar{\Sigma}', \quad (22) $$

where $n$ denotes the phase volume fraction $l$, $v$, or $nc$. Substituting the above equations, Eqs. (21) and (22), alongside the source terms, Eqs. (1) and (11), leads to the following set of equations that are solved:

$$ \begin{aligned} & \frac{\partial \bar{\rho} \bar{\alpha}_l}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_l}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ D_{\alpha_l} \frac{\partial \bar{\alpha}_l}{\partial x_j} \right] + \dot{m} \left( \frac{1}{\rho_l} - \alpha_l \left( \frac{1}{\rho_l} - \frac{1}{\rho_v} \right) \right), \\ & \frac{\partial \bar{\rho} \bar{\alpha}_v}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_v}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ D_{\alpha_v} \frac{\partial \bar{\alpha}_v}{\partial x_j} \right] - \dot{m} \left( \frac{1}{\rho_v} - \alpha_v \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right) \right), \quad (23) \\ & \frac{\partial \bar{\rho} \bar{\alpha}_{nc}}{\partial t} - \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_{nc}}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ D_{\alpha_{nc}} \frac{\partial \bar{\alpha}_{nc}}{\partial x_j} \right], \\ & \frac{\partial \bar{\rho} \bar{\Sigma}'}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\Sigma}'}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ D_{\Sigma'} \frac{\partial \bar{\Sigma}'}{\partial x_j} \right] + C_{\Sigma} \frac{\Sigma}{\tau_{\Sigma}} \left( 1 - \frac{\Sigma}{\Sigma_{eq}} \right), \quad (24) \end{aligned} $$

where $D_\alpha$ and $D_{\Sigma'}$ are the sub-grid turbulent diffusivity coefficients and are considered proportional to the sub-grid turbulent viscosity, $D_\alpha = D_{\Sigma'} = D_{\text{sgs}} = \nu_{\text{sgs}}/S_{\text{sgs}}$, where $S_{\text{sgs}}$ is the turbulent Schmidt number which is considered to be 0.9.¹⁸

## D. Interface capturing method (ICM)

Inspection of Eq. (23) highlights the diffusive characteristic of this model. The system is closed using a diffusive turbulent term, and an interface compression term seen in Eq. (1) is absent. Following the method of Anez et al.,²⁹ the ICM approach re-introduces this interface compressive term in the volume fraction $\alpha$ equations

$$ \begin{aligned} & \frac{\partial \bar{\rho} \bar{\alpha}_l}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_l}{\partial x_j} + \underbrace{\frac{\partial}{\partial x_j} C_{\alpha} (\alpha_l \alpha_v (U_l - U_v) + \alpha_l \alpha_{nc} (U_l - U_{nc}))}_{\text{ICM}} \\ & = \underbrace{\frac{\partial}{\partial x_j} \left[ (1 - C_{\alpha}) D_{\alpha} \frac{\partial \bar{\alpha}_l}{\partial x_j} \right]}_{\text{ELSA}} + \mathbf{S}_{\alpha}, \\ & \frac{\partial \bar{\rho} \bar{\alpha}_v}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_v}{\partial x_j} + \underbrace{\frac{\partial}{\partial x_j} C_{\alpha} (\alpha_v \alpha_l (U_v - U_l))}_{\text{ICM}} \\ & = \underbrace{\frac{\partial}{\partial x_j} \left[ (1 - C_{\alpha}) D_{\alpha} \frac{\partial \bar{\alpha}_v}{\partial x_j} \right]}_{\text{ELSA}} - \mathbf{S}_{\alpha}, \\ & \frac{\partial \bar{\rho} \bar{\alpha}_{nc}}{\partial t} + \frac{\partial \bar{\rho} \bar{u}_j \bar{\alpha}_{nc}}{\partial x_j} + \underbrace{\frac{\partial}{\partial x_j} C_{\alpha} (\alpha_{nc} \alpha_l (U_{nc} - U_l))}_{\text{ICM}} \\ & = \underbrace{\frac{\partial}{\partial x_j} \left[ (1 - C_{\alpha}) D_{\alpha} \frac{\partial \bar{\alpha}_{nc}}{\partial x_j} \right]}_{\text{ELSA}}, \end{aligned} \quad (25) $$
where C_α is the interface compression coefficient. Conventionally, this term determines whether the volume fraction interface is considered sharp, but here, it also determines where to apply the gradient law for the LES closure of the ELSA model. In this case, this coefficient is considered a dynamic binary switch. When set to zero, the surface is considered to be diffusive, and when set to one, the surface is considered to be sharp. The setting of this switching is determined but two criteria for IRQ. Once a threshold has been reached with either IRQ, then the surface is treated as either sharp or diffusive.

$$IRQ_{\Sigma} = \frac{\Sigma_{min}}{\Sigma}, \qquad (26)$$

$$IRQ_{\kappa} = \frac{1}{\Delta \cdot \kappa}. \qquad (27)$$

The first criterion IRQ_Σ is based on the local surface density. It is a ratio of the minimum surface density and the local surface density. In local regions of high surface density fluctuations, then this value will be low. Here, an adopted threshold of IRQ_Σ < 0.5 for C_α = 0 is adopted, with C_α = 1 when IRQ_Σ ≥ 0.5. The second criterion IRQ_κ is based on the cubic root of the cell volume and the liquid surface curvature κ. The more curvature there is locally, lower IRQ_κ is with a threshold of IRQ_κ < 1 for C_α = 0 was adopted here. If either one of these criteria is reached, then the setting for dilute region is considered. The interface compression term then switches as a binary step function between 0 and 1. The criteria in which the dilute and sharp surfaces have been clearly outlined which show that with little curvature and surface density locally then the SGS fluctuations are not considered and the interface is artificially compressed. This is advantageous when attempting to model wall attached cavitation as this allows for the sharp leading interface to be resolved. The nature in which the compressive coefficient C_α switches may become problematic in far field atomization as this could lead to an over-prediction of sharp surfaces. However, it is argued by such flow regions; then, more appropriate Lagrangian models can be used to model far-field dispersed sprays.²⁹

### E. Discretization

The solution procedure for the surface density transport is split into the decomposition of Eq. (15). First, the minimum surface density Σ_min is determined with the use of Eq. (14). Then, the additional transported surface density Σ' is solved using the following discretized form of Eq. (12)

$$\frac{\Sigma'^{,n+1} - \Sigma'^{,n}}{\Delta t} V_P + \sum_f F_f = \sum_f D_{sgs,f} S_f \cdot (\nabla \Sigma'^{,n+1})_f. \qquad (28)$$

The source term is treated separately with an explicit and implicit term. The source term linearization ensures that the system is diabolically dominant

$$S_{\Sigma} = C_{\Sigma} \frac{\Sigma}{\tau_{\Sigma}} \left(1 - \frac{\Sigma}{\Sigma_{eq}}\right) = \underbrace{\frac{C_{\tau}}{\tau_{sgs}} \frac{\Sigma}{\sigma_{eq}} (\Sigma_{eq} - \Sigma_{min})}_{\text{Explicit}} - \underbrace{C_{\tau} \tau_{sgs} \frac{\Sigma}{\sigma_{eq}} \Sigma'}_{\text{Implicit}}. \qquad (29)$$

### III. LIQUID JET

To demonstrate the model performance in terms of predicting the underlying liquid breakup processes occurring in injection systems, a DNS test case developed by Ménard et al.³² is presented first. This setup represents a liquid jet formation from a generic Diesel fuel injector that resembles typical experimental injection systems. The operating conditions of this system are not fully representative of realistic industrial injectors; however, this case can still be considered representative of high Weber and Reynolds number injectors where the jet momentum drives the liquid breakup.¹⁷ The turbulent diffusion closure handles the liquid dispersion well for this type of atomization. It should be pointed out that, for this case, only a brief discussion of the physics is included since the case is used mostly for numerical validation focusing on the comparison of our data with the available DNS data from previous studies to analyze the numerical sensitivity of the different model parameters.

The mesh adopted consists of a domain of 25 × 6 × 6 d_j, where d_j is the initial jet diameter and the center of this diameter is the point of origin of the domain, see Fig. 2. Three mesh resolutions were considered for this case: 320 × 40 × 40 (M_j), 240 × 30 × 30 (M_m), and 160 × 120 × 120 (M_c) cells. To induce turbulence directly at the inlet, a digital filter generator proposed by Klein et al.³³ was adopted. This boundary condition considers an artificial mesh that is beyond the domain which is mapped at the inlet.³³ A turbulence intensity of 5% and a length scale of L ≈ 10μm are used for this case to define the turbulent fluctuations within the artificial domain. The selection of these boundary conditions was based on previous studies.³⁴ Alongside variation of the mesh resolution, the critical Weber number We* was varied using the values of We* = 1, 3, and 6. The WALE turbulence model was adopted in this case.³⁵ The chosen probe locations are based on the locations used in previous studies at x/d_j = 5, 10, and 20 perpendicular to the center of the nozzle (see Fig. 2). The first two points represent the transition from the dense core to dilute flow. At the final probe, surface breakup is caused by droplet collision as well as the turbulent viscous effects.

Figure 3 is a visual comparison of the fully formed liquid jet surface between previous DNS simulation¹⁷ and VOF-ELSA model. The DNS model employed a coupled VOF and level-set method where the presented interface is defined by the level set function. The interface presented on the ELSA model is an iso-contour where α_l = 0.1 that

<figure><img src="image_14.png" /><figcaption>FIG. 2. Examples of (top) the computational domain and (bottom) the location of the probe locations.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 3. Comparison of the DNS level-set interface with the LES-$\Sigma$ transport interface defined as $\alpha_l = 0.1$. With the LES interface, a center slice of a color distribution of the surface density $\Sigma$ has been superimposed on the bottom half to show the dispersion of liquid structures in the dilute region.</figcaption></figure>

represents the threshold in which the flow is considered dilute. The LES filtering of the volume fraction transport produces an interface from the cell average of the phase indicator $\alpha_l$. The introduction of the ELSA model complements the volume fraction distribution by tracking the various sub-grid surface interactions with the dispersed droplets in the dilute phase observed by the surface density. The average surface density across the center of the domain can be seen in Fig. 4 with the iso-contours of $\alpha_l = 0.5$ (blue) and $\alpha_l = 0.1$ (green). Higher concentration of surface density can be observed in the transition between the dense liquid core and the dilute region. The surface density extends into the dilute region following the angle of the spray. The spray angle was measured to be approximately $18^\circ - 21^\circ$ based on the iso-contour of $\langle\alpha_l\rangle = 0.01$ which is in agreement with the Reitz and Bracco correlation value of $19.2^{\circ}$.³⁶ No liquid surface is observed deep within the liquid core (as expected) and starts being generated close to

<figure><img src="image_16.png" /><figcaption>FIG. 4. A color distribution of the mean surface density $\Sigma$. A blue line for $\alpha_l = 0.5$ and a green line for $\alpha_l = 0.1$ has been superimposed.</figcaption></figure>

$\alpha_l = 0.5$. The tracked surface density allows for quantitative extraction of comparable experimental quantities such as the SMD. That said, even though the surface density can give us a good indication of the SGS liquid structures, explicit information on the interface is not available such as ligament detachment from the main liquid core. However, given the computational expense of DNS simulations, the ELSA model still provides a good approximation of the liquid structures and quantitative analysis.

We start our analysis with the comparison of the liquid volume fraction profiles with the DNS data included for comparison. The distributions are sampled from the probe locations indicated in Fig. 2. The time averaged volume fraction profiles taken from these locations are presented in Fig. 5. The initial axial decay of the volume fraction represents the boundary of the dense liquid core with the LES profiles showing a shorten core compared to the DNS results regardless of the mesh resolution. The fastest decay can be seen with the coarse mesh, which is an indication that the early breakup can be attributed potentially to grid resolution. The volume fraction profiles show good agreement with DNS results further away from the nozzle. Similar results were presented in previous LES studies.²⁸,³¹

The probed distributions of the liquid surface density are presented in Fig. 6 with DNS also included for comparison. A critical Weber number of $We^* = 1$ was adopted for comparison of the various mesh resolutions. The differences between mean surface densities

<figure><img src="image_17.png" /><figcaption>FIG. 5. A comparison of the mesh resolutions of the liquid volume fraction profiles at the axial and radial probe locations of $x/d_j = 5$, $x/d_j = 1$, $x/d_j$. $We^* = 1$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 6. A comparison of the mesh resolutions of the liquid surface density profiles at the axial and radial probe locations of $x/d_j = 5$, $x/d_j = 1$, $x/d_j$, $We^* = 1$.</figcaption></figure>

from the variation of mesh refinement are consistent with the volume fraction distributions outlined previously. The early breakup of the dense liquid core can be observed at the axial profile, with the finer mesh producing slightly improved results. The early generation of surface density $\Sigma$ from our calculations compared to the DNS profile is also consistent with observations made with the mean volume fraction. The diffusive characteristics of the gradient-law LES closure hinder the models ability to produce a sharp profile especially around the liquid core. Here, second-order schemes were adopted for the flux convection, however, using a higher-order scheme such as the WENO scheme could produce a better profile for the liquid core.³⁷ Nonetheless, even with the early breakup of the liquid core, the maximum generated surface from all cases shows good agreement with experiment. The axial profile, after the peak surface density, follows the same decay as the

DNS data with differences observed between the mesh resolutions. The surface distribution radially, further away from the center, shows good agreement with a peak $\langle\Sigma\rangle$ found at $r/d_j = 0.5$. The profiles observed at $x/d_j = 10$ and $x/d_j = 20$ agree reasonably well with DNS with a slight over-estimation in the generated surface for all cases. These results overall highlight that although the model can predict reasonably well the surface generation and destruction mechanisms the models diffusive nature creates discrepancies.

Variation of the critical Weber number $We^*$ profiles can be found in Fig. 7, performed on the fine mesh. All cases follow a similar profile with surface generating and collapsing at similar locations; however, there are differences in the overall surface generated. Observation of all distributions shows that increasing $We^*$ leads to less surface generated. This provides for a better radial profile at $x/d_j = 5$ although here the

<figure><img src="image_19.png" /><figcaption>FIG. 7. A comparison of different critical Weber $We^*$ numbers of the liquid surface density profiles at the axial and radial probe locations of $x/d_j = 5$, $x/d_j = 1$, $x/d_j$. Cases were performed using the fine mesh $M_f$.</figcaption></figure>
dense liquid core is present in the DNS case (zero surface at $r/d_j = 0$) and has defused by this location for the LES results. All radial profiles demonstrate a convergence of the surface density away from the center showing that $We^*$ has little influence in this region. Sensitivities are shown to be isolated to the amount of surface density $\Sigma$ generated with a reasonable profile obtained with $We^* = 1.5$. However, the maximum amount of surface density was under-estimated in this case, thus an ideal value for $We^*$ would be in the range $1 < We^* < 1.5$. A value of $We^* = 1$ is adopted here throughout as the differences between the two cases are not significant where the ideal critical Weber number $We^*$ varies on a case-by-case basis.

#### IV. CAVITATING JET

The test case detailed above focused on investigating the numerical sensitivity of the surface density parameters for predicting external spray formation. To investigate the surface evolution process occurring at cavitating sprays, a simplified nozzle geometry of a high pressure injection system is adopted. This case is a scaled representation of a generic fuel (Diesel) injector present in IC engines featuring an enlarged orthogonal geometry. The operating conditions result in $Re$ numbers that are typical in a range of common industrial systems, and thus, the selected case is considered a suitable case to identify the effect that cavitation has on spray atomization. The WALE turbulence model was adopted in this case.³⁵ Flow contour visualization and some quantitative analysis of both the jet formation and cavitation dynamics are available from experiments.³ In the experimental setup, transparent material is used on the front and back of the injector which allows for visualization of the cavitation structures and their fluid surfaces, but it does not provide quantitative information about the surfaces within the nozzle. More experimental data are available for the jet formation

region. The cases that are explored correspond to a turbulent spray (without cavitation) and a cavitation spray. For these two cases, SMDs from the experiments are available. The droplet SMD measurements were performed using PDA method at set points around the spray topology. These measurements are taken at a specific period of time and then averaged over all appropriate points.

##### A. Setup

Figure 8 shows the locations in which experimental data are provided and droplets are sampled. The sampling locations for the numerical simulations have been extended further closer toward the exit of the nozzle for comparison with cavitating conditions. Thermophysical properties can be found in Table 1. The cases considered here maintain outlet pressure at atmospheric pressure with variation of the inlet pressure. For the turbulent region, which corresponds to flow conditions in which cavitation does not occur, the pressure difference is $P_{inj} = 0.16$ MPa, and for the cavitating region, corresponding to flow conditions where cavitation is present, the pressure difference is $P_{inj} = 0.3$.

##### B. Results and discussion

The top of Fig. 9 shows a comparison of the radial overall SMD profile at 70 mm away from nozzle exit between numerical and experimental data.³ The SMDs are in reasonably good agreement with experiment; however, the experimental profile gives droplet with almost constant size ($d_{32} \approx 49.6$ µm). Lower SMDs are observed from the simulations except for approximately 2 mm away from the center of the jet. This is largely consistent with the liquid jet test case presented in Secs. III and IV which demonstrated a slight over estimation of the

<figure><img src="image_20.png" /><figcaption>FIG. 8. (Top) The grid domain and (bottom) the probe locations in which the SMDs are sampled. The green line denotes sampled locations that are compared to experiment.</figcaption></figure>
<table><caption>TABLE I. Thermodynamic properties of diesel.³⁸</caption><tr><td>Vapor pressure (MPa)</td><td>0.00128</td></tr><tr><td>Liquid density (kg/m³)</td><td>830</td></tr><tr><td>Vapor density (kg/m³)</td><td>0.1361</td></tr><tr><td>Liquid viscosity (Ns/m²)</td><td>0.00223</td></tr><tr><td>Liquid-vapor surface tension</td><td>0.278</td></tr><tr><td>Liquid-air surface tension</td><td>0.0261</td></tr></table>

surface generation, and thus a tendency of the model to predict smaller droplets. This can be overcome by increasing the critical Weber number We* which, in turn, can reduce the over-estimation of surface generation. The location of 2 mm away from the center is roughly in the axial direction of the nozzle wall, and subsequently the side of the jet surface. Here, droplet collision and breakup are under-predicted leaving larger droplets. This could be attributed to the ICM criteria under detecting dilute regions due to an under-prediction of interface curvature upstream after nozzle exit. The bottom of Fig. 9 shows a histogram of the droplet distribution sampled radially 70 mm away from nozzle exit. Mean overall SMDs for simulation and experiments are overlaid. A peak is observed for droplets of approximately $d_{32} = 10$ µm though a mean overall SMD was observed at $d_{32} = 43.69$ µm due to the presence of larger droplets. The surface evolution that resulted in these far field droplet diameters was initiated within the nozzle from the onset of cavitation.

<figure><img src="image_21.png" /><figcaption>FIG. 9. Comparison between the overall SMDs reported in literature and the turbulent region. (Top) SMD radial profiles 70 mm away from nozzle exit. (Bottom) Histogram of the droplet distribution of the turbulent region sampled radially across 70 mm away from nozzle exit. The red black lines denote the measured overall SMDs for simulation and experiment, respectively.⁵</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIG. 10. Comparisons between experiment, three-phase cavitating ELSA-ICM model and equivalent three-phase cavitating VOF model of (top) the nozzle exit velocities, (middle) the respective cavitation numbers, and (bottom) the spray cone angles. The dashed line indicates the threshold in which the flow is considered to be in the cavitating region.</figcaption></figure>

Further experimental comparisons can be found in Fig. 10 that is used to investigate the wider impact these surfaces have to the spray. Results from an equivalent three-phase cavitating VOF solver without any surface density coupling and surface diffusion are included. The threshold in which the flow is considered to be in the cavitating region is indicated with a dashed line. Comparisons of the nozzle exit velocities at the top of Fig. 10 show good agreement with experiment, where velocities increase with an increase in injection pressure. The velocity values produced by the ELSA-ICM model were slightly closer to experimental results than the generic VOF model especially for the cavitating region. Subsequently, this leads to cavitation numbers that are slightly more representative to experiment as seen in the middle of Fig. 10. The potential for the flow to cavitate is indicated by cavitation numbers under 1, the turbulent and cavitation regions produced cavitation numbers of approximately 1.59 and 0.67, respectively. The ELSA-ICM solver can model the surface generation and destruction above grid and at the SGS where the flow is less impacted by the internal viscous effects from the cavitation.
The bottom of Fig. 10 shows a comparison of the measured spray cone angles. The ELSA-ICM model replicated the spray angle well. The generic VOF model generally under-predicted the spray cone angles. The surface dynamics at the SGS are transported with the ELSA-ICM framework which further impacts the subsequent spray formation. This is neglected within the generic VOF solver. The generation and destruction of cavitation structures cause the spray angle to widen, the turbulent and cavitating regions produced angles of approximately 16° and 29°, respectively. To get a better understanding of how the surface dynamics of cavitation impacts, the spray formation not only by widening the spray angle but also by reducing droplet diameters visual inspection of the flow field is required.

Figure 11 presents the surface density contour distribution across the center plane of the domain with experimental visualizations of the jet for comparison.³ A white iso-contour for $\alpha_l = 0.5$ and a green contour for $\alpha_l = 0.1$ have been superimposed to show the location of the liquid jet. The occurrence of cavitation causes a wider spray cone angle (in comparison to the turbulent jet) as seen in the bottom of Fig. 10. Observations of the surface density $\Sigma$ near the nozzle show that for the turbulent region the surface is concentrated in the region of the volume fraction interface. Significant momentum-driven breakup has yet to occur as well as sufficient turbulent effects to deform the surface. For the cavitating case, generated surface is present in the nozzle due to the generated wall attached cavitation. The cavitation break-off causes the deformation of the surface between the liquid and the vapor which, in turn, forms additional vapor structures that are transported downstream to the nozzle exit. By the exit of the nozzle, greater surface density is present due to the bubble breakup energy. This greater surface density at the gravitating case is maintained throughout the dilute

phase of the flow. We observe that the presence of cavitation allows for a longer dense liquid core of the jet. This can be attributed to the characteristic of wall attached cavitation which exists largely within the boundary layer of the nozzle. The dynamics of cloud break-off and collapse is contained within this region where by nozzle exit it is this region that interacts with the turbulent shear effects with the atmosphere. Without the convection of cavitation structures within the boundary layer, the turbulent shear stresses along the surface of the turbulent region are transferred to the center of the jet more rapidly causing the liquid core to diffuse earlier.

A typical characteristic of wall attached cavitation as in the case presented here is the presence of large vapor structures that break off and collapse downstream.³⁹ The surface energy of these structures is conserved and affects the subsequent atomization process.³ It is important, therefore, that the interfaces of these large structures to be tracked numerically and accurately. The advantage of adopting the ICM with the ELSA model is the ability to appropriately maintain a sharp interface for large structures. This is achieved by imposing artificial interface compression to the volume fraction transport, through a binary coefficient $C_\alpha$. Figure 12 shows a color distribution of the interface compression coefficient $C_\alpha$ alongside a color distribution of the surface density $\Sigma$ with iso-contours of the liquid volume interface at $\alpha_l = 0.5$ and $\alpha_l = 0.1$. For both cases, the majority of the spray after nozzle exit is treated by the ELSA framework with interface compression inactive. This shows that the viscous turbulent effects are promoting surface deformation of the main liquid jet. Interface compression is active in areas outside of the liquid core where the droplet distribution is dispersed. Here, the liquid volume fraction is very dilute giving a smooth interface definition. Furthermore, the dilute liquid volume fraction and

<figure><img src="image_23.png" /><figcaption>FIG. 11. Visual comparison between the turbulent and cavitating region. The color distribution shows the surface density $\Sigma$ across the center of the domain. Contours for $\alpha_l = 0.5$ (white) are included. The black and white pictures are experimental photographs.³</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**FIG. 12.** Combined color distributions of the interface compression coefficient $C_\alpha$ and surface density $\Sigma$ with the turbulent region top and the cavitation region bottom. Contours for $\alpha_l = 0.5$ (white) and $\alpha_l = 0.1$ (green) are included.</figcaption></figure>

smooth interface curvature lead to a low $\Sigma_{eq}$ which limits the surface generation. Thus, in these areas, the minimum surface density $\Sigma_{min}$ because of the presence of the liquid dominates the overall definition for surface density $\Sigma$. Cavitation generates liquid surfaces within the nozzle which is absent for the turbulent region. Interface compression is active almost throughout the cavitation structures. This demonstrates that the cavitation structures within the nozzle are larger than the grid scale with a sharp interface. Further downstream can be seen that interface compression becomes slightly more inactive as the cavitation structures evolve more into bubbly mixtures. The sharp interfaces of the large cavity structures travel to nozzle exit and interacts with the initial interface of the liquid jet. This generates greater interface curvature of the main jet and widens the distribution of surface density along the interface.

The influence of cavitation on the liquid surface dynamics at the nozzle exit and the widening of the spray cone angle is also investigated here. Figure 13 shows the evolution of the in-nozzle flow as cavitation develops over time. By $t=1.5$ ms, cavitation structures have formed, and by $t=2$ ms, wall attached cavities have formed with some cloud break-off already occurring. The liquid interface after nozzle exit is largely laminar with little surface generated beyond the volume fraction interface. Some of the cavity structures have fully collapsed without any bubbly structures yet reaching nozzle exit by $t=3$ ms. The greater surface deformation at nozzle exit can be observed due to the energy that is conserved from cavity collapse, and the spray angle is increased. A consistent reentrant jet along the nozzle wall underneath the cavitation is present. Once $t=3.5$ ms is reached, small cavity structures and bubbly mixtures exit the nozzle which increases the amount of surface density $\Sigma$ further away from the liquid volume fraction interface. These generated surfaces are then transported downstream contributing to further droplet break-off and collision.

Greater surface density in the spray formation will lead to a decrease in mean diameters, and this was reported by Suh et al.³ Figure 14 is a comparison of the overall SMD values downstream from the nozzle between the turbulent and cavitating regions. Left is the axial profile. Initially close to the nozzle, it can be seen that the cavitating region consists of larger droplets. Thereafter however, smaller SMDs are observed in the cavitating region with this region showing smaller SMDs overall. The mean difference in SMD between the two regions was 2.48 µm which is relatively larger but in reasonable agreement with experiment with a reported difference of approximately 1 µm.³ The right of Fig. 14 is the radial distribution with the profile sampled 30 mm away from the nozzle exit. Both profiles follow the same as Fig. 9 where closer to the center of the spray a higher overall SMD is observed and then reduces to a lower consistent SMD further away from the center. There is less deviation of SMD for the cavitating region. The turbulent region fluctuates slightly moving away from the center before reducing relatively rapidly further away from the center. Without cavitation, the atomization/break-off for the turbulent region is driven by the momentum of the liquid and the turbulent shear

<figure><img src="image_25.png" /><figcaption>**FIG. 13.** The evolution of the internal nozzle cavitation and the near nozzle surface density over time. The color distribution shows the surface density $\Sigma$ across the center of the domain. Contours for $\alpha_l = 0.5$ (white) and $\alpha_l = 0.1$ (green) are included.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>FIG. 14. Comparison of SMDs between the turbulent and cavitating regions. (Left) Axial profile and (right) is the radial profile sampled 30 mm away from nozzle exit.</figcaption></figure>

stresses. Subsequently larger droplets are formed toward the center of the jet with little other influence other than the turbulent interaction with the atmosphere. Small droplets are seen further away from the center of the jet as small ligaments break off from the surface from

shear effects. The presence of cavitation introduces bubble breakup energy that is conserved and affects the surface breakup downstream of nozzle exit. This and with small cavity structures leaving the nozzle increases the amount of surface at nozzle exit, see Fig. 11, thus the

<figure><img src="image_27.png" /><figcaption>FIG. 15. Histograms showing the evolution of the droplet diameter distribution over the distance away from the nozzle exit. The left column shows the distributions for the cavitating region, and the right presents the turbulent region.</figcaption></figure>
ICM criteria considers this region as consistently dilute. Furthermore, wall attached cavitation, which is observed here, is inherently a boundary layer phenomenon which reducing the width of the subsequent liquid core. Though close to nozzle exit, the extent of the liquid core is slightly longer due to the higher exit velocity, further downstream the stored bubble breakup energy diverges throughout the liquid core which enhances the atomization with greater droplet collision events.

The overall SMDs show that cavitation affects the atomization of the liquid producing smaller droplets on average than flows without cavitation present. In order to understand better how the atomization has been affected an inspection of the droplet diameter distributions is needed. Figure 15 is array of histograms showing the droplet distributions are the sampled axial locations corresponding to the left of Fig. 14. Initially, both regions have similar distribution of droplet diameters even so the cavitating region contains some larger droplets that are absent from the turbulent region. Further downstream from nozzle exits; however, larger droplets develop progressively for the turbulent region. This is not the case for the cavitating region where the maximum SMD grows slightly before reducing gradually further downstream. A common observation of the distributions is that all cases have reasonably equitable peak SMDs at around $d_{32} \approx 10 \ \mu m$; thus, the main contributor to the differences in overall SMD between the two regions is the presence of larger droplets for the turbulent region. These larger droplets do not persist with cavitation as the greater droplet collision break up these structures.

## V. CONCLUSIONS

The introduction of an SGS surface density quantity which is based on the ELSA-ICM approach into a three phase cavitating solver was made. This allowed for greater insight into the internal and near-nozzle dynamics of high pressure spray injectors. The transport of the generated surface was validated against DNS data using a range of grid resolutions. ELSA-ICM simulations were performed using a generic orthogonal nozzle with results compared with experiment. The main findings from this chapter can be summarized as follows:

* Validation showed that the surface density is slightly over predicted. This is due to the under-prediction of the dense liquid core where this structure becomes diffusive more rapidly compared to DNS results. The grid resolution variation demonstrated that this was largely grid dependent, and this diffusive trend becomes less pronounced.

* Increasing the critical Weber number $We^*$ reduced the amount of surface density generated. It was shown that an optimum value is somewhere in the region of $1 < We^* < 1.5$ with a value of $We^* = 1$ was used throughout the subsequent investigations.

* The measured spray cone angles agree well with experiment. The presence of cavitation caused a wider spray cone angle on the jet formation due to the greater curvature and surface density at nozzle exit.

* Smaller droplet size distributions were observed with cavitation considered; this is consistent throughout. This is due to the surface energy of collapsed cavity structures being conserved and transported to nozzle exit, increasing local interface curvature and surface density. The change in overall droplet size in the presence of cavitation compared to systems without cavitation was shown to be slightly larger than experiment. More work is needed to understand the numerical sensitivity of surface density

parameters such as the critical Weber number $We^*$ on cavitating systems.

* The integration of the ICM with the ELSA approach allows for the liquid interface to be modeled as either dilute or sharp suitably depending on local surface conditions. The model was shown to largely consider the cavitation surface structures within the nozzle as sharp interfaces, with the jet formation surfaces considered as dilute. This shows that the cavitation interface is modeled appropriately as wall attached cavitation typically consistent of relatively large cavity structures.

## ACKNOWLEDGMENTS

The authors would like to acknowledge funding by the UK Engineering and Physical Science Research Council support through the Grant Nos. EP/N509607/1 and EP/S001824/1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Paul McGinn:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal). **Giovanni Tretola:** Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Konstantina Vogiatzaki:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (lead); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

¹C. Crua and M. R. Heikal, "Time-resolved fuel injector flow characterisation based on 3D laser doppler vibrometry," *Meas. Sci. Technol.* **25**, 125301 (2014).

²B. A. Sforzo, K. E. Matusik, C. F. Powell, A. L. Kastengren, S. Daly, S. Skeen, E. Cenker, L. M. Pickett, C. Crua, and J. Manin, "Fuel nozzle geometry effects on cavitation and spray behavior at diesel engine conditions," in *Proceedings of the 10th International Symposium on Cavitation (CAV2018)* (ASME, 2018).

³H. K. Suh and C. S. Lee, "Effect of cavitation in nozzle orifice on the diesel fuel atomization characteristics," *Int. J. Heat Fluid Flow* **29**, 1001–1009 (2008).

⁴A. Sou, S. Hosokawa, and A. Tomiyama, "Effects of cavitation in a nozzle on liquid jet atomization," *Int. J. Heat Mass Transfer* **50**, 3575–3582 (2007).

⁵A. Ahmed, "LES of atomization and cavitation for fuel injectors," Ph.D. thesis (Université de Rouen Normandie, 2019).

⁶P. McGinn, D. Pearce, Y. Hardalupas, A. M. K. P. Taylor, and K. Vogiatzaki, "Cavitation bubble cloud break-off mechanisms at micro-channels," *Fluids* **6**, 215 (2021).

⁷S. Jahangir, W. Hogendoorn, and C. Poelma, "Dynamics of partial cavitation in an axisymmetric converging-diverging nozzle," *Int. J. Multiphase Flow* **106**, 34–45 (2018).
<footnote n="8">⁸A. Peters, H. Sagar, U. Lantermann, and O. el Moctar, "Numerical modelling and prediction of cavitation erosion," *Wear* **338–339**, 189–201 (2015).</footnote>

<footnote n="9">⁹M. Ylönen, *Cavitation Erosion Monitoring by Acoustic Emission*, Tampere University Dissertations – Tampereen yliopiston väitöskirjat (Tampere University, 2020).</footnote>

<footnote n="10">¹⁰A. Sou, B. Bicer, and A. Tomiyama, "Numerical simulation of incipient cavitation flow in a nozzle of fuel injector," *Comput. Fluids* **103**, 42–48 (2014).</footnote>

<footnote n="11">¹¹F. Örley, T. Trummler, S. Hickel, M. S. Mihatsch, S. J. Schmidt, and N. A. Adams, "Large-eddy simulation of cavitating nozzle flow and primary jet break-up," *Phys. Fluids* **27**, 086101 (2015).</footnote>

<footnote n="12">¹²M. G. Mithun, P. Koukouvinis, and M. Gavaises, "Numerical simulation of cavi- tation and atomization using a fully compressible three-phase model," *Phys. Rev. Fluids* **3**, 064304 (2018).</footnote>

<footnote n="13">¹³F. P. Kärrholm, H. Weller, and N. Nordin, "Modelling injector flow including cavitation effects for diesel applications," in 5th Joint ASME/JSME Fluids Engineering Conference, 2007.</footnote>

<footnote n="14">¹⁴C. Yu, Y. Wang, C. Huang, T. Du, C. Xu, and J. Huang, "Experimental and numerical investigation on cloud cavitating flow around an axisymmetric pro- jectile near the wall with emphasis on the analysis of local cavity shedding," *Ocean Eng.* **140**, 377–387 (2017).</footnote>

<footnote n="15">¹⁵F. Giussani, F. Piscaglia, G. Saez-Mischlich, and J. Hèlie, "A three-phase VOF solver for the simulation of in-nozzle cavitation effects on liquid atomization," *J. Comput. Phys.* **406**, 109068 (2020).</footnote>

<footnote n="16">¹⁶A. Vallet, A. A. Burluka, and R. Borghi, "Development of an Eulerian model for the atomization of a liquid jet," *Atomization Spray* **11**, 24 (2001).</footnote>

<footnote n="17">¹⁷R. Lebas, T. Menard, P. Beau, A. Berlemont, and F. Demoulin, "Numerical sim- ulation of primary break-up and atomization: DNS and modelling study," *Int. J. Multiphase Flow* **35**, 247–260 (2009).</footnote>

<footnote n="18">¹⁸G. Tretola, K. Vogiatzaki, and S. Navarro-Martinez, "Sub-grid effects in atomisation process using stochastic fields," in *ICLASS 2018*, 14th Triennial International Conference on Liquid Atomization and Spray Systems, 2018.</footnote>

<footnote n="19">¹⁹K. G. Lyras, S. Dembele, and J. X. Wen, "Numerical simulation of flashing jets atomisation using a unified approach," *Int. J. Multiphase Flow* **113**, 45–58 (2019).</footnote>

<footnote n="20">²⁰P. McGinn, G. Tretola, and K. Vogiatzaki, "Unified modeling of cavitating sprays using a three-component volume of fluid method accounting for phase change and phase miscibility," *Phys. Fluids* **34**, 082108 (2022).</footnote>

<footnote n="21">²¹M. S. Plesset and R. B. Chapman, "Collapse of an initially spherical vapour cav- ity in the neighbourhood of a solid boundary," *J. Fluid Mech.* **47**, 283–290 (1971).</footnote>

<footnote n="22">²²P. Zwart, A. G. Gerber, and T. Belamri, "A two-phase flow model for predicting cavitation dynamics," in Fifth International Conference on Multiphase Flow, 2004.</footnote>

<footnote n="23">²³J. Sauer and G. H. Schnerr, "Unsteady cavitating flow—a new cavitation method based on a modified front capturing method and bubble dynamics," in 2000 ASME Fluids Engineering Summer Conference, 2000.</footnote>

<footnote n="24">²⁴C. Morel, "On the surface equations in two-phase flows and reacting single-phase flows," *Int. J. Multiphase Flow* **33**, 1045–1073 (2007).</footnote>

<footnote n="25">²⁵D. Lhuillier, C. Morel, and J.-M. Delhaye, "Bilan d'aire interfaciale dans un mélange diphasique: Approche locale vs approche particulaire," *C. R. Acad. Sci., Ser. IIB - Mech.-Phys., Astron.* **328**, 143–149 (2000).</footnote>

<footnote n="26">²⁶J.-M. Delhaye, "Some issues related to the modeling of interfacial areas in gas-liquid flows. I. The conceptual issues," *C. R. Acad. Sci. - Ser. IIB - Mech.* **329**, 397–410 (2001).</footnote>

<footnote n="27">²⁷O. Séro-Guillaume and N. Rimbert, "On thermodynamic closures for two-phase flow with interfacial area concentration transport equation," *Int. J. Multiphase Flow* **31**, 897–920 (2005).</footnote>

<footnote n="28">²⁸S. Navarro-Martinez, "Large eddy simulation of spray atomization with a prob- ability density function method," *Int. J. Multiphase Flow* **63**, 11–22 (2014).</footnote>

<footnote n="29">²⁹J. Anez, A. Ahmed, N. Hecht, B. Duret, J. Reveillon, and F. Demoulin, "Eulerian–Lagrangian spray atomization model coupled with interface capturing method for diesel injectors," *Int. J. Multiphase Flow* **113**, 325 (2019).</footnote>

<footnote n="30">³⁰S. Navarro-Martinez, G. Tretola, M. R. Yosri, R. L. Gordon, and K. Vogiatzaki, "An investigation on the impact of small-scale models in gasoline direct injection sprays (ECN Spray G)," *Int. J. Engine Res.* **21**, 217–225 (2020).</footnote>

<footnote n="31">³¹J. Chesnel, J. Reveillon, T. Ménard, and F. X. Demoulin, "Large eddy simulation of liquid jet atomization," *Atomization Spray* **21**, 711–736 (2011).</footnote>

<footnote n="32">³²T. Ménard, S. Tanguy, and A. Berlemont, "Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet," *Int. J. Multiphase Flow* **33**, 510–524 (2007).</footnote>

<footnote n="33">³³M. Klein, A. Sadiki, and J. Janicka, "A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations," *J. Comput. Phys.* **186**, 652–665 (2003).</footnote>

<footnote n="34">³⁴G. Tretola, K. Vogiatzaki, and S. Navarro-Martinez, "Implementation of a probabilistic surface density volume of fluid approach for spray atomisation," *Comput. Fluids* **230**, 105121 (2021).</footnote>

<footnote n="35">³⁵F. Ducros, N. Franck, and T. Poinsot, "Wall-adapting local eddy-viscosity mod- els for simulations in complex geometries," in *Numerical Methods for Fluid Dynamics VI* (1998).</footnote>

<footnote n="36">³⁶R. Reitz and F. Bracco, "Mechanism of atomization of a liquid jet," *Phys. Fluids* **25**, 1730–1742 (1982).</footnote>

<footnote n="37">³⁷J. W. Gärtner, A. Kronenburg, and T. Martin, "Efficient WENO library for OpenFOAM," *SoftwareX* **12**, 100611 (2020).</footnote>

<footnote n="38">³⁸European Committee for Standardization, EN590:2009 *Automotive fuels – Diesel – Requirements and Test Methods* (European Committee for Standardization, 2009).</footnote>

<footnote n="39">³⁹J.-P. Franc, "Physics and control of cavitation," in *Design and Analysis of High Speed Pumps* (NATO, 2006), pp. 2-1-2-36.</footnote>