Romain Janodet
email: romain.janodet@coria.fr

Carlos Guillamón

Vincent Moureau

Renaud Mercier

Ghislain Lartigue

Pierre Bénard

Thibaut Ménard

Alain Berlemont

A massively parallel accurate conservative level set algorithm for simulating turbulent atomization on adaptive unstructured grids

Keywords: Multiphase flows, Atomization, Unstructured grids, Incompressible flow LES, Conservative level set, Adaptive mesh refinement

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multiphase flows are ubiquitous in nature and in industrial systems, and understanding the formation of sprays is crucial for designing complex injection technologies. More specifically, turbulent primary atomization is of first importance in aeronautical combustors, into which fuel is injected in liquid form, goes under atomization, evaporation, turbulent mixing and eventually combustion. All these physical processes are strongly dependent, and thus improving numerical predictions of pollutant formation requires deep numerical insights into the atomization process. This process designates the evolution of large liquid structures into small droplets and thus governs the liquid droplet distribution in the combustor and the subsequent evaporation. Realistic industrial injector geometries for atomization simulations are very complex, including sharp edges and narrow corners, which hinders the meshing task with cartesian grids: the use of unstructured meshes is thus well suited for this type of simulations [START_REF] Leparoux | Primary atomization simulation applied to a jet in crossflow aeronautical injector with dynamic mesh adaptation[END_REF]. Although the literature remains limited about unstructured atomization computations, previous studies on static meshes can be found, such as [START_REF] Jofre | A 3D Volume-of-Fluid advection method based on cell-vertex velocities for unstructured meshes[END_REF][START_REF] Jofre | Parallel load balancing strategy for Volume-of-Fluid methods on 3D unstructured meshes[END_REF][START_REF] Ivey | Accurate interface normal and curvature estimates on three-dimensional unstructured non-convex polyhedral meshes[END_REF][START_REF] Ivey | Conservative and bounded volume-of-fluid advection on unstructured grids[END_REF][START_REF] Evrard | Estimation of curvature from volume fractions using parabolic reconstruction on two-dimensional unstructured meshes[END_REF][START_REF] Patel | A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids[END_REF][START_REF] Xie | A conservative solver for surface-tension-driven multiphase flows on collocated unstructured grids[END_REF][START_REF] Maric | Unstructured un-split geometrical Volume-of-Fluid methods -A review[END_REF] for Volume Of Fluid (VOF) interface capturing techniques [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flows[END_REF], or [START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF][START_REF] Kees | A conservative level set method suitable for variable-order approximations and unstructured meshes[END_REF][START_REF] Balcazar | A finite-volume/level-set method for simulating two-phase flows on unstructured grids[END_REF][START_REF] Bøckmann | A gradient augmented level set method for unstructured grids[END_REF][START_REF] Lin | A volume-conserving balanced-force level set method on unstructured meshes using a control volume finite element formulation[END_REF] for Level Set methods (LS) [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF][START_REF] Sethian | Level Set Methods for Fluid Interfaces[END_REF].

Despite many years of study, high-performance atomization simulations remain a challenge: the wide range of involved time and space scales leads to important calculation costs. To tackle this issue, the use of Adaptive Mesh Refinement (AMR) is particularly helpful, allowing the allocation of computational resources to the interface dynamics. A popular technique for solving incompressible two-phase flows on adaptive cartesian grids is tree-based AMR (quad/octree), as in the work of Popinet [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. Regarding unstructured meshes, Dai and Schmidt [START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF] proposed a moving-mesh algorithm to simulate free-surface flows with large deformations on adaptive tetrahedral grids. Both studies are proposed in the context of VOF methods, and the meshes are typically refined based on the local curvature of the interface. For classical Level Set methods on adaptive tetrahedral grids, Morgan and Waltz [START_REF] Morgan | 3D level set methods for evolving fronts on tetrahedral meshes with adaptive mesh refinement[END_REF] introduced a method to treat evolving fronts with specified velocity, while Rodriguez et al. presented a parallel two-phase flow solver using a finite-element formulation in [START_REF] Rodriguez | A parallel adaptive mesh method for the numerical simulation of multiphase flows[END_REF], where the adaptive mesh is refined based on the distance to the interface.

In this work, the Accurate Conservative Level Set (ACLS) method is used to capture the interface on 2D and 3D unstructured grids: this method accurately predicts the interface dynamics while conserving liquid mass [START_REF] Olsson | A conservative level set method for two phase flow[END_REF][START_REF] Olsson | A conservative level set method for two phase flow II[END_REF][START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF]. ACLS on unstructured grids has been proposed before in [START_REF] Kees | A conservative level set method suitable for variable-order approximations and unstructured meshes[END_REF][START_REF] Balcazar | A finite-volume/level-set method for simulating two-phase flows on unstructured grids[END_REF] on static meshes, and by Antepara et al. [START_REF] Antepara | Tetrahedral adaptive mesh refinement for two-phase flows using conservative level-set method[END_REF] for adaptive tetrahedral grids. In this latter, non-conforming AMR with an octree data structure is used. In the present work, massively distributed adaptive conforming unstructured grids are employed, and adapted using both metric-based and skewness-based quality criteria. We extend these previous unstructured ACLS studies by including an adapted and improved reinitialization method, first introduced by Chiodi and Desjardins for cartesian grids [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] and that has never been implemented to date for unstructured grids. In particular, interface-merging properties of this reinitialization equation are discussed. Moreover, a narrow band around the liquid-gas interface is built to compute interface features only in the area of interest [START_REF] Adalsteinsson | A Fast Level Set Method for Propagating Interfaces[END_REF], allowing significant computational savings.

This article is organized as follows: the various parts of the procedure implemented in the YALES2 unstructured low-Mach number code [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF] are first presented, namely the ACLS interface-capturing technique, the fractional-step integration of Navier-Stokes equations with explicit treatment of jump conditions, and the dynamic mesh adaptation strategy. The accuracy, spatial convergence, and computational efficiency of the unstructured ACLS method are assessed on classical interface transport test cases. The behavior of the overall algorithm in the case of canonical capillary-driven two-phase flows is then verified. Finally, two LES simulations in complex geometries are performed on adaptive tetrahedral meshes: the break-up of a turbulent liquid jet from a low-pressure compound nozzle [START_REF] Grout | Fractal analysis of atomizing liquid flows[END_REF][START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF] and of a non-reactive high-pressure liquid jet in crossflow are shown at various resolutions, demonstrating the accessibility of jet instabilities with 3D adaptive unstructured meshes. Results are validated against experiments performed in the same conditions [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF][START_REF] Ragucci | Trajectory and momentum coherence breakdown of a liquid jet in highdensity air cross-flow[END_REF].

Mathematical formulation 2.1 Accurate Conservative Level Set framework

In the Accurate Conservative Level Set (ACLS) framework [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF], the liquid-gas interface Γ is represented using a hyperbolic tangent profile:

ψ(x, t) = 1 2 tanh φ(x, t) 2ε + 1 (1)
where the parameter ε sets the thickness of the profile, and φ(x, t) = ±|x(t)x Γ (t)| is the signed-distance function.

Using ψ, the interface is located at the iso-level 1/2:

Γ(t) = {x ∈ R 3 | ψ(x, t) = 1/2} (2)
as depicted in the one-dimensional example of Fig. 1. In this figure, a characteristic profile thickness δ ψ is used, computed based on the maximum gradient of ψ:

δ ψ = 1 max |∇ψ| = 4ε (3)
Assuming the flow velocity field u is divergence free, the scalar ψ is advected by the fluid: Interface capturing on a one-dimensional example using the ACLS method. x Γ = 0; δ ψ = 4ε. and then reshaped using the reinitialization equation of [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF]:

∂ψ ∂t + ∇ • (ψu) = 0 (4)
∂ψ ∂τ = ∇ • ( ε(∇ψ • n)n Diffusion -ψ(1 -ψ)n Resharpening ) ( 5
)
where τ is a pseudo-time, and n is the interface normal:

n = ∇φ |∇φ| (6)
A new form of the reinitialization equation was recently proposed by Chiodi et al. in [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF]:

∂ψ ∂τ = ∇ •   1 4 cosh 2 φmap 2ε (|∇φ map • n| -1) n   (7)
where φ map = ε ln (ψ/(1ψ)) is an analytical signed-distance function, mapped for ψ ∈]0; 1[. More details regarding the history of conservative level-set reinitialization equations can be found in [START_REF] Mirjalili | Interface-capturing methods for two-phase flows: An overview and recent developments[END_REF].

The signed-distance function φ is reconstructed at nodes in the narrow band around the interface using a Geometric-Projection Marker Method (GPMM), firstly introduced in [START_REF] Janodet | An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows[END_REF], to estimate the smallest distance to the interface. The mean curvature κ is computed directly from φ using Goldman's formula [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF]:

κ = Tr (H(φ)) -∇φ T |∇φ| • H (φ) • ∇φ |∇φ| |∇φ| (8)
where H(φ) is the Hessian matrix of the signed-distance function.

The ACLS method is not strictly conservative for two reasons. First, the approximate and exact enclosed volumes of liquid are respectively:

V l = V ψ dV and V l = V H(ψ -1/2) dV ( 9
)
where H is the Heaviside function (0 in gas, 1 in liquid). In the ACLS method, the approximate volume V l is conserved, but differs from the exact volume V l due to interface curvature, i.e. in the ubiquitous case of non-straight interfaces.

A second reason for losing liquid mass on irregular meshes is the variation of the mesh resolution from fine to coarse, which leads to the smoothing of ψ, and may imply a loss of the ψ = 1/2 iso-surface. In this case, the interface is not detected anymore and gas can be set instead of liquid.

Finally, the reinitialization process on non-homogeneous meshes may also induce errors in the transport speed: the profile thickness parameter ε is mesh dependant, typically set to ∆x/2. Spatial variations of this parameter induce changes in the profile thickness, which in turn may induce errors on the transport velocity of the ψ = 1/2 iso-surface. Because of all these numerical errors on non-homogeneous meshes, special attention is needed when coupling the ACLS procedure with dynamic mesh adaptation, so that the interface never encounters strong cell-size gradients.

Navier-Stokes equations

The momentum equation for an incompressible Newtonian fluid is written in each of the phases as:

∂u ∂t + ∇ • (u ⊗ u) = - 1 ρ ∇p + ∇ • τ + g (10)
with ∇ • u = 0. The flow velocity is assumed continuous in both normal and tangential directions across the liquid/gas interface Γ, and hence there is no velocity jump: [u] Γ = 0. However, pressure and density are discontinuous across Γ [START_REF] Delhaye | Jump conditions and entropy sources in two-phase systems. Local instant formulation[END_REF]. Surface tension σ is embedded in the pressure jump [p] Γ , and, without Marangoni forces, the pressure jump is written as:

[p] Γ = p l,Γ -p g,Γ = σκ Γ + 2[µ] Γ n t • ∇u • n (11)
It is numerically treated using the Ghost-Fluid Method (GFM) [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF]. In Eq. [START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF], n is the interface normal from Eq. ( 6), κ Γ is the total interface curvature, p l,Γ and p g,Γ are the liquid and gas pressure at the interface, respectively. The material properties are reconstructed as:

ρ(x, t) = ρ g + (ρ l -ρ g )H(ψ(x, t) -1/2) and µ(x, t) = µ g + (µ l -µ g )ψ(x, t) (12)
for the density and the dynamic viscosity, respectively. The linear solver used to solve the Poisson equation for pressure is the Deflated Preconditioned Conjugate Gradient of [START_REF] Malandain | Optimization of the deflated Conjugate Gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF].

3 Numerical formulation for unstructured grids 3.1 Accurate Conservative Level Set algorithm

Narrow-band ACLS

In order to reduce computational costs, interface features such as signed distance, normals, curvature and reinitialization fluxes are computed in a restricted region around the interface, called "narrow band" [START_REF] Adalsteinsson | A Fast Level Set Method for Propagating Interfaces[END_REF]. Each node i belonging to the narrow band is identified with a signed integer b i (b > 0 in liquid and b < 0 in gas), depending on the node layer it is part of. The flagging of the narrow band is illustrated in Fig. 2

Signed-distance function reconstruction: the GPMM algorithm

To compute the interface normals (Eq. ( 6)) and curvature (Eq. ( 8)), a reconstructed signed-distance function φ is used. As detailed in [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF], using this smooth function instead of ψ for the computation of normal avoids spurious orientation of this normal vector, which is thus appropriate for the reinitialization equation. This distance reconstruction is challenging on unstructured meshes, and very few massively parallel algorithms are able to do it. In the work of Dapogny and Frey [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], the signed-distance function is reconstructed on unstructured meshes from a triangulation of the interface using the unsteady Eikonal equation. In particular, the reduction of triangulation-induced numerical errors when refining the mesh close to the interface is discussed and quantified.

In this work, the interface is also approximated by a triangulation before reconstructing the distance, but it has been chosen to proceed geometrically using a pair-based Fast-Marching Method (FMM) for its computation from the explicit surface. The methodology is the following: (i) a preliminary sorting of the node pairs in the narrow band is performed to be able to sweep from the close-interface band layers to the ones further; (ii) a segmentation (2D) or triangulation (3D) of the interface within grid cells is done; (iii) a Fast-Marching Method (FMM) is then performed along the edges to build a collection of the closest segments (2D) or triangles (3D) at each node of the band; (iv) the signed-distance function is computed from an exact geometric projection at each node of the band. All these steps constitute the Geometric-Projection Marker Method (GPMM), first introduced in [START_REF] Janodet | An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows[END_REF] and presented here in more details.

Sorting of the node pairs In order to propagate data efficiently from the first band level |b| = 1 to the furthest |b| = |b| max , sorting the pairs of nodes layer by layer in the narrow band is necessary. To do so, a pair-based band-layer index B ik is defined as follows:

B ik = 1 if b i b k = -1 min (|b i |, |b k |) + 1 if b i b k = -1 (13)
with 1 ≤ B ik ≤ |b| max + 1. This index is depicted in Fig. 2 (b).

Determination of intersection points

It is then needed to compute the position x P of the intersection points P between the interface Γ and the edges ik crossing the interface (B ik = 1). Given the conservative level-set functions ψ i and ψ k (Eq. ( 1)) at the nodes i and k belonging to a crossing pair ik, one can invert these functions to build two distance functions to the interface φ inv,i and φ inv,k :

φ inv,i = 2ε i atanh (2ψ i -1) φ inv,k = 2ε k atanh (2ψ k -1) (14)
The level-set position θ on a pair ik crossing the interface is then calculated as:

θ = φ inv,i φ inv,i -φ inv,k (15)
which means that θ = 0 and θ = 1 correspond to Γ at nodes i and k, respectively. Thus, if the pair of nodes ik verifies the condition (ψ i -0.5)(ψ k -0.5) < 0, the position of the intersection point is eventually computed by interpolating the nodal coordinates as:

x P = x i + θ(x k -x i ) = x i + θ∆x ik (16)
Segmentation/Triangulation of the interface In order to compute the distance to an explicit interface, noted Γ * in the following, a segmentation in 2D or a triangulation in 3D of the implicit interface Γ is first performed. This step is illustrated in 2D in Fig. 3. The concept of marker stems from this procedure. A marker M is an array containing the coordinates x K (M ) = (x K (M ), y K (M ), z K (M )) of its associated intersection points P K (M ) seen in Fig. 3. Hence, for 2D triangular grids, each interface segment is defined by a set of two intersection points (K = 1, 2) stored in a marker. For 3D tetrahedral meshes, there can be three or four intersection points depending on how the interface intersects a tetrahedron. For a three-point intersection, the interface within the cell is defined by a single triangle. In the case of a four-point intersection, the resulting quadrilateral is divided into four triangles to ease the storing of the surface elements, leading in both cases to three points stored in one marker in 3D (K = 1, ..., 3). Each node of the narrow band knows a maximum fixed number of markers N M defined by the user (I = 1, ..., N M ), which is usually set to 3 in 2D and 10 in 3D. Indeed, as studied in [START_REF] Janodet | An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows[END_REF], propagating several markers significantly improves the accuracy and spatial convergence of the method. A set of markers known by a node i of the narrow band is a marker list. Each marker M , which is a node-based array M i , also contains the absolute distance d(M i ) from the node i to the associated surface element, and two identificators id 1,2 (M i ). id 1 is the processor number and id 2 is the mesh element number in which the intersection points are. A 2D nodal marker list is presented in Table 1.

Table 1: A 2D marker list for a node i of the narrow band: 2 identificators id 1,2 , 4 coordinates

(K = 1, 2), 1 distance d. M 1,i ... M I,i ... M N M ,i id 1 (M 1,i ) ... id 1 (M I,i ) ... id 1 (M N M ,i ) id 2 (M 1,i ) ... id 2 (M I,i ) ... id 2 (M N M ,i ) x 1 (M 1,i ) ... x 1 (M I,i ) ... x 1 (M N M ,i ) y 1 (M 1,i ) ... y 1 (M I,i ) ... y 1 (M N M ,i ) x 2 (M 1,i ) ... x 2 (M I,i ) ... x 2 (M N M ,i ) y 2 (M 1,i ) ... y 2 (M I,i ) ... y 2 (M N M ,i ) d(M 1,i ) ... d(M I,i ) ... d(M N M ,i )
The markers will then be sorted by distance for each node i, within the marker list so that:

d(M 1,i ) < ... < d(M I,i ) < ... < d(M N M ,i ) (17)
Sorting of markers and distance computation When a marker is evaluated to update and fill marker lists, a check is first performed based on the two marker identificators to identify whether the marker is already present in the marker list or not. If not, the marker distance needs to be estimated to place the new marker in the list. Two cases are then considered based on the nodal band-layer index b i , and a threshold |b| exact = 2: (i) in the first two band levels, i.e. if |b i | ≤ |b| exact , the distance d is computed by performing an exact geometric projection; (ii) for the furthest band levels, i.e. if |b| exact < |b i | ≤ |b| max , the distance is initially approximated in order to speed-up the marker-sorting procedure. We focus on 2D triangular grids in the following for sake of simplicity.

In case (i), the geometric projection algorithm is used, illustrated in Fig. 4. If the projection is inside the element (Fig. 4 (a)), the distance to the segment [P 1 P 2 ] stored in marker M I,i at node i is computed geometrically as:

d(M I,i ) = (x i -x P1 ) • n P1P2 n P1P2 (18)
where n P1P2 is the normal vector to the interface segment [P 1 P 2 ]. Otherwise, if the projection is outside the element (Fig. 4 (b)), the distance d is the distance between the node i and the closest intersection point P K :

d(M I,i ) = min(|x i -x P1 |, |x i -x P2 |) (19)
For case (ii), the distance is simply approximated for the sorting as:

d(M I,i ) = min   x i -x P1 , x i -x P2 , (x i -x P1 ) + (x i -x P2 ) 2   (20)
Once the marker distance is estimated, it is then needed to sort the marker by distance and thus include the new one in the nodal marker list. This procedure is presented in Algorithm 1.

Algorithm 1 Update of a marker list at a node i of the narrow band. The new marker to be included in the list is M * ,i .

1: if d(M * ,i ) ≤ d(M 1,i ) then 2: for I = N M -1, 1, -1 do 3: M I+1,i = M I,i 4:
end for 5: end for 18: end if Optimization and parallelization The present algorithm is iterative for parallelism purposes. YALES2 parallelization uses a double-domain decomposition: within a processor, grid elements are assembled in groups. To optimize the algorithm, it is necessary to check local changes, within a cell group or a processor.

M 1,i = M * ,i 6: else if d(M * ,i ) ≥ d(M N M ,i ) then 7:
In addition to the marker list array, another array which contains marker flags F(M ) (0 or 1) is created before propagating data (i.e. after treating the first band level) in the narrow band. F i (M ) = 1 means that the marker M is new in the list at node i: e.g. this flag would be turned to 1 right after lines 5 and 14 of Algorithm 1. This operation allows to have updated information on the status of markers during data propagation, and thus avoids unnecessary repetition of computational operations (e.g. redundant updates and sorting of markers).

Another optimization is to store the squared distance d 2 (M I,i ) in marker lists instead of the distance d(M I,i ), to avoid many calculations of square roots, which have seen to be very expensive due to the large number of distance computations in the present algorithm. Indeed, this minor change does not affect the comparison and sorting of markers.

Final nodal signed-distance computation in the narrow band Once all markers are sorted, M 1,i contains the smallest squared distance d 2 (M 1,i ). The signed-distance is thus simply computed at node i as

φ i = ± d 2 (M 1,i ) if |b i | ≤ |b| exact ,
otherwise for the nodes where the distance d in markers has been only approximated for the sorting (|b| exact < |b i | ≤ |b| max ), the distance φ i is eventually computed using the projection algorithm presented above (Fig. 4) onto the closest surface element, known thanks to the marker sorting.

Finite-volume framework for unstructured grids: spatial discretization and notations

YALES2 is a pair-based finite-volume low-Mach solver with node-centered control volumes [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF]. These control volumes V i are based on the dual mesh, as shown in Fig. 5. The spatial discretization is ensured by 4th-order central finite-volume schemes, written on the dual mesh constructed from edge and element centers (dual cells connect centroids of cells to centroids of faces).

In this framework, second and fourth-order gradients of a quantity ζ are computed at a given node i using the gradient theorem as:

∇ζ| i = 1 V i k∈Ni ζ ik dA ik (21)
Figure 5: Finite-volume pair-based discretization in YALES2, using node-centered control volumes V i .

∇ζ| 2nd i = 1 V i k∈Ni ζ i + ζ k 2 dA ik ( 22
)
∇ζ| 4th i = 1 V i k∈Ni ζ i + ζ k 2 + ∇ζ| 2nd i -∇ζ| 2nd k 6 • ∆x ik dA ik ( 23
)
with N i the set of neighbors of node i, and ∆x ik the direction of the pair ik. It is useful to define the unit pair direction e ik = ∆x ik /∆x ik , which will be used in the following sections. More details about the high-order finite-volume framework used can be found in [START_REF] Bernard | A framework to perform high-order deconvolution for finite-volume method on simplicial meshes[END_REF].

Using the divergence theorem, the divergence of a flux F writes:

∇ • F| i = 1 V i k∈Ni F ik • dA ik (24)
with F ik the pair-based flux computed as in the gradient to either second or fourth order.

Interface normals and curvature computation

The interface normals and curvature are computed in the narrow band around the interface, using the reconstructed signed-distance function φ from the GPMM algorithm presented above.

The interface normals are pair-based, and computed as:

n ik = ∇φ GPMM | 4th ik = (∇φ GPMM | 4th i + ∇φ GPMM | 4th k )/2 (25)
where ∇φ GPMM | 4th i or k are fourth-order nodal gradients, computed using Eq. ( 23). The interface curvature κ Γ is computed in two steps: (i) the curvature is first estimated at nodes (κ i ) in the narrow band using Eq. ( 8); (ii) for the node pairs crossing the interface (b i b k = -1), a linear interpolation of the nodal curvature radii R i = 1/κ i is performed to obtain κ Γ .

Computation of curvature at nodes The local mean curvature at a node i is given by discretizing Eq. ( 8):

κ i = Tr (H i (φ GPMM )) - ∇φ T GPMM | 2nd i ∇φGPMM| 2nd i • H i (φ GPMM ) • ∇φGPMM| 2nd i ∇φGPMM| 2nd i ∇φ GPMM | 2nd i ( 26
)
where H i (φ GPMM ) is the discretized Hessian operator at node i applied to the reconstructed signed-distance function φ GPMM . Two discretizations of the Hessian operator are studied in this work: compact and noncompact. The compact Hessian computation necessitates only neighboring nodes (i.e. a compact stencil), and is written as:

H C i (φ GPMM ) = 1 V i k∈Ni     φ GPMM,k -φ GPMM,i ∆x ik e ik Aligned with ik + ∇φ GPMM | 2nd ik -∇φ GPMM | 2nd ik • e ik e ik Orthogonal to ik     ⊗ dA ik (27)
whereas the non-compact formulation is simply obtained by applying two consecutive second-order gradients, leading to a non-compact stencil:

H NC i (φ GPMM ) = ∇ ∇φ GPMM | 2nd 2nd i ( 28
)
Interpolation for crossing pairs Once the nodal curvatures are computed, a linear interpolation of the curvature radii is performed for the pairs ik crossing the interface Γ, resulting in a pair-based interface curvature:

κ Γ,ik = 1 θR i + (1 -θ)R k = κ i κ k θκ k + (1 -θ)κ i (29)
with θ the interpolation weight (0 < θ < 1), defined in Eq. ( 15), representing the position of the interface on the crossing pair ik.

Discretization of the reinitialization equation

The reformulated reinitialization (Eq. ( 7)) is discretized using Eq. ( 24):

ψ n * +1 i -ψ n * i ∆τ = 1 V i k∈Ni     1 4 cosh 2 φ n * map,ik 2ε ik ∇φ map | n * ik • n ik -n ik • n ik     n ik • dA ik ( 30
)
where φ map,ik = (φ map,i + φ map,k )/2 and ε ik = (ε i + ε k )/2 are linearly interpolated interface distance and profile thickness, respectively. ∇φ map | ik is computed using a compact term in the direction of the node pair as in Eq. ( 27): [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF] where ∇φ map | 2nd i or k are second-order nodal gradients computed using Eq. ( 22). The interface normals n ik are computed with Eq. ( 25). It has been chosen to use non-unit normals, to avoid spurious gradient orientations when several interfaces are close to each other, and thus improve the front-merging properties of the implemented reinitialization.

∇φ map | ik = φ map,k -φ map,i ∆x ik e ik + ∇φ map | 2nd i + ∇φ map | 2nd k 2 - ∇φ map | 2nd i + ∇φ map | 2nd k 2 • e ik e ik
Special attention has to be paid in regions of the narrow band where ψ is close to 0 and 1. In these regions, φ map is difficult to estimate and two improvements in the reinitialization algorithm are introduced: (i) in order to ensure a correct computation of ∇φ map | ik , which requires several nodal values of φ map , an under/overshoot tag is attributed to all pairs with at least one node concerned by ψ falling outside ]0, 1[; (ii) the resharpening term is removed so that the reinitialization becomes purely diffusive. Using the original isotropic diffusion term of Olsson and Kreiss [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] in this case, the reinitialization then takes the form:

ψ n * +1 i -ψ n * i ∆τ = 1 V i k∈Ni ε ik ∇ψ| n * ik • dA ik (32)
with ∇ψ| n * ik discretized as in Eq. [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF]. Once the overshoots or undershoots of ψ disappear, the reinitialization turns back to its main form given by Eq. [START_REF] Grout | Fractal analysis of atomizing liquid flows[END_REF].

Usually, 3 steps of reinitialization are performed and the pseudo-timestep of each reinitialization step is simply determined by a Fourier condition on the diffusive part of the reinitialization. This stability condition leads to:

∆τ = min pairs ik Fo.∆x 2 ik ε ik ( 33
)
where ∆x ik stands for the length of the pair ik, and Fo is the Fourier number set to 0.5.

Transport of the ACLS function -Time integration

The ACLS function ψ is advected using Eq. ( 4). In YALES2, a 4th-order discretization is performed for convective terms, leading to the following semi-discrete equation for a node i:

∂ψ i ∂t + ∇ • (ψu) i = 0 (34)
Then by applying Eq. ( 24):

∂ψ i ∂t + 1 V i k∈Ni ψ ik u ik • dA ik = 0 ( 35
)
Figure 6: Full ACLS procedure: the blue boxes represent computations in the narrow band.

Using a 4th-order spatial discretization as in Section 3.1.3, Eq. ( 35) yields:

∂ψ i ∂t + 1 V i k∈Ni ψ i + ψ k 2 + ∇ψ| 2nd i -∇ψ| 2nd k 6 • ∆x ik u i + u k 2 + ∇u| 2nd i -∇u| 2nd k 6 • ∆x ik •dA ik = 0 (36)
Time integration is performed using the TFV4A scheme of [START_REF] Kraushaar | Application of the compressible and low-Mach number approaches to Large-Eddy Simulation of turbulent flows in aero-engines[END_REF], which is a 4th-order, Runge-Kutta-like, time-marching scheme. The timestep ∆t is computed from the classical explicit CFL stability criterion:

∆t = min pairs ik CFL|dA ik |∆x ik |u ik • dA ik | ( 37
)
where ∆x ik is the length of the pair ik. The full ACLS procedure is depicted in Fig. 6.

Low-Mach flow solver

Ghost-Fluid Method

The Ghost-Fluid Method [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF][START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF] is a robust and accurate way to explicitly deal with jump conditions. Let us consider the jump at a liquid-gas interface in one dimension of a physical quantity ζ, written

[ζ] Γ = [ζ](x Γ ) = ζ l (x Γ ) -ζ g (x Γ ) = ζ l,Γ -ζ g,Γ .
The liquid and gas quantities in each phase can be extended in the other phase by continuity:

ζ l (x) = ζ l,Γ + (x -x Γ ) ∂ζ ∂x l,Γ + O(x -x Γ ) 2 and ζ g (x) = ζ g,Γ + (x -x Γ ) ∂ζ ∂x g,Γ + O(x -x Γ ) 2 (38)
which leads to the expression of the jump everywhere in the fluid domain:

[ζ](x) = [ζ] Γ + (x -x Γ ) ∂ζ ∂x Γ + O(x -x Γ ) 2 (39)
The extension to 3D is straightforward:

[ζ](x) = [ζ] Γ + (x -x Γ ) • [∇ζ] Γ + O(x -x Γ ) 2 (40)
Fig. 7 illustrates the principle of the Ghost-Fluid Method. In the following, this methodology is used to discretize the pressure gradient in the numerical framework of this study: pair-based finite-volume schemes. Discretization of the pressure gradient A first-order finite-difference approximation of the pressure term gives:

1 ρ ∂p ∂x l,i+ 1 2 = 1 ρ l p l,i+1 -p l,i ∆x + O(∆x) (41)
Using the pressure jump at x i , the liquid pressure is p l,i = [p] i + p g,i and can be replaced as follows:

1 ρ ∂p ∂x l,i+ 1 2 = 1 ρ l p l,i+1 -[p] i -p g,i ∆x + O(∆x) (42)
The expression for the pressure jump [p] i is given by Eq. ( 39):

[p] i = [p] Γ + (x i -x Γ ) ∂p ∂x Γ + O(x i -x Γ ) 2 (43)
[p] Γ is a known quantity (Eq. ( 11)) but [∂p/∂x] Γ is unknown a priori. With the assumption that the velocity is continuous across the interface, the pressure term in the momentum equation Eq. ( 10) should also be continuous:

1 ρ ∂p ∂x Γ = 0 ⇐⇒ 1 ρ l ∂p ∂x l,Γ = 1 ρ g ∂p ∂x g,Γ (44)
After many developments, which can be found in [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF], the final form of the discretized pressure gradient for finite differences is obtained, for Γ ∈ [x i ; x i+1 ]:

1 ρ ∂p ∂x l,i+ 1 2 = 1 ρ ∂p ∂x g,i+ 1 2 = 1 ρ * p l,i+1 -p g,i ∆x - 1 ρ * [p] Γ ∆x + O(∆x) (45)
with ρ * = θρ g + (1θ)ρ l a modified density, computed using the interpolation weight defined in Eq. [START_REF] Lin | A volume-conserving balanced-force level set method on unstructured meshes using a control volume finite element formulation[END_REF].

Implementation in the Poisson solver YALES2 is a pair-based finite-volume solver, as a result implementation of finite-difference schemes can be performed on the pairs. The expression of the discretized pressure term of the momentum equation written in Eq. ( 45) is integrated in the Poisson solver of YALES2 to account for interfacial jump conditions.

To solve the incompressible Navier-Stokes equations, a projection method is employed [START_REF] Chorin | The numerical solution of the Navier-Stokes equations for an incompressible fluid[END_REF]. The Poisson equation, used to enforce the solenoidal nature of the velocity field u n+1 , is:

∇ • 1 ρ ∇p = 1 ∆t ∇ • u * (46)
with u * the predicted velocity, estimated from u n (see next section 3.2.2). Eq. ( 46) can be discretized as presented in Section 3.1.3:

k∈Ni 1 ρ ∇p ik • dA ik = 1 ∆t k∈Ni u * ik • dA ik (47)
in which the pressure gradient is computed at pairs. A linear system of the form Ap = b that includes the jump conditions needs to be solved, and is written as:

k∈Ni 1 ρ * p k -p i ∆x ik e ik • dA ik Ap = 1 ∆t k∈Ni u * ik • dA ik + k∈Ni 1 ρ * [p] Γ ∆x ik e ik • dA ik b ( 48
)
The second part of the RHS b, which contains the pressure jump, is added for all pairs ik crossing the interface. The linear system is solved using the Deflated Preconditioned Conjugate Gradient of [START_REF] Malandain | Optimization of the deflated Conjugate Gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF]. The correction step is written as:

u n+1 -u * ∆t = - 1 ρ ∇p i (49)
in which the pressure gradient is calculated at nodes. For a pair of nodes which does not cross the interface (i.e. within one incompressible phase, liquid or gas), the pressure gradient is simply discretized as in Section 3.1.3, with a constant density ρ corresponding to the phase. For a crossing pair ik, the jump conditions need to be accounted for:

1 ρ ∇p 2nd i = 1 ρ i V i k∈Ni p i + p k + [p] k 2 dA ik (50)
where the pressure jump is expressed from the developments above using the Ghost-Fluid Method:

[p] k = ρ i ρ * [p] Γ + 1 - ρ i ρ * (p i -p k ) (51)

Velocity prediction with momentum correction

In order to avoid velocity mixing between liquid and gas, and thus limit spurious kinetic energy transfers at the interface, Desjardins and Moureau introduced a momentum correction at the interface in [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF]. In this section, we only consider the convection term in Eq. ( 10) for convenience. The prediction step writes:

u * -u n ∆t = -∇ • u n+1/2 ⊗ u n+1/2 (52)
which can be rewritten in terms of momentum:

ρ * * u * -ρ n u n ∆t = -∇ • ρ n+1/2 u n+1/2 ⊗ u n+1/2 (53)
where ρ * * is a predicted density computed from the following step:

ρ * * -ρ n ∆t = -∇ • ρ n+1/2 u n+1/2 (54)
Eq. ( 54) can be discretized in space as:

ρ * * i -ρ n i ∆t = - 1 V i k∈Ni ρ n+1/2 ik u ik • dA ik (55)
with ρ n+1/2 ik computed using a simple upwinded scheme to ensure the boundedness of the density flux at the interface:

ρ n+1/2 ik =      (1 -U ik ) ρ n i + ρ n k 2 + U ik ρ n i if u ik • dA ik ≥ 0 (1 -U ik ) ρ n i + ρ n k 2 + U ik ρ n k if u ik • dA ik < 0 (56)
where the pair-based limiter U is computed as:

U ik = max |b| upw -min (|b i |, |b k |) + 1 |b| upw , 0 (57)
with |b| upw = 5, meaning that the first five narrow-band levels from the interface are concerned by the upwinding.

A similar numerical strategy is applied for the conservative flux of Eq. ( 53):

ρ * * i u * i -ρ n i u n i ∆t = - 1 V i k∈Ni ρu n+1/2 ik u ik • dA ik (58)
with ρu n+1/2 ik calculated to 4th-order as:

ρu n+1/2 ik =        (1 -U ik ) ρ n i u n i + ρ n k u n k 2 + U ik ρ n i u n i + ρ n+1/2 ik ∇u| 2nd i -∇u| 2nd k 6 • ∆x ik if u ik • dA ik ≥ 0 (1 -U ik ) ρ n i u n i + ρ n k u n k 2 + U ik ρ n k u n k + ρ n+1/2 ik ∇u| 2nd i -∇u| 2nd k 6 • ∆x ik if u ik • dA ik < 0 (59)

Time integration -Full solution procedure

As in Section 3.1.6, time integration is performed using the TFV4A scheme of [START_REF] Kraushaar | Application of the compressible and low-Mach number approaches to Large-Eddy Simulation of turbulent flows in aero-engines[END_REF], which is a 4th-order, Runge-Kutta-like, time-marching scheme. The timestep ∆t is computed from the explicit capillary constraint:

∆t = min crossing pairs ik STN 2 ∆x 2 ik min(ρ g , ρ l ) σ|κ Γ,ik | (60)
where STN is the surface tension number, usually set to 0.5. The full solution procedure is here summarized:

(i) Using the unstructured ACLS algorithm, advance the interface implicitly through the ψ field with Eq. ( 36) using u n , then reinitialize ψ to reform the hyperbolic tangent profile with Eq. ( 30).

(ii) Set the new narrow band around the interface (Section 3.1.1).

(iii) Compute interface features in the narrow band: geometric signed distance φ GPMM (Section 3.1.2), normals n (Eq. ( 25)), curvature κ Γ (Eq. ( 26) then Eq. ( 29)).

(iv) Predict the velocity field u * from u n , using the momentum correction at the interface of [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF] adapted to unstructured meshes (Section 3.2.2).

(v) Project the velocity field by solving the Poisson equation, in which interface jump conditions are treated with the Ghost-Fluid Method (Eq. ( 48)).

(vi) Correct the velocity, with the pressure gradient that includes the jump conditions using the Ghost-Fluid Method (Eq. ( 49)), eventually leading to u n+1 so that ∇ • u n+1 = 0.

(vii) Grid adaptation (see next section 3.3).

Isotropic dynamic mesh adaptation strategy

To resolve the small physical scales at the interface on tetrahedral grids at a moderate cost, an isotropic Adaptive Mesh Refinement (AMR) technique is used [START_REF] Bénard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF], employing the MMG library [START_REF]Mmg Platform -Robust, Open-source, and Multidisciplinary Software for Remeshing[END_REF]. Special attention is needed when coupling AMR with ACLS to avoid liquid mass loss, so that the interface never encounters cell-size gradients. Thus, the mesh is refined around the interface based on the distance D = |φ| to it, so that Γ always stays in a protected region of constant cell size ∆x min (Fig. 8). It is important to stress that this metric evolution ensures that the remeshing process always occurs upwind and downwind of the interface but never at the interface itself. This strategy avoids any interpolation error at the interface location. The width of this region is 2N p ∆x min , with N p a user-defined parameter, usually set between 6 and 12. Away from this area, the metric evolution follows a linear law until it reaches the maximum cell size ∆x init , allowing only smooth and controlled cell-size gradients in the domain. The adaptation process is triggered automatically when the interface approaches the edge of the protected region, leading to considerable computational savings in remeshing. A 2D/3D static droplet of radius R d is initialized on homogeneous triangular/tetrahedral grids to assess the accuracy of the GPMM procedure, both in 2D and 3D, for the signed-distance reconstruction on irregular meshes. Errors are calculated as follows:

L 2 (φ) = 1 N i∈|b| N i∈|b| i=1 (φ th,i -φ GPMM,i ) 2 and L ∞ (φ) = max i∈|b| φ th,i -φ GPMM,i (61)
where |b| is the absolute narrow-band level, φ th,i is the theoretical distance at node i to the circular/spherical interface, and N i∈|b| is the total number of nodes in the corresponding absolute narrow-band level. The calculations of L 2 and L ∞ are performed for nodes belonging to edges crossed by the interface (|b| = 1) and in the sixth band level (|b| = 6, far from the interface). As homogeneous triangles/tetrahedra are used to tessellate the domain, the mesh cell size is computed as ∆x = √ V i in 2D, and as ∆x = 3 √ V i in 3D, where V i is the control volume based on the dual mesh around a node i. Results are shown in Fig. 9 and Fig. 10 for the 2D and 3D cases, respectively.

It can be seen in Fig. 9 (a) that in the first band level (|b| = 1), i.e. close to the interface, the GPMM procedure converges with second-order accuracy (L 2 and L ∞ ) with low error values, which is satisfactory for 2), a sufficient amount of data is available at a node far from the interface to compute the projected distance to the closest surface element. The advantage of the GPMM algorithm, which is an exact geometric projection, is that the numerical errors are only induced by segmentation (2D) or triangulation (3D): a perfect line (2D) or plane (3D) on an unstructured mesh lead to a distance computed within machine accuracy.

On the same test case, the accuracy of the interface curvature κ Γ (Eq. ( 29)) is evaluated in 2D and 3D using the following error norms:

L 2 (κ) = 1 κ Γ,th 1 N cp Ncp ik=1 (κ Γ,th -κ Γ,ik ) 2 and L ∞ (κ) = 1 κ Γ,th max cp κ Γ,th -κ Γ,ik (62)
where N cp stands for the total number of node pairs crossing the interface Γ, κ Γ,th is the constant theoretical curvature of the circle/sphere, and κ Γ,ik is the pair-based interface curvature, for a node pair ik among the crossing pairs, obtained with Eq. [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF]. Results are shown in Fig. 11, for the two different discretizations of the Hessian operator H (φ GPMM ) in Eq. ( 26). The non-compact discretization consists in applying two consecutive second-order gradients to the reconstructed signed-distance function φ GPMM (Eq. ( 28)), leading to a noncompact stencil, wheareas a smaller stencil is used in the compact formulation: only direct-neighbor points are used to compute the Hessian (Eq. ( 27)). Slight convergence of the L 2 -norm is observed for the 2D and 3D interface curvatures computed using the non-compact Hessian matrix. However, using a compact stencil leads to divergence of the curvature for both L 2 and L ∞ norms, in 2D and 3D. As a consequence, the non-compact formulation will be used in the following.

Interface transport tests

In order to assess the accuracy, spatial convergence and robustness of the ACLS procedure, classical interface transport tests are performed.

Static sphere on static tetrahedral grids

A 3D static sphere is initialized on a homogeneous tetrahedral mesh. The two reinitialization equations Eq. ( 5) [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] and Eq. ( 7) [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] are compared at two grid resolutions: 6 and 12 points in the sphere radius, and each of them is performed 50, 000 times throughout the runs. As homogeneous tetrahedra are used to tessellate the domain, the mesh cell size is computed as ∆x = 3 √ V i , where V i is the control volume based on the dual mesh around a node i. Results are presented in Fig. 12. As presented in [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] for cartesian meshes, it is observed that the reformulated reinitialization Eq. ( 7) induces much less artificial deformation of the interface at high levels of reinitialization, here on tetrahedral meshes. Thus, the conclusion of [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] can be extended to unstructured meshes.

Rotation of a notched disk on static triangular grids

In order to assess the robustness of the methods presented above for the transport of interfaces with topology that includes sharp edges and corners, the classical notched disk rotation case, introduced by S.T. Zalesak in [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF], is first run on homogeneous triangular grids. A two-dimensional disk of radius R = 0.15, with a notch of width w = 0.05, is initialized at (x, y) = (0, 0.25) in a [-0.5; 0.5] 2 square. Given the angular velocity Ω = 1rad/s and the origin of the domain as the center of rotation (x 0 , y 0 ) = (0, 0), the velocity field is:

u(x, y) = -Ω(y -y 0 ) Ω(x -x 0 ) = -y x (63)
Thus, the disk performs one rotation about the origin in T = 2π/Ω = 2πs. In the present paragraph, 10 rotations are done on triangular meshes with a CFL of 0.48, which in theory should not change the shape of the slotted disk. Three reinitialization steps are done at each temporal iteration and ε = ∆x/2. The narrow band is of width |b| max = 12 to ensure that |∇ψ| = 0 at the edge of the band and thus better preserve mass. Final interface contours after 10 rotations are presented in Fig. 13 for various mesh resolutions. The results are superimposed over a highly-resolved simulation (400 points in the disk radius). Good mesh convergence properties for interface dynamics are demonstrated: at the highest resolution shown (100 points in the disk radius), the initial shape is well retrieved. Fig. 14 shows two scalar fields, namely the signed narrow-band levels b and the reconstructed signed-distance function using the GPMM algorithm φ GPMM , after 10 rotations for R/∆x = 50. These fields are seen to be smooth and well computed. Fig. 15 depicts the node normal vector field around the interface in the same computational conditions: as detailed in Section 3.1.2, computing the normals from a signed-distance function avoids spurious normal orientation.

To quantitatively assess the accuracy of the method, we use a methodology inspired by [START_REF] Owkes | A mass and momentum conserving unsplit semi-Lagrangian framework for simulating multiphase flows[END_REF] to compute boundedness, shape, and mass errors at the end of the runs t = T :

E bound (T ) = max -min i=1...Ncv (ψ i (T )V i ), max i=1...Ncv ((ψ i (T ) -1)V i ) (64) E shape (T ) = Ncv i=1 ψ i (T ) -ψ i (0) V i (65)
ψ i (T )V i - Ncv i=1 ψ i (0)V i Ncv i=1 ψ i (0)V i (66)
These errors allow the estimation of the largest overshoot/undershoot in the domain, the discrepancies in the final shape compared to the initial shape, and the mass loss throughout the run, respectively. The computational efficiency is measured through a reduced computational time RCT:

RCT = WCT × N cores N iter × N cv (67)
with WCT the Wall Clock Time. The results are presented in Table 2. Although machine accuracy is not reached, the level of errors in mass remains very low for all resolutions. Moreover, neither undershoot nor overshoot are present at the end of the simulation (10 rotations) for half of the resolutions: some can be appear during these runs with an extremely small amplitude but are naturally corrected over time, as the interface topology becomes smoother. The shape error is plotted in Fig. 16: the method exhibits between first and second-order convergence on two-dimensional triangular meshes for this test case. The parallel efficiency of the conservative level set solver is demonstrated through the RCT values, which remain low for all calculations. 10 0 E shape (t = T )

Figure 16: Shape error E shape after 10 rotations of the 2D notched disk (t = T ), for various mesh resolutions: from 12.5 to 400 points in the disk radius. The dashed and solid lines are first and second-order convergence, respectively.

Circle deformation in a vortex on static triangular grids

A two-dimensional deformation of a disk in a vortex is then run [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF]. This test is used to evaluate the behaviour of the algorithms at resolution limit, which is useful for two-phase flow simulations, as thin ligaments are present.

A disk of radius R = 0.15 is initially centered at (x, y) = (0, 0.25) in a 2D box of dimensions [-0.5; 0.5] 2 . The disk is deformed by a swirling velocity field:

u(x, y, t) = -2 sin 2 (πx) sin(πy) cos(πy) cos(πt/T ) 2 sin 2 (πy) sin(πx) cos(πx) cos(πt/T ) (

The disk is first stretched for 0 ≤ t ≤ T /2, and then unstretched for T /2 ≤ t ≤ T . At t = T , the disk is theoretically back to its original circular shape. The total simulation time is chosen as T = 8s. The initial CFL number is 0.64 for all performed simulations of this paragraph. Three reinitialization steps are done at each temporal iteration and ε = ∆x/2. The narrow band is of width |b| max = 12. Maximum and final interface deformation, at t = T /2 and t = T respectively, are presented in Fig. 17 for various mesh resolutions. The results are superimposed over a highly-resolved simulation (400 points in the disk radius). Good mesh convergence properties are demonstrated: at the highest resolution shown (100 points in the disk radius), the original circular shape is perfectly retrieved. It should also be observed that at resolution limit, no locking phenomenon on the mesh nor artificial surface tension occur: the tail breaks up into small circles.

Using the same methodology than in the previous section, boundedness, shape, and mass errors (Eq. ( 64) to (66)) are computed at the end of the runs (t = T ). The results are presented in Table 3, with the reduced computational times of the runs (Eq. ( 67)). Again, the levels of error in bound and mass remain very low. Good mass conservation is ensured even at coarse grid resolutions. Convergence of the geometric error is displayed in Fig. 18: the method exhibits between first and second-order convergence on two-dimensional triangular meshes.

Sphere deformation in a vortex on adaptive tetrahedral grids

The previous test case can be extended to 3D [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF]. A sphere of radius R = 0. The simulation time is chosen as T = 3s. At t = T , i.e. after one full stretch/unstretch cycle, the initial sphere is supposed to be recovered. At maximum deformation (t = T /2), a thin sheet is formed and breaks up numerically into small spheres if the mesh is not fine enough. Adaptive tetrahedral meshes are used for this test case, dynamically and locally refined during the runs around the interface location, with a metric ∆x min . Here, it has been chosen to change the background mesh from one simulation to another so that ∆x max /∆x min = 5 for all runs. Larger ratios will be considered in the next sections of this article. For all calculations of this section, the refined region is of half-width 10∆x min (N p = 10, see Fig. 8), and AMR is triggered each time the interface is closer than 5 cells from the border of this fine region. The narrow-band width |b| max is set to 10, and the initial CFL number is 0.64. Interface contours are seen in Fig. 19 and contours along with the tetrahedral cell-size distribution in Fig. 20: good convergence is observed, and the final shape is almost retrieved for R/∆x = 60. The results are in accordance with the ones obtained for the two-dimensional case.

The computational efficiency, measured through the reduced computational time RCT (Eq. ( 67)), is presented in Table 4. The increase in computational cost remains low when increasing the resolution, and is mainly due to adaptive mesh refinement. The proposed ACLS procedure is scalable. The damping of a two-dimensional surface wave is then studied. Two superposed immiscible viscous fluids at rest are seperated by a flat interface, initially perturbed by a small-amplitude sinusoidal capillary wave. The initial signed-distance function is written as:

φ 0 (x, y) = -y -A 0 cos 2πx λ ( 70
)
in a two-dimensional [-λ/2; λ/2] 2 domain which is periodic in the x-direction, and with slip walls for top and bottom boundaries. The initial wave amplitude is chosen as A 0 = λ/100, with λ = 2π the wavelength of the perturbation. The initial perturbation is damped over time due to viscous dissipation. The results are compared to analytical solutions from Prosperetti's theory [START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF], which quantifies the viscous damping for vanishingly small amplitudes and infinite domains as an initial-value problem. The theoretical solution for two fluids with equal kinematic viscosity ν and λ = 2π is:

A th (t) = 4(1 -4β)ν 2 8(1 -4β)ν 2 + ω 2 0 A 0 erfc √ νt + 4 i=1 z i Z i A 0 ω 2 0 z 2 i -ν exp (z 2 i -ν)t erfc z i √ t (71)
where z i are the four roots -two pairs of conjugated complex numbers -of the following fourth-order polynomial equation:

z 4 -4β √ νz 3 + 2(1 -6β)νz 2 + 4(1 -3β)ν 3/2 z + (1 -4β)ν 2 + ω 2 0 = 0 ( 72
)
and

Z i = 4
j=1,j =i (z jz i ). In the case of equal density ρ, β = 1/4 and the normal-mode oscillation frequency ω 0 is computed from the dispersion relation as follows:

ω 2 0 = σ 2ρ ( 73
)
The non-dimensional time is t * = ω 0 t. The calculations are performed on three homogeneous triangular grids (λ/∆x = 30; 100; 200) up to T * = ω 0 T = 25, which corresponds to approximately 4 oscillations. The unity density ratio case is investigated, gravity is neglected, and 3 steps of reinitialization are performed at each timestep. The physical and numerical parameters of the runs are summarized in Table 5. The temporal evolution of the wave amplitude over time is depicted in Fig. 21 (a). The damping is correctly reproduced from 100 points in the wavelength λ in the case of equal densities. Fig. 21 (b) shows the instantaneous relative amplitude error |A(t) -A th (t)|/A 0 and confirm the convergence properties. To quantify further the differences between theory and numerical simulations, the following L 2 norm is computed, as in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]:

L 2 = 1 λ 1 T * T * t * =0 (A(t * ) -A th (t * )) 2 dt * (74)
where A and A th are the computed and theoretical wave amplitudes at a given instant, respectively. The L 2 -norm values and the computational costs of the runs are presented in Table 6. Again, the results show fair convergence and accordance with theoretical data in the case of unity density ratio. It validates the overall algorithm robustness and physical meaning in the case of surface tension/viscosity interactions on unstructured grids.

Spurious currents on static tetrahedral grids

This elementary test is used to assess the numerical errors induced by the surface tension computation procedure. For a static droplet, it aims to show the capability of a method to predict Laplace's law, i.e. the balance between surface tension and pressure forces. For a droplet in equilibrium, the theoretical velocity is 0 and thus nonphysical velocities indicate the presence of spurious (or parasitic) currents. Exact balance (to within machine accuracy) has been shown first on cartesian grids in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. Most past studies have been conducted with low Laplace numbers. In this work, we study spurious currents to check the robustness and assess the numerical errors of the algorithms in physical conditions close to the ones found in the range of the targeted practical applications, i.e. realistic liquid injection systems. To this mean, a high Laplace number is investigated: La = 2, 000, 000.

To study their magnitude over time, the maximum Capillary number is used: Ca max = µ l u max /σ. Time is normalized using a visco-capillary timescale: t * = t/τ , with τ = µ l D/σ. A three-dimensional droplet of diameter D = 2mm is initialized in a unit centimetric square box. The simulations are performed on three static tetrahedral grids of resolutions ∆x = 200µm, 100µm, 50µm, up to t * = 10, 000. Three reinitialization steps are performed at each iteration, and the explicit capillary time-step constraint STN in Eq. ( 60) is set to 0.71 for all resolutions. The physical and numerical parameters are presented in Table 7.

Table 7: Physical and numerical parameters for the measurements of spurious currents. ) are seen for all resolutions, and a significant reduction of the spurious current magnitude is observed for D/∆x = 40. The spatial convergence of the spurious currents, which dominate the solution in such physical conditions, is displayed, and the presented procedure is demonstrated to be highly robust on unstructured grids. µ l /µ g = 50.

ρ l /ρ g µ l /µ g La = ρ l σD/µ 2 l D/∆x

Rayleigh-Plateau instability

In this section, following the studies of [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Dai | Adaptive tetrahedral meshing in free-surface flow[END_REF][START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF][START_REF] Menard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet[END_REF], numerical predictions of the capillary instability are performed on 3D static unstructured meshes and validated against Weber's linear theory [START_REF] Weber | Zum Zerfall eines Flüssigkeitsstrahles[END_REF]. A static column of water of density ρ l , dynamic viscosity µ l , and surface tension σ in quiescent air is considered. This column is infinite in the x-direction and its radius is r 0 . A small axisymmetric perturbation wave is imposed at the column surface. The solution is of the form:

r(x, t) = r 0 + r 0 exp (ωt + ikx) (75)
where is the perturbation, assumed to be small, k = 2π/λ is the wavenumber, and ω is the growth rate. Using the non-dimensional wavenumber ξ = kr 0 and growth rate ω 0 = σ/(2ρ l r 3 0 ), the equation for the growth rate in Weber's theory applied to static jets, for ξ < 1, is given by:

ω 2 + 3µ l ξ 2 ρ l r 2 0 ω = ω 2 0 ξ 2 (1 -ξ 2 ) (76)
The mathematical developments leading to this result can be found in [START_REF] Sterling | The instability of capillary jets[END_REF].

The simulations are performed in a 3D, periodic in the x-direction, cylindrical domain of length [0; λ] and radius [0; 3r 0 ]. The signed-distance function is initialized as:

φ 0 (x, y, z) = r 0 -r 0 cos (kx) -y 2 + z 2 (77)
where r 0 = (1/3) × 10 -3 m, and the initial velocity is set to 0. In order to keep the interface free from spurious signals, the initial perturbation is set to be of the same order of magnitude than the homogeneous tetrahedralcell size: r 0 ∼ ∆x, here = 0.06. The runs are performed at various ξ values, ranging between 0 and 1 (unstable oscillatory solutions), and using an effective resolution of r 0 /∆x = 6. The definition of the cases can be found in Table 8. For all cases, the narrow-band width is |b| max = 12 and 3 steps of reinitialization are 9. The Laplace number is La = ρ l r 0 σ/µ 2 l = 18, 530. The interface contours for ξ = 0.52 (Case 3) can be seen in Fig. 23: implicit interface topology changes are correctly predicted. The relative deformation based on the maximum radius is defined at x = λ/2 as: Table 9: Physical parameters for the Rayleigh-Plateau instability simulations. ρ l /ρ g µ l /µ g La 829 62 18,530

Figure 23: Rayleigh-Plateau instability on a homogeneous tetrahedral mesh: interface contours at various non-dimensional times t/t σ . La = 18, 530. ξ = kr 0 = 0.52. r 0 /∆x = 6.

t σ = 1/ω 0 . r max -r 0 r 0 = exp (ωt) (78)
The characteristic time is t σ = 1/ω 0 = 1.013ms. Using t σ , Eq. ( 78) becomes:

ln r max -r 0 r 0 = ω ω 0 t t σ (79)
This term is linear in time and the derivative is the non-dimensional growth rate. However, due to the use of fully 3D tetrahedral meshes, several causes of inaccuracy might deteriorate the quality of the computed growth rate: (i) the jet can be slightly off-centered in the radial direction, (ii) the maximum amplitude can be strongly shifted from the middle of the jet in the longitudinal direction. While problem (ii) is solved by choosing an initial perturbation of the same order of magnitude than the cell size, problem (i) requires a dedicated methodology.

To post-process the results, the jet is discretized into rings of chosen thickness δ = 80µm, and area A. For each ring i of each temporal solution n, the mean radius is computed from the mean center:

r n i = 1 A A (y -y) 2 + (z -z) 2 dA ( 80
)
where y = (1/A) A y dA and z = (1/A) A z dA. The maximum mean radius r n max is then computed for each temporal solution n. The temporal evolution of the relative deformation for ξ = 0.52 (Case 3) is plotted in Fig. 24 (a). Despite a non-linear phase at the beginning of the calculation, remaining until the disturbance has grown quite large, the non-dimensional growth rate ω/ω 0 is correctly predicted. To quantify the growth rate only in the linear part, a temporal mean is performed from t/t σ = 4 to the breakup. The standard deviation SD of the growth-rate signal is also computed within the range considered:

SD = 1 N -1 N -1 n=1 ω ω 0 n - ω ω 0 mean 2 ( 81
)
where N is the total number of temporal solutions. The theoretical and numerical growth rates for all computed non-dimensional wavenumbers ξ are presented in Table 10 and depicted in the dispersion diagram of Fig. 24 (b). Larger standard deviations are obtained in the second and sixth cases (ξ = 0.26 and ξ = 0.90 respectively): the linear regime is more oscillatory, less clear to identify, in these two conducted simulations.

3D droplet collision using adaptive mesh refinement

This section aims to validate the unstructured discretization of the reformulated reinitialization equation in the case of a interface merging scenario. This is a non-trivial case: a bad prediction of the resharpening term leads to the locking of the interface on the grid, and as a consequence keeps the two interfaces of the two droplets from The physical parameters involved in this process are the Weber number We = ρ l u 2 r D/σ, with u r the relative velocity of the two droplets; the Ohnesorge number Oh = µ l / √ σρ l D; the non-dimensional impact parameter X, which is 0 for a head-on collision; and the droplet diameter ratio ∆. The values are presented in Table 11 and correspond to the experiment of [START_REF] Ashgriz | Coalescence and separation in binary collisions of liquid drops[END_REF], in the reflexive separation regime. The narrow-band width is |b| max = 10, and the minimum cell size ∆x min = 6µm is imposed in a refined region of half-width 10∆x min (N p = 10, see Fig. 8). The mesh adaptation process is triggered automatically when the interface is 5 cells away from the border of the refined region, which leads to massive computational savings. Three steps of reinitialization are performed at each iteration.

Instantaneous interface positions and meshes are presented in Fig. 25, and show good interface quality and accuracy, even after a very large number of iterations (∆t ∼ 0.01µs). Despite the presence of an entrapped bubble due to the lack of perfect symmetry with the unstructured grid, implicit topology changes are shown, and fair overall dynamics is displayed compared to the experimental results of [START_REF] Ashgriz | Coalescence and separation in binary collisions of liquid drops[END_REF]. The computational efficiency is presented in Table 12. This study follows the experiments of [START_REF] Grout | Fractal analysis of atomizing liquid flows[END_REF] and the simulations of [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF]. A single low-pressure compound nozzle is considered: three shifted cylinders are superimposed, as depicted in Fig. 26, from which water is injected in quiescent air at Reynolds number Re = 3653, based on the discharge orifice diameter d inj = 180µm, and liquid Weber number We l = 1061. The shift of the disks induces non-axial velocity components, and a double-vortex flow is formed, as illustrated in Fig. 27. In these physical conditions, primary atomization is mainly driven by capillary effects such as ligament breakup and jet forming by bubble burst, and enhanced by liquid turbulence. The issued water jet is slow, and hence aerodynamic forces are negligible. As the injector geometry strongly conditions the primary breakup process, the unsteady resolution of the internal flow is important, as mentioned in [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF].

In this work, the whole flow is computed by the mean of one Large-Eddy Simulation (LES), from the internal flow to the primary atomization process. The LES turbulence model is the dynamic Smagorinsky model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]. Adaptive tetrahedral meshes are used and two resolutions are investigated: ∆x min = 5µm and ∆x min = 2.5µm. For the two simulations, the refined region is of half-width 10∆x min (N p = 10, see Fig. 8), and AMR is triggered each time the interface is 5 cells away from the edge of this region. The jet morphology and angle are studied, and compared to the ones of [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF].

Internal flow

As presented above, the adaptive mesh refinement strategy is based on interface displacement, and thus the mesh size inside the injector, where only liquid is present, is the same for the two conducted calculations: the y + field represented at the walls is shown in Fig. 28. Converged temporal statistics of the flow velocity projected in the nozzle exit plane are shown in Fig. 29: the double-vortex secondary flow of the discharged disk (Disk 3 of Fig. 26) is clearly displayed.

External flow

A comparison of interface topologies and vorticity-dominated regions at the two simulated resolutions are shown in Fig. 30. A clear difference is seen in the appearance of jet instabilities: a closeup on some capillary effects, especially the Rayleigh-Plateau instability is shown for the 2.5µm case in Fig. 31. It depicts the temporal evolution of a ligament, circled in red in Fig. 31 (a), which is pinched off in Fig. 31 (b) until breakup (c) to produce two droplets. These outcome drops later collide and merge. These observations are simply not accessible with ∆x min = 5µm. The jet topology in the 2.5µm case is close to the one obtained by a DNS simulation at 1.44µm, performed with a symmetry plane in [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF]. This similarity is encouraging to target an even higher resolution. Fig. 32 shows the mesh cell-size distribution for both computed resolutions. The background mesh is kept identical for both, and it does not affect the quality of the solution: only smooth cell-size gradients are present in the domain, even with a higher cell-size ratio ∆x max /∆x min = 22.8 in the 2.5µm case.

Eventually, the spray angle α is estimated at various distances z/d inj away from the nozzle exit for d inj /∆x min = 72, in order to compare with the simulations and experiments of [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF]. An image processing technique has been employed on 250 frames, corresponding to 250µs of physical time, of the projected shadow of the spray to obtain a single mean image: Fig. 33 (a). The angle values are presented in Fig. 33 (b): the angles are rather overestimated, especially close to the nozzle, compared to the experimental values presented in [START_REF] Cousin | Primary breakup simulation of a liquid jet discharged by a lowpressure compound nozzle[END_REF], which lie between 20 • and 25 • . This could be explained by the lack of atomization close to the nozzle in our numerical simulations. Nevertheless the overall behavior is correct, despite a more significant decrease of the angle with the distance from the nozzle exit.

Computational performances

The computational cost and performances of the two conducted simulations on adaptive tetrahedral grids are shown in Table 13. Good scalability of the unstructured ACLS/AMR procedure is displayed. The cost of the simulations is moderate: in comparison, the same run at d inj /∆x min = 72 on a static grid requires 8192 cores and 1.62B cells. Hence, the present method allows huge computational savings for simulations on tetrahedral grids. a turbulent boundary layer along a flat plate. Its thickness is calculated following a 1/7th power law. The outer part of the velocity profile follows a parabolic law. The injected mass flow rate is the same one as reported in the experiments. The initial mesh is the same for all simulations: it consists of 32M tetrahedral elements (N cells =32M). The number of cores N cores allocated to each simulation ensures a constant ratio N cells /N cores between 100,000 and 150,000. For all simulations, the refined region is of half-width 10∆x min (N p = 10, see Fig. 8), and AMR is triggered each time the interface is 4 cells away from the edge of this region. Each case will yield different meshes as function of the operating condition and mesh resolution, see Table 16. Therefore, the simulations will have different costs, see Table 17. Physical running times simulated depend on the computational resources available.

Vertical trajectory and experimental validation

Numerical simulations are validated with the experimental correlation obtained by [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF] for the vertical jet penetration. It has been obtained by testing experimentally several operating conditions, and has a standard deviation of value 0.81. The correlation corresponds to the trajectory of the jet's windward side, and is given by:

z d inj = 1.57q 0.36 ln 1 + 3.81 x d inj (83)
This trajectory is valid up to the streamwise position x/d inj = 22 and depends solely on the q factor, as the authors state that the influence of the Weber number on the trajectory is negligible. Other experimental studies, such as the one by Ragucci [START_REF] Ragucci | Trajectory and momentum coherence breakdown of a liquid jet in highdensity air cross-flow[END_REF], provide experimental correlations with a dependence on We. While the role of We on the trajectory is still an open question, the results of our simulations are only compared to the ones of [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF] (Eq. ( 83)).

The resulting mean trajectories are shown in Fig. 35. The experimental correlation is given by the black solid line, and the shadowed area denotes its confidence region. The red and blue lines are the numerical mean trajectories for the coarse and fine mesh resolutions, d inj /∆x min = 22.5 and d inj /∆x min = 45, respectively. They are obtained by averaging the instantaneous trajectories of the windward side. Fig. 35 (a) shows the first operating point from Table 14. It can be seen that the numerical penetration in the near-nozzle region is slightly underestimated for both mesh resolutions. However, as the jet moves downstream its penetration fits into the experimental limits. The effect of mesh resolution on the trajectory is not very significant for this operating point. The second operating point, Fig. 35 (b), shows a larger difference in trajectory with the minimum cell size. For the coarsest resolution d inj /∆x min = 22.5, the trajectory is underestimated and does not fall within the confidence interval of the correlation, although its tendency is similar to the experimental results. For the finest resolution d inj /∆x min = 45, the trajectory follows very closely the experimental correlation. This discrepancy with respect to the other trajectories is due to liquid separation at the nozzle's convergent exit section in the fine cases. The effective area of the jet is reduced when injected to the plenum and, due to mass conservation, its velocity increases. Consequently, the jet penetrates further than the one of the coarse mesh, where liquid separation is not observed. Fig. 36 shows clearly this phenomenon. This difference is attributed to the resolution of the detached boundary layer within the nozzle by the fine mesh, which is not resolved by the coarse one. For the operating point at lower We g , liquid separation is obtained for both mesh resolutions. This can be explained by a lower liquid Reynolds number for this operating point, therefore the boundary layer is thicker and can be captured by coarse mesh. These observations, where solutions with liquid detachment match their trajectories to the experimental one, suggest that this phenomenon was actually occurring in the experiments of [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF].

The accuracy of the mean numerical trajectories can be quantitatively assessed by defining a L 2 error as in Eq. (84):

L 2 = 1 N N i=1 z d inj exp,i - z d inj num,i 2 (84)
where N is the total number of sample points along the trajectory coordinate x/d inj . The values for each case are shown in Table 15. As also observed in Fig. 35, the two resolutions for the operating point at low We g have very similar errors. For the high We g simulations, the coarse mesh shows a larger error, while the fine mesh yields the lowest error of all the simulations.

Jet topology and breakup

The effect of mesh resolution (see Table 16) is shown in the instantaneous snapshots of Fig. 36 for the operating point We g = 1470. The liquid-gas interface and the mesh at the plane y = 0 are displayed. Both pictures show the column and surface breakup features characteristic of jet in crossflow atomization [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF]. The interface cell size has a direct effect onto breakup and atomization. Both meshes capture different column breakup modes: the coarse mesh shows Kelvin-Helmholtz instabilities developing in the windward side of the jet, while the fine mesh displays Rayleigh-Taylor instability waves. This directly affects primary atomization, as it can be seen from the different liquid ligaments and structures that result from the breakup of the liquid dense core. The detachment of the liquid jet inside the nozzle observed for the finest resolution also affects the development of surface instabilities that will eventually contribute to breakup. Subsequent secondary atomization is also strongly affected by the resolution, as the finer mesh can capture smaller droplets than the coarse mesh (see Table 16).

Regarding the operating point, Fig. 37 shows instantaneous snapshots of the liquid-gas interface for the finest case (d inj /∆x min = 45) for both values of We g . The same time instant than in Fig. 36 is displayed. By comparing the operating points, it is seen that both cases present similar column breakup features following Rayleigh-Taylor instabilities. The simulation with We g = 830 shows that ligaments formed during primary atomization stretch and form lobes more pronounced than for We g = 1470. Furthermore, the number of wavelengths until breakup is lower for We g = 1470 than for We g = 830. This causes an earlier rupture of the liquid ligaments.

Computational performances

Table 16 shows the number of cores, tetrahedral elements and droplets generated at t = 0.3ms for all simulations. By comparing the mesh resolutions for each operating point, it is observed that six times more droplets are generated when the minimum cell size is reduced twice. Thus, as more droplets are resolved, more cores are needed to keep the ratio N cells /N cores constant for optimal performances: when the resolution increases twice, the number of elements rises near tenfold. When the operating points are compared at identical mesh resolutions, 50% more droplets are present in the highest gaseous-Weber case. This is due to a higher crossflow velocity, which increases the number of droplets generated by surface breakup.

The computational performances and costs are shown in Table 17. The Reduced Computational Times (RCT) are calculated according to Eq. ( 67). In all cases, the time spent in the ACLS procedure is almost constant and moderate despite the large liquid-gas surface to be resolved, demonstrating the high scalability of the interface-capturing algorithm. Most of the total computational time is spent in the AMR routine: indeed, the jet velocity is high and hence the AMR is automatically triggered much more frequently than in the lowpressure jet simulations of Section 5.1. In general, jet in crossflow simulations are more expensive than the Figure 36: Lateral view of meshes and interface contours near the injector at time instant t = 0.3ms. q = 6.

We g = 1470.

Figure 37: Instantaneous front (left) and lateral (right) views at time instant t = 0.3ms for the two operating conditions. d inj /∆x min = 45. q = 6. We g = 830 (top), We g = 1470 (bottom).

low-pressure jet ones (see Table 13 for comparison): while this latter jet is governed by capillary forces and aerodynamic effects are negligible, the jet in crossflow atomization is driven by aerodynamics and inertia (the limiting timestep of the simulations is always due to the CFL number).

It can be concluded that most of the computational resources allocated to the jet in crossflow cases are used for solving the AMR routine, making these computations highly expensive for cases where very small droplets are generated and a higher resolution of the atomization is required. A promising solution to this problem is coupling the presented methodology to a Lagrangian solver that could identify spherical droplets and convert them to Lagrangian Point Particles (LPP). This Eulerian-Lagrangian two-way coupling would allow spherical droplets to move following a ballistic motion and break up according to secondary atomization models. In this way, droplets would not be resolved with the Eulerian mesh by the ACLS/AMR methodology, saving huge computational resources.

Conclusion

An efficient and robust overall procedure for the computation of two-phase flows on unstructured meshes has been presented. Spatial convergence and accuracy of the interface capturing technique have been first demonstrated through well-known interface transport tests, both in 2D and 3D, and confirmed by the results obtained on canonical two-phase flow examples. Capillary-driven flows are correctly computed on triangular and tetrahedral grids, despite the challenges that poses the use of these irregular meshes. The accuracy and convergence of the numerical framework presented in this paper can be further improved by using the higher-order finitevolume schemes of [START_REF] Bernard | A framework to perform high-order deconvolution for finite-volume method on simplicial meshes[END_REF], which have proven to upgrade the operator discretization on unstructured grids and hence the interface curvature computation. Adaptive mesh refinement allows significant computational savings, especially in the case of complex turbulent 3D flows: highly-resolved flows with large density ratios in complex geometries using unstructured meshes are now accessible. The results obtained in the low-pressure compound nozzle case demonstrate the potential of the method to access a deep numerical insight of jet instabilities and internal flow dynamics with 3D unstructured meshes. This is a major progress for the computation of realistic industrial flows. To further increase the computational gain, especially in the liquid jet in crossflow simulations where many droplets need to be resolved, implementing an Eulerian-Lagrangian two-way coupling to treat small droplets as lagrangian point particles (LPP) would allow local coarsening of the mesh. This coupling would thus result in much cheaper simulations at equivalent resolutions. Eventually, this multiscale approach would pave the way to both extremely-resolved liquid jet cores using unstructured meshes and accurate droplet statistics far from the injector.

Figure 1 :

 1 Figure1: Interface capturing on a one-dimensional example using the ACLS method. x Γ = 0; δ ψ = 4ε.

  Figure 2: (a) Flagging of the narrow band around the interface Γ, here ik is a pair of nodes crossing the interface; (b) Pair-based band-layer index B ik for data propagation in the Geometric-Projection Marker Method (GPMM). methodology, an edge ik crossing the interface verifies b i b k = -1. In practice, the band width is set between 8 and 12 cells on each side of the interface: 8 ≤ |b| max ≤ 12.

Figure 3 :

 3 Figure 3: Segmentation of the interface in 2D. P K (M I ) is the K-th intersection point associated with the I-th marker M I .

Figure 4 :

 4 Figure 4: Projection algorithm to compute the exact geometric absolute distance d to the triangulated interface Γ * .

return 8: else 9 :

 9 for I = 1, N M -1 do 10: if d(M * ,i ) < d(M I+1,i ) then 11: for I = N M -1, 1 + I, -1 do 12: M I +1,i = M I ,i

Figure 7 :

 7 Figure 7: Ghost-Fluid Method: extension of the jump of a physical quantity ζ.

Figure 8 :

 8 Figure 8: Adaptive mesh refinement strategy: refinement based on the distance D to the interface.

Figure 9 :Figure 10 :

 910 Figure 9: L 2 and L ∞ of the reconstructed signed-distance function φ with the 2D GPMM algorithm (Geometric-Projection Marker Method), in first (a) and sixth (b) narrow-band levels |b|. The dashed and solid lines are first and second-order convergence, respectively.

Figure 11 :

 11 Figure 11: L 2 and L ∞ norms of the 2D (a) and 3D (b) interface curvature κ Γ , considering two spatial discretizations of the Hessian operator: with compact and non-compact stencils.

Figure 12 :

 12 Figure 12: 3D static sphere transport test on two tetrahedral grids: comparison between Eq. (5) (original ACLS reinitialization) and Eq. (7) (reformulated reinitialization) after 50,000 iterations of reinitialization.

Figure 13 :

 13 Figure 13: Interface contours after 10 rotations for various mesh resolutions. The thin line represents a simulation performed with 400 points in the disk radius.

Figure 14 :

 14 Figure 14: Scalar fields after 10 rotations with R/∆x = 50. (a) Signed narrow-band levels b; (b) Reconstructed signed-distance function φ GPMM . Black contours represent the interface.

Figure 15 :

 15 Figure 15: Node normal vector field n i = ∇φ GPMM

Figure 17 :

 17 Figure 17: Interface contours at maximum deformation (t = T /2; top) and final shape (t = T ; bottom) for various mesh resolutions. The thin line represents a simulation performed with 400 points in the disk radius.

Figure 18 :

 18 Figure18: Shape error E shape at the end of the 2D circle deformation in a vortex test (t = T ), for various mesh resolutions: from 12.5 to 400 points in the disk radius. The dashed and solid lines are first and second-order convergence, respectively.

Figure 19 :

 19 Figure 19: Interface contours at maximum stretch (t = T /2; top) and final shape (t = T ; bottom) for various mesh resolutions.

Figure 20 :

 20 Figure 20: Tetrahedral cell size distribution at maximum stretch (t = T /2) for various mesh resolutions.

  15 is initialized at (x, y, z) = (-0.15, -0.15, -0.15) in a [-0.5, 0.5] 3 three-dimensional box. The advection velocity is specified as: u(x, y, z, t) =   2 sin 2 (πx) sin(2πy) sin(2πz) cos(πt/T ) sin(2πx) sin 2 (πy) sin(2πz) cos(πt/T ) sin(2πx) sin(2πy) sin 2 (πz) cos(πt/T )

Figure 21 :

 21 Figure 21: (a) Temporal evolution of the dimensionless wave amplitude A/λ over non-dimensional time t * = ω 0 t; (b) Instantaneous relative amplitude error |A -A th |/A 0 .

  are presented in Fig.22. Reasonably low errors (Ca max ∼ 10 -3

= 10 D/∆x = 20 D/∆x = 40 Figure 22 :

 10204022 Figure22: Spurious currents around a static droplet on tetrahedral meshes. La = 2, 000, 000. ρ l /ρ g = 1000.µ l /µ g = 50.

Figure 24 :

 24 Figure 24: Rayleigh-Plateau instability. La = 18, 530. r 0 /∆x = 6. (a) Temporal evolution of the relative deformation based on maximum radius for ξ = kr 0 = 0.52, t σ = 1/ω 0 . (b) Dispersion diagram; the error bars represent the standard deviation SD of the growth-rate signal ±SD (see Eq. (81)).

Table 11 :Figure 25 :

 1125 Figure 25: Water/Air droplet collision. We = 23; Oh = 0.0047; X=0; ∆ = 1; D/∆x min = 133.

Figure 26 :

 26 Figure 26: Compound nozzle geometry: triple-disk injector.

Figure 27 :

 27 Figure 27: Illustration of the double-vortex internal flow in the cavity disk (Disk 2 in Fig. 26): velocity magnitude along with streamlines.

Figure 28 :

 28 Figure 28: y + distribution represented at boundaries. d inj = 180µm. Re = 3653. We l = 1061.

Figure 29 :

 29 Figure 29: Velocity statistics: (a) Mean and (b) RMS velocity projected in the exit section. d inj = 180µm. Re = 3653.

Figure 30 :

 30 Figure 30: Interface contours and Q-criterion of the water jet atomizing from a low-pressure compound nozzle in quiescent air at the two investigated resolutions: (a) ∆x min = 5µm and (b) ∆x min = 2.5µm.

Figure 31 :

 31 Figure 31: Closeup on some capillary effects during the atomization process of a water jet in quiescent air from a low-pressure compound nozzle. d inj = 180µm. Re = 3653. We l = 1061. ∆x min = 2.5µm.

Figure 32 :

 32 Figure 32: Mid-plane colored by cell size. d inj = 180µm. (a) ∆x min = 5µm; (b) ∆x min = 2.5µm.

Figure 33 :

 33 Figure 33: Estimation of the spray angle α at different distances z/d inj from the nozzle exit. d inj /∆x min = 72.

6 , 5 d 6 , 5 d

 6565 We g = 830 Experimental correlation d inj /∆x min = 22.inj /∆x min = We g = 1470 Experimental correlation d inj /∆x min = 22.inj /∆x min = 45

Figure 35 :

 35 Figure 35: Numerical mean jet trajectories for (a) We g = 830 and (b) We g = 1470 compared with the experimental correlation from Becker [32]. The grey area represents the standard deviation.

Table 2 :

 2 Error norms for each mesh resolution at the end of the notched disk test, and computational cost of the runs.

	R/∆x	E bound (T )	E shape (T )	E mass (T )	mass loss [%] N cores RCT [µs]
	12.5 25 50	4.562 × 10 -10 1.360 × 10 -2 1.410 × 10 -7 2.406 × 10 -4 0 7.544 × 10 -3 1.158 × 10 -7 1.986 × 10 -4 1.119 × 10	1 1	3.3 2.9

-11

2.226 × 10 -3 1.736 × 10 -9 2.980 × 10 -6 2 3.5 100 0 8.427 × 10 -4 1.874 × 10 -8 3.219 × 10 -5 8 4.9 200 1.274 × 10 -11 3.434 × 10 -4 3.473 × 10 -9 5.965 × 10 -6 28 5.8 400 0 1.224 × 10 -4 3.076 × 10 -9 5.284 × 10 -6 112 5.8

Table 3 :

 3 Error norms at the end of the circle deformation test, and computational cost of the runs.

	R/∆x	E bound (T )	E shape (T )	E mass (T )	mass loss [%] N cores RCT [µs]
	12.5 25 50 100	2.128 × 10 -7 3.642 × 10 -2 1.790 × 10 -9 2.518 × 10 -6 1.212 × 10 -7 8.452 × 10 -3 4.637 × 10 -8 6.551 × 10 -5 5.266 × 10 -8 1.970 × 10 -3 3.223 × 10 -8 4.558 × 10 -5 2.012 × 10	1 1 2	4.7 3.8 3.7

-8

6.016 × 10 -4 1.053 × 10 -8 1.490 × 10 -5 8 5.2 200 1.917 × 10 -9 2.719 × 10 -4 3.155 × 10 -9 4.464 × 10 -6 28 6.1 400 8.623 × 10 -10 1.289 × 10 -4 1.165 × 10 -9 1.649 × 10 -6 112 6.4

Table 4 :

 4 Computational costs of the sphere deformation simulations R/∆x min N cores N cells,ini /N cores RCT [µs] RCT ACLS [µs] RCT AMR [µs]

	15	10	93,361	153.6	77.2	73.6
	30	28	122,645	210.6	86.3	120.3
	60	224	71,239	372.4	99.4	267.2
	4.3 Capillary-driven flows
	4.3.1 Damping of a capillary wave on static triangular grids

Table 5 :

 5 Physical and numerical parameters for the viscous damping of a capillary wave case.ρ l /ρ g µ l /µ g La = ρλσ/µ 2

	λ/∆x	∆t *

Table 6 :

 6 Error levels of the L 2 norm (Eq. (74)) and reduced computational time (RCT; Eq. (67)) for each resolution.

	λ/∆x	L 2	N cores N cells /N cores RCT [µs]
	30 100 200	1.528 × 10 -4 3.250 × 10 -5 1.647 × 10 -5	1 1 4	2028 22,476 22,520	23.2 13.3 20.2

Table 8 :

 8 Cases definition for the Rayleigh-Plateau instability simulations.

	Case λ [mm] λ/r 0 ξ = kr 0 N cells
	1	16	48	0.13	2.5M
	2	8	24	0.26	1.3M
	3	4	12	0.52	637K
	4	3.3	10	0.63	525K
	5	2.8	8.5	0.74	453K
	6	2.3	7	0.90	372K

performed at each timestep. Classical water/air parameters are used, and are summarized in Table

Table 10 :

 10 Numerically predicted mean growth rates versus theoretical growth rates from Weber's theory, with r 0 /∆x = 6.

	Case 1	ξ 0.13	ω/ω 0 | mean ω/ω 0 | th Relative error Standard deviation SD 0.109 0.130 -16.2% 0.011
	2	0.26	0.273	0.252	+8.3%	0.059
	3	0.52	0.452	0.442	+2.3%	0.020
	4	0.63	0.469	0.483	-2.9%	0.019
	5	0.74	0.470	0.489	-3.9%	0.034
	6	0.90	0.355	0.383	-7.3%	0.048
	merging. In order to demonstrate that the presented procedure allows implicit topology changes, a head-on
	collision of two equal-size droplets of water in air is computed.

Table 12 :

 12 Computational performances of the water/air head-on droplet collision simulation. D/∆x min ∆x max /∆x min N cells,ini N cores RCT AMR [µs] RCT [µs]

	133	17	36M	448	23.8	239.4
	5 Application to injection systems
	5.1 LES of a low-pressure water jet from a compound nozzle

Table 13 :

 13 Computational costs and performances of the low-pressure compound injector simulations on adaptive grids. The simulations have been conducted on 2nd Gen AMD EPYC processors. d inj /∆x min N cores N cells RCT [µs] RCT ACLS [µs] RCT AMR [µs] RCT Poisson [µs] CPU time [h]/Phys. time [ms]

	36	600	48M	342.3	99.7	137.2	50.7	52,174
	72	1400 141M	610.0	142.0	208.6	199.2	915,531

Table 15 :

 15 Error levels of the L 2 norm (Eq. (84)) for jet in crossflow simulations We g d inj /∆x min L 2

	830	22.5	1.079
	830	45	0.661
	1470	22.5	1.376
	1470	45	0.442

Table 16 :

 16 Number of computational cores, mesh cells and droplets generated for jet in crossflow simulations. All values correspond to time instant t = 0.3ms.We g d inj /∆x min N cores N cells N drops

	830	22.5	560	78M	769
	830	45	5120 734M	3961
	1470	22.5	800	117M	1174
	1470	45	5960 876M	6615

Table 17 :

 17 Computational performances for jet in crossflow simulations at time instant t = 0.3ms. We g d inj /∆x min RCT [µs] RCT ACLS [µs] RCT AMR [µs] RCT Poisson [µs] CPU time [h]/Phys. time [ms]

	830	22.5	836.1	88.9	586.8	97.4	76,800
	830	45	1612.0	116.8	1151.2	167.2	922,100
	1470	22.5	833.4	100.4	503.9	152.6	105,600
	1470	45	2158.3	113.2	1660.8	190.8	1,200,000

Acknowledgements

Financial support for the conducted research presented in this article has been provided by the French National Association for Research and Technology ANRT and Safran SA under the CIFRE convention No. 2018/0002, and by the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 765998 in the project ANNULIGhT. Computer resources have been provided by CRIANN, France, under the allocation 2012006, and by GENCI, France, under the allocations A0072B06880 and A0092B11072. The authors would also like to acknowledge their academic partners at LEGI laboratory, Grenoble, France, for the fruitful discussions on this topic.

LES of a high-pressure kerosene jet in crossflow

In order to test the capabilities of the algorithm for resolving atomization in aeronautical injection systems, a non-reactive liquid kerosene jet in crossflow (JICF) injection has been simulated. Liquid is injected through a nozzle into a plenum where a stream of air flows in the perpendicular direction to the liquid (the crossflow). Then, the jet bends towards the air direction and atomizes due to the aerodynamic interaction. The chosen configuration is the experimental test-bench by Becker and Hassa [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF], which consists of a kerosene JICF injected into a high-pressure environment representative of the ambient conditions found within aeronautical gas turbines. The physics of the JICF is governed by the two following parameters: the kinetic energy ratio between the liquid and gaseous phases q and the gaseous Weber number We g . They are defined as follows:

where ρ is the density, u is the bulk velocity of the flow, d inj is the liquid nozzle diameter (0.45mm in Fig. 34) and σ is the surface tension coefficient. The subindexes l and g denote liquid and gaseous phases respectively.

Numerical setup

Fig. 34 (a) shows the numerical setup replicating the experiment from [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF]. It consists of a plenum of dimensions (L x , L y , L z ) = (270mm, 40mm, 25mm) and a tapered nozzle for liquid injection. A close-up view of the nozzle is shown in Fig. 34 (b). The injected liquid is JET-A kerosene with density ρ l = 795kg.m -3 , dynamic viscosity µ l = 1.5 × 10 -3 kg.m -1 .s -1 and surface tension σ = 22 × 10 -3 kg.s -2 . The gaseous crossflow is pressurized air at 6bar with density ρ g = 7.21kg.m -3 and dynamic viscosity µ g = 1.82 × 10 -5 kg.m -1 .s -1 . The density and dynamic viscosity ratios are respectively ρ l /ρ g = 110 and µ l /µ g = 83. The Ohnesorge number is Oh = µ l / ρ l σd inj = 0.017. Two operating points studied experimentally by [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF] are simulated. They are given by the dimensionless numbers q and We g shown in Table 14. Two interface resolutions are considered for both operating points: ∆x min = 20µm and 10µm, that can be expressed relative to the injection diameter as d inj /∆x min = 22.5 and 45 respectively. Therefore, a total of four JICF cases are computed. These simulations are performed with LES where a dynamic Smagorinsky closure is used for modeling the unresolved turbulence scales [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]. For kerosene injection, a Poiseuille profile is prescribed at the liquid inlet. The mean velocity of this profile is equal to the bulk velocity u l given in Table 14. Regarding the gaseous phase, the experiments from [START_REF] Becker | Breakup and atomization of a kerosene jet in crossflow at elevated pressure[END_REF] report a gaseous boundary layer thickness of 5mm right upstream the injector. As the numerical domain is smaller than the experimental test bench, a velocity profile with near-wall and outer regions is specified at the gaseous inlet. The near-wall region is modeled considering that the gaseous phase evolves in the streamwise direction as