# An Euler-Lagrange coupling strategy for pressurized liquid jet atomization in Large Eddy Simulation

Thibault Gioud, Thomas Laroche, Thomas Schmitt, Bénédicte Cuenot, Odier Nicolas

► To cite this version:

Thibault Gioud, Thomas Laroche, Thomas Schmitt, Bénédicte Cuenot, Odier Nicolas. An Euler-Lagrange coupling strategy for pressurized liquid jet atomization in Large Eddy Simulation. Computers and Fluids, 2025, E&S, 294, pp.106625. 10.1016/j.compfluid.2025.106625. hal-05045837

HAL Id: hal-05045837

https://hal.science/hal-05045837v1

Submitted on 24 Apr 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Euler-Lagrange coupling strategy for pressurized liquid jet atomisation in Large Eddy Simulation

Thibault Gioud<sup>a</sup>, Thomas Laroche<sup>a</sup>, Thomas Schmitt<sup>b</sup>, Bénédicte Cuenot<sup>a</sup>, Odier Nicolas<sup>a</sup>

<sup>a</sup>CFD Team, CERFACS, 31057 Toulouse, France

<sup>b</sup>Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91190, Gif-sur-Yvette, France.

## Abstract

The atomisation process is crucial for numerous industrial and engineering applications, ranging from fuel injection systems to spray coating technologies and remains a challenge for numerical prediction. This study presents a coupled Euler/Lagrange methodology to model the atomisation processes, at affordable numerical costs. An Eulerian framework is used to resolve the liquid core and the largest liquid structures, while a Lagrangian approach tracks the smallest droplets and their interactions, addressing transport, secondary breakup, and evaporation phenomena. The Eulerian two-phase flow method considered here is a multi-fluid diffuse interface approach that assumes equilibrium of temperature, pressure, velocity, and Gibbs potentials at the liquid-gas interface. An originality of this paper is the transfer of liquid from the under resolved Eulerian liquid structures to the Lagrangian formalism when user-defined geometrical criteria are met. A validation of the coupling strategy is performed on a jet in cross flow configuration. Although the spatial distribution of the droplet volumetric flux is not exactly retrieved, other key characteristics of the spray, such as the Sauter Mean Diameter (SMD) and the outer trajectory, show very encouraging results. Finally, this paper shows the importance of taking into account the liquid core in the simulation, rather than directly injecting already atomized droplets.

**Keywords:** Two-phase flows, Multi-fluid method, diffuse interface, Atomisation, Coupling strategy, Euler-Lagrange

Email addresses: gioud@cerfacs.fr (Thibault Gioud), thomas.laroche@safrangroup.com (Thomas Laroche), thomas.schmitt@centralesupelec.fr (Thomas Schmitt), cuenot@cerfacs.fr (Bénédicte Cuenot), odier@cerfacs.fr (Odier Nicolas)
# 1. Introduction

In liquid rocket engines, atomization is a dominant process that drives the flame structure and should be accurately predicted. Liquid core atomization creates a dense cloud of small droplets, maximizing liquid-gas surface exchange and enhancing combustion efficiency [47]. Atomization process involves complex phenomena across a large range of spatio-temporal scales [3]. The liquid breakup is driven by three main phenomena, the capillary forces, the inertial forces and the viscosity force [24, 87]. Moreover, two types of liquid breakup can be distinguished: the breakup of the liquid core, which is referred to as primary breakup, and the further breakup of these atomized structures into smaller droplets, which is referred to as secondary breakup [35, 52]. A clear understanding of these phenomena is yet to be reached, and still requires both experimental and numerical investigations.

Extensive experimental research ([20, 28, 30, 44, 45, 53]) has been conducted on liquid jets or sheets atomization for the understanding of basic involved phenomena, and evaluation of properties of the resulting spray. Most of these studies focus on the liquid core destabilisation and the primary atomization processes. Detailed reviews are provided in [4, 22]. The present study focuses on a Jet In Cross Flow (JICF) configuration that is a widely studied atomization test case ([8, 26, 80, 89]), involving the disintegration of a liquid jet core by a gaseous cross flow. The work of Broumand & Birouk [15] provides a comprehensive review on this subject. These investigations have allowed the proposition of jet trajectory correlations as a function of global jet properties, and the measurement of droplets size distributions. Nevertheless, the mechanisms that control spray formation and evolution are far from being understood, due to the lack of experimental data on local jet properties. Indeed, the speed and small scales at which the processes occur and the opaque cloud created by the smallest droplets prevent from accurate measurements. This lack of knowledge is highly unfortunate, since these atomization processes play a major role in industrial fuel injectors, either in aeronautical engines [47], or in subcritical liquid rocket engines [55].

Computational Fluid Dynamics (CFD) can be a helpful tool to complement experimental measurement, and to help physical understanding of processes involved in these industrial injectors. On one hand, numerical methods have been developed to allow interface capturing, involving accurate numerical schemes to handle complex topologies and leading to significant computational costs due
to the multiscale nature of atomization. These "sharp interface methods" mainly consist of Volume of Fluid (Hirt & Nichols [40]) and Level-Set methods (Sussman et al. [84]), possibly coupled (Sussman & Puckett [83]). These techniques, often relying on incompressible set of equations, have allowed Direct Numerical Simulations (DNS) of jet atomization ([5, 58, 79, 92]), assuming that all the physics is predicted on the mesh, providing insightful physical investigation, at the price of a significant computational cost. In order to alleviate this computational cost, Adaptive Mesh Refinement (AMR) technique have been developed, for both structured ([1, 48, 67]) and unstructured meshes [42]. Nevertheless, the cost of such DNS remains significant, prohibiting their daily use in the industry.

On the other hand, an affordable computational method for numerical spray prediction consists in Lagrangian Particle Tracking (LPT), assuming point particle. Droplets are assumed to be spherical and rigid punctual spheres, not discretized on the mesh, and may account for the subgrid interactions between the liquid droplet and the gaseous phase. Such approaches allow the use of affordable meshes, and have led to computations of industrial flows such as Diesel injection ([34], [2]), aircraft engines injection ([21, 41, 59, 73, 85]), and sub-critical rocket engines injection ([11, 12, 68]). If very interesting results can be obtained for spray prediction in various conditions, at a low computational cost, LPT method unfortunately fails to describe liquid injection. Instead, LPT injects already formed droplets, and therefore requires the statistical knowledge of droplet size distributions, and possibly size-velocity correlations to be injected. Such data can be provided through challenging experimental measurement or DNS, but are in practice rarely known. Moreover, this liquid injection strategy imposes a frozen spray at injection, without interaction with the flow.

In order to better describe liquid injection and avoiding the tuning of the initial Lagrangian spray characteristics, coupling between Eulerian and Lagrangian approaches have been proposed. The liquid core and the largest liquid structures are resolved using an Eulerian two-phase flow method, while the LPT method is used to address the under-resolved droplets and secondary atomization processes. Herrmann [39] first proposed a parallel Eulerian/Lagrangian multi-scale coupling procedure for two-phase flows, using a Level-Set method. Several other coupling algorithms ([38, 49, 90]) have been proposed for the Eulerian-Lagrangian droplet transfer. For example, [38] relied on a Connected Component Labelling (CCL) algorithm to couple the VOF method with LPT model.
The global concept of the Euler-Lagrange coupling procedure involves three key steps: (1) identification of candidate liquid structures for atomization, (2) assessment of structure properties against predefined criteria, and (3) transfer of liquid structures satisfying the criteria from the Eulerian to the Lagrangian formalism. Typically, spatial resolution criteria are established based on geometric considerations to ensure the transfer of under-resolved liquid structures only. However, the criteria for determining atomized structures can vary among authors. For instance, in [27, 39] the transfer criterion is based on the volume of liquid structures, while in [92], the criterion is based on the number of cells in the liquid structures. Additionally, many coupling algorithms use a sphericity criterion to ensure the transfer of only spherical structures. Furthermore, in some cases, a physical criterion such as the Weber number is used, in addition with geometric criteria, as in [27]. A common method for assessing coupling strategies involves simulating Jet in Cross Flow (JICF) configurations, as in [17, 88, 91] among others. To the authors knowledge, these studies allow the prediction of the external jet trajectory, but mostly rely on incompressible frameworks, and no assessment for the spray statistics (SMD) is proposed for these studies.

The objective of this paper is to propose a framework intended for Large Eddy Simulations (LES) of industrial injectors, while allowing affordable computational costs for a direct use within design offices. We furthermore intend our methodology to rely on a compressible set of equations, to comply with a potential pressurized, reactive context. The authors are aware of some very recent works [10, 50, 82] proposing promising results for VOF methods coupled to Adaptive Mesh Refinement. In this work, we instead propose to couple a compressible multi-fluid Diffuse Interface Method (DIM) with a LPT algorithm, to predict both the main liquid core and the subsequent spray. The use of a DIM methods intends to avoid the use of adaptive mesh refinement techniques, which might induce extra CPU costs and challenges in an unstructured context. The multi-fluid Diffuse Interface Method (DIM) relies on a 3-equation transport system recently developed within the work of Pelletier et al. [61], and derived from Saurel & Abgrall [71]. The non-ideal real-gas thermodynamics is handled through the Peng & Robinson [62] (PR) equation of state. This method has been used in realistic reactive liquid rocket engine configurations [32, 33, 75, 76], providing encouraging results at several pressure levels. The proposed algorithm identifies closed liquid structures relying on a threshold value of the liquid volume fraction ($\alpha_{liq}$). Spherical structures are transferred into a single Lagrangian particle, while elongated structures are transferred into multiple Lagrangian
particles. Allowing the transfer of elongated structures facilitates the transfer of larger Eulerian liquid structures, which tend to be located in closer proximity to the injector, thereby reducing the overall mesh load. Finally, a secondary atomization model relying on the stochastic methodology proposed in [2, 34] transforms the primary Lagrangian droplets into secondary droplets.

The paper is structured as follows. Section 2 presents the multi-fluid model, the considered thermodynamic closure and the Lagrangian formalism description. The implementation of the surface tension numerical model is presented in Section 3. Finally, Section 4 presents the Euler-Lagrange transition algorithm. Finally, Section 5 presents a comprehensive comparison between experimental and numerical results for a Jet-in-Crossflow (JICF) configuration.

## 2. Governing equations and numerical models

The numerical simulations are performed with the massively parallel, explicit Large-Eddy Simulation (LES) solver AVBP developed at CERFACS [31, 77]. The numerical scheme is a two-step Taylor Galerkin method [18], 3rd order accurate in space and time. To mitigate numerical oscillations and improve numerical stability during the transient regime, a density-based selective sensor [75] is used. This sensor introduces artificial viscosity in regions characterized by steep thermophysical gradients. Once the flow is established, the artificial viscosity model is configured to be activated exclusively for the most severe phenomena. Unresolved scales are modelled thanks to the WALE Subgrid Scale (SGS) model [56]. Constant turbulent Schmidt and Prandtl numbers (fixed at 0.6 for both) are used for the sub-grid species and heat diffusion terms.

### 2.1. System of equations

Assuming equilibrium of pressure, velocity, temperature and chemical potential between the liquid and the gas, the 3-equations [6, 61] homogeneous equilibrium for the Favre-filtered fully compressible Navier-Stokes equations yields to Eqs. (1):

$$ \frac{\partial \bar{\rho} \tilde{Y}_k}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{Y}_k \tilde{\mathbf{u}}) = - \nabla \bar{\mathcal{J}}_k - \nabla \mathcal{J}_k^t + \bar{\omega}_k \quad \text{with} \quad k = [1, \dots, N_s] \qquad (1a) $$

$$ \frac{\partial \bar{\rho} \tilde{\mathbf{u}}}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{\mathbf{u}} \otimes \tilde{\mathbf{u}} + \bar{P} \underline{\mathbf{I}}) = \nabla (\mu \underline{\tau}) + \nabla (\mu \underline{\tau}^t) \qquad (1b) $$

$$ \frac{\partial \bar{\rho} \tilde{e}_{tot}}{\partial t} + \nabla \cdot (\bar{\rho} e_{tot} + \bar{P}) \tilde{\mathbf{u}} = \nabla (\mu \underline{\tau} \cdot \tilde{\mathbf{u}}) - \nabla \underline{\mathbf{q}} - \nabla q^t + \bar{\omega}_T \qquad (1c) $$
where $\bar{\tau}$ and $\bar{\tau}$ are the spatial and mass-weighted Favre spatial filters, respectively. The density, velocity and pressure are respectively represented by $\rho$, $\mathbf{u}$ and $P$. The overall $N_s$ (gas + liquid) species mass fractions are $Y_k = m_k/m$. The viscous stress tensor is noted $\underline{\tau}$ and the total energy is $e_{tot} = e_c + e_s$, with $e_c = \mathbf{u}^2/2$ the kinetic energy and $e_s$ the sensible energy. The heat fluxes are denoted by $q$. The subgrid scale fluxes, represented by the subgrid scale stress tensor $\underline{\tau}^t$, SGS energy fluxes $q^t$ and SGS species fluxes $J_k^t$, are considered in the formulation. Furthermore, the species reaction rate and heat release rate are denoted as $\dot{\omega}_k$ and $\dot{\omega}_T$, respectively. In this approach, the liquid phase and the gas phase are likely to coexist, and the mixture density and energy are defined in Eq. (2):

$$ \rho = \alpha_l \rho_l^{SAT} + (1 - \alpha_l) \rho_g^{SAT} \qquad (2a) $$

$$ \rho e_s = \alpha_l \rho_l^{SAT} e_{s,l}^{SAT} + (1 - \alpha_l) \rho_g^{SAT} e_{s,g}^{SAT} \qquad (2b) $$

where $\alpha_l$ is the liquid volume fraction defined in Eq. (3) and the superscript SAT represents the properties at saturation. Finally, subscripts $l$ and $g$ respectively denote the liquid and the gas phase.

$$ \alpha_l = \frac{\rho - \rho_g^{SAT}}{\rho_l^{SAT} - \rho_g^{SAT}} \qquad (3) $$

## 2.2. Thermodynamic closure

In the present work, the thermodynamic closure relies on the cubic equation of state proposed by Peng & Robinson [62] (PR) and defined in Eq. (4):

$$ P(\rho, T, \underline{Y}) = \frac{\rho \bar{r} T}{1 - \bar{b} \rho} - \frac{\bar{a}(T) \rho^2}{1 + 2 \bar{b} \rho - \bar{b}^2 \rho^2} \qquad (4) $$

where $\bar{r} = R/\bar{W}$ with $\bar{W}$ being the mixture molar mass, $\bar{a}$ and $\bar{b}$ the attractive interaction coefficient and the mixture covolume. Those coefficients are computed following the Van der Waals mixing laws detailed in Poling et al. [66].

The determination of the two-phase thermodynamic state is achieved thanks to the minimization of a $N_s$-variable function, where $N_s$ is the number of chemical species, in order to minimize the Gibbs energy of the whole system. Some authors have relied on such an approach for large eddy simulation
of Diesel [54] or rocket engines [86] applications. The exact computation of the thermodynamic equilibrium is however expensive, particularly as the number of species increases. To circumvent this issue, Pelletier et al. [61] have proposed an approximate two-phase equilibrium formulation for multi-components systems which keeps the same complexity regardless of the number of species involved in the interface region. This approximate equilibrium assumes that the liquid and gaseous phases have the same chemical composition, i.e., for any species k, $Y_{k,l} = Y_{k,v} = Y_k$. This leads to the equilibrium assumptions defined in Eq. (5).

$$P_l^{SAT} = P_v^{SAT} \qquad (5a)$$

$$T_l^{SAT} = T_v^{SAT} \qquad (5b)$$

$$g_l^{SAT} = g_v^{SAT} \qquad (5c)$$

where g is the partial Gibbs energy. Finally, $\rho_l^{SAT}$ and $\rho_v^{SAT}$ are obtained thanks to a tabulation of the saturation curves as presented in Pelletier et al. [61].

## 2.3. Lagrangian formalism

The Lagrangian droplets are characterized by their coordinates $X_p$ and transported at each temporal iteration thanks to their velocities $U_p$ and directions vectors. In the context of massively parallel computation, the approach proposed in Haselbacher et al. [37] is considered to locate particles in the mesh. The unsteady forces of Basset [7], Boussinesq [13] and Oseen [57] are neglected, as proposed in Elghobashi & Truesdell [25] for high density ratio ($\rho_{ptcl}/\rho_{gas} \gg 1$). Finally, considering a dilute spray regime and small, rigid and spherical particles, the Lagrangian equations for droplets motion and energy are Eqs. (6–9):

$$\frac{d\mathbf{X}_{ptcl}}{dt} = \mathbf{U}_{ptcl} \qquad (6)$$

$$\frac{dm_{ptcl} \mathbf{U}_{ptcl}}{dt} = \mathbf{F}_{ptcl} \qquad (7)$$

$$\frac{dm_{ptcl}}{dt} = \dot{m}_{ptcl} \qquad (8)$$

$$\frac{dC_{p,l}T_{ptcl}}{dt} = \frac{1}{m_{ptcl}}(-\phi_g^c + \dot{m}_{ptcl}L_v(T_{ptcl})) \qquad (9)$$
with $C_{ptcl,l}$ the liquid heat capacity and $\mathbf{X}_{ptcl}$, $\mathbf{U}_{ptcl}$, $m_{ptcl}$, $T_{ptcl}$ the droplet position, velocity, mass and temperature, respectively. The sum of all external forces acting on the droplet is $\mathbf{F}_{ptcl}$, $\dot{m}_{ptcl}$ is the evaporation rate, $L_v$ the latent heat of evaporation and $\phi_g^c$ the gaseous conductive heat flux at the droplet surface. In the present Lagrangian solver the liquid density and surface tension coefficient are assumed to be constant. Gravity is neglected in our configuration and the external forces $\mathbf{F}_{ptcl}$ (Eq. (7)) are reduced to a drag force (Eq. (10)) which assumes a stationary flow in the particle vicinity:

$$ \mathbf{F}_{drag} = \frac{1}{2} \rho_g C_{drag} S_{ptcl} ||\mathbf{U}_g - \mathbf{U}_{ptcl}|| (\mathbf{U}_g - \mathbf{U}_{ptcl}) \quad (10) $$

where $S_p$ is the droplet surface area and $C_{drag}$ the drag force coefficient. The drag force coefficient is modelled through the Schiller & Nauman [74] empirical correlation whereas evaporation is eventually computed through the Spalding [81] approach. Such approaches have the advantage of being adapted to massively parallel solvers, easy to implement and robust. Note that the drag force relies on the Eulerian gaseous velocity $\mathbf{U}_g$, which itself accounts for the turbulent dissipation provided by the WALE model.

Secondary atomization is ensured by the secondary atomization model FAST, implemented in AVBP by Senoner [78], which builds upon the prior contributions of Gorokhovski [34] and Apte et al. [2]. Notably, FAST has demonstrated promising results in previous realistic aeronautical configurations, as highlighted in the works of Treleaven et al. [85] and Carmona et al. [16]. Details of the implemented model can be found in these recent works. The triggering of the secondary atomization model relies on two key parameters, namely the break-up time $t_{bu}$ and the critical radius $r_{crit}$, as defined in Eqs. (11b) and (11a):

$$ t_{bu} = B \frac{2r_{ptcl}}{u_{rel}} \sqrt{\frac{\rho_l}{\rho_g}} \quad (11a) $$

$$ r_{crit} = \frac{We_{crit}\sigma}{2\rho_g u_{rel}^2} \quad (11b) $$

where $B = \sqrt{3}$, $r_{ptcl}$ is the particle radius, $\sigma$ is the surface tension coefficient and $u_{rel}$ represents the relative velocity between the Lagrangian particle and the gaseous Eulerian field. The critical Weber number $We_{crit}$ only exhibits minor variations for usual spray characteristics, typically around the
value of 12. Pilch & Erdman [64] proposed a correlation to determine the critical Weber number (Eq. (12)):

$$ We_{crit} = 12 (1 + 1.077Oh^{1.6}) \qquad (12) $$

with,

$$ Oh = \frac{\mu_l}{\sqrt{\rho_l 2 r_{ptcl} \sigma}} \qquad (13) $$

A Lagrangian particle is atomized thanks to the secondary atomization model if both conditions of Eq. (14) are satisfied:

$$ r_{ptcl} > r_{crit} \qquad (14a) $$

$$ t_{ptcl} > t_{bu} \qquad (14b) $$

where $t_{ptcl}$ is the lifetime of the droplet. Finally, the relative velocity is identified as a critical parameter, determining the smallest droplet diameter, which is likely to vary depending on the local flow conditions.

### 3. Surface tension implementation

The surface tension is modeled using the incompressible formulation proposed by Brackbill et al. [14]. In order to handle sharp density gradients which cannot be explicitly computed on the numerical mesh, Brackbill et al. [14] proposed to apply surface tension over a transition zone of approximately one cell width, using a color function. In that respect, the sharp discontinuity is replaced by a zone of continuous transition between liquid and gas. The resulting volume force called the Continuum Surface Force (CSF) mimics the effect of the surface force, and follows Eq. (15):

$$ F_{ST}(x) = \sigma\kappa(x)\nabla C \qquad (15) $$

where $\sigma$ is the surface tension coefficient, $\kappa$ the curvature, and $C$ a given color function. In the context of diffuse interface models, the mass or volume fraction is often used as color functions, as
in [43, 46, 63, 72]. The present work proposes to rely on the volume fraction $\alpha_l$ as the color function to account for surface tension forces in the 3-equation model, following Eq. (16):

$$ \mathbf{F}_{ST}(x) = -\sigma\kappa(x)\nabla\alpha_l \qquad (16) $$

with,

$$ \kappa(x) = -\nabla \cdot \left( \frac{\nabla \alpha_l}{||\nabla \alpha_l||} \right) \qquad (17) $$

Finally, Perigaud & Saurel [63] extended the incompressible CSF formulation to a compressible context. This paper adopts their non-conservative formulation due to its ease of implementation. The final system given in Eq. (18) is written as:

$$ \frac{\partial \rho Y_k}{\partial t} + \nabla \cdot (\rho Y_k \mathbf{U}) = 0 \quad \text{with} \quad k = [1, \dots, N_s] \qquad (18a) $$

$$ \frac{\partial \rho \mathbf{U}}{\partial t} + \nabla \cdot (\rho \mathbf{U} \otimes \mathbf{U} + P \underline{\mathbf{I}}) = \mathbf{F}_{ST} \qquad (18b) $$

$$ \frac{\partial \rho e_{tot}}{\partial t} + \nabla \cdot (\rho e_{tot} + P) \mathbf{U} = \mathbf{U} \mathbf{F}_{ST} \qquad (18c) $$

The surface tension implementation has been previously validated on various academic cases, as presented in [32, 60].

## 4. Euler-Lagrange transition model

The considered methodology in this paper is sketched in Fig. 1. Once atomized from the liquid core, liquid structures that are considered detached are transformed into Lagrangian particles, if satisfying some given geometrical criteria. Liquid structures are considered detached when their surfaces form a closed boundary. The algorithm first determines the surface topology of liquid structures, evaluates their volume, and replaces the Eulerian liquid region by Eulerian gas. Lagrangian particles are finally introduced in order to account for the exactly same volume of liquid.

### 4.1. Eulerian liquid structure surface identification

The surface determination of Eulerian liquid structures follows the methodology proposed in Rochette et al. [70] for flame front detection. Their method relies on the use of fictive Lagrangian
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the overall workflow of the atomization algorithm (left) illustrated on a three-dimensional co-axial jet (right)</figcaption></figure>

particles as scout tracers to fetch and retrieve the front location. With the help of the Haselbacher algorithm [36] to locate particles, this technique was found the most efficient in the context of parallel unstructured mesh computation. The surface identification steps are depicted Fig. 2. A
<figure><img src="image_3.png" /><figcaption>Figure 2: Surface identification steps, using fictive particles as scout tracers. Blue shaded field represents the $\alpha_l$ field. Dark blue solid line represents the iso-contour for $\alpha_l = \alpha_l^{thresh}$. Red circles represent scout tracers. Grey squares represent flagged interface cells.</figcaption></figure>

threshold value of liquid volume fraction $\alpha_l^{thresh}$ defines what should be considered as liquid or gas in this diffuse interface context, as shown in Fig. 2a (dark blue line). All cells crossed by the iso-contour $\alpha_l = \alpha_l^{thresh}$ are flagged as interface cells (grey cells in Fig. 2b). The condition for these cells to be flagged as interface cells is that at least one of its node fulfills the condition $\alpha_l > \alpha_l^{thresh}$ (treated as liquid) and at least one does not (treated as gas), as illustrated in Fig. 3.

Once the interfacial cells are located, tracers are initialized in these cells, as shown in Fig. 2c. These Lagrangian tracers transport a number of properties, namely the initial interfacial cell ID, the partition number, and their initial position. The scout particles are sent from their initial position into the domain following the interface normal direction, as illustrated in Fig. 2c. These scout particles travel in the domain until they reach an opposite interfacial cell (Fig. 2d) or exit the domain. The distance between the initial position and the opposite interface cell is then stored, and associated to the scout tracer. This distance allows to measure a characteristic size for the Eulerian cluster.
<figure><img src="image_4.png" /><figcaption>Figure 3: Example of tetrahedral interfacial cell crossed by the iso-surface $\alpha_l = \alpha_l^{\text{thresh}}$</figcaption></figure>

These identified interface cells are then associated to a unique surface ID, characteristic of the considered liquid structure. To do so, a connectivity matrix is constructed, in which each interface node is linked to its surrounding interface cells. Fig. 4 illustrates interface cells and corresponding nodes, leading to the connectivity matrix Eq. (19).

$$ \left( \begin{array}{c|ccc} \text{Node id} & \multicolumn{3}{c}{\text{Surrounding interface cells}} \\ \hline \textcolor{red}{n_{c_1}} & c_1 & c_2 & \cdot \\ \textcolor{red}{n_{c_2}} & c_2 & c_3 & c_4 \\ \textcolor{red}{n_{c_3}} & c_3 & c_4 & \cdot \\ \textcolor{red}{n_{c_4}} & c_4 & c_5 & \cdot \end{array} \right) \qquad (19) $$

The dimensions of the matrix are $[n_{n_{S_i}} \times n_{cell}]$, where $n_{n_{S_i}}$ is the number of nodes in the cluster and $n_{cell}$ the number of cells in the partition. This connectivity matrix is then reduced, comparing each line to the following one. If at least one cell is shared between these two lines, the list of cells associated to the higher node number is copied to the line associated to the lower node number, and discarded. This allows to reduce the connectivity matrix to a single line, Eq. (20) for our example.
$$ (n_{c_1} \ c_1 \ c_2 \ c_3 \ c_4 \ c_5) \qquad (20) $$

This single line corresponds to a single node ID ($n_{c_1}$), which can be associated to the entire Eulerian identified surface structure.

<figure><img src="image_5.png" /><figcaption>Figure 4: Scheme of the clusterization process. Grey cells represent interfacial cells, blue cells represent liquid cells, white cells represent gaseous cells and orange dots represent interfacial nodes. The dark blue line represents the iso-contour $\alpha_l = \alpha_l^{thresh}$. $nc_i$ is the smallest ID node belonging to cell $c_i$.</figcaption></figure>

## 4.2. Dealing with parallel computation

In the context of massively parallel solver, the identification node number may not be unique. Hence, uniqueness of such identification number is ensured according to Eq. (21) :

$$ ID_{struc} = n_{part}N_{node}^{max} + ID_{local} \qquad (21) $$

where $ID_{struc}$ is the single liquid structure identification number, $n_{part}$ the number of partitions, $N_{node}^{max}$ the maximal number of nodes in partitions and $ID_{local}$ the local identification number, e.g $n_{c_1}$ in Eq. (20). Since a liquid structure can spread over several partitions, the highest $ID_{struc}$ is given to the partition shared nodes, then to the whole interfacial nodes. The different steps
<figure><img src="image_6.png" /><figcaption>Figure 5: Definition of the global ID of a detached structure</figcaption></figure>

described above are depicted Fig. 5.

### 4.3. Eulerian liquid structure volume identification

The surface clusterization described above encloses a given volume of liquid material, that needs to be identified as part of the liquid structure. Similarly as for the surface clusterization, the use of scout tracers is once again considered here. They are initialized in previously identified interface cells, and transported into the liquid structure following the interface normal direction, while carrying the global interface identification number (Fig. 6a). When reaching another interface, its identification number is compared to the one carried by the scout tracer (Fig. 6b). If similar, the scout particle is sent back, and provides the identification number to all cells crossed (Fig. 6c). The scout particle is deleted otherwise. The isolated liquid structure is finally defined as a group of cells sharing the same identification number (Fig. 6d). Only closed structures, for which interface cells are continuously adjacent, are transformed into a Lagrangian particle.

### 4.4. Shape analysis

After the Eulerian liquid structure identification, a shape analysis is conducted. This allows to differentiate rather spherical droplets from elongated ligaments, which will require different treatments. The eccentricity computation first requires the knowledge of the structure barycenter, defined following Eq. (22). This barycenter position will serve for the average structure length definition, which will be compared to the structure width.
<figure><img src="image_7.png" /><figcaption>Figure 6: Volume clusterization. The red dot represents a scout tracer. The grey and red cells represents the interfacial cells and the inner liquid structure cells, respectively.</figcaption></figure>

$$ (\vec{x}_{\text{barycenter}})_i = \frac{\sum_{k=1}^{N_{\text{cell}_i}} (\alpha_l)_k \vec{x}_k}{\sum_{k=1}^{N_{\text{cell}_i}} (\alpha_l)_k} \qquad (22) $$

where $N_{cell_i}$ is the number of cells belonging to the liquid structure $S_i$, $\vec{x}_k$ the center position of cell $k$, and $(\alpha_l)_k$ the liquid volume fraction in cell $k$.

In some specific situations, the structure appears to be so much distorted that the geometrical barycenter (Eq. (22)) is located outside the liquid structure. In such cases, the closest liquid cell to the geometrical barycenter is considered. Finally, the liquid structure characteristic length is computed following Eq. (23):

$$ L_{\text{max}} = ||\vec{x}_{\text{furthest}} - \vec{x}_{\text{barycenter}}|| \qquad (23) $$

with $\vec{x}_{furthest}$ the position of the furthest cell from the barycenter position (Fig. 7). The average liquid structure width is computed from the widths $\Gamma_k$ evaluated thanks to the scout tracers during
the surface structure identification (Section 4.1, Fig. 2), following Eq. (24) :

$$ \Gamma_{S_i} = \frac{\sum_{k=1}^{N_{cell_i}} (\Gamma_k)_{S_i}}{N_{cell_i}} \qquad (24) $$

Therefore, the computation of liquid structure eccentricity, which allows the differentiation between droplets and stretched ligaments, is computed following Eq. (25)

$$ \varepsilon_{S_i} = \sqrt{1 - \frac{\left(\frac{\Gamma_{S_i}}{2}\right)^2}{(L_{max}^2)_{S_i}}} \qquad (25) $$

An eccentricity threshold $\varepsilon_{thresh}$ is chosen, below which the liquid structure is considered as a droplet and converted into a single Lagrangian particle. Otherwise, the liquid structure is considered to be a ligament. Based on a stretching factor $\Lambda_{S_i}$ defined Eq. (26), liquid structures are themselves differentiated between ligaments, transformed into several Lagrangian droplets, and large structures that are kept in the Eulerian formalism.

$$ \Lambda_{S_i} = \frac{2(L_{max})_{S_i}}{\Gamma_{S_i}} \qquad (26) $$

For a stretching factor higher than a user-defined threshold $\Lambda_{thresh}$, positions and properties of injected Lagrangian droplets are to be determined. Single Lagrangian droplets are injected at the position $\mathbf{x}_{barycenter}$. The associated velocities and temperatures correspond to the averaged Eulerian liquid velocity and averaged Eulerian temperature (Eq. (27)).

$$ \Phi_{S_i} = \frac{\sum_{k=1}^{N_{cell_i}} \Phi_k}{N_{cell_i}} \qquad (27) $$

where $\Phi_{S_i}$ is the average intensive quantity $\Phi$ for structure $S_i$ and $\Phi_k$ the cell value. When condition $\varepsilon < \varepsilon_{thresh}$ is not verified, the liquid structure is replaced by $N_D$ Lagrangian particles of diameter $d_D$ equally distributed along the main direction of the structure, characterized by the $\mathbf{r}$ vector defined in Eq. (28). In a similar way as for the single droplet, a unique velocity and temperature is computed for the $N_D$ Lagrangian droplets injected.

$$ \mathbf{r} = \mathbf{x}_{furthest} - \mathbf{x}_{barycenter} \qquad (28) $$

The $N_D$ droplets (Eq. (33)) are located along the direction of $\mathbf{r}$, with a distance $\delta_D$ between two droplets following Eq. (29):
$$ \delta_D = \frac{2(L_{max})S_i}{N_D - 1} \qquad (29) $$

Eq. (29) determines the distance between each Lagrangian droplet within the same ligament. It corresponds to the droplet diameter plus the inter-droplet spacing. Then, the position $\mathbf{x}_k$ of the $k^{th}$ droplet is computed as follows:

$$ \mathbf{x}_k = \mathbf{x}_{\text{barycenter}} + \mathbf{r} - \mathbf{r}_k \quad \forall k \in \{1, \dots, N_D\} \qquad (30) $$

with,

$$ \mathbf{r}_k = (k-1)\delta_D \frac{\mathbf{r}}{\|\mathbf{r}\|} \qquad (31) $$

The Lagrangian droplets must then be located within the mesh, along with the identification numbers of the injection cells. For the injection of multiple particles (ligaments: $\epsilon > \epsilon_{thresh}$), determining the IDs of the injection cells would require looping over the entire domain, significantly increasing the computational cost. Furthermore, the massively parallel context of AVBP would make this task even more complicated. To address this issue, the Haselbacher algorithm [36] is considered once again. The droplet is therefore created at its final position $\mathbf{x}_k$ (Fig. 7).

<figure><img src="image_8.png" /><figcaption>Figure 7: Sketch of a liquid structure transferred from the Diffuse Interface Method (DIM) Eulerian formalism to the Lagrangian formalism. $D_k$ is the distance $|| \mathbf{x}_k - \mathbf{x}_{\text{barycenter}} ||$</figcaption></figure>
The location at $x_k \neq x_{barycenter}$ may require to cross a partition limit. As transport among the partitions in the domain is already embedded in the Lagrangian solver, very few additional treatments are required. The determination of particle location is illustrated in Fig. 7.

### 4.5. Mass, Momentum and Energy Transfer

#### 4.5.1. Mass transfer

As illustrated Fig. 7, the detached liquid structure is replaced by $N_D$ Lagrangian droplets. A first estimation of the mass of each droplet $m_{D_{approx}}$ is computed using Eq. (32).

$$m_{D_{approx}} = \frac{4}{3}\pi \left(\frac{a_{S_i}}{2}\right)^3 \rho_l \quad (32)$$

with $a_{S_i}$ the estimated diameter of the particle that corresponds to the mean liquid structure width $\Gamma_{S_i}$ (Eq. (24)) and $\rho_l$ the imposed liquid Lagrangian particle density, which is assumed constant in the present LPT approach. Considering the effective ligament mass $m_{S_i}$, the number of droplet $N_D$ of mass $m_{D_{approx}}$ is provided Eq. (33).

$$N_D = \left[ \frac{m_{S_i}}{m_{D_{approx}}} \right] \quad (33)$$

with $\lfloor \rfloor$ a floor function and $m_{S_i}$ defined as,

$$m_{S_i} = \sum_{k=1}^{N_n} \Delta m_{n_k} = \Delta \rho \cdot V_{n_k} \quad \forall k \in \{1, \dots, N_n\} \quad (34)$$

where $V_n$ represents the nodal volume of node $k$, $N_n$ the node number inside the candidate liquid structure and $\Delta \rho$ denotes the difference between the liquid density and the one of the replacing gas. Naturally, there is a residual liquid mass noted $m_r$, defined Eq. (35).

$$m_r = m_{S_i} - N_D m_{D_{approx}} \quad (35)$$

A part of this residual mass is given to each droplet equally, allowing to compute the final mass $m_D$ of each droplet using Eq. (36).

$$m_D = m_{D_{approx}} + \frac{m_r}{N_D} \quad (36)$$

Finally, the Lagrangian droplets diameter $d_D$ of a liquid structure is computed using Eq. (37).
$$d_D = 2 \sqrt[3]{\frac{3m_D}{4\pi\rho_l}} \qquad (37)$$

The geometric properties of each Lagrangian particle to be injected are therefore known.

Using Eqs. (32–37), the liquid mass is entirely transferred to the Lagrangian formalism ensuring that the liquid mass is conserved. As the transition between the two solvers transforms an Eulerian volume of liquid $V_D$ into a point source particle, Eulerian empty spaces need to be filled with gas. The density of the replacing gas $\rho_{g_{input}}$ must be determined.

In the context of multi-component mixtures, the replacing gas can be pure fuel, pure oxydizer or a mixture of both. In this work, the last option is retained. The idea is that the droplet is surrounded by its own vapor as well as a part of the surrounding gas. Therefore, the empty spaces created by the transition algorithm is replaced by a gas mixture equivalent to the liquid structure surrounding mixture. From the algorithmic point of view, 3 steps are needed. The first step is to determine which nodes are considered as belonging to the droplet’s surrounding gas. The second step is to get the droplet’s surrounding gas averaged pressure, temperature and mass fraction. Finally, $\rho_{g_{input}}$ is computed using Eq. (38):

$$\rho_{g_{input}} = \rho_{vap}^{SAT}(P_{mean}, T_{mean}, Y_{k,mean}) \qquad (38)$$

where the supercript SAT represents the saturation conditions. For the first step, the considered nodes are those belonging to a cell crossed by the user defined iso-contour $\alpha_l = \alpha_l^{thresh}$, as shown Fig. 8. Then the average pressure, temperature and mass fraction are computed thanks to Eq. (27). To sum up, the liquid mass is totally transferred to the Lagrangian formalism and is replaced in the Eulerian multi-fluid formalism by a gaseous mixture whose physical properties correspond to those of the gaseous surrounding mixture. Note that adding an Eulerian mass of gas is required due to the point source approximation of the Lagrangian solver.

### 4.5.2. Momentum transfer

As a consequence of the liquid transfer between the Eulerian and the Lagrangian formalism, the Eulerian momentum must be updated in the Eulerian conservation equations (Eq. (18)), following Eq. (39).
<figure><img src="image_9.png" /><figcaption>Figure 8: Sketch of the retained nodes used to compute the replacing gas properties</figcaption></figure>

$$ (\rho \mathbf{u})^{\text{new}} = (\rho \mathbf{u})^{\text{old}} - \Delta \rho \frac{(\rho \mathbf{u})^{\text{old}}}{\rho^{\text{old}}} \qquad (39) $$

with,

$$ \Delta \rho = \rho^{\text{old}} - \rho_{g_{\text{input}}} \qquad (40) $$

where the superscripts old and new respectively denote the properties before and after the transfer of Eulerian liquid structures into the Lagrangian formalism.

Finally, all created Lagrangian particles from a liquid structure have a unique velocity that corresponds to the considered Eulerian liquid structure mean velocity.

### 4.5.3. Energy transfer

The total Eulerian energy transported by the numerical solver is reminded in Eq. (41):

$$ \rho e_{\text{tot}} = \rho e_s + \rho e_c \qquad (41) $$
where $e_s$ represents the sensible energy while $\rho e_c = \frac{1}{2}\rho||u||^2$ represents the kinetic energy. First, the focus is put on the variation of $e_s$. Since the predicted Eulerian state is single phase and the pressure $P_{mean}$ and density $\rho_{g_{input}}$ are known, the temperature $T_{mean}$ of the new gaseous state is computed by solving the cubic equation of state (Eq. (4)) using an iterative method. Then, the sensible energy of the gas is computed using Eq. (42):

$$e_s^{new} = e_s(\rho_{g_{input}}, T_{mean}) \qquad (42)$$

In the Lagrangian framework, the thermodynamic knowledge of the droplet lies in its temperature $T_{drop}$. In the present work, $T_{drop} = T_{mean}$ because the surrounding gas and liquid temperature values are generally close at the breakup instant.

Finally, the kinetic energy in the Eulerian multi-fluid formalism is written as presented in Eq. (43):

$$e_c^{new} = \frac{1}{2}\rho_{g_{input}}(u^2 + v^2 + w^2) \qquad (43)$$

### 4.6. Flowchart

A flowchart illustrating the proposed transition algorithm is shown in Fig. 9, with the two geometrical tests (A) and (B).

## 5. Validation of Euler-Lagrange transition algorithm: Application to a Jet in Cross flow configuration.

### 5.1. Operating conditions

The Euler-Lagrange transition algorithm presented in the previous sections is assessed on a jet in cross flow (JICF) configuration. Test case is based on the experimental study of Song et al. [80]. The computational domain is sketched Fig. 10. The injector diameter is $d_0 = 0.508$ mm and is located 30 mm downstream of the air inlet. Experimental measurements are conducted in planes P1 and P2, depicted in Figure 10. The jet trajectory is measured in plane P1, while the droplet size distribution is measured in plane P2. The operating conditions are summarized in Tab. 1. The momentum flux ratio $q$ is defined as $q = (\rho_l u_l^2)/(\rho_g u_g^2)$. The aerodynamic Weber number is defined
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Flowchart of the transition algorithm</figcaption></figure>
as $We_a = (\rho_g u_g^2 d_0)/\sigma$ and the liquid Reynolds number is $Re_l = (\rho_l u_l d_0)/\mu_l$. When the liquid is transferred to the Lagrangian formalism, a constant particle density ($\rho_{ptcl} = \rho_l$) is assumed.

<table><caption>Table 1: Flow properties</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Momentum flux ratio</td><td>q</td><td>10</td><td>—</td></tr><tr><td>Aerodynamic Weber</td><td>We<sub>a</sub></td><td>500</td><td>—</td></tr><tr><td>Liquid Reynolds</td><td>Re<sub>l</sub></td><td>300</td><td>—</td></tr><tr><td>Surface tension coefficient</td><td>σ</td><td>0.0234</td><td>N.m<sup>-1</sup></td></tr><tr><td>Liquid Density</td><td>ρ<sub>l</sub></td><td>830</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Gas Density</td><td>ρ<sub>g</sub></td><td>8</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Liquid Temperature</td><td>T<sub>l</sub></td><td>294</td><td>K</td></tr><tr><td>Gas Temperature</td><td>T<sub>g</sub></td><td>294</td><td>K</td></tr><tr><td>Pressure</td><td>P<sub>bulk</sub></td><td>6.89</td><td>bar</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Figure 10: Computational domain for Jet In Cross Flow (JICF) case: geometrical quantities</figcaption></figure>
## 5.2. Boundary conditions

The outlet boundary condition consists in a characteristic non-reflecting boundary condition [65], imposing static pressure $P_s = 6.89$ bar, corresponding to the experimental bulk pressure. Both inlet boundary conditions use a characteristic non-reflecting inlet [19] formulation. For the air inlet, mass flow rate and temperature are imposed, while for the liquid inlet density and velocity are imposed. Adiabatic law-of-the-wall is imposed on walls, except within the liquid injector for which slip conditions are imposed, due to the limited number of cells in the injector. The normalized wall distance $y^+$ lies between 45 and 350 along the bottom flat plate, from the gaseous injection to the outlet, as illustrated in Fig. 11.

<figure><img src="image_13.png" /><figcaption>Figure 11: Bottom plate $Y^+$ field</figcaption></figure>

## 5.3. Comparison with experimental data

Comparisons between the simulation and the experiment focus on the jet trajectory measured in plane P2, and on the droplet size distribution measured in plane P1.

### 5.3.1. Mesh sensitivity

To assess the mesh sensitivity of the whole process consisting in an Eulerian liquid injection, transition to Lagrangian particles, and secondary atomization, 4 simulations are performed on different mesh resolutions. These meshes respectively consist of 30 millions, 57 millions, 79 millions and 87 millions tetrahedral cells. The mesh refinement mainly targets the jet liquid core vicinity, while few refinements are conducted in the jet plume region, as shown Fig. 12. Details regarding mesh
properties are provided in Tab. 2. This study aims at investigating the influence of refining the liquid core and primary droplets surrounding flow on the droplet distribution and the jet trajectory. The computational time required by these meshes is 10 000 hours (M0), 16 000 hours (M1), 58 000 hours (M2), and 88 000 hours (M3) per convective time (1.5 ms).

<table><caption>Table 2: Cell characteristic length ($\Delta_x$) for the 4 meshes at various locations.</caption><thead><tr><th>Mesh</th><th>Injector</th><th>Refined liquid region</th><th>Plume</th></tr></thead><tbody><tr><td>M0</td><td>0.07 mm</td><td>0.1 mm</td><td>0.145 mm</td></tr><tr><td>M1</td><td>0.07 mm</td><td>0.083 mm</td><td>0.13 mm</td></tr><tr><td>M2</td><td>0.036 mm</td><td>0.053 mm</td><td>0.11 mm</td></tr><tr><td>M3</td><td>0.032 mm</td><td>0.04 mm</td><td>0.11 mm</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Figure 12: Center plane meshes, focus on the near injection and plume region</figcaption></figure>

### 5.3.2. Jet topology: a brief overview

Snapshots of M1 and M3 results are depicted in Fig. 13, where the liquid core is visualized as an Eulerian iso-contour with a threshold value of $\alpha_l = 0.5$. Lagrangian particles are visually depicted as turquoise spheres, representing their physical volume and position. The largest particles correspond to the primary droplets generated by the transition algorithm, while the smallest particles represent the droplets resulting from the secondary atomization algorithm. A more coherent liquid core is observed in the M3 case, in contrast to M1. Furthermore, the M3 case exhibits a greater number of primary droplets with more dispersed initial diameters. According to previous section,
<figure><img src="image_15.png" /><figcaption>Figure 13: Instantaneous snapshots of liquid core (iso-surface $\alpha_l = 0.5$) and Lagrangian particles for M1 ((a)) and M3 ((b))</figcaption></figure>

the two droplet populations are also observed in the M3 case. In the following sections, a mesh sensitivity study is performed to investigate jet trajectory and droplet distribution in the spray. Then, a quantitative analyse is conducted on both Eulerian and Lagrangian fields.

### 5.3.3. Mesh sensitivity on jet trajectory, $D_{32}$ and droplet size distribution

The position of the maximum Mie scattering intensity, corresponding to the maximum droplet volumetric flux, is reported in Song et al. [80]. It is compared to the numerical temporally- and spatially- averaged droplet volumetric fluxes $V_f$, shown in Fig. 14, with $V_f$ defined in Eq. (44):

$$V_f = \frac{\overline{\rho_{ptcl} u_{ptcl} \alpha_l}}{\langle \rho_{inj} u_{inj} \rangle} \qquad (44)$$

with $\rho$ the droplet density, $u$ the axial velocity and subscripts ptcl and inj corresponding to liquid injection (Euler) and Lagrangian particles properties. The temporal averaging operator is repre-
<figure><img src="image_16.png" /><figcaption>Figure 14: Volumetric flux field $V_f$ displayed on the P2 plane in shaded blue for the M0, M1, M2 and M3 meshes. The black curve represents the experimental maximum intensity mie-scattering</figcaption></figure>

sented as < . >, while the spatial averaging operator is denoted by . . The considered averaging
times are 10 ms (6.5 convective times) on M0 and M1 and 8 ms (5.3 convective times) on M2 and
M3. The spatial average is computed by azimuthally averaging 11 *Oxy* planes in the jet diameter,
with *z* ∈ [-0.25, +0.25] mm.

The maximum volumetric flux $V_f$, shown in Fig. 14 reaches higher vertical positions with increasing
mesh resolution, but the position is still under-estimated. Also note that field of $V_f$ becomes more
and more homogeneous with mesh refinement, indicating that more primary droplets are created
thanks to the Euler-Lagrange transition algorithm, having smaller diameters, and evolving into a
more homogeneous spray through secondary atomization processes. Finally, zoomed-in regions in
the liquid injection region evidence an increase of droplets injected with a dominant vertical tra-
jectory. Note that the M3 mesh exhibits closer agreement with the measurements.

The position of the maximum volumetric flux is not a common data reported in the literature,
contrarily to the outer jet trajectory, for which several correlations exist ([9, 29, 51, 69]), in studies
considering the same fuel and similar flow conditions. These correlations are reported in Fig. 15,
compared to the droplet number fields for all meshes. Note that Lubarsky et al. [51] reported two
<figure><img src="image_17.png" /><figcaption>Figure 15: Droplet number field over the P2 plane for M0, M1, M2, M3 mesh refinements. Markers represent various experimental correlations</figcaption></figure>

correlations, one for the mean outer trajectory position $y_{mean}$, the other corresponding to maximum spray penetration ($y_{max} = y_{mean} + 2.8y_{rms}$). The droplet number is derived by superimposing all droplets within the range of $z_{drop} = [-0.25, 0.25]$ mm from instantaneous solutions over an averaging time of 6.5 ms.

First, it can be observed that the outer trajectory for the coarsest mesh (Fig. 15) is significantly underestimated when compared to the correlations. This can be attributed to the initial trajectories which, as previously shown in Fig. 14, follow a nearly horizontal direction. The outer trajectory for M2 to M3 are similar and are in reasonable agreement with the experimental correlations, as depicted in Fig. 15.

The comparison between numerical prediction of the outer jet trajectory and experimental results has been reported by several groups relying on sharp interface methods in recent numerical works
([39, 48]), with a general good agreement. However, the numerical prediction of the Sauter Mean Diameter (SMD), or droplet size distribution is rarely assessed. Comparison between numerical results and experimental ones for the SMD and droplet sizes distribution is conducted in this subsection. The SMD is defined in Eq. (45):

$$D_{32} = \frac{\sum_{i=1}^{N_{drop}} D_i^3}{\sum_{i=1}^{N_{drop}} D_i^2} \qquad (45)$$

A correlation for the $D_{32}$ diameter at plane P1 and the flow properties is proposed in [80] and reported in Fig. 16. The numerical result of the present work is reported with blue circle. An overestimation of $D_{32}$ on the mesh M0 can be evidenced, whereas an excellent agreement is found for other mesh refinements. If different typical sizes of primary droplets are generated thanks to the Euler-Lagrange transition algorithm depending on the mesh refinement, the secondary atomization model allows to produce a final diameter $D_{32}$ which is almsot no sensitive to the mesh refinement.

Going further, the droplet size distribution is shown in Fig. 17. The horizontal axis corresponds to the droplet diameter, and the vertical axis to the volume fraction of each droplet size to the total volume of droplets. Experimental data are represented with black circles, whereas the present simulations results are represented in colored symbols. The numerical PDFs are in reasonable agreement with the experimental data, although an underestimation of the peak diameter is observed for all meshes. Additionally, too many large droplets ($D > 110$ µm) are still present at this axial location compared to the experimental data. A low mesh sensitivity is also found for the prediction of the volumetric droplet size distribution.

### 5.3.4. Effect of mesh resolution on Eulerian fields

This section investigates several Eulerian fields in order to enhance our understanding of jet trajectory. The time-averaged density fields around the liquid injection for all mesh refinements are presented Fig. 18. First, it can be observed that increasing mesh refinement leads to the sharpening of density gradients along x-axis, due to a lower horizontal entrainment of the Eulerian liquid core. Furthermore, only the M2 (Fig. 18c) and M3 (Fig. 18d) cases exhibit a curved liquid core, indicating that the jet is coherent enough to be deflected by the gas flow. The axial expansion of the interface is significantly reduced in the M3 case. To complete this analysis, the averaged vertical
<figure><img src="image_18.png" /><figcaption>Figure 16: Evolution of SMD with respect to flow properties. The black curve corresponds to the experimental correlation obtained by [80], while the black circles represent the experimental data obtained at various locations and under different flow conditions. Circular marker represents M0 configuration, cross-shaped marker represents reference M1 configuration, and square-shaped marker represents M2 configuration. The SMD values of the present work cases are obtained at a single operating point ($We_a = 500, q = 10$) and located at P1.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17: Volumetric droplet distribution at the P1 plane location for different mesh refinements, where dark circles represent the experimental data and white circle the DIM injection results.</figcaption></figure>

liquid velocity $V_l$, defined in Eq. (46), is depicted in Fig. 19 for all mesh refinements.
<figure><img src="image_20.png" /><figcaption>Figure 18: Eulerian density fields in the vicinity of the injection region for M0 ((a)), M1 ((b)), M2 ((c)) and M3 ((d)) mesh refinements.</figcaption></figure>

$$V_l = \frac{\overline{\alpha_l v}}{\overline{\alpha_l}} \qquad (46)$$

The presence of a negative liquid velocity region downstream of the liquid core can be evidenced. This shows that the liquid jet exhibits an unsteady vertical motion directed toward the bottom wall, as illustrated in Fig. 20.

The axial velocity $u$ and corresponding root mean square velocity $u_{rms}$ fields are respectively shown in Figs. 21 and 22. The white line in Fig. 21 represents the iso-contour corresponding to $u = 0$ m/s. Focusing on the near injection region, axial velocity fields indicate that the recirculation bubble expands and higher pressure drop is observed as the cell size decreases, suggesting a more coherent liquid core. The $u_{rms}$ fields depicted in Fig. 22 show that the thickness of the perturbed zone and the intensity of $u_{rms}$ increase as the mesh is refined. The lower gradients observed in the x-direction
<figure><img src="image_21.png" /><figcaption>Figure 19: Eulerian vertical liquid velocity fields in the vicinity of the injection region for M0, M1, M2 and M3 mesh refinements.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 20: Iso-contour of $\alpha_l = 0.5$ at 4 different instants illustrating the liquid impingement against the bottom wall for M3 case.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 21: Eulerian axial velocity field in the vicinity of the injection region for M0, M1, M2 and M3 mesh refinements. The iso-contour of $u = 0$ m/s is represented by a white line.</figcaption></figure>

are attributed to a lower diffusion of the mesh.

### 5.3.5. Effect of mesh resolution on Lagrangian fields

The effect of mesh resolution on the Lagrangian properties is investigated in this subsection. The evolution of the Lagrangian field colored by the Stokes number is depicted in Fig. 23, with the Stokes number defined in Eq. (47):

$$ \mathrm{St} = \frac{1}{18} \frac{\rho_{ptcl} d_{ptcl}^2 ||\vec{U}_g||}{\mu_g d_0} \qquad (47) $$

with $\rho_{ptcl}$ the Lagrangian particle density, $d_{ptcl}$ the Lagrangian particle diameter and $\mu_g$ the surrounding gas viscosity. Darkest regions (Fig. 23) indicate areas where droplets follow their own motion due to inertia, while lighter regions correspond to areas where the droplets motion is mainly driven by the gaseous flow. As previously mentioned, the largest primary droplets are generated on the coarser meshes with mainly streamwise trajectories. This results in higher Stokes numbers near the injector, which persist further along the streamwise direction, as depicted in Fig. 23. Overall, mesh refinement results in lower gradients within the Stokes fields, suggesting more gradual breakup processes.
<figure><img src="image_24.png" /><figcaption>Figure 22: Eulerian root mean square velocity field $u_{rms}$ for M0, M1, M2 and M3 mesh refinements.</figcaption></figure>

To complete this analysis, fields of the relative velocity, defined in Eq. (48), are shown for all meshes in Fig. 24.

$$|u_{rel}| = ||\vec{U}_{ptcl} - \vec{U}_{g}|| \qquad (48)$$

The relative velocity is a key parameter for the secondary atomization model that drives the resulting daughter particles diameter, the particle breakup time and the critical radius, following Eqs. (11a) and (11b). The main idea is that higher relative velocities increase the probability of secondary atomization processes to occur [34]. Furthermore, higher relative velocities result in a higher ratio between the diameter of the atomized mother particle diameter and the resulting daughter particle diameters. Since the primary droplets are larger on coarser meshes, they require a longer distance to reach the gas velocity, resulting in longer region of high relative velocities than on refined meshes. Moreover the relative velocity field becomes more homogeneous with mesh refinement, indicating more gradual breakup processes.

## 6. Pure Lagrangian injection

To show the benefit of the proposed methodology, a pure Lagrangian injection is conducted, as in [2, 23]. A similar mesh as M1 (Fig. 12) is considered, without the injector region. Instead, Lagrangian blobs with a constant diameter $d_{lag}$ are injected at the injector location. Two different
<figure><img src="image_25.png" /><figcaption>Figure 23: Lagrangian Stokes number (St) field for M0 ((a)), M1 ((b)), M2 ((c)) and M3 ((d)) mesh refinements.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 24: Absolute relative velocity $|u_{rel}|$ for M0 (a), M1 (b), M2 (c) and M3 (d) mesh refinements</figcaption></figure>
diameters for the initial blob are considered, $d_{lag} = 0.5$ mm and 0.25 mm, respectively corresponding to the case LAG1 and LAG2. The maximum droplet volumetric flux for both cases is shown in Fig. 25. Interestingly, an overall better agreement can be found for the jet trajectory with this pure Lagrangian injection, compared to the proposed methodology. Furthermore, the initial Lagrangian blob diameter significantly impacts the spray dispersion around the average jet trajectory. Such an effect is also noted in Fig. 26, where the effect on the initial diameter on the near field of the injector is significant. This primary analysis shows that an accurate tuning of the initial blob diameter may allow to provide an accurate jet trajectory.

<figure><img src="image_27.png" /><figcaption>Figure 25: Pure Lagrangian injection: Volumetric flux for LAG1 ((a)) and LAG2 ((b))</figcaption></figure>

However, Fig. 27 showing the resulting D32 obtained thanks to these pure Lagrangian injections, and compared with the Eulerian injection, experimental measurements and correlations, evidences the fact that a significant D32 underestimation is obtained considering only a pure Lagrangian injection. Similarly, the droplet size distribution cannot be retrieved with this pure Lagrangian injection, as shown in Fig. 28, and thus, whatever the initial blob diameter considered. Six instantaneous snapshots of the axial velocity and Lagrangian particles are shown in Fig. 29. The droplet trajectory is mainly vertical after injection. Furthermore, since Lagrangian particles have no volume, no stagnation point occurs as the crossflow interacts with the jet, leading to a misprediction of the jet wake. A more homogeneous Eulerian velocity field is crossed by the particles, which induces a low dispersion in the droplet size distribution.

Despite the fact that some improvements may be obtained by tedious tuning of the initial blob distribution size on the resulting jet trajectory compared to the proposed approach, no real im-
<figure><img src="image_28.png" /><figcaption>Figure 26: Droplet number density field over the P2 plane for LAG1 (top) and LAG2 cases (bottom). Markers represent various experimental correlations</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 27: Evolution of SMD with respect to flow properties. The black curve corresponds to the experimental correlation obtained by [80], while the black circles represent the experimental data obtained at various locations and under different flow conditions. The SMD values obtained using the Diffuse Interface Method (DIM) for the liquid injection on M2 mesh, at a single operating point (We = 500, $q = 10$) located at P1 are represented by the blue ring, while the cross and square denote the SMD values for the two pure Lagrangian cases.</figcaption></figure>

provements can be obtained on the final SMD and droplet size distribution. This evidences the need to account for the liquid core.
<figure><img src="image_30.png" /><figcaption>Figure 28: Volumetric droplet distribution at the P1 plane location with the experimental data represented by dark circles, Eulerian liquid injection results by white circles, and pure Lagrangian injections results by white crosses and squares.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 29: Plane cuts of axial velocity field at six different instantaneous times for the LAG1 case.</figcaption></figure>

## 7. Conclusion

A framework for spray prediction in a Large-Eddy Simulation context is presented, relying on a coupling between a compressible multi-fluid (MF) method and a Lagrangian Particle Tracking (LPT) approach. Large scale interfacial features including the destabilized liquid core are resolved thanks to a diffuse interface method [61], accounting for surface tension, implemented in a compressible, multi-species parallel solver. An algorithm dedicated to transfer closed liquid structures into a Lagrangian description is then proposed. An eccentricity criterion allows to differentiate
droplets, transferred into a single Lagrangian particle, and liquid ligaments, transformed into sev-
eral Lagrangian particles, if satisfying a stretching condition. The proposed algorithm is well-suited
for massively parallel and unstructured simulations. A secondary atomization model, relying on
a stochastic method, allows to transform the primary Lagrangian droplets into the final spray,
for an affordable mesh discretization. This framework is validated on a Jet-In-CrossFlow (JICF)
configuration, with a focus on the jet trajectory, the SMD and the droplet size distribution. A
very good agreement is found regarding the SMD on a rather coarse mesh, as well as a satisfying
droplet size distribution. The jet trajectory is however mispredicted on this coarse mesh. Improving
the mesh resolution allows a better liquid core prediction, and thus a better liquid jet trajectory.
Interestingly, both the Sauter Mean Diameter and the droplet size distribution downstream of the
liquid jet are not sensitive to mesh resolution. A comparison with a pure Lagrangian simulation is
finally conducted. This highlights the need to account for some liquid core for an accurate jet wake
prediction, and thus an accurate SMD and droplet size distribution. These encouraging results
demonstrate the robustness of the proposed framework to estimate the spray characteristics in a
3D pressurized configuration, at affordable computational cost.

Acknowledgments

This work has been funded by CNES and ArianeGroup and was performed using HPC from CERFACS and benefited from GENCI-CINES HPC ressources (Grant 2024 - A0152B10157). Nicolas Odier acknoweldges the French National Research Agency for the funding of the project MADIRAN (grant ANR-21-CE50-0043)
References

[1] Agbaglah, G., Delaux, S., Fuster, D., Hoepffner, J., Josserand, C., Popinet, S., Ray, P., Scardovelli, R., & Zaleski, S. (2011). Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method. *Comptes Rendus Mécanique*, **339**, 194–207.

[2] Apte, S., Gorokhovski, M., & Moin, P. (2003). LES of atomizing spray with stochastic modeling of secondary breakup. *International Journal of Multiphase Flow*, **29**, 1503–1522. URL: https://doi.org/10.1016/S0301-9322(03)00111-3. doi:10.1016/S0301-9322(03)00111-3.

[3] Ashgriz, N. (2011). *Handbook of atomization and sprays: theory and applications*. Springer Science & Business Media.

[4] Ashgriz, N. (2011). *Handbook of Atomization and Sprays: Theory and Applications*.

[5] Asuri Mukundan, A., Ménard, T., Brändle de Motta, J. C., & Berlemont, A. (2022). Detailed numerical simulations of primary atomization of airblasted liquid sheet. *International Journal of Multiphase Flow*, **147**, 103848.

[6] Barret, M., Faucher, E., & Herard, J.-M. (2002). Schemes to compute unsteady flashing flows. *AIAA Journal*, **40**, 905–913. URL: https://doi.org/10.2514/2.1727. doi:10.2514/2.1727.

[7] Basset, A. B. (1888). *A treatise on hydrodynamics: with numerous examples* volume 2. Bell and Company.

[8] Becker, J., & Hassa, C. (2002). Breakup and atomization of a kerosene jet in crossflow at elevated pressure. *Atomization and Sprays*, **12**, 49–67.

[9] Bellofiore, A., Cavaliere, A., & Ragucci, R. (2007). Air density effect on the atomization of liquid jets in crossflow. *Combustion Science and Technology*, **179**, 319–342. URL: https://doi.org/10.1080/00102200600809563. doi:10.1080/00102200600809563.

[10] Bhatia, B., Johny, T., & De, A. (2023). Understanding the liquid jet break-up in various regimes at elevated pressure using a compressible VOF-LPT coupled framework. *International Journal of Multiphase Flow*, **159**, 1–55. doi:10.1016/j.ijmultiphaseflow.2022.104303.

[11] Blanchard, S. (2021). *Multi-physics large-eddy simulation of methane oxy-combustion in liquid rocket engines*. Phd thesis Ecole doctorale MEGEP - Dynamique des fluides - CERFACS.
[12] Blanchard, S., & Cuenot, B. (2019). Large eddy simulation of a single-injector lox/gch4 com-
bustion chamber using analytically reduced chemistry. In *DFG-TRR-40 Summer Program*
*2019* 1 (pp. 69–86).

[13] Boussinesq, J. (1903). *Théorie analytique de la chaleur: mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière* volume 2. Gauthier-Villars.

[14] Brackbill, J., Kothe, D., & Zemach, C. (1992). A continuum method for modeling surface tension. *Journal of Computational Physics*, *100*, 335–354. URL: https://doi.org/10.1016/0021-9991(92)90240-y. doi:10.1016/0021-9991(92)90240-y.

[15] Broumand, M., & Birouk, M. (2016). Liquid jet in a subsonic gaseous crossflow: Recent progress and remaining challenges. *Progress in Energy and Combustion Science*, *57*, 1–29.

[16] Carmona, J., Treleaven, N., Odier, N., & Cuenot, B. (2023). A Lagrangian simulation methodology for Large Eddy Simulations of prefilming airblast injectors. *Journal of Propulsion and Power*.

[17] Chéron, V., Brändle de Motta, J. C., Ménard, T., Poux, A., & Berlemont, A. (2023). A coupled eulerian interface capturing and lagrangian particle method for multiscale simulation. *Computers & Fluids*, *256*, 105843.

[18] Colin, O., & Rudgyard, M. (2000). Development of high-order taylor-galerkin schemes for LES. *Journal of Computational Physics*, *162*, 338–371. URL: https://doi.org/10.1006/jcph.2000.6538. doi:10.1006/jcph.2000.6538.

[19] Daviller, G., Oztarlik, G., & Poinsot, T. (2019). A generalized non-reflecting inlet boundary condition for steady and forced compressible flows with injection of vortical and acoustic waves. *Computers & Fluids*, *190*, 503–513. URL: https://doi.org/10.1016/j.compfluid.2019.06.027. doi:10.1016/j.compfluid.2019.06.027.

[20] descamps, M., Matas, J.-P., & Cartellier, A. (2008). Gas-liquid atomisation: gas phase char-
acteristics by piv measurements and spatial evolution of the spray, .

[21] Domingo-Alvarez, P., Bénard, P., Moureau, V., Lartigue, G., & Grisch, F. (2020). Impact of spray droplet distribution on the performances of a kerosene lean/premixed injector. *Flow, Turbulence and Combustion*, *104*, 421–450.
[22] Dumouchel, C. (2008). On the experimental investigation on primary atomization of liquid streams. *Experiments in Fluids*, **45**, 371–422.

[23] Eckel, G., Rachner, M., Clercq, P. L., & Aigner, M. (2016). Semi-empirical model for the unsteady shear breakup of liquid jets in cross-flow. *Atomization and Sprays*, **26**, 687–712. URL: https://doi.org/10.1615/atomizspr.2015011185. doi:10.1615/atomizspr.2015011185.

[24] Eggers, J., & Villermaux, E. (2008). Physics of liquid jets. *Reports on progress in physics*, **71**, 036601.

[25] Elghobashi, S., & Truesdell, G. C. (1992). Direct simulation of particle dispersion in a decaying isotropic turbulence. *Journal of Fluid Mechanics*, **242**, 655–700. URL: https://doi.org/10.1017/s0022112092002532. doi:10.1017/s0022112092002532.

[26] Eslamian, M., Amighi, A., & Ashgriz, N. (2014). Atomization of liquid jet in high-pressure and high-temperature subsonic crossflow. *AIAA Journal*, **52**, 1374–1385.

[27] Evrard, F., Denner, F., & van Wachem, B. (2019). A multi-scale approach to simulate atomization processes. *International Journal of Multiphase Flow*, **119**, 194–216.

[28] Farago, Z., & Chigier, N. (1992). Morphological classification of disintegration of round liquid jets in a coaxial air stream. *Atomization and Sprays*, **2**, 137–153.

[29] Freitag, S., & Hassa, C. (2008). Spray characteristics of a kerosene jet in cross flow of air at elevated pressure. In *Proceedings ILASS-Europe* (pp. 8–10).

[30] Fuster, D., Matas, J.-P., Marty, S., Popinet, S., Hoepffner, J., Cartellier, A., & Zaleski, S. (2013). Instability regimes in the primary breakup region of planar coflowing sheets. *Journal of Fluid Mechanics*, **736**, 150–176.

[31] Gicquel, L. Y., Gourdain, N., Boussuge, J.-F., Deniau, H., Staffelbach, G., Wolf, P., & Poinsot, T. (2011). High performance parallel computing of flows in complex geometries. *Comptes Rendus Mécanique*, **339**, 104–124. URL: https://doi.org/10.1016/j.crme.2010.11.006. doi:10.1016/j.crme.2010.11.006.

[32] GIOUD, T. (2024). *Large Eddy Simulations for atomization in subcritical liquid rocket engines*. Ph.D. thesis Université de Toulouse.
[33] Gioud, T., Odier, N., Cuenot, B., Schmitt, T., Saucereau, D., & Martin-Benito, M. (2022). Injection modelling in LOX/GCH4 rocket engines with a diffuse interface method. In *Space Propulsion Conference - 09-13 May*.

[34] Gorokhovski, M. (2001). The stochastic lagrangian model of drop breakup in the computation of liquid sprays. *Atomization and Sprays, 11*, 15. URL: https://doi.org/10.1615/atomizspr.v11.i5.20. doi:10.1615/atomizspr.v11.i5.20.

[35] Guildenbecher, D. R., López-Rivera, C., & Sojka, P. E. (2009). Secondary atomization. *Experiments in Fluids*, *46*, 371-402.

[36] Haselbacher, A., Najjar, F., & Ferry, J. (2007). An efficient and robust particle-localization algorithm for unstructured grids. *Journal of Computational Physics*, *225*, 2198-2213.

[37] Haselbacher, A., Najjar, F., & Ferry, J. P. (2007). An efficient and robust particle-localization algorithm for unstructured grids. *Journal of Computational Physics*, *225*, 2198-2213.

[38] Heinrich, M., & Schwarze, R. (2020). 3d-coupling of volume-of-fluid and lagrangian particle tracking for spray atomization simulation in openfoam. *SoftwareX*, *11*, 100483.

[39] Herrmann, M. (2010). Detailed Numerical Simulations of the Primary Atomization of a Turbulent Liquid Jet in Crossflow. *Journal of Engineering for Gas Turbines and Power*, *132*, 061506.

[40] Hirt, C., & Nichols, B. (1981). Volume of fluid (vof) method for the dynamics of free boundaries. *Journal of Computational Physics*, *39*, 201-225. doi:https://doi.org/10.1016/0021-9991(81)90145-5.

[41] Jaegle, F., Senoner, J.-M., García, M., Bismes, F., Lecourt, R., Cuenot, B., & Poinsot, T. (2011). Eulerian and lagrangian spray simulations of an aeronautical multipoint injector. *Proceedings of the Combustion Institute*, *33*, 2099-2107.

[42] Janodet, R., Guillamón, C., Moureau, V., Mercier, R., Lartigue, G., Bénard, P., Ménard, T., & Berlemont, A. (2022). A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids. *Journal of Computational Physics*, *458*. URL: https://hal.science/hal-03024186.
[43] Laroche, T., Odier, N., Schmitt, T., Pelletier, M., & Cuenot, B. (2021). A diffuse interface method with real-gas thermodynamic equilibrium closure applied to capillary problems. *International Conference on Liquid Atomization and Spray Systems (ICLASS)*, 1. URL: https://doi.org/10.2218/iclass.2021.5996. doi:10.2218/iclass.2021.5996.

[44] Lasheras, J. C., & Hopfinger, E. J. (2000). Liquid jet instability and atomization in a coaxial gas stream. *Annual Review of Fluid Mechanics*, 32, 275-308.

[45] LASHERAS, J. C., VILLERMAUX, E., & HOPFINGER, E. J. (1998). Break-up and atomization of a round water jet by a high-speed annular air jet. *Journal of Fluid Mechanics*, 357, 351-379. doi:10.1017/S0022112097008070.

[46] Le Touze, C., Dorey, L.-H., Rutard, N., & Murrone, A. (2020). A compressible two-phase flow framework for large eddy simulations of liquid-propellant rocket engines. *Applied Mathematical Modelling*, 84, 265-286. URL: https://doi.org/10.1016/j.apm.2020.03.028. doi:10.1016/j.apm.2020.03.028.

[47] Lefebvre, A., & McDonell, V. (2017). *Atomization and Sprays (2nd ed.)*. CRC Press.

[48] Li, X., & Soteriou, M. C. (2018). Detailed numerical simulation of liquid jet atomization in crossflow of increasing density. *International Journal of Multiphase Flow*, 104, 214-232. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2018.02.016. doi:10.1016/j.ijmultiphaseflow.2018.02.016.

[49] Ling, Y., Zaleski, S., & Scardovelli, R. (2015). Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model. *International Journal of Multiphase Flow*, 76, 122-143.

[50] Liu, C., Xi, X., Liu, H., Li, W., Jia, M., & Li, R. (2024). Multi-Scale methodology of breakup and atomization for liquid jets in crossflow. *Applied Thermal Engineering*, 248. doi:10.1016/j.applthermaleng.2024.123309.

[51] Lubarsky, E., Shcherbik, D., Bibik, O., Gopala, Y., & Zinn, B. (2012). Fuel jet in cross flow—experimental study of spray characteristics. *Advanced Fluid Dynamics*, 59.

[52] Luo, K., Shao, C., Chai, M., & Fan, J. (2019). Level set method for atomization and evaporation simulations. *Progress in Energy and Combustion Science*, 73, 65-94.
[53] Matas, J.-P., Marty, S., & Cartellier, A. (2011). Experimental and analytical study of the shear instability of a gas-liquid mixing layer. *Physics of Fluids, 23*.

[54] Matheis, J., & Hickel, S. (2018). Multi-component vapor-liquid equilibrium model for LES of high-pressure fuel injection and application to ECN Spray A. *International Journal of Multi-phase Flow, 99*, 294–311. URL: https://doi.org/10.1016/j.ijmultiphaseflow.2017.11.001. doi:10.1016/j.ijmultiphaseflow.2017.11.001.

[55] Mishra, D. (2017). *Fundamentals of rocket propulsion*. CRC Press.

[56] Nicoud, F., & Ducros, F. (1999). Subgrid-scale stress modelling based on the square of the velocity gradient tensor. *Flow, Turbulence and Combustion, 62*, 183–200. URL: https://doi.org/10.1023/a:1009995426001. doi:10.1023/a:1009995426001.

[57] Oseen, C. W. (1927). *Neuere methoden und ergebnisse in der hydrodynamik* volume 1. Akademische Verlagsgesellschaft MBH.

[58] Pairetti, C., Villiers, R., & Zaleski, S. (2021). On shear layer atomization within closed channels: Numerical simulations of a cough-replicating experiment. *Computers & Fluids, 231*, 105125.

[59] Patel, N., Kirtaş, M., Sankaran, V., & Menon, S. (2007). Simulation of spray combustion in a lean-direct injection combustor. *Proceedings of the Combustion Institute, 31*, 2327–2334.

[60] Péden, B., Carmona, J., Boivin, P., Schmitt, T., Benedicte, C., & Odier, N. (2023). Numerical assessment of diffuse-interface method for air-assisted liquid sheet simulation. *Computers & Fluids, 266*, 106022.

[61] Pelletier, M., Schmitt, T., & Ducruix, S. (2020). A multifluid taylor-galerkin methodology for the simulation of compressible multicomponent separate two-phase flows from subcritical to supercritical states. *Computers & Fluids, 206*, 104588. URL: https://doi.org/10.1016/j.compfluid.2020.104588. doi:10.1016/j.compfluid.2020.104588.

[62] Peng, D.-Y., & Robinson, D. B. (1976). A new two-constant equation of state. *Industrial & Engineering Chemistry Fundamentals, 15*, 59–64. URL: https://doi.org/10.1021/i160057a011. doi:10.1021/i160057a011.
[63] Perigaud, G., & Saurel, R. (2005). A compressible flow model with capillary effects. *Journal of Computational Physics*, **209**, 139–178. URL: https://doi.org/10.1016/j.jcp.2005.03.018. doi:10.1016/j.jcp.2005.03.018.

[64] Pilch, M., & Erdman, C. (1987). Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop. *International journal of multiphase flow*, **13**, 741–757.

[65] Poinsot, T., & Lele, S. (1992). Boundary conditions for direct simulations of compressible viscous flows. *Journal of Computational Physics*, **101**, 104–129. URL: https://doi.org/10.1016/0021-9991(92)90046-2. doi:10.1016/0021-9991(92)90046-2.

[66] Poling, B., Prausnitz, J., & O'Connell, J. (2001). *The Properties of Gases and Liquids*. McGraw-Hill Education.

[67] Popinet, S. (2009). An accurate adaptive solver for surface-tension-driven interfacial flows. *Journal of Computational Physics*, **228**, 5838–5866.

[68] Potier, L., Duchaine, F., Cuenot, B., Saucereau, D., & Pichillou, J. (2022). Prediction of wall heat fluxes in a rocket engine with conjugate heat transfer based on large-eddy simulation. *Entropy*, **24**.

[69] Ragucci, R., Bellofiore, A., & Cavaliere, A. (2007). Breakup and breakdown of bent kerosene jets in gas turbine conditions. *Proceedings of the Combustion Institute*, **31**, 2231–2238. URL: https://doi.org/10.1016/j.proci.2006.07.204. doi:10.1016/j.proci.2006.07.204.

[70] Rochette, B., Riber, E., Cuenot, B., & Vermorel, O. (2020). A generic and self-adapting method for flame detection and thickening in the thickened flame model. *Combustion and Flame*, **212**, 448–458. URL: https://doi.org/10.1016/j.combustflame.2019.11.015. doi:10.1016/j.combustflame.2019.11.015.

[71] Saurel, R., & Abgrall, R. (1999). A multiphase godunov method for compressible multifluid and multiphase flows. *Journal of Computational Physics*, **150**, 425–467. URL: https://doi.org/10.1006/jcph.1999.6187. doi:10.1006/jcph.1999.6187.
[72] Saurel, R., & Pantano, C. (2018). Diffuse-interface capturing methods for compressible two-phase flows. *Annual Review of Fluid Mechanics*, *50*, 105–130. URL: https://doi.org/10.1146/annurev-fluid-122316-050109. doi:10.1146/annurev-fluid-122316-050109.

[73] Schiavo, E. L., Laera, D., Riber, E., Gicquel, L., & Poinsot, T. (2021). On the impact of fuel injection angle in euler-lagrange large eddy simulations of swirling spray flames exhibiting thermoacoustic instabilities. *Combustion and Flame*, *227*, 359–370.

[74] Schiller, L., & Nauman, A. (1935). A drag coefficient correlation. *VDI Zeitung*, *77*, 318–320.

[75] Schmitt, T. (2020). Large-eddy simulations of the mascotte test cases operating at supercritical pressure. *Flow, Turbulence and Combustion*, *105*, 159–189.

[76] Schmitt, T., & Ducruix, S. (2023). Evaluation of large-eddy simulation coupled with an homogeneous equilibrium model for the prediction of coaxial cryogenic flames under subcritical conditions. *Aerospace*, *10*, 98. URL: https://doi.org/10.3390/aerospace10020098. doi:10.3390/aerospace10020098.

[77] Schönfeld, T., & Rudgyard, M. (1999). Steady and unsteady flow simulations using the hybrid flow solver AVBP. *AIAA Journal*, *37*, 1378–1385. URL: https://doi.org/10.2514/2.636. doi:10.2514/2.636.

[78] Senoner, J.-M. (2010). *Simulation aux grandes échelles de l' écoulement diphasique évaporant dans un brûleur aéronautique par une approche Euler-Lagrange*. Ph.D. thesis Institut National Polytechnique de Toulouse - Dynamique des Fluides.

[79] Shinjo, J., & Umemura, A. (2010). Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation. *International Journal of Multiphase Flow*, *36*, 513–532.

[80] Song, J., Cain, C. C., & Lee, J. G. (2014). Liquid jets in subsonic air crossflow at elevated pressure. *Journal of Engineering for Gas Turbines and Power*, *137*. URL: https://doi.org/10.1115/1.4028565. doi:10.1115/1.4028565.

[81] Spalding, D. (1953). The combustion of liquid fuels. *Symposium (International) on Combustion*, 4, 847–864. URL: https://doi.org/10.1016/s0082-0784(53)80110-4. doi:10.1016/s0082-0784(53)80110-4.
[82] Sun, Y., Li, Y., Dreßler, L., Nishad, K., & Sadiki, A. (2024). Multiscale numerical modeling of a complete spray evolution including breakup of liquid jet injection in gaseous cross flow. *International Journal of Multiphase Flow*, 170, doi:10.1016/j.ijmultiphaseflow.2023.104655.

[83] Sussman, M., & Puckett, E. G. (2000). A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. *Journal of Computational Physics*, 162, 301–337. doi:10.1006/jcph.2000.6537.

[84] Sussman, M., Smereka, P., & Osher, S. (1994). A level set approach for computing solutions to incompressible two-phase flow. *Journal of Computational Physics*, 114, 146–159. doi:https://doi.org/10.1006/jcph.1994.1155.

[85] Treleaven, N., Laera, D., Carmona, J., Odier, N., Gentil, Y., Dombard, J., Daviller, G., Gicquel, L., & Poinsot, T. (2022). Coupling of combustion simulation with atomisation and filming models for les in swirled spray flames. *Flow Turbulence and Combustion*, 109, 759–789.

[86] Tudisco, P., & Menon, S. (2019). Numerical investigations of phase-separation during multi-component mixing at super-critical conditions. *Flow, Turbulence and Combus- tion*, 104, 693–724. URL: https://doi.org/10.1007/s10494-019-00101-4. doi:10.1007/s10494-019-00101-4.

[87] Villermaux, E. (2007). Fragmentation. *Annu. Rev. Fluid Mech.*, 39, 419–446.

[88] Wen, J., Hu, Y., Nakanishi, A., & Kurose, R. (2020). Atomization and evaporation process of liquid fuel jets in crossflows: A numerical study using Eulerian/Lagrangian method. *International Journal of Multiphase Flow*, 129, 103331.

[89] Wu, P.-K., Kirkendall, K. A., Fuller, R. P., & Nejad, A. S. (1997). Breakup processes of liquid jets in subsonic crossflows. *Journal of Propulsion and Power*, 13, 64–73.

[90] Yu, H., Goldsworthy, L., Ghiji, M., Brandner, P., & Garaniya, V. (2017). A parallel volume of fluid-lagrangian parcel tracking coupling procedure for diesel spray modelling. *Computers & Fluids*, 150, 46–65.

[91] Zhou, W., Chen, B., Zhu, Q., Rao, S., & Xu, X. (2023). Numerical simulation of angled-injected liquid jet breakup in supersonic crossflow by a hybrid VOF-LPT method. *International Journal of Multiphase Flow*, 166.
[92] Zuzio, D., Estivalèzes, J.-L., & DiPierro, B. (2016). An improved multiscale eulerian–lagrangian method for simulation of atomization process. *Computers & Fluids*, **176**, 285–301.