Pierre Milbeo

Matthieu Simon

Claude Didierjean

Emmanuel Wenger

Emmanuel Aubert

Jean Martinez

Muriel Amblard

Monique Calmès

Baptiste Legrand

A bicyclic unit reversal to stabilize the 12/14-helix in mixed homochiral oligoureas

published or not. The documents may come

A bicyclic unit reversal to stabilize the 12/14-helix in mixed homochiral oligoureas Pierre Milbeo a , Matthieu Simon a , Claude Didierjean b , Emmanuel Wenger b , Emmanuel Aubert b , Jean Martinez a , Muriel Amblard a , Monique Calmès a* , Baptiste Legrand a* a IBMM, Univ Montpellier, CNRS, ENSCM, Montpellier, France b Université de Lorraine, CNRS, CRM2, Nancy, France

The insertion of cyclic building blocks in oligoureas to stabilize or modulate the properties of the 12/14-helix was often fruitless. We herein propose a fully compatible highly constrained building block that could be incorporated into oligoureas to develop highly stable and functionnal oligoureas helices.

Peptidomimetic foldamers, differing from α-peptides by the addition of at least one atom in the peptide backbone, can mimic the structural and functional properties of biopolymers. They are interesting candidates for biological applications, as catalysts and in the field of materials. 1 The development of a wide range of -, -or δ-amino acids and the design of oligomers with heterogeneous backbones result in a high diversity of helical structures. [START_REF] Horne | As selected reviews[END_REF] Besides, the use of peptidomimetic -building blocks owning an additional heteroatom in the backbone led to the construction of oligomers such as oligoureas. 3 They exhibit particular folding properties mainly related to the ability of the urea link to form bifurcated hydrogen bonds. 2c,3e-I The group of Guichard has first shown that they adopt a 12/14-helix stabilized by hydrogen bonding network forming 12-and 14-membered pseudocycles. 3e-i More recently, they have demonstrated that helical structures are also obtained from hybrid peptides combining urea/amide or urea/carbamate bonds. 4 Furthermore, robust helical structures can be achieved by conformational restriction of the building blocks using cyclic units as mainly described in -and -peptide foldamers. 2b,5 Nevertheless, incorporation of appropriate constrained cyclic residues in oligoureas is challenging and only few examples are reported in the literature. Our group has first demonstrated that short homo-oligourea sequences based on Bicyclic Amino Carbamoyl units (BAC units, derived from the (S)-2aminobicyclo[2.2.2]octane carboxylic acid, [(S)-ABOC]) 6 (Fig. 1), lead to a highly stable 12/14-helix both in solution and in the solid state. 7 No trans/cis isomerization urea bond was detected in solution for the (S)-BAC homo-oligoureas whereas it occurs for oligoureas derived from acyclic β 3 -amino acids. 3h Then, it has been reported that replacement of the central residue of an acyclic oligourea pentamer by a chiral cis-1,2diaminocyclohexane residue slightly enhances helix stabilization, while the use of the corresponding trans residue compromises the helical folding. 8,[START_REF] Wechsel | The use of a chiral cis-1,2-diaminocyclohexane unit has been also reported as a screw-sense controller of helical conformation of achiral acyclic homo-oligoureas (hAib derivatives)[END_REF] Later, we reported on the synthesis of mixed oligoureas alternating (S)-and (R)-BAC units with acyclic β 3 -amino acid derivatives, i.e. (S)-N-(2aminopropyl)carbamoyl (APC) and (S)-N-(2-amino-3phenylpropyl)carbamoyl (APPC) units, in order to elaborate robust and versatile oligoureas (Fig. 1a,b). 10 Unfortunately, extensive structural analyses revealed that the homochiral S/S sequences did not fold, while the corresponding heterochiral R/S sequences led to a distorted left-handed 12/14-helix with solvent molecules intercalated in its bifurcated hydrogen-bond network.

In this piece of work, we have modified the substitution pattern of the bicycle by reversing its orientation within oligoureas using the chiral (S) -1,2-diaminobicyclo [2.2.2] octane motif [(S)-DABO] recently developed in our group (Fig. 1a). 11 Therefore, the synthesis of the mixed oligoureas can be reversed, allowing the permutation of the bridgehead position and the 2-position in the final molecule, compared to previously synthesized oligoureas based on the ABOC residue (Fig. 1b). In this orientation, we hypothesized that the reverse Bicyclic Amino Carbamoyl unit [(S)-rBAC] derived from the (S)-DABO could share similar torsion angle values and signs in the mixed (S)-rBAC/(S)-AAC oligoureas inducing canonical mixed oligourea helices. Thus, we report here the synthesis and the extensive structural characterization of oligoureas alternating the (S)-rBAC and (S)-AAC residues. We tentatively rationalized the folding behaviours of this class of oligoureas incorporating the (S)-BAC, (R)-BAC and (S)-rBAC bicyclic units using density functional theory (DFT) calculations.

We used stepwise solution-phase method to prepare four oligoureas 1-4 (Fig. 1c, Scheme S1) of different lengths alternating (S)-rBAC and (S)-AAC residues. The first coupling was performed by adding an equimolar amount of H-(S)-DABO-Boc (Fig. 1a) to the succinimidyl carbamate of H--Ala-OBn (SuO-CO--Ala-OBn) in the presence of N,N-diisopropylethylamine (DIEA). Chain elongation was achieved by successive Boc removal, trifluoroacetic acid (TFA) salt neutralization, and urea bond formation using alternatively acyclic and cyclic residues (see Scheme S1 in the SI for full details). Succinimidyl carbamate derivatives Boc-APC-OSu and Boc-APPC-OSu (Fig. 1d), 12 were used to introduce acyclic residues. The (S)-rBAC unit was correctly incorporated using freshly activated H-(S)-DABO-Boc with the DSC reagent under optimized conditions that limited formation of the by-product (i.e. cyclic biuret, Fig. 1d) by limiting the activation time, and by repeating the coupling step until completion of the reaction (detailed procedure available in SI).

We first assessed the conformational preferences of the newly synthesized oligoureas 1-4 recording their CD spectra in CH3OH at 20°C. All compounds shared the characteristic signature of the oligoureas right-handed 12/14-helix with a strong positive maximum around 205 nm (Fig. 2a). These results indicated that homochiral mixed oligoureas alternating (S)-rBAC and (S)-AAC motifs might adopt right-handed 12/14-helices in solution. Such hypothesis was stressed using NMR spectroscopy in CD3OH and DMSO-d6. We first assigned nearly all resonances of the compounds 2-4 combining COSY, TOCSY and ROESY spectra (Tables S1-7 and Fig. S4-11). The 3 J(NH, α CH) coupling constants of the (S)-rBAC and (S)-AAC residues in the oligomers were characteristic of an antiperiplanar arrangement with average values of 9.5 ± 0.7 Hz and 9.5 ± 0.6 Hz in CD3OH and DMSO-d6, respectively (Table S8). Moreover, we measured large chemical shift differences between the diastereotopic α CH2 protons of the central acyclic residues APC2 and APPC4, i.e. Δδ ( α CH2) > 1.0 ppm (Table S9) and found a dozen of typical i,i+2 and i,i+3 medium and long range NOE cross-peaks, which are fully compatible with the right-handed 12/14-helical fold (Fig. 2b) in both solvents. Such evidences indicated that (S)-rBAC/(S)-AAC oligoureas adopted a stable 12/14-helix in CD3OH in agreement with the CD data, but also in the chaotropic DMSO-d6. When we used a complete set of NOEs as distance restraints for the NMR structure calculation of hexamer 3, the resulting 15 lowest-energy conformations exhibited a welldefined right-handed 12/14-helix (rmsd value of 0.44 Å on the backbone heavy atoms) stabilized C=O(i)⋅⋅⋅HN(i+1) and C=O(i)⋅⋅⋅HN′(i+2) hydrogen bonds (Fig. 2b). Importantly, highresolution crystals structures tetramer 2 and octamer 4 suitable for X-ray diffraction (XRD) were obtained by slow evaporation of methanol solutions. Both asymmetric units consisted of one oligourea, which adopted a canonical right-handed 12/14-helix with 2.5 residues per turn and a pitch around 5.1 Å (Fig. 2c 6°, 1 = 57 ± 5° and 2 = 72 ± 8° when the last residue torsion angle values were omitted. Such values were very close to those reported for the canonical oligourea 12/14-helix composed of acyclic  3 -amino acid derivatives 3h . The backbone structures of the canonical (S)-AAC octamer (ccdc75017) and the mixed (S)-rBAC/(S)-AAC oligourea 4 were well superimposed as shown in Fig. 3a (the RMSD value of the backbone heavy atoms was 0.86 Å). The main structural deviations between compounds 2 and 4 occurred at both terminal parts, which are key areas for the crystal packing. The extremities linked the oligomers in a head to tail fashion as previously observed in oligourea crystal structures 3h, 7 (Fig. S15).

We then compared the relative stability of the (S)-rBAC/(S)-AAC oligourea helix to the (S)-AAC and (S)-BAC homooligoureas, and heterochiral mixed (R)-BAC/(S)-AAC helices. As previously described, 3f,7,10 the signal intensity at 205 nm (θ205) increased with the oligomer's length and only a slight and reversible decrease was observed when the temperature rose to 55°C (Fig. S3). The weak thermal loss of structure was fully reversible and the oligomer's stability increased with their lengths. Few NOEs could be detected on the ROESY spectrum of 1 indicating that it should be disordered or highly flexible. In agreement, θ205 was about twice less weak for 1 than for the (S)-BAC dimer (θ205 = 10×10 3 and 17.8×10 3 deg.cm².mol -1 , respectively) that was able to initiate a 12/14 turn with only two units 7 . Nevertheless, it is worth noting that compound 3 and the (S)-BAC hexamer shared comparable maxima values, θ205~35×10 3 deg.cm².mol -1 (while a maximum at 27×10 3 deg.cm².mol -1 was reported for an acyclic oligourea hexamer) 10 . Thus, the (S)-rBAC/(S)-AAC 1:1 pattern shared a comparable stability with the (S)-BAC helix exclusively composed of bicyclic units. The bicycle bridgehead steric hindrance prevented the trans/cis isomerization of the adjacent urea bond. We detected only N'H(i)-Hβ(i+1) NOEs assigned to trans/cis isomerization of the urea bonds preceding the AAC residues.

Finally, we tentatively rationalized why the 12/14-helix could fully accommodate the r(BAC) residue in homochiral oligoureas thanks to DFT calculations and a close analysis of the optimized structures. First, it is worth noting that the BAC residue imposed a steric unfavourable (-)-synclinal conformation to the (S)-acyclic residues with a θ1 backbone torsion angle average value of -56°, while a (+)-synclinal arrangement around the main chain C-C bonds was typically observed for the β 3 , γ 4 -peptides and the acyclic oligourea canonical helix. 3e-I Homochiral 12/14-helices models with similar parameters were built, i.e. (S)-BAC/(S)-AAC and (S)-rBAC/(S)-AAC, and optimized at the B3LYP-D3/6-31G (d,p) level of theory (detailed protocol available in SI). Calculations converged toward stable helices in vacuum and the structural deviations were very slight between the initial and the optimized helices (the torsion angles did not vary significantly, Fig. S16-17). Interestingly, the theoretical study was in good agreement with the experimental data. The (S)-rBAC/(S)-AAC was the most stable helix (Table S11). The disordered (S)-BAC/(S)-AAC oligourea had significantly higher energy when it was optimized as a 12/14-helix single conformer (ΔG = +31.05 kJ.mol -1 ). A close inspection of the models showed that various steric effects involved the AAC side chains (Fig. 3). In the (S)-rBAC/(S)-AAC helix, these latter were in optimal equatorial position while they were virtually projected in axial within the heterochiral oligourea, generating strong steric clashes with the adjacent backbone carbonyls. In both (S) and (R)-BAC/(S)-AAC, the side chains were projected in inadequate intermediate positions, which prevented the folding of the homochiral (S)-BAC/(S)-AAC oligourea, but could be accommodated in the (R)-BAC/(S)-AAC helix. 10 In this latter, the local crowding was offset by the insertion of solvent molecules in its bifurcated hydrogen bonds, which significantly bent its axis while the helical structure of compounds 2 and 4 was ideally straight (Fig. 3).

In summary, a high content of rBAC (1:1 pattern) could be incorporated without altering the 12/14-helix structure while increasing its stability compared to oligourea helices composed solely of highly constrained bicyclic residues. Furthermore, the rBAC incorporation locked the urea bond beside the bicycle bridgehead in trans configuration. However, the minimal content of bicyclic units to be incorporated in acyclic oligoureas for a maximum stabilization of the 12/14-helix remains to be determined. Interestingly, such incorporation in functional oligoureas would be valuable to tune their properties.

The authors thank the "Laboratoire de Mesure Physique" (LMP) of the university of Montpellier, the SCBIM and "Université de Lorraine" for NMR spectroscopic and XRD facilities. The EXPLOR mesocentre is also thanked for providing access to computing facility (project 2019CPMXX0984/wbg13). There are no conflicts to declare.

Fig 1 .

 1 Fig 1. (a) Monomeric units and torsion angles nomenclature. (b) General oligourea sequences based on (S)-BAC (in blue), (R)-BAC (in red) and (S)-AAC residues, (S)-rBAC (in green) and (S)-AAC motifs. (c) Sequences of (S)-rBAC-containing oligoureas 1-4. (d) Chemical structures of succinimidyl carbamate derivatives of β 3 -amino acids (Boc-APC-OSu and Boc-APPC-OSu), and cyclic byproducts derivatives (cyclic urea, i.e. allophanate) and cyclic biuret).

Fig. 2

 2 Fig. 2 (a) CD spectra of 1 (plain line), 2 (dashed line), 3 (dash-dotted line) and 4 (dotted line) recorded at 20°C. The absolute values of the CD positive maxima (in deg.cm 2 .dmol -1 , λ = 205 nm) of the homo-(S)-BAC and (S)-rBAC/(S)-AAC oligoureas are indicated in the right corner. (b) Characteristic non-sequential backbone/backbone and backbone/side chain NOE cross-peaks detected both in CD3OH and DMSO-d6. Superimposition of the 15 lowest energy NMR solution structures of hexamer 3. (c) Crystal structures of 2 and 4. The (S)-AAC residues are in yellow. The hydrogen bonds are shown in yellow dashed lines.

Fig 3 .

 3 Fig 3. (a) Superimposition of the crystal structures of the (S)-AAC (ccdc75017, in yellow) and mixed (S)-rBAC/(S)-AAC oligourea octamers (4, in purple) (left panel). Crystal structure of the mixed heterochiral (R)-BAC/(S)-AAC oligourea octamer (right panel). (b) Schematic representation of the homochiral 12/14-helix models optimized by DFT. Steric clashes due to non-ideal configuration and atoms promiscuity was indicated as red arrows.


César I Pairetti

Santiago Márquez Damián

Norberto M Nigro

Stéphane Popinet

Stéphane Zaleski

Mesh resolution effects on DNS-VOF of primary atomization

Keywords: Atomization, Volume of Fluid (VOF), Adaptive Mesh Refinement (AMR), Basilisk

In this work, we use Direct Numerical Simulation (DNS) through a Volume of Fluid (VOF) solver with Adaptive Mesh Refinement (AMR) to analyze the atomization of a pulsating round liquid jet with a narrow length scale range configuration. We propose three grid sizes based on characteristic scales we estimate from the deformation and fragmentation processes of this problem. We compute drop statistics and general spray features for each case. We found that mesh resolution affects the atomization rate and the probability density function of droplets size and position, not only for the under-resolved drops but all liquid structures. The two simulations with higher grid refinement presented volume-weighted distribution with minor differences. Therefore, we propose assessing the accuracy of atomization simulations based on the volume fraction of under-resolved structures.

Introduction

Primary atomization is a complex phenomenon of great interest because, as described by [START_REF] Lefebvre | Atomization and sprays[END_REF], it is present in several natural processes (e. g. sea-waves, waterfalls) and technical applications (e. g. internal combustion engines, spray coating, spray drying, disease propagation).

Atomization is a multiscale phenomenon whose characteristic length scales depend mainly on the turbulent intensity and the acting interface forces. As [START_REF] Gorokhovski | Modeling primary atomization[END_REF] discussed in their paper, this is one of the features that make it difficult to apply Direct Numerical Simulations (DNS) to study primary atomization problems. In particular, high-speed injection of round liquid jets presents atomization at high Reynolds and Weber numbers which produce a wide range of flow length-scales. This feature increases the computational cost of DNS, as discussed by [START_REF] Chesnel | Subgrid analysis of liquid jet atomization[END_REF], and also makes it more difficult to study the problem experimentally, as noted by [START_REF] Dumouchel | On the experimental investigation on primary atomization of liquid streams[END_REF].

Nevertheless, the growth of computational power and the optimization of numerical methods during the first two decades of the 21st century have driven some research on this topic. Some of the goals from these works are improving the understanding of fragmentation mechanisms, especially in the dense spray region, and developing low-cost models with direct technological applications, such as injector nozzle design. In the following paragraph, we summarize some research projects that illustrate the current state of the art on primary atomization DNS.

On one of the first fuel injection DNS analysis, [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] simulated a diesel jet at moderate speeds, applying advanced interface representation techniques (CLSVOF) and a turbulent boundary condition based on the works of [START_REF] Klein | A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[END_REF].

From these results, [START_REF] Chesnel | Subgrid analysis of liquid jet atomization[END_REF] estimated the effect of momentum transport phenomena at the sub-grid scales and proposed algebraic interface models in the context of an ELSA (Eulerian-Lagrangian Spray Atomization), developed and validated by [START_REF] Lebas | Numerical simulation of primary break-up and atomization: Dns and modelling study[END_REF] and [START_REF] Anez | Liquid atomization modeling in openfoam (r)[END_REF] among others.

A few years later, [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF] studied different atomization regimes with moderate Reynolds and Weber numbers using a second-order level-set method. In their analysis, they describe in detail the jet tip deformation and fragmentation process and how they interact with the vortex dynamics. Although the boundary conditions used in these analyses do not accurately represent the internal flow of a fuel injector, these simulations capture in great detail the growth of hydrodynamic instabilities on the surface of the jet. Recently, [START_REF] Shinjo | Fluid dynamic and autoignition characteristics of early fuel sprays using hybrid atomization les[END_REF] proposed atomization models using DNS results to develop a Sub-Grid Scale (SGS) breakup model, based on ligament pinch-off, to perform Large Eddy Simulations (LES) of diesel combustion.

In the particular case of fuel injection, validation through experimental data is still a big challenge. Even if there are benchmark injection experiments, such as the spray A developed by the [START_REF]Engine Combustion Network, Engine Combustion Network[END_REF], they measure only the disperse spray zone during periods of a few miliseconds. In contrast, [START_REF] Salvador | Determination of critical operating and geometrical parameters in diesel injectors through one dimensional modelling, design of experiments and an analysis of variance[END_REF] studied the spray A case in a low-velocity regime using the Volume of Fluid (VOF) method on uniform cartesian meshes. They reproduced the first 20 µs of injection, reaching the breakup length and the early development stages of the dense spray.

Regarding more fundamental studies, [START_REF] Jarrahbashi | Vorticity dynamics for transient highpressure liquid injection[END_REF] and [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF] analyzed the deformation of cylindrical jets for a broad spectrum of flow regimes. They describe the vortex dynamics and how it affects the instabilities developed along with the jet core, analyzing the interaction between turbulence and fragmentation mechanisms. Using similar methods [START_REF] Zandian | Planar liquid jet: Early deformation and atomization cascades[END_REF] described the breakup process on planar jets and how the fragmentation modes develop for a different set of parameters. Even if these analyses characterize the growth of interface instabilities in detail, they do not give further insights on how this process impacts on the final drop size distribution.

Besides these difficulties, DNS is still used to study early atomization stages on fuel injection problems. [START_REF] Ling | Direct numerical simulation of an atomizing biodiesel jet: Impact of fuel properties on atomization characteristics[END_REF] and [START_REF] Zhang | Effect of fuel viscosity on the atomization of diesel and biodiesel jets[END_REF][START_REF] Zhang | Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under nonevaporative operating conditions[END_REF] applied the VOF technique with Adaptive Mesh Refinement (AMR) to analyze the primary atomization of diesel, biodiesel, and gasoline on injection problems at moderate speeds. Through a complete statistical characterization of the cases, these works describe the overall evolution of the fuel spray. These simulations employed the Basilisk C library, by [START_REF] Popinet | The basilisk code[END_REF]; we obtained the results presented in this manuscript using an updated version of these numerical tools.

In this work, we employ the Basilisk VOF-AMR solver to analyze the primary atomization of a pulsating liquid jet, using different grid sizes based on the characteristic length scales of the flow. Section 2 describes the mathematical model and numerical techniques implemented on the Basilisk solver. Section 3 defines the atomization problem designed for this study and the criteria employed on the grid convergence analysis. Section 4 describes the general features of the injection process and presents the drops statistics, comparing the results for the three refinement levels. Finally, section 5 summarizes the conclusions and gives some insight regarding future works.

The main goal of the present analysis is to find quantitative descriptions for low grid resolution effects on primary atomization DNS. The following questions serve as a guide for that objective:

• Is it possible to determine a priori the grid resolution needed to reproduce the physical behavior of primary atomization?

• Is it possible to assess the numerical accuracy of a simulation without running a mesh convergence analysis?

• What should be the threshold to determine which liquid structures are resolved or not? How to measure the impact of the under-resolved drops?

• How does poor mesh resolution affect the global properties of a spray, such as the atomized volume or the interface area density?

Mathematical model and numerical method

The Momentum-Conserving VOF (MCVOF) solver described in this section is based on the compressible solver by [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF], considering incompressible flow. These schemes are described in detail by [START_REF] Pairetti | Bag mode breakup simulations of a single liquid droplet[END_REF] and, more recently, by [START_REF] Zhang | Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under nonevaporative operating conditions[END_REF]. The three articles aforementioned present validation cases to test the consistency and robustness of the momentum-conserving advection schemes and the discrete balance of surface tension forces. In particular, [START_REF] Zhang | Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under nonevaporative operating conditions[END_REF] also report numerical results for the spray G benchmark case. The general layout of the model and numerical methods are given next.

Incompressible Multiphase Flow model

The mass and momentum conservation equations for incompressible and isothermal flow can be expressed as:

∇ • u = 0 (1) ∂ρ u ∂t + ∇ • (ρ u u) = -∇p + ∇ • (2µD) + f σ (2)
where u( x, t) is the velocity field and p( x, t) is the pressure field. Tensor D is defined as

1 2 ∇ u + (∇ u) T .
f σ = σκ n s δ s (3)
where the surface tension coefficient σ is considered constant. The force only acts at the free surface, hence the Dirac function δ s , and also depends on the interface shape, particularly on its curvature κ and normal n s .

In the context of the one-fluid formulation for multiphase flows, the color function c( x, t) gives the phase spatial distribution on the domain, taking value c( x p , t) = 1 if the point p is filled with one phase and c = 0 otherwise.

The interface is then located in the discontinuity surface of the c function, defining δ s , n s , and κ. The color function transport is given by the following Equation:

∂c ∂t + ∇ • (c u) = c ∇ • ( u). (4)
where the right-hand side is equal to zero due to Equation (1).

Numerical method

The numerical implementation used in this work represents the freesurface using the Piecewise Linear Interface Capturing (PLIC) VOF method by [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] and [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF]. In the Finite Volume Method (FVM) context, the mean value of the c function on a cell is:

f Ω = 1 ∆Ω Ω c(x, t) dV, ( 5
)
where ∆Ω is the volume of the cell Ω. Then f Ω is the volume fraction of the main phase (c = 1) in the cell. The mixture properties on the cell can then be computed by arithmetic means:

ρ Ω = f ρ l + (1 -f )ρ g µ Ω = f µ l + (1 -f )µ g (6)
From this point on, we drop the cell-based subscript for the sake of clarity.

The following equations are valid on each cell and must be solved for all of them. In this context, the approximate projection method by [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] is used to solve coupling between Equations ( 1) and ( 2), considering that the velocity ( u) is staggered in time VOF (f ) and pressure (p), which is noted next by the superscript indicating the time-step for each variable. The discrete equations can then be expressed as in [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]:

f n+ 1 2 -f n-1 2 ∆t + ∇ • ( u n f n ) = c c ∇ • u n (7) ρ u * -ρ u n ∆t + ∇ • ρ n+ 1 2 u n u n = ∇ • [µ n+ 1 2 (D n + D * )] + (σκδ s n s ) n+ 1 2 (8) ∇ • ∆t ρ n+ 1 2 ∇p n+ 1 2 = ∇ • u * (9)
u n+1 = u * - ∆t ρ n+ 1 2 ∇p n+ 1 2 (10)
Here c c is the contraction function used in the split volume fraction advection scheme from [START_REF] Weymouth | Conservative volume-of-fluid method for free-surface simulations on cartesian-grids[END_REF], c c = 1 if f > 0.5 and c c = 0 otherwise. The advection terms in Equations ( 7) and ( 8) involve VOF and momentum fluxes respectively. In this context:

∇ • ( u n f n ) = f,i F f,i ∆Ω F f,i = f a ( u f,i S f ), (11)
∇• ρ n+ 1 2 u n u n = f,i F ρ u f,i ∆Ω F ρ u f,i = [(ρ l u) f f a + (ρ g u) f (1 -f a )] ( u f,i S f ), (12)
where F f,i and F ρ u f,i are the VOF and momentum fluxes through the face i, computed from the face fraction f a = V a /(u f,i S f ) represented in Figure 1: the left side shows a cell where the light gray region is the main phase volume and the dark gray region is the main phase volume transported to the neighbor cell. The face-reconstruction of the velocity, u f , and momentum for each phase, (ρ l u) f and (ρ g u) f , are computed by the second-order upwind BCG scheme, by [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF], using a minmod slope limiter. The right side of Figure 1 shows that the VOF function f implicitly defines the momentum fields that must be advected following Equation 12 to ensure consistency between mass and momentum transport.

We use the semi-implicit Crank-Nicholson scheme to compute the diffusive flux due to the viscous term in Equation [START_REF] Anez | Liquid atomization modeling in openfoam (r)[END_REF]. We discretize the surface tension force term at the face with the same scheme employed to compute the pressure gradient. This ensures a well-balanced formulation that reduces spurious currents, as explained by [START_REF] Popinet | Numerical models of surface tension[END_REF]. The interface curvature is computed using second-order stencils based on height functions computed by an analytical formulation from [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. Equations ( 9) and ( 10) are the projection steps that will ensure mass conservation for the velocity field at the step (n + 1).

The solver procedure is as follows: MCVOF method Given the fields at a time-step n: f n-1/2 , u n : 1. Compute face fluxes from equations ( 11) and ( 12)

2. Integrate VOF function by Equation ( 7) to get f n+1/2

3. Integrate the momentum Equation ( 8), adding the advection fluxes and the right hand-side terms, to get the velocity approximation u * 4. Solve Equation ( 9) to compute pressure p n+1/2

5. Correct u * by Equation [START_REF] Shinjo | Fluid dynamic and autoignition characteristics of early fuel sprays using hybrid atomization les[END_REF] to compute the next step velocity u n+1

In the simulations reported, we applied Adaptive Mesh Refinement (AMR) in the octree grid Basilisk implementation. In this context, the whole domain is the root cell, with refinement level 0 and side L 0 ; adding a refinement level to a cell consist of dividing it into eight children with a grid size equal to half the size of its parent. Given the maximum refinement level (L), the smallest grid size ∆ = L 0 /2 L . The refinement criterion is based on bounding a wavelet-based error estimation. Further details on this technique theoretical basis can be found in the works by [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF], [START_REF] Popinet | A quadtree-adaptive multigrid solver for the serre-greennaghdi equations[END_REF], and [START_REF] Van Hooft | Towards adaptive grids for atmospheric boundary-layer simulations[END_REF].

Simulation setup

We analyze the atomization of a circular jet injected on a gas-filled cubical chamber. This problem is solved using AMR with three grid refinement levels and the same refinement criteria. The physical properties of each phase are based on low-velocity diesel injection, following [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF], using a pulsating boundary condition for the inlet, as done by [START_REF] Ling | Direct numerical simulation of an atomizing biodiesel jet: Impact of fuel properties on atomization characteristics[END_REF][START_REF] Zhang | Effect of fuel viscosity on the atomization of diesel and biodiesel jets[END_REF]. This configuration has a narrow length-scale range and a simple deterministic injection velocity profile independent of grid refinement. The same adaptive refinement criteria are applied in all the simulations.

We bound the wavelets error estimates of the volume fraction f and the velocity field u, defining as maximum thresholds for the local error estimations E(f ) < 10 -5 and E( u) < 1 m/s respectively.

Problem definition and boundary conditions

The jet of diameter D inj is injected in a cube with a length of 30 D inj , as represented in Figure 2. The boundary condition on the injection face (x = 0) imposes the no-slip condition everywhere except on the liquid section

( y 2 + z 2 < D inj /2)
where the injection velocity is:

u x (t) = U inj [1 + ζ sin (f 2πt)] f = St U inj D inj . ( 13
)
This periodic perturbation has been employed in the works of [START_REF] Ling | Direct numerical simulation of an atomizing biodiesel jet: Impact of fuel properties on atomization characteristics[END_REF] to accelerate the atomization process near the nozzle. Moreover, this deterministic boundary condition produces simple hydrodynamic instabilities that are well resolved for all the grid size employed. It is worth to mention that [START_REF] Jarrahbashi | Vorticity dynamics for transient highpressure liquid injection[END_REF][START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF] reproduced early spray formation dynamics on simulations of axisymmetric jets with a single frequency perturbation, showing that this level of simplification still preserves the behavior of the problem of interest. All the simulations were performed using the same boundary condition with ζ = 0.1 and St = 1.3. On the remaining cube sides, we allow free outflow:

∂ n u Γ = 0 and p Γ = 0.
We use the phase properties of an academic diesel injection problem, following the work [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF], summarized in Table 1. This configuration presents a relatively small ratio between macro and micro length-scales, reducing the cost of resolving all the relevant scales for atomization phenomena.

D inj [µ m] U inj [m/s] ρ g [kg/m 3 ] ρ l [kg/m 3 ] µ g [kg/(m s)] µ l [kg/(m s)] σ[N/m] 100 100 25 696 10 -5 1.2 × 10 -3 0.06
In this context, the present mesh convergence study considers three maximum refinement levels, named M1, M2, and M3. The mesh step for each case will be defined by the end of the next section, based on the flow characteristic length scales.

Dimensionless groups and characteristic scales

The dimensionless groups relevant to a round liquid jet atomization problem are:

We l = ρ l U 2 inj D inj σ Re l = ρ l U inj D inj µ l ρ * = ρ l ρ g µ * = µ l µ g (14)
where ρ l and ρ g are the densities of each phase, in this case, liquid and gas respectively. The viscosity coefficients are noted by µ l and µ g . σ is the surface tension coefficient, U inj is the injection mean velocity and D inj is the nozzle diameter. The Reynolds and Weber numbers can be computed for both phases. They also can be used to express the Ohnesorge number

(Oh l = √ We l /Re l ).
In order to verify the accuracy of the incompressibility hypothesis, the Mach number (Ma g = U inj / √ kRT ) should be lower than 0.3.

The characteristic time based on the problem kinematics is

t c = (D inj /U inj ),
which in this case is exactly equal to 1 µs. We analyze several characteristic lengths, described in detail in Appendix .1, to define the grid resolutions for our simulations. We summarize the values of these length scales (l ) on Table 3. We briefly describe the relevance of each length-scale next.

The Kolmogorov scales in each phase (η g and η l ) define the requirements for resolving all the turbulent structures: according to [START_REF] Pope | Turbulent flows[END_REF], resolving all scales by pure DNS is only possible if ∆ < 2η. The Hinze scale (ζ) represents the biggest drop that does not suffer breakup by turbulent effects, as described by [START_REF] Ling | Multiscale simulation of atomization with small droplets represented by a lagrangian point-particle model[END_REF].

Given the effects of the shear stresses on the jet surface, Kelvin-Helmholtz Instability (KHI) dominates the deformation stage. Therefore, we consider the critical wavelength for KHI, λ KH c , based on linear stability theory considering a simplified shear flow on the jet surface. This is a conservative reference for the smallest unstable interface perturbations.

Finally, we consider [START_REF] Senecal | Modeling high-speed viscous liquid sheet atomization[END_REF] estimations of the relevant length scales for sheets and ligaments breakup: λ s and D lig . The short wavelength λ s is reported by these authors to be independent of the jet geometry, as it is also discussed by [START_REF] Reitz | Mechanism of atomization of a liquid jet[END_REF] who use a similar estimation for round jets. Regarding the diameter D lig , we use the same mass conservation relation to estimate the size of ligament detached from the sheet developed from the jet tip, measuring the film thickness locally on previous simulations.

It is important to recall that all these estimations are conservative, based on simplified geometries and flow conditions, aiming to consider the most challenging scenarios. Based on these length scales, we choose the mesh step for simulation M1 to be ∆ 1 = 1.46µm corresponding to a maximum octree-refinement level of 11. M2 and M3 will use level 12 and 13 respectively, resulting in ∆ 2 = 0.73µm and ∆ 3 = 0.37µm.

In this context, the M1 simulation will not be able to resolve shortwavelength unstable hydrodynamic instabilities. M2 will barely represent λ KH c with two cells and M3 will even capture stable interface perturbation.

Regarding fragmentation, the M1 grid resolution will not be able to represent ligaments predicted by [START_REF] Senecal | Modeling high-speed viscous liquid sheet atomization[END_REF] formulation, as of ∆ 1 > 0.5D lig . The M2 resolution will be able to reproduce ligaments with D > 0.54 D lig and M3 with D > 0.26 D lig .

We run the cases M1 and M2 until they reached 25 µs of physical time.

These simulations required approximately 3 500 and 50 000 CPU-hours respectively. At t = 25µs, mesh M1 had 87 237 340 cells and mesh M2 had 414 813 435 cells. We run the M3 simulation until t = 11.2µs, which took 218 000 CPU-hours and finished with 602 257 783 cells. We performed all the simulations on the machine Irene administrated by TGCC.

Results and discussion

In this section, we describe the physical problem and the fragmentation mechanism and we compare statistics for the three simulations to analyze the grid resolution effects. We first define the parameters to characterize the spray and individual droplets. The Sauter Mean Diameter (SMD) of a spray is defined as:

SM D = 6 N i V i N i A i (15)
where V d and A d are the volume and area of the "d" droplet respectively.

N is the total amount of drops in the spray. The SMD is equivalent to the diameter of a spherical drop with the same surface area density than the whole spray; smaller SMD imply higher area to volume ratio.

For single droplets, diameter and shape will be measured in terms of equivalent volume diameter (D 30 ) and sphericity (ψ) which can be defined as:

D 30 = 6 V d π (1/3) D 32 = 6 V d A d ψ = D 32 D 30 (16)
It should be noticed that D 32 is equivalent to the SMD of a single drop, and it is always smaller than D 30 ; except in spherical shapes, for which ψ = 1. The graphs in subfigures (e) and (f) show the evolution of the ratio be-

General overview

tween droplet volume V drops (t) = i N drops V i and injected volume V inj (t) = t Q(τ ) dτ .
The linear scale graph, subfigure (e), shows that at the beginning of injection (t < 5 µs) the total atomized volume is less than 1%. Between 5 µs and 10 µs the atomization rate increases significantly and at t = 20µs, around 10% of the injected volume is atomized. The graph with the logarithmic scale, subfigure (f), allows us to observe the evolution of the volume fraction in drops smaller than 3µm. After t = 10 µs, they represent less than 0.5% of the total injected volume for the three cases.

The last graph of Figure 3, subfigure (g), represents the spray SMD. Before 5 µs the amount of drops is so small that each new liquid structure affects the SMD significantly; it is also worth noticing that, as seen in subfigure (f), mostly small drops are produced. The slope of the SMD graph changes significantly near 8 µs. At the same instant, the atomization rate increases, as seen in subfigure (e). This change can only be possible due to the production of large drops, given that the small drops fraction remains almost constant, as can be noted in subfigure (f).

These three graphs show that small drops are produced from early stages of atomization, in contrast to large drops and ligaments that need more time to develop from shear-induced instabilities. Regarding grid resolution effects, we see that M1 presents a higher atomization rate, subfigure (e), and a lower volume fraction of small drops, subfigure (f). Remarkably, M1 and M2 present good agreement on the spray SMD after t = 10 µs, despite the aforementioned differences.

Figure 4 compares the evolution of the core's D 32 and spray SMD for the three simulations. M2 and M3 present good agreement (less than 2% difference) for both drops and core behavior, even at early atomization stages.

It is worth noticing that, despite this agreement, the liquid core of M1 has a much lower D 32 and, therefore, a higher surface density. From this result, we can claim that the M1 generates drops and surface area faster than the finer simulations. A better understanding of the relation between the jet evolution analyzed in the frames (a-d) in Figure 3 and the grid resolution effect on the spray characteristics requires a more detailed observation of the breakup mechanisms. These processes have been studied by [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF] and [START_REF] Ling | Direct numerical simulation of an atomizing biodiesel jet: Impact of fuel properties on atomization characteristics[END_REF], among others.

From this starting point, we highlight some aspects of the fragmentation process in Figure 5 to study how grid resolution can affect them. In this figure, the black circle marks the formation of holes in a region of a liquid sheet with velocities between 0 and 75m/s, denoting the presence of shear-induced instabilities. On the other hand, the detached filament marked by the red ellipse is practically at rest but its diameter is considerably reduced; the high aspect ratio and deformation resulting from aerodynamic effects will produce nonlinear capillary breakup, even if the shear stresses are not enough to induce fragmentation.

The previous analysis is also valid for the rest of the jet. The growth of instabilities in all the jet surface can be seen in Figure 6, where the liquid fraction on the middle plane is represented in black. This figure shows that perturbations in the rest of the jet surface produce thin structures as well.

We can also observe in this frame the loss of mass near the axis. Even if Helmholtz Instability (KHI) on the axis direction. The KHI growth rapidly generates thin liquid sheets that suffer fragmentation before reaching x = 3D inj , as we show in the red box detail; after this point, the instability has lost its axial symmetry due to the turbulent behavior of the mixing layer.

These dynamics follow the generation of transverse instabilities described by [START_REF] Jarrahbashi | Vorticity dynamics for transient highpressure liquid injection[END_REF] and [START_REF] Jarrahbashi | Early spray development at high gas density: hole, ligament and bridge formations[END_REF]. As the holes expand, they form ligaments and drops that detach from the core and impact on the upwind instabilities on the jet surface, increasing its fragmentation rate.

In this context, regarding grid resolution effects on the fragmentation process, it is worth to remark that both the sheet perforation and the ligament pinch-off occur when two different interfaces meet in the same cell; from that instant on, the grid cannot resolve the coherent structure. This fact could explain why coarse grids present faster fragmentation and, therefore, higher atomization rate.

Another spray feature that might be affected by poor mesh resolution is the spatial distribution of the drops, partially described in Figure 7.

The left column shows the complete jet for each simulation, colored by axial velocity (U x ) and the right column shows a two-way histogram where each pixel color represents the volume fraction of the drops aligned with the ray parallel to the axis and located at that particular (Y, Z) position.

The side views show several effects on the drop positions. First, in finer grids more drops are generated in the near nozzle region. Another effect of mesh refinement is the radial expansion of the spray: the M3 simulation presents the furthest drops. In the three cases, the drops detached from the jet tip are almost still, which means that they will not travel downward significantly. We can also notice that the M1 simulation has larger liquid structures near the jet tip; most of them present a high degree of deformation.

This fact will will be addressed quantitatively in the next section.

The scatter histograms show that, as we use higher refinement levels, mass distribution around the axis gets more even and covers a larger radius.

Increasing mesh resolution also increases the entrainment rate of drops in the mixing layer: the cases M2 and M3 present drops almost R inj /2.

Effects of mesh-resolution on drop statistics

In this section, we analyze the mesh resolution effects on drop size and shape. We first compare the droplets D 30 distributions for the times t = [START_REF] Klein | A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[END_REF][START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF][START_REF]Engine Combustion Network, Engine Combustion Network[END_REF] µs in Figure 8 and then focus on the spray at the latter time, fitting the distributions in Figure 9.

At all times, the most frequent diameter is 2 ∆; smaller diameters cannot be resolved, hence their presence is minimal. For example, M3 presents almost 10.000 drops in the range D 30 ∈ [0.5 : 1]µm while M2 produces less than 1000 drops in the same range. These graphs also show that big droplets (D 30 > 12 µm) take longer periods of time to develop, as there are almost no droplets in that range at 6µs, but all the simulations show more than 10 drops bigger than 10 µm at t = 11 µs.

The histograms can be compared quantitatively by, for example, fitting them using a log-normal distribution,

N (d) = C d exp - (ln d -μ) 2 2 σ2 . ( 17
)
Figure 9 shows the fitted histogram (weighted by D 30 ) using a logarithmic scale on both axes. This representation has been used by [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: a numerical closeup[END_REF] to ease the comparison between distributions. We also report the mean (μ) and standard deviation (σ) in this Figure . It is worth to notice that the mean value for the M1 simulation is of the same order that the grid size, while M3 could resolve 395 its mean diameter with almost three cells. If we now consider the volume-weighted PDF in Figure 10, the volume fraction of the small droplets is negligible: drops with ∆ D 30 < 2 µm represent less than 2% of the atomized volume. Moreover, even if there are just a few large liquid structures, which will eventually break up, these are much more significant in terms of volume fraction. For longer simulation times, the number of large drops will be higher, and the effect of small droplets will be less significant. These facts explain why the SMD in Figure 3 are considerably This graph shows that, at t = 11 µs, the M1 simulation predicts that the atomized liquid volume will be 30% higher than the one predicted by the other simulations (about 4.5% of the injected volume), which is consistent with the (e) and (g) graphs in Figure 3. At this instant, M2 and M3 present good agreement for all the diameter ranges. Both simulations predict that less than 4.5% of the injected liquid will atomize in drops.

Figure 12 shows the sphericity distribution normalized by the number of drops and liquid volume.

We should recall that lower ψ values correspond to higher deformation.

For example, a cylinder with a 5:1 length-to-radius ratio has ψ ≈ 0.69. This structure is unstable and will eventually suffer fragmentation or turn spherical due to surface tension. Therefore, the number of low ψ ligaments and sheets will only be significant if their production rate is equal or higher than their fragmentation rate. We can see in the graph that highly deformed structures (ψ ≤ 0.7) represent almost 30% of the drop count for M1 and more than 70% of the atomized volume. In contrast, M2 predicts that around 31% of the atomized volume will be highly deformed, much closer to the M3 prediction of 25%. Another relevant observation on these results is that more than 45%, both in volume and drop count, of the atomized liquid structures, are almost spherical (ψ ≥ 0.9) for M2 and M3, while only 10 % of the atomized volume has this characteristic in M1. To summarize, coarse simulations cannot resolve the stable structures and have a higher production rate of skewed ligaments, the ψ distribution is shifted to lower values. In contrast, more than 60% of the drops on M3 present a stable shape and will probably not suffer further fragmentation.

All these results concern a single instant on spray development. As Figure 3 shows, t = 11µs is just by the end of the initial transient development of jet topology. To get a better comparison of the three cases during this process, Figure 13 presents the temporal evolution of the accumulated drops volume fraction grouped in ranges of 2µm

Final remarks and future work

We presented Detailed Numerical Simulations of a pulsating round liquid jet for a regime based on diesel injection at low velocities using the VOF technique with AMR. Based on these results, we described the breakup mechanisms and discussed how poor grid resolution accelerates fragmentation by merging interfaces that meet in a single cell. This effect increases the atomization rate of the jet and over-predicts its surface density.

We then analyzed the impact that grid resolution has on the statistical features of the spray generated by atomization. We found that it affects not only the capability of generating small droplets but also the frequency of ligament detachment. As a consequence, a simulation with coarser mesh will have less small drops, more large liquid structures, and probably a higher atomized volume. We found that comparing the Probability Distribution Function weighted by volume highlights both issues. The cumulative PDF, presented in Figure 11, is a useful tool to evaluate the atomization error due to mesh resolution and its extension to Figure 13 shows the trend of this error through time.

For the analyzed primary atomization problem, with Re l = 5800 and We l = 11600, simulations with a grid resolution of ∆ = D inj /137 (M2) and ∆ = D inj /274 (M3) present drop size distributions with similar volume fraction in the small diameter range (D 30 < 4 µm), which are of order 10% by the end of the initial tip development process. The agreement through all the diameter ranges (less than 2%, as seen in Figure 13) indicates that using higher refinement levels will affect a small volume fraction of the drops, mainly the under-resolved structures. The prediction of shape, based on sphericity, presented higher discrepancies: M2 stable (ψ > 0.7) volume fraction was 5% higher than the M3 prediction. Results from the M1 simulation differ more than 40%, which impairs the capability of this simulation to predict, for example, the stability of drops in the dense spray region. For the particular application of designing a near nozzle fuel injection DPM, for example, the small discrepancies between M2 and M3 would not make significant differences, but using M1 results would be unacceptable.

We also found that the increase of the atomization rate due to poor mesh resolution affects the position of the drops: faster fragmentation of liquid sheets from the jet tip reduces the radial expansion of the spray, as shown in Figure 7.

Based on these results, we briefly discuss each of the questions at the end of section 1.

Regarding the determination a priori of the grid size required to perform DNS of primary atomization: we have proposed a set of length scales, defined by the problem dimensionless groups, that give the minimal refinement requirements to reproduce turbulence, deformation and breakup phenomena. If the grid cannot resolve any of these scales, especially those related to deformation and breakup, the simulation will not capture the breakup mechanisms properly. Nevertheless, using grid sizes even smaller than these reference lengths does not guarantee that the numerical error will be negligible, that will also depend on the accuracy requirements.

Regarding the assessment of numerical accuracy without running a mesh convergence analysis, we believe that a reliable parameter to quantify the impact that under-resolved structures have on the overall result is the volume fraction of the smallest relevant drops. This criterion would be valid for most of the applications that care about areas and volumes and not on drop count or minimal drop diameter.

The threshold to classify under-resolved structures is also an open question that depends on the particular application. From a geometric perspective, drops with less than 2 cells per diameter cannot be represented with interface capturing methods; moreover, the numerical error of VOF-PLIC representation will probably be significant for drops with less than 8 cells per diameter. From a practical point of view, we need that drops preserve their physical behavior during processes such as deformation, fragmentation, and translation. In the context of atomization, if drops with 8 cells per di-ameter are present in a stable regime and if this stability is preserved in the numerical results we could consider that these drops are resolved enough for our purpose, even if their shape is not accurately represented.

Regarding how poor mesh resolution affects the global properties of the simulated sprays: we found that it affects the shape and position of both small and large liquid structures; hence, under-resolved simulations predict denser sprays with a lower fraction of stable drops.

It is worth to notice that the simulated time did not reach a fully developed breakup length, neither produced a statistically stationary spray. One of our future goals will be extending these simulations to analyze the behavior after the jet tip detachment. Moreover, we will also study the impact of turbulence on atomization. Another future use of these results is to develop Sub-Grid Scale models to perform Large Eddy Simulations in Basilisk.

Afterward, we hope to simulate fuel injection cases closer to practical applications (e. g. spray A) and analyze how the flow regime (mainly We l and Re l ) impacts on the mesh resolution required to bound numerical breakup.

context, [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] use linear stability theory to analyze the growth of a periodic perturbation (e. g. sine function) on a planar interface separating two different fluids on layers which have a relative velocity U . They give, considering inviscid flow, the following dispersion relation:

ω(k) = ρ g ρ l (∆U ) 2 (ρ g + ρ l ) 2 k 2 - σ (ρ g + ρ l ) k 3 1/2 -→ λ c = 2π (ρ g + ρ l )σ ρ g ρ l (∆U ) 2 . ( .2)
where ω is the wave frequency, which depends on the flow parameters and on the perturbation wave number k = 2π/λ. The wave will only be unstable if ω(k) has a positive imaginary component. Hence, the critical value λ c is the minimal wavelength for which the disturbance grows and therefore the perturbations with λ < λ c will remain stable.

On the other hand, fragmentation mechanisms related to capillarity are also affected by shortwave instabilities, as discussed by [START_REF] Rangel | The linear and nonlinear shear instability of a fluid sheet[END_REF][START_REF] Senecal | Modeling high-speed viscous liquid sheet atomization[END_REF][START_REF] Dumouchel | Multi-scale analysis of simulated capillary instability[END_REF]. As the sheets and filaments reduce their thickness and speed, the perturbations evolve as Rayleigh-Taylor (sheets) or Plateau-Rayleigh (ligaments) instabilities. In both cases, the breakup occurs by a reduction of the section until the opposite surfaces of the liquid structure come into contact. In particular, [START_REF] Senecal | Modeling high-speed viscous liquid sheet atomization[END_REF] approximation of the wavelength with the maximum growth rate based on linear analysis of liquid sheets on gas streams:

k s = ρ g (U ) 2 2σ → λ s = 4σπ ρ g (U ) 2 . ( .3)
From this length, they also estimate the diameter of the ligaments as a function of the liquid sheet initial thickness (e 0 ), from which they also give an approximation for the size of drops produce from ligament breakup,

D lig = 8e 0 k s → D drop = 1.88 D lig (1 + 3 Oh l ) 1/6 . ( .4)
The initial thickness of the sheet will, of course, depend on the jet geometry and flow conditions. A priori estimations could be done based on, for example, KHI wavelength. Nevertheless, the values reported on Table 3 consider e 0 as the thickness of sheets still attached to the jet core, as can be seen in Figure .14 which shows a side view of the VOF field on the middle plane. It can be noted that as the sheet moves away from the tip its thickness is reduced, forming films of less than 2 µm thick.

  The properties of the flow (ρ and µ) are density and viscosity. The last term on the right-hand side on the Navier-Stokes Equation 2 represents the surface tension force:

Figure 1 :

 1 Figure 1: Flux computation example, considering a single square cell. Left: VOF flux computation on the left face, where normal velocity is u f , of a cell with interface normal n s . Right: schematics on how the density field affects momentum transport.

Figure 2 :

 2 Figure 2: Isometric view:domain and detail of initial evolution (t = 0.2 µs) with adaptive refinement

Figure 3

 3 Figure 3 shows a side view of the interface at regular time intervals of 5 µs for the M2 simulation, at subfigures (a-d), and the evolution of some spray parameters for the three simulations. The snapshots in the figure let us grasp the general aspects of the core deformation and atomization. The frames (a) and (b) show that the tip deformation requires more than 5 µs to produce the first liquid sheets. Between 5 and 10 µs, the first group of drops and short ligaments detaches from the tip. This ring of liquid structures expands radially and moves on the axial direction much slower than the rest of the jet, as can be seen in frames (c) and (d). Between 15 and 20µs a second, denser, group has formed. It is worth to notice that at this time, the drops in the first ring are almost spherical.

Figure 3 :

 3 Figure 3: Side view of M2 simulation (a-d). Temporal evolution of drop volume fraction in linear (e) and logarithmic (f) scales, to compare the total atomized volume and the volume of the small drops respectively. Evolution of the spray Sauter Mean Diameter (SMD) for the three mesh resolutions (g).

Figure 4 :

 4 Figure 4: Core and drops Sauter Mean Diameter (SMD) for the three mesh resolutions.

Figure 5 :

 5 Figure 5: M2 simulation: side view of the jet at t = (14.4, 14.6, 14.8, 15.0)µs, colored by axial velocity (u x ). The top row shows the jet's core and detached drops, where a hole formation and expansion process is highlighted with a circle. The middle row shows the evolution of the jet core, filtering the drops out; the oval here highlights the detachment of two ligaments due to capillary effects. The bottom row shows a zoomed image of the head region highlighted by a rectangle in the top row, showing the hole dynamics with more detail.

Figure 6 :

 6 Figure 6: f (top) and u x (bottom) fields at plane z = 0µm for t = 20µs of M2 simulation. The box on the top-left corner contains an isometric zoomed view detail of the near nozzle region. The box on the right shows an isometric view of the jet head.

Figure 6

 6 Figure 6 also shows the growth of the long-wavelength perturbation in the near nozzle region (x ∈ [0, 2D inj ]), where shear stresses induce a Kelvin

Figure 7 :

 7 Figure 7: Mesh convergence: side view colored by velocity (left) and scatter histograms showing radial mass distribution on drops (right) for t = 11 µs. Simulations M1 (top), M2 (middle) and M3 (bottom).

Figure 8 :

 8 Figure 8: Drop size (D 30 ) histograms for the three simulations (M1, M2, M3) at three different times t = (6, 9, 11) µs. Bin width is 0.5 µm.

Figure 9 :Figure 10 : 20 Figure 11

 9102011 Figure 9: Log-normal PDF fitting for range [1-30]µm. Vertical lines at D 30 = 2∆. Bin widths based on power of 2, starting from [0.125-0.25]µm

Figure 11 :

 11 Figure 11: Cumulative distribution weighted by volume at t = 11 µs. Vertical lines placed atD 30 = 2∆.

Figure 12 :

 12 Figure 12: Sphericity distribution, weighted by drop number and atomized volume.

Figure 13 :

 13 Figure 13: Evolution in time of the cumulative distribution functions of the three simulations M1(red), M2(green), and M3(blue). The drops are grouped in ranges from 0 to 20 µm with a width of 2 µm. The first group of drops (D 30 ∈ (0, 2)µm) is gray colored to highlight that it contains the under-resolved drops. The right side of each graph shows the cumulative volume fractions at 12 µs.

Figure . 14 :

 14 Figure .14: Sheet thickness measurement over VOF field with on the M3 simulation.

Table 1 :

 1 Physical properties of the system

Table 2 :

 2 Dimensionless description of the case

	ρ *	µ *	Re l	Re g	We l We g	Oh l	Ma g
	27.8 83.3 5800 25000 11600 430 0.0186 0.175

Table 3 :

 3 Characteristic length scales.

	l	η g	η l	ζ	λ KH c	λ s D lig
	[µm] 0.078 0.131 9.71 1.56 3.01 2.69

  . In this representation, for example, we see that the M2 simulation (green) has only 15% of the liquid volume in drops with D 30 < 2µm at t = 1µs, and we see that this fraction decreases with time, reaching 2% at 12µs. These graphs show that, for all the cases, liquid structures with D 30 > 20µm are formed after 8 µs. The M1 case presents The biggest difference between these results in the ranges of D 30 < 10 µm is 2% of the total volume. Taking the M3 result as a reference, M1 underestimates the volume fraction of D 30 < 10 µm by

	both cases these fractions are approximately 2% by t = 11µs. For the M1
	simulation, the volume fraction of drops in this range is less than 0.5% after
	2 µs.
	Even if the size distribution is quite similar for M2 (∆ = 0.73µm) and
	M3 (∆ = 0.37µm) grids, it is worth noticing that trends are not the same
	and the breakup length has not yet been reached. Moreover, at t = 10µs
	for example, the fraction on the range (D 30 > 20µm) shows the biggest
	discrepancy between M2 and M3. At the final sampling time (t = 12µs)
	the volume fraction on the lowest range (D 30 < 2µs) is almost negligible on
	both M2 and M3 cases. approximately 15%.
	the earliest generation of big ligaments and the highest volume fraction of
	large structures (e. g. drops with D 30 > 20 µm represent 40% of the volume
	at t = 9.4µs).
	Regarding small droplets, we can notice that M2 and M3 present similar
	behavior after 9µs: the volume fraction on that range decreases slowly; in

Appendices

Appendix .1. Estimation of flow characteristics scales

Turbulence is one of the main features of atomization that produce multiscale behavior. As described by [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF], the Kolmogorov and Hinze scales are relevant on atomization. [START_REF] Kolmogorov | A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF], based on self-similar turbulence theory, gives an estimate of the length scale of viscous dissipation phenomena. Following [START_REF] Pope | Turbulent flows[END_REF], this value can be estimated from the relevant macro scale (D inj ) and the Reynolds number. [START_REF] Hinze | Turbulence. an introduction to its mechanism and theory[END_REF], on the other hand, defines the maximum diameter of drops that can resist pressure fluctuations without breaking up. In this context,

where η is the Kolmogorov scale, is the turbulent energy dissipation rate and ν is the relevant kinematic viscosity. The Hinze scale ζ is, as aforementioned, a measure on the biggest droplets that can resist turbulence fluctuations. It should be noted that in the problems where a mixing layer develops, the average droplet size is usually considerably less than ζ, as explained by [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF].

On the other hand, jet deformation and fragmentation mechanisms involve the growth of hydrodynamic instabilities. On primary atomization problems, the effects of tangential relative velocity to the interface, related to Kelvin-Helmholtz Instability (KHI), are particularly important. In this