# Direct numerical simulation of a subcritical coaxial injection in fiber regime using sharp interface reconstruction

Jean-Christophe Hoarau, Luc-Henry Dorey, Davide Zuzio, Florian Granger, Jean-Luc Estivalezes

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

## ► To cite this version:

Jean-Christophe Hoarau, Luc-Henry Dorey, Davide Zuzio, Florian Granger, Jean-Luc Estivalezes. Direct numerical simulation of a subcritical coaxial injection in fiber regime using sharp interface reconstruction. International Journal of Multiphase Flow, 2024, 180, pp.104974. 10.1016/j.ijmultiphaseflow.2024.104974. hal-04212765v3

HAL Id: hal-04212765

https://hal.science/hal-04212765v3

Submitted on 2 Sep 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
# Direct numerical simulation of a subcritical coaxial injection in fiber regime using sharp interface reconstruction.

J.-C. Hoarau<sup>a</sup>, L.-H. Dorey<sup>a,*</sup>, D. Zuzio<sup>b</sup>, F. Granger<sup>a</sup>, J.-L. Estivalezès<sup>b</sup>

## Abstract

The numerical simulation of space launchers combustion chambers is a topic of increasing interest, as it may help the development of safer and more efficient designs. Understanding fuel injection is a particularly severe challenge. The liquid oxygen is injected by a round orifice surrounded by an annular gaseous stream of fuel, leading in subcritical conditions to a two-phase assisted atomization process. The result of this process is a very dense and polydisperse two-phase flow, which strongly influences the behavior of the chamber. Experimental investigation of this flow is difficult due to the axisymmetric geometry and the dense characteristic of the spray. Neither RANS nor Large Eddy Simulation (LES) possess reliable models able to reproduce the smallest scales of atomization, one of the reasons being the lack of relevant experimental data. Therefore, this work aims to provide detailed information on the atomization process using Direct Numerical Simulation. This paper presents a Direct Numerical Simulation (DNS) of a coaxial liquid-gas assisted atomization in the typical fiber regime encountered in cryogenic injectors. This study aims to better understand the evolution of liquid topology and extract relevant information that may help develop larger-scale models. The flow was first analyzed in terms of topology statistical data, using a dedicated detection and classification algorithm that could characterize the individual liquid structures. These include the central liquid core, the ligament created during primary atomization, and the spherical droplet obtained at the end of the atomization process. Subsequently, a more global statistical topology indicator was investigated, namely the interface area density distribution. This quantity is used in larger-scale RANS or LES models to predict the smallest scales of atomization. Therefore, understanding its behavior in a realistic case is of utmost importance. The interface area density distribution was correlated to the global jet behavior and the liquid topology data obtained by the detection algorithm. The results showed, in particular, a strong correlation between the initial increase of liquid-gas interface area density with the generation of ligaments and between the continuous decrease of the interface area density during droplet formation and stabilization.

**Keywords:** Primary atomization, coaxial injector, liquid rocket engine, DNS, textural atomization

## 1. Introduction

Liquid rocket engines rely on high-energy assisted atomization to produce the best oxidant-fuel mixture for optimal combustion. When sufficiently low temperature and pressure conditions are encountered in the chamber, typical injection conditions are liquid oxygen and gaseous fuel. These conditions are called *subcritical* and are typically encountered during engine ignition, for permanent regimes of low-thrust engines, and during low-thrust phases of more powerful throttleable engines. In these conditions, the typical coaxial injection

produces a round liquid oxygen jet surrounded by an annular, high-speed gaseous fuel (hydrogen or methane) flow. The interaction between these flows triggers two-phase instabilities, which lead to liquid breakup and subsequent atomization.

Two main stages can be distinguished in the assisted atomization process [34]. The first occurs near the orifice, where the liquid jet becomes subject to longitudinal instabilities [49, 40, 41, 42], which are the results of the shearing effect triggered by the interaction between the liquid and gas flows. The instabilities involved in this phase are the Kelvin-Helmholtz, induced by the interfacial shearing, and the Rayleigh-Taylor, caused by the heavy liquid acceleration [42]. Both sheared liquid

<footnote><sup>a</sup>DMPE, ONERA, Université Paris-Saclay, 91120 Palaiseau - France</footnote>

<footnote><sup>b</sup>DMPE, ONERA, Université de Toulouse, 31000 Toulouse - France</footnote>

<footnote>*Corresponding author
Email address: luc-henry.dorey@onera.fr (L.-H. Dorey)</footnote>
Nomenclature

<table><tr><td>(●)</td><td>Temporal average</td><td>ρ<sub>l</sub></td><td>Liquid density</td></tr><tr><td>⟨●⟩<sub>θ</sub></td><td>Azimuthal average</td><td>σ</td><td>Surface tension</td></tr><tr><td>J</td><td>Momentum Ratio</td><td>Σ</td><td>Liquid-gas interface area density</td></tr><tr><td>Re<sub>g</sub></td><td>Gas Reynolds number</td><td>VOF</td><td>Volume of fluid</td></tr><tr><td>We<sub>g</sub></td><td>Gas Weber number</td><td>φ</td><td>Level-set</td></tr><tr><td>Re<sub>l</sub></td><td>Liquid Reynolds number</td><td>δ<sub>g</sub></td><td>Vorticity length</td></tr><tr><td>We<sub>l</sub></td><td>Liquid Weber number</td><td>d<sub>30</sub></td><td>Liquid structure characteristic diameter</td></tr><tr><td>U<sub>g</sub></td><td>Gas injection Velocity</td><td>δ<sub>lip</sub></td><td>Lip length</td></tr><tr><td>U<sub>l</sub></td><td>Liquid injection Velocity</td><td>D<sub>l</sub></td><td>Liquid injector diameter</td></tr><tr><td>ρ<sub>g</sub></td><td>Gas density</td><td>h<sub>g</sub></td><td>Gas injection gap</td></tr></table>

round jets at low Weber numbers and sheets main bodies present a similar flapping motion in this zone, as reported by [24], [43] and [28]. This zone usually extends up to several diameters from the injector, where the main breakup point is reached. The liquid forms a population of individual structures, in particular elongated streamwise ligaments. All these phenomena make up primary atomization. The secondary atomization occurs when the released liquid structures are, in turn, affected by the shearing gas flow, in a continuous fragmentation in smaller and smaller structures up to the final spray, under the Rayleigh-Plateau instability [30, 9, 16]. Coaxial jet atomization regimes can be classified by function of the aerodynamic Weber number. Farago and Chigier [14] described five regimes: axisymmetric Rayleigh breakup, non-axisymmetric Rayleigh breakup, membrane breakup, fiber breakup, and superpulsating breakup, in function of the gaseous Weber number and the liquid Reynolds number. In [31], the momentum flux ratio $M$ was added as a relevant non-dimensional number to be considered, in addition to Weber and Reynolds numbers. More recently, four breakup regimes were distinguished by [36]: first-wind induced breakup regime ($We_a \approx 9$), second-wind induced breakup regime ($We_a \approx 50$), bag breakup regime ($We_a \approx 150$) and the fiber-type atomization ($We_a \approx 500$). This last regime is commonly used in liquid rocket engines to obtain the fastest oxygen atomization. In this regime, the liquid core stretches under the gaseous stresses, periodically releasing large liquid blobs in the stream. The instabilities growing on the liquid core develop in many streamwise thin ligaments,

which subsequently undergo Rayleigh-Plateau breakup. This mechanism generates a complex two-phase flow within a few jet diameters from the injection point, in which the liquid core coexists with large liquid structures, ligaments, and initial droplets.

The assisted atomization is believed to play a significant role in the global behavior of the chamber. In particular, it may participate, together with the thermal and acoustic behavior of the thrust chamber, in the development of combustion instabilities. The prediction and control of combustion instabilities is still a very challenging task, as several complex mechanisms interact. The experimental approach is a meaningful way to study the instability onset and development better. However, full-scale engine tests are very expensive, and the measurement systems face severe restrictions in the operating conditions of a combustion chamber. Atomization poses the particular challenge of a dense liquid-gas zone investigation: even with available optical access, the liquid topology is hard to characterize near the injector. Peripheral droplet populations can be detected and measured, but little information about their formation mechanisms exists.

Direct Numerical Simulation (DNS) can provide helpful information about primary atomization. Appropriate numerical methods have been developed in the last decades to simulate complex two-phase flows involving atomization. The Level-Set approach and its variant, the Conservative Level-Set (CLS), were employed at first to simulate low density ratio diesel jets and liquid jet in cross-flow (LJICF) [10, 23]. The Volume-of-Fluid (VOF) approach was used to simulate
the break-up of a liquid jet by a high-speed coaxial gas jet [54, 37]. The two approaches were coupled by [52] in the Coupled Level-Set/Volume of Fluid method, and used to simulate turbulent jet atomization [44, 51, 35]. CLS and CLSVOF were used to simulate the primary atomization of a gas-sheared liquid sheet by [46] and [61], albeit with reduced density ratios. In [1, 38, 26] DNS of air-assisted thick film atomization were performed. The results were analyzed regarding interfacial wave development, turbulence and two-phase mixing layer development, and two-phase flow turbulence.

Numerical simulation of a coaxial air-water jet was performed in [58] using a CLSVOF approach, making use of a modified extrapolated liquid velocity field to ensure proper liquid kinetic energy preservation and a LES (Large Eddy Simulation) approach for both liquid and gas flow. The considered air-water configuration at ambient pressure was taken from [6]. The OpenFOAM VOF approach was used in [45] to perform simulations of the primary breakup of a high-viscosity glycerol-water mixture using an external mixing twin fluid nozzle. Compressible LES were performed to capture the morphology of the primary breakup and the important flow field characteristics. The VOF approach was used to track the unsteady evolution and breakup of the liquid jet. This work reported good agreement concerning the liquid jet's breakup frequency, velocity fields, and morphology. However, the low-Weber regime was unfavorable for atomization, and the liquid jet showed mostly large-scale deformations induced by pulsating and flapping instabilities. In [60], an extensive study of the primary atomization behavior of high-viscous liquid (glycerol-water mixture) jets was performed again with the OpenFOAM VOF approach. Different gas pressures were investigated. The simulations showed a pulsating mode instability of liquid jet, along with a fiber-type disintegration. This study focused on the very dense zone of primary atomization in terms of liquid core length and kinetic energy distribution. However, the analysis did not extend to the atomized part of the liquid jet. Moreover, in that case, the viscosity ratio $\mu_l/\mu_g$ is very high (around 5500), giving a liquid Reynolds number of the order of 35, far from injection conditions of liquid rocket engines. It should be noted that none of these LES approaches for atomization simulations took into account subgrid-scale modeling of surface tension forces. Indeed, as they work with a one-fluid approach for two-phase flow, this LES modeling must be understood as dealing only with the velocity field.

Besides interface capturing methods described above, other so-called multi-fluid approaches are often used to simulate multiphase flows in compressible flows, con-

trary to the previous ones used primarily for incompressible flows. Multi-fluid approaches can be classified in the family of diffuse interface methods since the interface is no longer captured or tracked as a discontinuity between two non miscible phases, but rather represented in each computational cell by the volume fraction of each phase, which is transported by the local velocity field. It should be noted that although this method is less accurate for the same mesh size than an interface capturing method, it is easily generalized to the case of unstructured meshes. Indeed, this method was successfully employed by [4, 5] to catch the primary atomization of an airblast atomizer at room conditions in the framework of large-scale LES-like simulations.

To simulate, in an industrial context, the two-phase flow from the injection system to the combustion of the final spray, it is not realistic to envisage direct simulations. On the other hand, a multi-fluid approach for the primary atomization coupled with a dispersed phase approach for the secondary atomization seems to be a good candidate. Using this strategy, full multi-scale LES simulations of a subcritical cryogenic LOX-GH2 fiber-regime jet were performed in [19] and [32] by explicitly coupling a multi-fluid diffuse interface method (for the dense zone or primary atomization zone) with an Eulerian dispersed phase solver (for the remaining dispersed zone). However, in the transition zone between dense and dispersed ones, the size and velocity of droplets were estimated as a function of the injected propellant properties and imposed as a boundary condition and not deduced from the dense zone. One way to overcome this shortcoming is to use an ELSA-type model initially developed in a RANS (Reynolds-averaged Navier-Stokes) context [55, 56] to apply it here in an LES context. As in turbulence RANS modeling, an additional quantity is followed in time, namely the liquid-gas interface area per unit of volume, or interface area density. This quantity can statistically represent the dispersion of the liquid in the local volume of control. The general form of the corresponding equation is a transport-diffusion, with multiple source terms representing production by turbulent mixing between liquid and gas, mean shear stress, break-up, and destruction by coalescence or evaporation. Despite the several propositions of source terms closure in RANS [25, 33, 12] and LES formalism [7, 2, 47], it remains challenging to understand the actual behavior of this quantity: the validation is often realized by comparing the final outcome of atomization, i.e. the final droplet size, with experiments. Indeed, especially for the coaxial jet configuration, the experimental investigation is strongly limited by the optical density of the produced
cloud of liquid structures, mainly along the jet cen- terline: superposition of the liquid structures hampers in this region the visual detection techniques so that the only available data involve the more dispersed tiny droplets population in the jet periphery.

That's why the present study aims to bring new re- sults to help understand the atomization phenomena and validate interface area density modeling for coaxial at- omization in fiber-type regime. Beyond this paper, we want to improve the accuracy of the numerical method- ology proposed by Le Touze et al. [32] to simulate cryo- genic flames of rocket engines. The approach is to im- prove atomization modeling with the help of interface area density transport and a source term specific to un- steady coaxial injection. To formulate this new model and validate interface area density estimation, we real- ized in this paper a DNS of the primary atomization in the fiber-type regime, that brings new insights into this phenomena along with interface area density reference results for validation.

## 2. Numerical tools

ONERA has been developing its in-house two-phase flow DNS solver, DYJEAT, since 2004 [8]. DYJEAT was successfully employed and validated in many com- plex two-phase flow simulations. In [61] both single- droplet behavior and a practical case of assisted atom- ization were investigated. An air-sheared liquid sheet was simulated, and the obtained droplets average diam- eter was found close to the measured one from the ex- perience of [13]. The work in [62] showed good capa- bilities of the code to reproduce the behavior of a rising bubble in water, as well as the deformation of a liquid jet in cross-flow with two different air-liquid momentum ratios. In [57] numerical simulations of liquid drops im- pacts on both dry and wet walls were successfully com- pared to experimental literature results. Finally, in the work of [29] the behavior of a water droplet impact- ing a rigid fiber was investigated, and shown to be very similar to the results from the experiences from [27], and this for different regimes ("splitting", "single drop", "capturing").

### 2.1. Physical model

DYJEAT model involves the direct numerical resolu- tion of the incompressible Navier-Stokes equations. A one-fluid sharp interface method is used to simulate two non-miscible phases. The interface is tracked in space and time to define the local properties of the fluid by a Level-Set function, coupled to a Volume-of-Fluid ap- proach to ensure mass conservation. The liquid and gas

properties will be referred to with the subscripts .l and .g, respectively. The two-phase incompressible Navier- Stokes equations read:

$$ \frac{\partial \phi}{\partial t} + U \nabla(\phi) = 0 \qquad (1) $$

$$ \frac{\partial C}{\partial t} + \nabla \cdot (UC) = -C \nabla \cdot U \qquad (2) $$

$$ \nabla \cdot U = 0 \qquad (3) $$

$$ \frac{\partial \rho U}{\partial t} + \nabla \cdot (\rho U \otimes U) = f + \nabla \cdot T \qquad (4) $$

where U is the velocity vector, ρ the density and f the external body forces such as gravity. φ is the Level-Set function, C the liquid volume fraction in each cell. The stress tensor is modeled as follows:

$$ T = -pI + D \qquad (5) $$

$$ D = \mu (\nabla U + (\nabla U)^T) \qquad (6) $$

Here p represents the hydrodynamic pressure, μ the dynamic viscosity, and D the viscous stress tensor. Jump conditions across the interface are set as follows: The normal velocity is conserved across the interface $U_l \cdot n = U_g \cdot n$; the surface tension effect and velocity no-slip condition balance viscous and pressure effects:

$$ [p] - n \cdot [2\nabla U] \cdot n = \sigma\kappa \qquad (7) $$

$$ t \cdot [2\nabla U] \cdot n = 0 \qquad (8) $$

where σ is the surface tension, κ is the local curvature, t and n are the tangential and normal unit vectors oriented from gas to liquid. The symbolic [.] represents the jump across the interface: $(.)_l - (.)_g$.

### 2.2. Numerical method

Level-Set and Volume-of-Fluid are widely spread in the literature for two-phase flow simulations. The Level-Set function allows for accurate computation of interface properties at the expense of mass conservation issues, whereas the VOF method ensures mass conser- vation. For these reasons, the Coupled Level Set-VOF method has been developed by [52] and [44]. This ap- proach uses both Level-Set and VOF functions to track the liquid, providing an accurate and reliable method to simulate two-phase flows. The Level-Set function, φ, represents the signed distance between any points in the domain and the interface. It is a scalar field driven by the scalar transport equation (1). The zero-level curve of φ provides the interface location. Interface properties
are calculated by :

$$n = \frac{\nabla \phi}{\|\nabla \phi\|} \qquad (9)$$

$$\kappa = \nabla \cdot \left( \frac{\nabla \phi}{\|\nabla \phi\|} \right) \qquad (10)$$

The Level-Set loses its properties of distance function during its advection. That loss is minimized by the redistance algorithm [53] which updates the Level-Set fields till the $||\nabla\phi|| = 1$ is discretely satisfied while keeping the zero level curve untouched. A reinitialization step is performed at each time step.

The color function C in the VOF method is defined for an incompressible flow as the volume fraction of liquid in each cell as:

$$C(x_i) = \frac{V_{li}}{V_i} \qquad (11)$$

$$V_i = V_{li} + V_{gi} \qquad (12)$$

with i the cell index and $V_i$ its volume. The color function is a scalar field driven by a conservation equation reading:

$$\frac{\partial C}{\partial t} + \nabla \cdot (CU) = 0 \qquad (13)$$

The interface is approximated in each cell by a plane, the so-called PLIC reconstruction from [59]. The normal vector components are inferred from the surrounding Level-Set values, whereas the plane itself is shifted to ensure that the bounded volume corresponds to the VOF value. The VOF PLIC allows for geometrical flux computation [52], requiring a directional splitting algorithm. Once the VOF is advected, the Level-Set function is locally shifted to be consistent with the local VOF value, thus ensuring mass conservation. The freeware VofTools library from [39] has been used to deal with the geometrical computations involving interface reconstruction.

Momentum conservation considerations led the authors to develop the corresponding consistent algorithm, called Conservative Momentum scheme, in [62]. The CMOM approach involves using the PLIC fluxes to compute momentum fluxes and solve a conservative form of the momentum equation. Using the mass fluxes ($\rho U$) in both mass and momentum equation ensures consistency between mass and momentum and global preservation of each phase energy. The resulting scheme dramatically improves the solver's robustness and accuracy when dealing with high flows at density ratio and significant phase velocity differences. The

Ghost-Fluid method [15] allows subgrid accurate computation of the pressure and stress tensor jumps. The continuous solution is extended on the two sides of the interface, and a corresponding source term is added in the numerical discretization.

## 2.3. Liquid structures detection

As described above, the massive dispersion of various multi-scale liquid objects in atomization processes makes studying their distribution very complex. To obtain statistics on the distribution, shape and size of the liquid structures generated in a spray, a detection algorithm was developed in the DYJEAT platform [3]. This new tool, called AlgoDetect, uses an algorithm inspired by the connected component labeling (CCL) method [23, 20]. Similar to the method proposed in [22], the AlgoDetect methodology is improved by the use of PLIC reconstruction to distinguish coherent structures separated by a distance of the order of one grid cell. The main difference with the CCL algorithm is that the AlgoDetect algorithm does not use dynamic lists containing the cells of coherent structures. A summary of the algorithm is given in the following paragraph.

Firstly, each cell containing liquid is labeled with a unique tag. Then, using VOF and PLIC reconstruction in these cells, the adjacent cells connected to the same coherent structures are identified. Finally, for each liquid-containing cell, the minimum tag of all connected adjacent cells is retained to assign the same tag in all cells of a coherent structure. The minimization procedure is iterated until any change is detected. Gaseous cells were previously identified with a higher value than all liquid cells to take into account PLIC reconstruction errors in the process of identifying connected adjacent cells. In this way, label minimization never retains the label of a gas cell and correctly separates each coherent liquid structure. After all this, using MPI communications, structures passing through several cores are connected using a procedure similar to that used in [23, 21, 48].

AlgoDetect can, therefore, detect, classify, and characterize individual liquid structures in a parallel computing environment. The detection and characterization are performed at runtime so that it is not necessary to store a considerable amount of time solutions, and the temporal resolution is arbitrarily chosen in function of the time step.

The detection provides information about each liquid structure, such as the position of the center of mass, its extent in each direction, its surface area and its volume. The output is consequently light in terms of storage. This tool allows the code to analyze complete 3D fields
without significantly impacting the storage needs. The CPU cost of detection strongly depends on the number of droplets and their spatial distribution. For the case investigated in this paper, the detection time is between 8 and 10 times the iteration time. Consequently, to avoid impacting execution time excessively, detection is only performed every 250 iterations to limit additional computing resources. Despite these substantial costs, detection remains highly effective and many of the analyses presented in this paper would have been prohibitively expensive without this tool.

The algorithm is built to perform structure detection and classification. This means that once the raw data have been collected, the liquid shape is analyzed and the structure is tagged in one of the following categories:

**Droplet:** This category contains all the structures that can be considered spherical droplets, which are unlikely to undergo further breakup unless submitted to extreme shearing;

**Ligament:** These particular structures are the typical elongated structures generated in many assisted atomization regimes, particularly the fiber regime, which is the focus of the current work. Ligaments are highly unstable structures as they rapidly undergo Rayleigh-Plateau instability and transform into several child droplets;

**Liquid core:** This is the main body, attached to the liquid injection hole.

Criteria used to classify the structures are based on physical data and are described later in the paper. It is worth noting that, in this work, the generic term *liquid structure* refers to any of the detected isolated liquid structures, including droplets and ligaments, minus the liquid core.

### 3. Simulated configuration

This section details the physical and numerical configurations investigated in this paper. This configuration was chosen to reproduce the case simulated by Rutard et al. [50]. The case reproduces an experimental case used in the CORIA laboratory [18, 17]. A schematic representation of the coaxial injector geometry considered in the simulation is shown in Figure 1.

The computational domain consists of a rectangular box of sizes [$L_x \times L_y \times L_z$] $\approx$ [18.6 $\times$ 5.8 $\times$ 5.8] $\times$ $D_l$, with $L_y$ and $L_z$ the width and depth of the domain and $L_x$ the domain length plus the injector length. The injector is modeled using the immersed interface method

<figure><img src="image_3.png" /><figcaption>Figure 1: Scheme of the injector geometry</figcaption></figure>

[29] and the mesh is a uniform cartesian grid for all simulations presented in this paper. Grid resolution effect is discussed in the following section.

The fluids are injected in the lower $x$ face of the cube by a Dirichlet inflow condition on the velocity, VOF, and $\phi$ fields. The liquid is introduced in the central part of the injector and the gas in the outer ring, around the liquid jet. Non-turbulent uniform velocity profiles are imposed as boundary conditions, and the portion of the injector of length $L_{inj} \approx 1.8 D_l$ helps establish realistic injection conditions in the gas and liquid. On the backplate of the injector, a no-slip boundary condition is imposed, and on all other faces, outflow boundary conditions are considered. Finally, the time step is set so the CFL number does not exceed 0.5 in all the domains.

It is important to notice the liquid-gas density ratio has been reduced to 100 to improve the simulation robustness. To this end, the gas density, viscosity and surface tension were adjusted. Adimensionalized numbers were preserved at same values of [50]: liquid Reynolds number $Re_l \approx 10^4$, gas Weber number $We_g = 490$ and momentum ratio $J \approx 4$. The operating point remains in the fiber regime by referring to the regime map given by [31].

### 4. Mesh refinement effect on liquid structure resolution.

The spatial resolution required to perform a direct numerical simulation is a crucial point in the setup of a massive simulation like the present DNS. In single-phase flow, it is well established that the Kolmogorov
micro-scale evaluates the universality of turbulence diffusion scales so that mesh convergence can be considered achieved once the mesh size can capture this scale. However, in multi-phase flow, this topic is still an open question. In this case, the unknown is the size of the smallest physical droplet created during the atomization process. On the other hand, it can be argued that capturing these smallest scales would probably be unnecessary if the associated liquid mass is a negligible fraction of the total liquid mass.

In the present work, a preliminary assessment of the encountered droplet sizes was made. This estimation was performed using an empirical correlation available in the literature, namely, the average ligament "size" $\langle D_0 \rangle$ given by Marmottant and Villermaux [42]:

$$ \langle D_0 \rangle \approx 0.23\lambda_T, \quad (14) $$

with:

$$ \lambda_T = 2.8\delta_g We_{\delta_g}^{-1/3} \left(\frac{\rho_l}{\rho_g}\right)^{1/3}, \quad (15) $$

$$ \delta_g = 5.6h_g Re_g^{-1/2}, \quad (16) $$

where the Reynolds and Weber numbers defined from gas quantities are:

$$ Re_g = \frac{\rho_g h_g u_g}{\mu_g}, \quad (17) $$

$$ \text{and } We_{\delta_g} = \frac{\rho_g u_g^2 \delta_g}{\sigma}. \quad (18) $$

Using the physical characteristics chosen for the present simulation, the resulting non-dimensional numbers would be $Re_g \approx 10^4$ and $We_{\delta_g} = 7.8$, resulting in a size of the atomized ligaments around $\langle D_0 \rangle \approx 0.55 \delta_{lip}$. As a trade-off between computational resources, spatial resolution and domain dimensions, we chose a mesh size of $\Delta = 0.17 \delta_{lip}$ that discretizes the main structures of the primary atomization process with three points. The grid composed of 2.8 billion cells is distributed as $[N_x \times N_y \times N_z] = [3072 \times 960 \times 960]$.

On the other hand, it is difficult to assess the actual number of grid points needed to resolve a given liquid structure accurately. Consequently, a mesh refinement study was carried out to compare the scales' characteristics resolved with each mesh level. Two more simulations were performed using reduced cell sizes. Besides the main simulation using a cell size of $\Delta_{fine} = 0.17 \delta_{lip}$, a second one employed a cell size of $\Delta_{medium} = 0.255 \delta_{lip}$, while a third one used a cell size of $\Delta_{coarse} = 0.34 \delta_{lip}$. To quantitatively compare the

performances of the three meshes, the global distribution density functions (pdf) of the diameter $d_{30} = \sqrt{V_d}$, where $V_d$ is the droplet volume, were extracted from the three simulations. The probability density functions were built using a Gaussian kernel estimation with the following function:

$$ \text{pdf}_{h_s}(x) = \frac{1}{nh_s} \sum_{i=1}^{n} K\left(\frac{x - x_i}{h_s}\right), \quad (19) $$

where $n$ is the number of droplets considered, $h_s$ the smoothing parameter and $K$ the Gaussian kernel function. In the estimation, the smoothing parameter of each density is imposed at $h_s = 0.0075 \delta_{lip}$ resulting in a good trade-off between the smoothing and the resolution of the pdf for the three cases considered. The results are plotted in Figure 2.

The inherent normalization of pdf that keeps their integral equal to 1 makes it difficult to compare pdf estimated with different mesh resolutions. Figure 2 a) compares the pdfs from the three mesh configurations using directly the equation (19). The shapes of the distributions are similar, with a peak value near $2\Delta$ where $\Delta$ is the cell size of the case considered. However, the absolute values of the pdfs are different: this difference is due to the number of droplets changing drastically between cases. Therefore, considering that a coarse and a fine mesh should simulate the same structures as long as they are resolved on both meshes, Figure 2 b) compares the pdfs normalized by the ratio of the droplet number of the finest mesh $n_{fine}$ over the number $n$ of droplets of the mesh considered.

This last comparison clearly shows that the medium and coarse case distributions fit well up to the resolution limit of $2\Delta$ once normalized. Hence, we can conclude that the medium and coarse cases reproduce the droplet diameter distribution well up to this value. Transposing this conclusion to the fine case means that the proposed DNS simulation can correctly describe the droplet diameter distribution up to the peak of the distribution corresponding to $2\Delta_{fine}$ i.e. $0.34 \delta_{lip}$. This conclusion is coherent with the typical Nyquist criterion reasoning that says it must be at least three points i.e. two cells in each direction to resolve a droplet correctly.

To strengthen this argumentation, the pdfs in Figure 3 were made with droplets bigger than the resolution threshold of the coarser grid $2\Delta_{coarse}$ for the three meshes. It is important to notice that no scaling was performed on those pdfs, only a selection of the liquid structures. The three pdfs superpose well using this selection method. To conclude this mesh influence study, the diameter of the smallest droplets resolved on
<figure><img src="image_4.png" /><figcaption>Figure 2: Probability density functions obtained on different mesh resolutions. a) pdf estimation using the formula (19), b) pdf with a scaling factor ($n/n_{fine}$).</figcaption></figure>

the mesh has been estimated. The following analyses of simulation results presented in this paper were performed on liquid structures larger than this threshold, while the smaller ones were discarded as numerical artifacts.

## 5. Recording time-averaged quantities

To extract valuable results of this DNS and calibrate models used in future LES calculations, time-averaged quantities are recorded for the three velocity components, $VOF$, $\Sigma$, and their fluctuations. The duration of the simulation's transient phase (or transient time)

<figure><img src="image_5.png" /><figcaption>Figure 3: Probability density functions with different mesh resolutions considering droplets bigger than $2\Delta_{coarse}$.</figcaption></figure>

had to be evaluated as the time to start recording time-averaged quantities. One can estimate the transient time using a characteristic speed in the flow. We proposed to use a velocity weighted by momentum to take the differences in fluid momentum into account:

$$V_w = \frac{U_g + JU_l}{1 + J}, \quad (20)$$

in which $J$ is the momentum ratio of the gas and liquid. It gives $V_w = 5.78 \text{ m/s}$ and the corresponding flow-through time is evaluated as $\tau_w = L_x/V_w = 0.018 \text{ s}$. The transient time is then considered to be in $[\tau_w; 2\tau_w]$.

Two quantities have been monitored to characterize the transient phase. The first one is the continuous liquid core length, whose evolution is presented in Figure 4. The solid vertical line, plotted at $t = 2\tau_w$, indicates the time from which the liquid core oscillates around its mean value and seems roughly established. The second quantity, presented in Figure 5, is the evolution of the number of small droplets, ie. with a characteristic length $\sqrt{L_x L_y L_z} \in [0, 5\Delta]$ (where $L_x, L_y, L_z$ are respectively the length of the structures in direction x, y and z). As small droplets are distributed all over the domain due to the whole atomization process, their number evolution gives a good overview of the overall flow behavior. At $t = 2\tau_w$, the number of small droplets is nearly stabilized, which confirms the end of the transient phase.

Consequently, the temporal mean values have been recorded from $t = 2\tau_w$ for $2\tau_w$, ie. between $t = 0.036 \text{ s}$ and $t = 0.073 \text{ s}$. Even if the calculation will be continued in future works, time-averaged quantities seemed sufficiently converged for this duration.
<figure><img src="image_6.png" /><figcaption>Figure 4: Continuous liquid core length evolution and the estimated transient time (solid vertical line).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of the total number of droplets with $\sqrt{L_x L_y L_z} \in [0, 5\Delta]$ and the estimated transient time (solid vertical line).</figcaption></figure>

# 6. Analysis of the liquid flow topology

## 6.1. Global description of the flow solution

The configuration of an assisted atomization jet in the fiber regime exhibits a complex dynamic with multi-scale structures. A global visualization of the simulated flow is proposed in Figure 6 in which the iso-contour $\phi = 0$ is displayed in yellow, representing the liquid interface. In addition, the transverse plane $z = 0$ is colored in grey-scale with the velocity magnitude to materialize the gas injection in the chamber. The dark-grey left-most plan is the iso-contour $\phi_{sol} = 0$ representing the solid surface of the injector. In the center of this plan, the beginning of the gas injector is noticeable, and the liquid interface covers the liquid part of the injector.

The flow structure is complex, with many intricate multi-size structures spread throughout the domain. The first feature to be described is the liquid core, number 1 in Figure 6. In this zone, the liquid core develops streamwise instabilities, which deform the interface. The second one is the location of the primary atomization, with the development of many ligaments of different shapes and sizes. The third one presents a polydispersion of large detached liquid ligaments and droplets, which assume a spherical shape under the effect of surface tension. These three zones are investigated in detail in the following sections. The analysis will be based on the detection and classification algorithm results, which allows a quantitative description of the local flow in terms of topology. Figure 7 shows a classification result on a snapshot of the liquid field. The liquid core, in gray, extends from the injector to about one-third of

the domain, with an increasingly complex shape due to the formation of ligaments. The detached ligaments, in red, are found in consequences from the periphery of the main jet and travel downstream under the convective flow. Droplets are mainly generated from the breakup of the ligaments and become predominant far from the injector.

## 6.2. Destabilization of the liquid core

The spray caused by the atomization originates in the liquid core oscillations generated by the fast surrounding gas jet. This section analyses some essential features of the liquid core to understand the genesis of the atomization process in zone 1. In Figure 6, longitudinal instabilities appear as soon as the liquid jet leaves the injector. In the close instantaneous view in Figure 8 a) lots of early destabilization waves are visible on the liquid core generated by the rapid gas flow. Those waves seem to be caused by the vortex shedding building up in the gas downstream of the injector outlet. In this area, two different structure types are noticeable. The first ones are small ripples on the surface only visible near the injection, identified in Figure 8 b) by red circles. The second ones seem more energetic, but they are tricky to observe in 3D visualizations of the liquid core. Those structures are only visible on slices of the interface, as in Figure 8 b), identified by blue circles.

In the following, a frequency analysis of the interface displacement is developed to characterize the structures just described further. For three different locations downstream from the backplate ($x/D_l = 0.5, x/D_l = 1$
<figure><img src="image_8.png" /><figcaption>Figure 6: Global visualization of an iso-contour of the interface in yellow together with a slice of velocity magnitude in black and white in the middle plane of the domain. Three zones can be distinguished in the flow: 1) the destabilization zone, 2) the atomization zone, and 3) the stabilization zone.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Visualisation of the core, the ligaments and spherical droplets detected in a flow snapshot.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: a) Focus on the destabilization wave on the liquid core. b) Instantaneous slice ($z = 0$) with the liquid phase in black. The red circles locate the tiny ripples, and the blue ones are the energetic structures strongly destabilizing the core.</figcaption></figure>

and $x/D_l = 1.5$), the position of the liquid core interface is recorded in the four cardinal directions ($Y^+ Y^- Z^+$ and $Z^-$) with respect to the central x-axis of the domain. At each x location, Fourier transform analysis is performed on each recorded signal, and the four spectra are averaged to bring out the principal structure frequencies. These averaged spectra are plotted at the three x locations in the figure 9.

This frequency analysis clearly shows a significant amplification of structures between 302 and 357 Hz. This spike frequency matches the energetic structures identified by the blue circles in Figure 8 b). Besides this significant spike, around 350 Hz, smaller peaks between 500 and 700 Hz are visible. These could correspond to the smaller structures identified using the red circles in Figure 8 b).

The origin of those specific frequencies is unclear and the genesis of the instabilities could be assessed using

<figure><img src="image_11.png" /><figcaption>Figure 9: Interface displacement spectra at three axial positions.</figcaption></figure>

the sensors' frequency content at well-chosen locations in the gas. Unfortunately, during the main simulation ie. $\Delta_{fine} = 0.17 \delta_{lip}$, the recording of high-frequency sensors was not yet implemented. However, the velocities have been recorded in some locations for the medium-resolution case. This information gives some clues in understanding the gas flow behavior and the origin of the frequency found on the interface displacement. At three x locations, the streamwise velocity has been recorded at a frequency of around 165 kHz by eight sensors shared around a circle whose radius is equal to $(D_l + h_g)/2 + \delta_{lip}$ the middle of the injector gas channel. Averaged spectra given in Figure 10 are obtained from these sensors at $x = 0.5 D_l$ and $x = 1 D_l$. They confirm the building up of specific frequencies in the gas, and the major amplification is observed for 309 Hz. This frequency perfectly corresponds to the large structures previously described in the interface displacement. The spectra also depict smaller amplification of frequencies between 350 Hz and 890 Hz. In addition, these structures are consistent with the small ripples described in the previous paragraphs. Those two observations indicate that the disturbances seen on the interface come from the turbulence and vortex shedding of the gas injection.

### 6.3. Droplet creation mechanism

To properly study the droplet creation mechanism, the central location of the atomization must be correctly located. As shown visually in Figure 6, zone 2 has the highest droplet concentration. To provide quantitative data and refine this area's location, some spatial statistics have been performed using the individual droplet
<figure><img src="image_12.png" /><figcaption>Figure 10: Streamwise velocity spectra on the medium case i.e. $\Delta_{medium} = 51 \ \mu m$.</figcaption></figure>

information supplied by the AlgoDetect tool. As analyzed in section 4, the droplets are correctly resolved if they have an equivalent diameter $d_{30} > 2\Delta$. The following analyses, therefore, exclude droplets smaller than this resolution criterion.

As the individual liquid structures detected by the algorithm are identified by their center of mass, the numerical domain was discretized in many macro-cells, each spanning several mesh cells in each direction. Therefore, the macro-cells constitute a regular Cartesian grid, each with a 800µm size. The valid liquid structures are then regrouped into the macro-cells according to the position of their center (clustering). This count is then averaged using 420 snapshots of $\phi$ and divided by a volume of a macro-cell $V_{MV}$, resulting in a density of liquid structures. Finally, the clustered distributions have been averaged in the azimuthal direction to account for the three-dimensional distribution of the liquid structures. In the following, all quantities averaged azimuthally will be noted with $\langle \bullet \rangle_{\theta}$.

The first and more straightforward analysis using this procedure counts all the valid detached liquid structures, i.e., those not excluded by the resolution criterion. The resulting distribution is presented in Figure 11, together with the same zone boundaries as in Figure 6. Droplets stuck on the backplate are not considered in the count because their behavior is not of interest here. The distribution shows maximum values located in the middle of zone 2, confirming the location of the dominant atomization process in this zone. An interesting remark can be made on the distribution by qualitatively comparing

the present results with Le Touze et al. [32]. The model used in [32] predicts a strong formation of fine droplets very close to the injection. On the contrary, zone 1 shows a large area where almost no droplet is observed. A destabilization time is necessary for the formation of the first droplets. However, this comparison should be taken with caution as significant differences exist between the conditions considered in [32] and the present simulation.

The surface instabilities that develop on the liquid core generate large pockets of liquid. With the action of the fast flow of gas, these large pockets expand to their point of rupture, where they become elongated structures called ligaments. The ligaments rapidly undergo Rayleigh-Plateau breakup, splitting themselves into many droplets. At the end of the breakup, the liquid has reached a stable state where all the resulting structures are spherical droplets, which do not undergo further atomization unless submitted to a renewed high-shearing effect. A visual example of a break-up chronology of a ligament just detached from a pocket is presented in Figure 12. Figure 12 a) shows a long ligament just detached from the core. Then, in the following Figures 12 b) and c), the structures generated by this ligament are tracked. These three figures illustrate the typical break-up mechanism and stabilization into spherical droplets occurring in the fiber regime.

To highlight the stabilization and quantify the characteristic size of the liquid structures across the domain, Figure 13 shows the local size distribution of the structures. The same methodology of clustering and azimuthal average (as in Figure 11) is employed here to make the Figure 13. However, in this figure, the numbers of structures are not simply added together, but their characteristic diameters are averaged in each cluster. Using a simple average of droplet diameters in each cluster is not a relevant option because it introduces a strong bias towards the smaller structures. The calculation of the average diameter over a cluster is therefore weighted using the volume of each structure to counteract this bias. The average characteristic diameter $D_v$ is then defined as follows:

$$D_v = \frac{\sum d_{30,p} V_p}{\sum V_p}, \quad (21)$$

where the sums consider each structure of a cluster, $V_p$ is the volume of a structure, and $d_{30,p} = \sqrt[3]{V_p}$ is a characteristic diameter often used in the literature to study atomization. These plots confirm the characteristics of the flow depicted in figure 7: a large number of smaller droplets are visible near the primary atomization zone, while a population of larger droplets, coming
<figure><img src="image_13.png" /><figcaption>Figure 11: Density of valid detached liquid structures (mm⁻³). Count gathered in all snapshots recorded in the DNS.</figcaption></figure>

from the stabilization of the larger axial ligaments, characterize the second half.

## 6.4. Ligament characteristics

The atomization process depicted in Figure 12 is directly due to the shear stress of the gas. It is well described in Dumouchel et al. [11] and defined as textural atomization. The creation of droplets seems to first spend most of the time through the formation of detached ligaments. Hence, their detection is a good indicator of the progress of this textural atomization. Therefore, detached ligament identification is made in post-processing using the data available after the AlgoDetect procedure. Two geometrical criteria are defined by equations (22) and (23) to detect ligaments structures:

$$ \frac{Max(L_x, L_y, L_z)}{Min(L_x, L_y, L_z)} > C_1, \quad (22) $$

$$ \text{and} \quad \frac{V}{\frac{\pi}{6} L_x L_y L_z} < C_2, \quad (23) $$

where $L_x$, $L_y$, and $L_z$ are the dimensions of the rectangular box containing the liquid structure, V its volume, and $C_1$ and $C_2$ constants wisely chosen to select the desired structures. In the following, those constants have been imposed at $C_1 = 1.6$ and $C_2 = 0.5$. The first criterion defined in equation (22) measures the aspect ratio of the bounding box of the structures. It, therefore, accurately detects elongated structures in the direction of the mesh, as shown in Figure 14 a). However, this criterion does not detect ligaments when they are diagonal to the mesh, as in Figure 14 b), because their bounding box is almost cubic. In this case, it is mandatory to have another criterion to detect these ligaments. The

second criterion defined in equation (23) compares the volume of the structure with the volume of an ellipsoid contained in the same box as the structure. These criteria are, therefore, complementary and allow the detection of most ligament structures, as shown in Figure 7. Spherical droplets are the remaining resolved liquid structures undetected by the criterion (22) and (23). Unresolved structures are liquid structures with a size under the resolution limit characterized in section 4.

Those criteria allow the detection of most ligaments. However, it is hard to predict the behavior of ligaments because of their various shapes and sizes. The early stage of atomization generates many tiny ligaments. These structures quickly transform into numerous tiny droplets. Downstream, ligament structures are still present, but they are significantly larger. Those larger ligaments seem to come from two different sources. The first ones, as early small ligaments, come from pockets of textural atomization [11] occurring later in the flow where the gas no longer has enough energy to disrupt them strongly. The second ones are structures coming directly from the center of the liquid core named structural atomization in [11]. Those two last types have a much longer stabilization time and do not produce many child structures or spherical droplets.

All those observations demonstrate the central role of the ligaments in the atomization process. The direction, velocity, and sizes of the final droplets seem to derive directly from the characteristics of ligaments. Moreover, their large ratio of surface over volume could strongly impact the averaged interface area density. The relation between ligaments and interface area density must be clarified because many authors use the interface area density directly to model atomization without consider-
<figure><img src="image_14.png" /><figcaption>Figure 12: Tracking of structures from a ligament detached from a liquid pocket in three subsequent snapshots.</figcaption></figure>

ing the presence of ligaments in the process. Therefore, Figure 7 perfectly highlights their possible small dimensions, and in some LES frameworks, those ligaments could easily be unresolved.

## 7. In-depth analysis of interface area evolution

The following analysis focuses on the spatial and temporal evolution of the interface area. In interface area density models, this quantity is used to statistically describe the dispersion of the local liquid mass, whose topology is unknown as the characteristic scale becomes smaller than the local mesh size. Using the resolved DNS fields, the spatial distribution of the calculated area can be extracted and correlated with the local physical phenomena occurring during primary atomization.

This knowledge can then be used to develop appropriate large-scale models.

### 7.1. Integrated interface area density evolution

The interface area density is computed in each cell by the ratio between the surface and the volume of the liquid :

$$ \Sigma = \frac{S_{\text{PLIC}}}{V_{\text{cell}}}, \qquad (24) $$

with $S_{\text{PLIC}}$ the area of the interface plane from the PLIC reconstruction (see section 2.2) and $V_{\text{cell}}$ the volume of the cell. The quantity of Equation (24) was computed in the whole simulation domain and averaged in time on all the time steps of the permanent regime. This field was then integrated in space along the azimuthal directions to obtain the axial profile, plotted in Figure 15. The field is also averaged in the azimuthal directions (in a similar way to Figure 11) to obtain the radial profiles, plotted in Figure 16.

The evolution of the axial profile of the interface area density integral (Figure 15) shows a steep growth near the injector nozzle, in particular in the highly unsteady ligament-forming zone 2 of Figure 6). After a peak located near $x/D_t = 2.5$, it slowly decreases with the distance to the injector. Radial profiles of interface area density are plotted in Figure 16 at the streamwise locations highlighted in Figure 6. The nearest position to the injector, $x = 0.5 D_t$, shows that the liquid-gas interface concentrates near the liquid core boundaries, forming a tall, thin peak. This can be easily explained by the low amplitude movement of the interface in this area, which concentrates $\Sigma$ in the mean field. Downstream, profiles become wider as the gas-liquid mixing layer thickness increases. Notably, the profiles of Figures 15 and 16 provide a suitable data set to evaluate and validate the models predicting the evolution of $\Sigma$. To the author's knowledge, no model has yet been able to reproduce the behavior shown by the DNS.

To understand the consequences of atomization on the interface area density, it is interesting to correlate the map of the identified liquid structures and their classification with the interface area density fields. The following section will provide more insight into the evolution of surface density and local fluid topology.

### 7.2. Mechanisms of interface area density production

Figure 17 shows the average interface area density field superposed to the iso-contours obtained from detecting and classifying the detached structures presented in section 6. Two sets of iso-contours are shown, the green representing the density of ligaments and the blue
<figure><img src="image_15.png" /><figcaption>Figure 13: Averaged structure size distribution recorded in all snapshots recorded in the DNS of the simulation adimensionalized by lip length.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Examples of ligaments in different configurations: a) aligned in the direction of the mesh, b) diagonal to the mesh directions.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Integral along the *x* axis of the time-averaged interface area density Σ.</figcaption></figure>

the density of spherical droplets detected as explained in section 6.4. The same azimuthal averaging procedure was performed for those two last quantities as for Figure 11.

Figure 17 confirms the atomization scenario proposed in the previous section. The ligaments need a characteristic time to form and then appear only at some distance from the injection plane, around one jet diameter. Indeed, an important ligament detachment from the liquid core occurs around the abscissa $x/D_l = 1$. Spherical droplets resulting from ligament rupture appear further downstream in the flow than ligaments. This observation confirms the scenario where atomization first passes over a ligament phase, and then creates droplets. Finally, the global shape of the spherical droplet contours in blue suggests a dispersion of these structures in space due to the expansion of the spray.

Figure 17 illustrates the complexity of interface area density interpretation. Indeed, the highest absolute values of Σ are located near the injection plane, where almost no individual structure exists. These high values are due to the low amplitude oscillations of the interface, which statistically concentrate at this location the instantaneous value of Σ. After this early stage, the oscillation grows larger, expanding the jet. At the same time, the global Σ increases from the generation of the ligaments, whose shape is characterized by high surface-to-volume ratios. Then, the Σ slowly decreases as the spherical structures gradually replace the ligaments. This observation is very interesting because an increase in the surface is the expected result of generat-
<figure><img src="image_18.png" /><figcaption>Figure 16: Radial profiles of spatial- and time-averaged interface area density $\langle\Sigma\rangle_{\theta}$.</figcaption></figure>

ing many small spherical droplets. However, this decrease is easily explained by the minimization of the surface area when the liquid structures take the form of spheres, which give the least possible surface-to-volume ratio for a given radius: the total surface still depends on the average droplet diameter.

To visualize the evolution of $\Sigma$ more precisely, Figure 18a) shows the evolution of the integral of $\Sigma$ along the $x$ axis together with integrated numbers of ligaments and spherical droplets in the total time average window. The integrated interface area density correlates well with the number of ligaments in function of the streamwise coordinate. The location of the peak of the interface area corresponds to the maximum number of detected detached ligaments. Downstream, as ligaments are quickly atomized and replaced by spherical droplets, the profile of $\Sigma$ follows the destruction of the ligaments and seems unaffected by the presence of the newly created spherical droplets: the surface carried by the latter is, therefore, small compared to the one carried by the ligaments at this location. To go further into the understanding of interface area creation and destruction phenomena, Figure 18b (resp. c) shows the evolution of the integral of the surface (resp. the volume) of the liquid core, the detached ligaments and the spherical droplets.

The surface evolution shows that the liquid core concentrates most of the interface area in the first $4 D_l$ (jet diameters) in the domain. The reason is that strong ripples and elongated structures are formed on the fringes of the liquid core. It takes some time for them to detach and transform into isolated filaments and, finally, into spherical droplets. The large number of detached ligaments and droplets in the $2 < x/D_l < 4$ region

possess, however, low volume and surface, these quantities increasing with the distance to the backplate. The first detached structures are smaller, with higher Weber numbers and lower characteristic times. Larger structures take more time to deform and break up under the shearing airflow. Downstream, starting from $x/D_l = 6$, the volumes of ligaments and spherical droplets increase to a plateau, the maximum surface of spherical droplets exceeding by a factor of two the ligaments one. The global interface area density continues to decrease steadily. This means that when the surface tension forces become dominant over inertial and viscous forces, the liquid structures assume a spherical form that globally decreases the total interface area density.

These graphs also confirm the evolution of the liquid structures deduced from Figure 12. The structures break up like ligaments and rapidly transform into many small spherical droplets. First, the surface created by the ligaments suddenly reaches a plateau around the abscissa $x/D_l = 2.5$. Moreover, this location corresponds to the peak of $\Sigma$ and ligament number, meaning the ligaments start disappearing at this location. Finally, this location also corresponds to a strong increase in spherical droplet number because the ligaments are massively transformed into droplets here.

Downstream the abscissa $x/D_l = 5$, the decrease of interface area density is still induced by stabilizing ligaments in spherical droplets, but another phenomenon seems to contribute to the decrease of $\Sigma$. As the flow progresses, fewer and fewer structures are present in the area while they become larger on average. The overall enlargement of the structures is visible in Figures 18a) and c) as the overall volume increases while the number of structures decreases, especially for ligament structures. Indeed, in zone 3, atomization creates larger structures coming directly from the core, with a low velocity, aligned with the $x$-axis, so these structures remain longer in the domain than tiny droplets. Dumouchel et al. [11] described this phenomenon as structural atomization.

It is worth noting that a deeper analysis of the interface area on the liquid core could allow us to understand better the interface area density spatial evolution in the domain. Some methodologies, proposed by [11], for example, are very promising for studying the evolution of the local topology of large structures like the liquid core.

### 7.3. Mechanisms of interface area density reduction

In the previous section, the results in Figures 18b) and c) showed a constant increase in the total volume of droplets with the distance from the injector. Conversely,
<figure><img src="image_19.png" /><figcaption>Figure 17: Half-cut of median plane colored by the spatial- and time-averaged interface area density using a logarithmic scale ($m^{-1}$), together with isolines of ligament density ($mm^{-3}$) in green and isolines of the density of spherical structures ($mm^{-3}$) in blue.</figcaption></figure>

the corresponding surface presented an initial increase up to a peak and subsequently a decline. These results may suggest that the detached liquid structures tend to stabilize and become spherical with the distance from the injection. This last section focuses on the stabilization achieved by the liquid structures and the characterization of their average size at the end of the domain.

The evolution of total volume structures in Figure 18 c) shows a plateau after $x/D_l = 7$, especially for the detached ligaments. At the same location, Figure 18 b) shows a decrease in the global surface. These two remarks combined show that liquid structures tend to stabilize into more spherical structures.

In the following, the quantitative evolution of the characteristic diameter $D_v$ will be evaluated with the help of profile extractions at specific locations plotted in Figure 19. The locations of the extractions are indicated using the vertical dashed lines in Figure 13.

Upstream profiles show no sign of atomization at the locations $x/D_l = 0.5$ and $x/D_l = 1$. Indeed, the only liquid droplets present after the core emphasis come from the recirculation zone. Only droplets with a characteristic size around $0.5 \delta_{lip}$ are small enough to be conducted by the recirculation flow. No droplets are directly created by the core in those locations, confirming the delay observed in the atomization process. Their average size below $0.5 \delta_{lip}$ is a direct consequence of the threshold of $2\Delta = 0.34 \delta_{lip}$ chosen to consider a liquid structure in the statistical analysis. The profile located at $x/D_l = 2$ is way more interesting regarding atomization. Located upstream in zone 2, it shows clearly the consequence

of atomization. The droplet created by the atomization in this location shows a constant characteristic diameter $D_v$ around $1.0 \delta_{lip}$. After the radius $r/D_l = 2$, the size of the liquid structures decreases, and the profile matches the first two profiles. This indicates the end of the predominance of atomization coming directly from the core. The droplet size of about $0.5 \delta_{lip}$ clearly indicates that the recirculation zone drives the structure size in this area. The following profile, located at $x/D_l = 4$, depicts the end of the strong atomization zone. This profile depicts two different phenomena. First, the average diameter increases in the spray passing from $1.0 \delta_{lip}$ in the profile $x/D_l = 2$ to $2.0 \delta_{lip}$ in this one. Secondly, as in the previous profile, the liquid structures' average size decreases but is much further away from the axis. Using the beginning of the decrease in size in the different profiles, it is possible to identify the extent of the spray generated by the primary atomization. Finally, the domain is not large enough to reach the end of the spray emphasis in this profile completely. This explains why the droplet size does not collapse with the three previous profiles at large radial distances from the axis. The profile at $x/D_l = 8$ shows the largest structures with a strong peak near the axis. The peak corresponds to large structures released from the liquid core by textural atomization. After that, a trend similar to all previous profiles is visible, with a decrease in size away from the axis. Finally, the last profile, located at $x/D_l = 10$, also shows a peak near the axis with a similar decrease trend away from the axis. The particular observation here is the overall decrease in size compared
<figure><img src="image_20.png" /><figcaption>a)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 18: Integral along the *x* axis of a) the time-averaged interface area density Σ, numbers of ligaments and spherical droplets, b) Surfaces of the liquid core, ligaments and spherical droplets, c) Volumes of the liquid core, ligaments and spherical droplets.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 19: Radial profiles of characteristic diameter *D*<sub>v</sub>.</figcaption></figure>

to the *x*/*D*<sub>*l*</sub> = 8 profile. The global stabilization of the structures could explain this decrease in size. Atomization would no longer play a determining role here, and the structures would become more and more spherical when their characteristic size decreases.

To conclude and summarize all the observations made in this section, an attempt will now be made to highlight the overall trends. First, the average structure size highly depends on the distance from the injector. It increases as you move away from the injection plane before stopping in the area where atomization no longer plays a significant role, and tends to decrease slightly at the end of the domain as the structures become more spherical. Secondly, textural atomization plays a significant role upstream, with tiny droplets created on the sides of the liquid core. It is replaced downstream by structural atomization in which large structures, mainly small liquid core parts, are concentrated on the axis. However, the transition between these two regimes does not appear sharp. In textural atomization, the droplet size seems to be relatively constant and only correlated to the distance from the injection; the decrease in size only occurs at the exit of the spray influence. In structural atomization, large structures seem to be more concentrated near the axis with a continuous decrease in structure size away from the axis.

**8. Conclusion**

To the best of the authors' knowledge, this paper presents the most resolved simulation to date of a round coaxial injection configuration in fiber atomization regime. A rapid gas flow is injected around a slow
liquid flow. The fiber regime's high-velocity difference between gas and liquid creates an intense atomization process. This paper aimed to give an in-depth insight into this complex phenomenon, with a detailed analysis of several aspects, and new data available for subgrid-scale model formulation and validation in the LES context. The focus was on the interface area, which is a driving quantity for atomization and evaporation when dealing with liquid fuels.

After the global description, a specific effort was made to justify the numerical refinement of the calculation and the physical time considered. The analysis was organized following the flow progression. First, a frequency analysis was performed on the liquid core to investigate the origin of the atomization process. Next, the new detection algorithm named AlgoDetect was used to analyze the main atomization, allowing the characterization of the size and shape of the liquid structures. The ligament structures are put into perspective in the atomization process, and the relation with the evolution of the liquid-gas interface area density was analyzed. Finally, the focus was on the droplets characteristic size in the spray and the stabilization of liquid structures of spherical shape. The flow analysis made it possible to specify the scenario occurring during atomization in the fiber regime. Close to the injector exit, a no-atomization zone was characterized by ripples visible on the liquid core, whose specific frequency is undoubtedly the primary origin of the whole atomization process of this regime. At some point in the liquid core, characterizing an atomization delay, the interface ripples become large enough to create large pockets of liquid together with elongated structures attached to the core. After that, the atomization of those large structures creates a population of elongated detached ligaments. Very quickly, these detached ligaments are also atomized into numerous spherical droplets. The resulting droplets are finally too small to be significantly affected by the atomization. This textural atomization was precisely located, with the number of detached ligament structures found as a good indicator of its progress. As for the interface area, it was observed that the liquid core involves the major part of the interface area during this process. The complete analysis of this liquid core, not discussed in this paper, is ongoing work.

In addition to the primary atomization scenario previously discussed, this article describes the liquid-gas interface area density in assisted atomization configuration. It was linked with the ligament and droplet volume and surface evolution to characterize the share of detached ligaments that induce interface area density compared to spherical droplets. The interface area density

profiles provided in this work will be essential for validating future LES simulations using subgrid-scale modeling of the interface area density.

This work's last significant assessment is about the droplet size evolution in the spray. A strong correlation links the creation's location and the structures' resulting characteristic size. The droplets created at an early atomization stage are much smaller than those created later in the flow. In addition, due to their small inertia, the smallest droplets quickly follow the gas expansion away from the axis.

The continuous growth of high-performance computing allowed the simulation described here, able to describe a flow as complex as coaxial atomization in fiber regime. Until recently, experimentation was the only way to gather data to validate the models used for intermediate-fidelity simulations like LES. However, the experimentation has some limitations. It is hard to characterize the atomization zone by visual observations because of the strong superposition of liquid structures in the field of view where the liquid is still dense. Size measurements can be made in the dispersed phase, but they only allow for an indirect validation of primary atomization modeling as they are located far from the atomization location. The DNS calculation presented in this article is a new and promising way to collect data to validate atomization models. For future studies, some improvements in post-processing should be made to help better understand the mechanism observed in atomization. The first point is tracking the evolution and displacement of individual structures. Structure tracking algorithms are already available in the DYJEAT framework and could easily be applied in this type of simulation, but their cost is currently prohibitive in this configuration. Furthermore, the resulting relation tree analysis between the structures is not trivial. This type of analysis could allow a better understanding of the correlation between the place of creation and the size of the droplets. To complete these analyses, the jet liquid core should also be characterized with in-depth advanced processing to understand processes behind primary atomization better.

## 9. Acknowledgment

For the purpose of Open Access, a CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

This work was performed using CCRT HPC resource (TOPAZE supercomputer).
[1] G. Agbaglah, R. Chiodi, and O. Desjardins. Numerical simulation of the initial destabilization of an air-blasted liquid layer. *J. Fluid Mech.*, 812:1024–1038, 2017.

[2] J. Anez, A. Ahmed, N. Hecht, B. Duret, J. Reveillon, and F.X. Demoulin. Eulerian-lagrangian spray atomization model coupled with interface capturing method for diesel injectors. *Int. J. Multiphase Flow*, 113:325–342, 2019.

[3] M. Averseng. Contribution à la modélisation de l'atomisation assistée par l'analyse de simulations haute fidélité. PhD thesis, ISAE Toulouse, 2022.

[4] G. Blanchard, P. Villedieu, and D. Zuzio. Numerical simulation of primary atomization of a sheared liquid sheet. part 1: model and numerical method description. In *ILASS Europe*, 2013.

[5] G. Blanchard, P. Villedieu, and D. Zuzio. Numerical simulation of primary atomization of a sheared liquid sheet. part 2: comparison with experimental results. In *ILASS Europe*, 2013.

[6] G. Charalampous, Y. Hardalupas, and A. Taylor. Structure of the continuous liquid jet core during coaxial air-blast atomization. *International Journal of Spray and Combustion Dynamics*, 1:389–415, 2009.

[7] J. Chesnel, J. Reveillon, T. Menard, and F.-X. Demoulin. Large eddy simulation of liquid jet atomization. *Atomization Sprays*, 21:711–736, 2011.

[8] F. Couderc. Development of a numerical code for the simulation of non-miscible fluid flow. Application to the air-assisted disintegration of a liquid jet. PhD thesis, Ecole nation supérieure de l'aeronautique et de l'espace, 2007.

[9] B. Dejean. Experimental study of liquid sheet atomization with and without prefilming zone for its modelling - Boundary condition influence. PhD thesis, Institut Supérieur de l'Aéronautique et de l'Espace (ISAE), 2015.

[10] O. Desjardins, J. McCaslin, M. Owkes, and P. Brady. Direct numerical and large-eddy simulation of primary atomization in complex geometries. *Atomization Sprays*, 23:1001–1048, 2013.

[11] C. Dumouchel, J.-B. Blaisot, F. Abuzahra, A. Sou, G. Godard, and S. Idlahcen. Analysis of a textural atomization process. *Exp. Fluids*, 60, 2019.

[12] B. Duret. Simulation numrique directe des coulements liquide-gaz avec vaporation : application l'atomisation. PhD thesis, INSA Rouen, 2013.

[13] Baptiste Djean, Pierre Berthoumieu, and Pierre Gajan. Experimental study on the influence of liquid and air boundary conditions on a planar air-blasted liquid sheet, part i: Liquid and air thicknesses. *International Journal of Multiphase Flow*, 79:202–213, 2016. ISSN 0301-9322. doi: https://doi.org/10.1016/j.ijmultiphaseflow.2015.09.002.

[14] Z. Farago and N. Chigier. Morphological classification of disintegration of round liquid jets in a coaxial air stream. *Atomization Sprays*, 1992.

[15] R. P Fedkiw, T. Aslam, B. Merriman, and S. Osher. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). *J. Comput. Phys.*, 152:457–492, 1999.

[16] V. Gutirrez Fernández. Experimental study of liquid sheet disintegration in a high pressure environment. Theses de doctorat, Institut Supérieur de l'Aéronautique et de l'Espace (ISAE), 2009.

[17] A. Ficuciello. Analysis of high frequency/high amplitude acoustic field effects on coaxial injection : application to liquid rocket engines. PhD thesis, Université de Rouen, 2018.

[18] A. Ficuciello, J.B. Blaisot, C. Richard, and F. Baillot. Investigation of air-assisted sprays submitted to high frequency transverse acoustic fields: Droplet clustering. *Phys. Fluids*, 29: 067103, 2017.

[19] P. Gaillard, C. Le Touze, L. Matuszewski, and A. Murrone. Numerical Simulation of Cryogenic Injection in Rocket Engine

Combustion Chambers. *Aerospace Lab*, page 16, 2016.

[20] L. He, X. Ren, Q. Gao, X. Zhao, B. Yao, and Y. Chao. The connected-component labeling problem: A review of state-of-the-art algorithms. *Pattern Recogn.*, 70:25–43, 2017.

[21] D. Hebert, D. Schmidt, D. Knaus, S. Philips, and P. Magari. Parallel VOF spray droplet identification in an unstructured grid. In *21th ILASS Conference*, 2008.

[22] K. Hendrickson, G. Weymouth, and D. Yue. Informed component label algorithm for robust identification of connected components with volume-of-fluid method. *Comput. Fluids*, 197: 104373, 2020.

[23] M. Herrmann. A parallel eulerian interface tracking/lagrangian point particle multi-scale coupling procedure. *J. Comput. Phys.*, 229:745–759, 2010.

[24] M. Hong. Atomisation et mélange dans les jets coaxiaux liquid-gaz. PhD thesis, INPG, 2003.

[25] S. Jay, F. Lacas, and S. Candel. Combined surface density concepts for dense spray combustion. *Combust. Flame*, 144:558–577, 2006.

[26] D. Jiang and Y. Ling. Destabilization of a planar liquid stream by a co-flowing turbulent gas stream. *Int. J. Multiphase Flow*, 122:103121, 2020.

[27] Sung-Gil Kim and Wonjung Kim. Drop impact on a fiber. *Physics of Fluids*, 28(4):042001, 04 2016. ISSN 1070-6631. doi: 10.1063/1.4945103. URL https://doi.org/10.1063/1.4945103.

[28] A. Kumar and S. Sahu. Liquid jet breakup unsteadiness in a coaxial air-blast atomizer. *International Journal of Spray and Combustion Dynamics*, 10:211–230, 2018.

[29] I. Lagrange, A. Orazzo, D. Zuzio, and J.-L. Estivalzes. Immersed interface method for the direct numerical simulation of air-blast primary atomization. In *ILASS America*, 2017.

[30] C. Larricq-Fourcade. Étude de la pulvérisation assistée en air d'une nappe liquide et influence d'un vent ionique sur les instabilités hydrodynamiques. Theses, ISAE, Université de Toulouse, 2006.

[31] J.-C. Lasheras and E.-J. Hopfinger. Liquid jet instability and atomization in a coaxial gas stream. *Annu. Rev. Fluid Mech.*, 32:275–308, 2000.

[32] C. Le Touze, LH. Dorey, N. Rutard, and A. Murrone. A compressible two-phase flow framework for Large Eddy Simulations of liquid-propellant rocket engines. *Appl. Math. Model.*, 84:265–286, 2020.

[33] R. Lebas, T. Menard, P.A. Beau, A. Berlemont, and F.X. Demoulin. Numerical simulation of primary break-up and atomization: DNS and modelling study. *Int. J. Multiphase Flow*, 35: 247–260, 2009.

[34] A. H. Lefebvre and V. G. McDonell. *Atomization And Sprays*. CRC Press, 2017.

[35] X. Li and M. C. Soteriou. Detailed numerical simulation of liquid jet atomization in crossflow of increasing density. *Int. J. Multiphase Flow*, 104:214–232, 2018.

[36] Y. Liang, L. C. Johansen, and M. Linne. Breakup of a laminar liquid jet by coaxial non-swirling and swirling air streams. *Phys. Fluids*, 34:093606, 2022.

[37] Y. Ling, S. Zaleski, and R. Scardovelli. Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model. *Int. J. Multiphase Flow*, 76:122–143, 2015.

[38] Y. Ling, D. Fuster, G. Tryggvason, and S. Zaleski. A two-phase mixing layer between parallel gas and liquid streams: multi-phase turbulence statistics and influence of interfacial instability. *J. Fluid Mech.*, 859:268–307, 2019.

[39] J. López and J. Hernández. Analytical and geometrical tools for 3d volume of fluid methods in general grids. *J. Comput. Phys.*,
227:5939–5948, 2008.

[40] A. Lozano, E. Calvo, J. A. García, and F. Barreras. Mode transitions in an oscillating liquid sheet. *Phys. Fluids*, 23:044103, 2011.

[41] A. Mansour and N. Chigier. Dynamic behavior of liquid sheets. *Phys. Fluids A*, 3:2971–2980, 1991.

[42] P. Marmottant and E. Villermaux. On spray formation. *J. Fluid Mech.*, 498:73–111, 2004.

[43] J.-P. Matas and A. H. Cartellier. Flapping instability of a liquid jet. In *7th International Conference on Multiphase Flow ICMF 2010*, page 413, 2010.

[44] 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.

[45] T. Müller, A. Snger, P. Habisreuther, T. Jakobs, D. Trimis, T. Kolb, and N. Zarzalis. Simulation of the primary breakup of a high-viscosity liquid jet by a coaxial annular gas flow. *Int. J. Multiphase Flow*, 87:212–228, 2016.

[46] N. Odier, G. Balarac, C. Corre, and V. Moureau. Numerical study of a flapping liquid sheet sheared by a high-speed stream. *Int. J. Multiphase Flow*, 77:196–208, 2015.

[47] L. Palanti, S. Puggelli, L. Langone, A. Andreini, J. Reveillon, B. Duret, and F.X. Demoulin. An attempt to predict spray characteristics at early stage of the atomization process by using surface density and curvature distribution. *Int. J. Multiphase Flow*, 147:103879, 2022.

[48] P. Pepiot and O. Desjardins. Numerical analysis of the dynamics of two- and three-dimensional fluidized bed reactors using an Euler-Lagrange approach. *Powder Technol.*, 220:104–121, 2012.

[49] F. Ben Rayana. *Contribution l'tude des instabilités interfaciales liquide-gaz en atomisation assiste et taille de gouttes*. PhD thesis, Grenoble INPG 2007, 2007.

[50] N. Rutard, L.H. Dorey, C. Le Touze, and S. Ducruix. Large-eddy simulation of an air-assisted liquid jet under a high-frequency transverse acoustic forcing. *Int. J. Multiphas. Flow*, 122:103144, 2020.

[51] J. Shinjo and A. Umemura. Surface instability and primary atomization characteristics of straight liquid jet sprays. *Int. J. Multiphase Flow*, 37:1294–1304, 2011.

[52] M. Sussman and E. G. Puckett. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. *J. Comput. Phys.*, 162:301–337, 2000.

[53] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressible two-phase flow. *J. Comput. Phys.*, 114:146–159, 1994.

[54] G. Tomar, D. Fuster, S. Zaleski, and S. Popinet. Multiscale simulations of primary atomization. *Comput. Fluids*, 39:1864–1874, 2010.

[55] A. Vallet. Contribution la modélisation de l'atomisation d'un jet liquide haute pression. PhD thesis, Université de Rouen Normandie, 1997.

[56] A. Vallet and R. Borghi. Modélisation eulerienne de l'atomisation d'un jet liquide. *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy*, 327:1015–1020, 1999.

[57] T. Xavier, D. Zuzio, M. Averseng, and J.-L. Estivalezes. Toward direct numerical simulation of high speed droplet impact. *Meccanica*, pages 1–15, 2019.

[58] F. Xiao, M. Dianat, and J.J. McGuirk. Les of turbulent liquid jet primary breakup in turbulent coaxial air flow. *Int. J. Multiphase Flow*, 60:103–118, 2014.

[59] D. L. Youngs. Time-dependent multi-material flow with large fluid distortion. *Numerical Methods for Fluid Dynamics*, 1982.

[60] F. Zhang, T. Zirwes, T. Müller, S. Wachter, T. Jakobs, P. Habisreuther, N. Zarzalis, D. Trimis, and T. Kolb. Effect of elevated pressure on air-assisted primary atomization of coaxial liquid jets: Basic research for entrained flow gasification. *Renewable Sustainable Energy Rev.*, 134:110411, 2020.

[61] D. Zuzio, J.-L. Estivàlèzes, and B. DiPierro. An improved multi-scale eulerian–lagrangian method for simulation of atomization process. *Computers and Fluids*, 176:285–301, 2018.

[62] D. Zuzio, A. Orazzo, J.-L. Estivàlèzes, and I. Lagrange. A new efficient momentum preserving level-set/VOF method for high density and momentum ratio incompressible two-phase flows. *J. Comput. Phys.*, 410:109342, 2020.