J.-L Estivalezes

W Aniszewski

F Auguste

Y Ling

L Osmar

J.-P Caltagirone

L Chirco

A Pedrono

S Popinet

A Berlemont

J Magnaudet

T Ménard

S Vincent

S Zaleski

A phase inversion benchmark for multiscale multiphase flows

Keywords: phase inversion flow, water-oil flow, atomization, multiphase flow benchmark, Volume Of Fluid Method, Level-Set Method

A series of benchmarks based on the physical situation of "phase inversion" between two immiscible liquids is presented. These benchmarks aim at progressing towards the direct numerical simulation of two-phase flows. Several CFD codes developed in French laboratories and using either Volume-of-Fluid or Level-Set interface tracking methods are used to provide physical solutions of the benchmarks, convergence studies and code comparisons. Two typical configurations are retained, with integral scale Reynolds numbers of 1.37 10 4 and 4.33 10 5 , respectively. The physics of the problem are probed through macroscopic quantities such as potential and kinetic energies, or enstrophy. In addition, scaling laws for the temporal decay of the kinetic energy are derived to check the physical relevance of the simulations. Finally the droplet size distribution is probed. Additional test problems are also reported to estimate the influence of viscous effects in the vicinity of the interface.

Introduction

The physical validation or verification of numerical codes devoted to multiphase flows and immiscible fluids (by immiscible we mean that interfaces remain sharp) is a major concern of modern CFD. A lot of work was carried out since the 1990s to validate interface-tracking algorithms in cases where the velocity field is given analytically [START_REF] Rider | Reconstructing volume tracking[END_REF][START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF]. Comparisons have also been achieved with available experiments or predictions of linear stability theories [START_REF]Test-cases for interface tracking methods[END_REF]. However, in these situations, the topology of the interface is simple and the flow problem involves a single scale and a single physical effect such as the continuity of stresses across the interface with the two-phase Poiseuille flow [START_REF] Benkenida | A method for the simulation of two-phase flows without interface reconstruction[END_REF], the Laplace law for a spherical or circular drop [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF], the oscillation modes of a free surface in the linear regime [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF], the head-on coalescence of two drops [START_REF] Ashgriz | Coalescence and separation in binary collisions of liquid drops[END_REF][START_REF] Tanguy | Application of a level set method for simulation of droplet collisions[END_REF], or the break-up of a low-velocity round jet due to the Rayleigh-Plateau instability [START_REF] Delteil | Numerical investigations in Rayleigh breakup of round liquid jets with VOF methods[END_REF] to mention just a few. None of these problems puts into play the interaction between fluctuating interfaces and a turbulent flow, which is a key point in many realistic cases involving a two-phase dynamics. Indeed, such multiphase interactions occur in many environmental and industrial applications. Let us just mention the atomization of a liquid jet in an engine [START_REF] Lasheras | Liquid jet instability and atomization in a coaxial gas stream[END_REF][START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF][START_REF] Eggers | Physics of liquid jets[END_REF], phase separators in chemical engineering processes, the ladle-based elaboration of steel, the boiling crisis in nuclear reactors [START_REF] Gilman | A self-consistent, physics-based boiling heat transfer modeling framework for use in computational fluid dynamics[END_REF], complex bubbly flows including the flow past Taylor bubbles in pipes [START_REF] Lu | Direct numerical simulations of multifluid flows in a vertical channel undergoing topology changes[END_REF], droplet or asteroid impacts [START_REF] Josserand | Drop impact on a solid surface[END_REF], wave impact on structures or the ubiquitous phenomenon of wave breaking [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF][START_REF] Mostert | High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplet production[END_REF][START_REF] Mostert | Inertial energy dissipation in shallow-water breaking waves[END_REF]. In most multiphase flows of environmental or engineering relevance, the multiscale character of the flow is a key issue to be handled by computational approaches. The interfacial length scales are associated on the one hand to large interface structures such as jets, films or large drops. On the other hand, small interfacial scales also exist, corresponding to a small-scale dispersed phase, i.e. small bubbles or droplets. These two widely distinct families of interfacial scales interact in a nonlinear way through ligament break-up or drop/bubble coalescence and the unsteady or turbulent nature of the carrier fluid motion plays a crucial role in these interactions. The definition of a reference benchmark for these multiscale sharp interface problems appears to be an important issue.

Homogeneous Isotropic Turbulence (HIT) often serves as a reference problem in the verification of singlephase DNS codes. In analogy, one would like to define a reference problem for the verification of CFD multiphase flow codes devoted to immiscible fluids. Our first idea was to try to define the reference problem on the basis of a well-documented experiment, such as atomization of liquid sheets or jets. However, the quality of the modeling and simulation of atomization is very sensitive to inflow conditions, in particular upstream turbulence, boundary layers or inlet pump oscillations [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Matas | Shear instability of an axisymmetric air-water coaxial jet[END_REF]. Although numerous experimental studies can be found in the literature [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF][START_REF] Lance | Turbulence in the liquid phase of a uniform bubbly air-water flow[END_REF][START_REF] Chanson | Air-water flow measurements with intrusive, phase-detection probes: Can we improve their interpretation[END_REF] they are plagued by two difficulties that make them unsuitable as references cases: either they have too simple flow structures (as for dam break on a dry bottom) or they are not sufficiently well documented or instrumented to provide all the information required for properly defining boundary and initial conditions for simulations. Moreover, many experiments lack instantaneous pressure, velocity and interface measurements (dam break on a wet bottom [START_REF] Kleefsman | Volume-of-fluid based simulation method for wave impact problems[END_REF]). As a consequence, we decided to define a synthetic reference case unrelated to a specific experiment. It is the phase inversion of an initial cube of light fluid inside a larger cube of a heavier fluid, both contained in a closed box (Figure 1). (Although an experiment in this configuration could in principle be considered, the setup resembles that of an oil and water phase separator or settling pond, or of a model Rayleigh-Taylor experiment, and Figure 2: Example of simulated phase inversion in a closed box [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF] illustrating the multiscale character of the interfacial flow. (The case is similar to our case 2 but with smaller surface tension.) with a 90 • -change in the gravity orientation, that of the dam-break problem [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF].) In this reference case, initial and boundary conditions are simple and well defined. We also simplify the problem somewhat by using a liquid-liquid density ratio (ρ 1 /ρ 2 = 9/10). There are two reasons for this simplification: 1) with the air-water density ratio, the air inertia is not large enough to defeat the liquid inertia and gravity (this is why the dam-break problem is not producing much atomization); and 2) the numerical issues are considerably simplified with moderate density ratios. (We note that a study of atomization with liquidliquid density ratios was published recently [START_REF] Constante-Amores | Direct numerical simulations of transient turbulent jets: vortex-interface interactions[END_REF]. The present study aims at the same kind of atomization but "in a box".) We wish to have complex flow features such as turbulence, interface rupture, coalescence, and multiple scales. Typically, the higher the Reynolds number Re and Weber number We the more interesting and complex will the flow be. Moreover, we note that in single-phase turbulence there is a large-Reynolds-number state where the range of scales is very wide and the statistical features of the flow on large and intermediate scales become relatively independent of Re, a kind of universality. Is there such an equivalent universal state at large Re and We in multiphase turbulence? This is not certain but there is clearly a benefit to reach an analogous state there. Such a fully-developed multiphase turbulence would be the relevant flow state with respect to common natural sciences or engineering situations. An example of the phase inversion flow at very high Re and We in which a broad range of interfacial scales coexist, multiple break-up and coalescence events take place, and large-scale sloshing motions are observed together with small dispersed droplets interacting with larger drops, is illustrated in Figure 2 (taken from [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF]).

However, with very high Re and We numbers, viscosity is considerably smaller than what is required to obtain fully-resolved simulations in multiphase as well as in single-phase flow and one cannot hope to perform truly a Direct Numerical Simulation (that is, a simulation with all length scales resolved and quantities such as dissipation and enstrophy converged). There is thus a dilemma or a balance to strike between the interest for "fully developed" multiphase turbulence and the desire to have a fully converged Direct Numerical Simulation for reference purposes. We chose the parameters given below in Table 3. In retrospect, once simulations are performed at great costs, it appeared clear (especially from the investigation of enstrophy convergence and particle size Probability Density Functions (PDF) as we shall see) that the current work inclines towards the high-Re side of the dilemma. We define two simulation cases in the phase inversion setup. Our simulations, especially in the second case (which are performed without explicit turbulent or subgrid droplet modeling), are finely-resolved implicit Large Eddy Simulations or Detailed Numerical Simulations. (Note that other combinations of parameters were used in [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Vincent | Test case n15: phase inversion in a closed box (pn, pe)[END_REF] and [START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF].) We also note that the second case presented in this work is a good illustration of an atomization process although at higher Re and We than previous works on atomization that got closer to the "converged" side of the dilemma [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF][START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF][START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF][START_REF] Tomar | Multiscale simulations of primary atomization[END_REF]. In a way, one may see the second case of the current problem as the study of "atomization in a box". An unexpected discovery of this study is that the non-convergence of certain flow characteristics is actually not only due to large Re and We numbers but also to the dependence of thin sheet breakup on the grid size, as shown in the analysis of droplets sizes in case 1 below.

Within the literature on high Re turbulent multiphase flow simulation one may distinguish two categories of works: 1) on diffusive interfaces and 2) on interfaces remaining sharp over time (immiscible fluids). The first class of works has investigated for example the interaction between Rayleigh-Taylor (RT) instabilities and turbulence [START_REF] Cook | The mixing transition in Rayleigh-Taylor instability[END_REF][START_REF] Ramaprabhu | A numerical study of the influence of initial perturbations on the turbulent Rayleigh-Taylor instability[END_REF][START_REF] Boffetta | The ultimate state of thermal convection in Rayleigh-Taylor turbulence[END_REF][START_REF] Abarzhi | Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids[END_REF] and a benchmark has been published in 2004 by Dimonte and co-workers [START_REF] Dimonte | Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids[END_REF]. This reference work involves physical phenomena similar to those occurring in the phase inversion configuration proposed here. Much attention has been paid to characterizing, with various models and numerical methods, the effects of initial perturbations of the interface, volume fraction profiles, self-similar bubble dynamics or mixing behavior. However, the major difference between reference [START_REF] Dimonte | Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids[END_REF] and the current work is the immiscible nature of the fluids and the discontinuous character of the interfaces (of course from a macroscopic point of view). It has been clearly demonstrated [START_REF] Abarzhi | Rayleigh-Taylor turbulent mixing of immiscible, miscible and stratified fluids[END_REF] that the diffuse character of the fluids has a major effect on the global dynamics of the turbulent two-phase flow in terms of spikes position for example, and that this effect increases over time. In diffuse modeling of RT flows, simulations mainly involve gas stratifications whereas our focus here is on immiscible fluids (with at least one liquid phase in the considered problem). Our benchmark thus belongs to the second category of problems. The phase inversion configuration has previously been investigated [START_REF] Vincent | Test case n15: phase inversion in a closed box (pn, pe)[END_REF] without the objective of analyzing the physical aspects of the multiscale character of the flow, nor comparing various interface tracking methods, but with the purpose of obtaining a priori estimates of subgrid terms that appear and must be modeled when LES is used in the framework of two-phase flows of immiscible fluids [START_REF] Aniszewski | A new approach to sub-grid surface tension for LES of two-phase flows[END_REF]. We note that for sharp interfaces and complex flows, RT simulations also exist [START_REF] Tryggvason | Numerical simulations of the Rayleigh-Taylor instability[END_REF][START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF], and intermediate modeling standing between diffuse and sharp interface approaches have also been proposed for simulating RT instabilities [START_REF] Jacqmin | Calculation of two-phase Navier-Stokes flows using phase-field modeling[END_REF]. The major concern was the model and numerical methods and there was no consideration of the turbulence characteristics of the flows. To our knowledge, neither benchmark nor detailed simulation works have been published about RT interaction with turbulence in the framework of sharp liquid-liquid interfaces.

The study reported here has two objectives. The first is to analyze convergence of the simulations with grid refinement, thus producing a reference set of results. Since no experimental measurements seem to exist for the phase inversion problem considered here, this first objective is a necessary step. The second objective is to compare the results of five distinct in-house DNS multiphase flow codes to the reference results. The results should highlight the capabilities but also the limitations of up-to-date numerical methods and interface models to tackle multiscale, multiphase, flow problems. The outline of the paper is as follows: the phase inversion setup is defined in the following Section 2. The codes and the numerical methods are described in Section 3. The reference simulations and results are in Section 4 while the results of all the codes are in Section 5, followed by concluding remarks in Section 6.

Definition of the phase inversion benchmark

We define the benchmark in this Section, following the motivations given in the previous section.

Geometry, initial and boundary conditions

As described in Figure 1, an initial cubic blob of light liquid, referred to as fluid 1, is placed in the bottom part of a cubic box filled with a heavier liquid, referred to as fluid 2. The size of the box is (H, H, H), while the size of the blob of light fluid is (H/2, H/2, H/2). All the outer walls are considered as free-slip impermeable walls, so that the normal velocity is zero and the tangential components obey a symmetry, i.e. shear-free, condition. A π/2 static contact angle is implicitly prescribed on the walls through the boundary condition imposed on the volume fraction (in the VOF approach) or the level-set approach. This is implemented by enforcing a null Neumann (symmetry) boundary condition n • ∇φ = 0, where φ is either the volume fraction or the level-set function. The gravitational acceleration is chosen as g = (0, -9.81, 0). Besides the simplicity of the initial and final configurations of the phase inversion problem, another of its advantages is that it provides the possibility to observe multiple coalescence and break-up events, although that introduces an additional complexity since in a simulation, coalescence is conditioned by numerical rather than physical properties.

Physical parameters of the test cases

Two test cases with different characteristics of the fluids are considered. The fluid properties may be expressed in terms of Re and We. In order to define these numbers, a velocity scale is required. We obtain it from energy considerations. In all of what follows, we use the characteristic or color function C that is defined with respect to fluid 1, so that C = 1 in fluid 1 and C = 0 in fluid 2. The potential energy of phase n (n = 1, 2) is

E p,n = Ω C n ρ n gy dV , (1)
where C n is the characteristic function of phase n, i.e. C 1 = C, C 2 = 1 -C, and g = g . The kinetic energy of phase n is

E k,n = 1 2 Ω C n ρ n u 2 dV , (2)
where u stands for the velocity field. Due to the simple topology of the interface in the initial and final stages, the potential energies are trivially found. At t = 0, E i p,1 = ρ

E k,1 + E k,2 ≤ E i p,1 + E i p,2 -E f p,1 -E f p,2 = 11 128 (ρ 2 -ρ 1 )gH 4 , (3)
and an upper bound on the velocity is ||u|| 2 < (2E k /ρ 1 ) 1/2 . (Note that the initial surface energy 3σH 2 /4 provides a small contribution to the total energy; however we neglect surface energy in this reasoning for simplicity.) Identifying the velocity with the upper bound, we obtain the velocity scale

U K = √ 11 8 ρ 2 -ρ 1 ρ 1 gH .
The dimensionless numbers are then defined based on U K and properties of fluid 2, namely

Re = ρ 2 HU K µ 2 , ( 4
) We = ρ 2 HU 2 K σ . (5)
As is customary in atomization problems, the Reynolds and Weber numbers may be combined in a way that eliminates the velocity scale, yielding the Laplace number La=σρ 2 H/µ 2 2 = Re 2 /We. To allow comparison with previous work on this setup, we may also adopt the same velocity scale as in previous studies by some of us [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Aniszewski | A new approach to sub-grid surface tension for LES of two-phase flows[END_REF], namely U g = (ρ 2 -ρ 1 )/ρ 1 gH/2, and a characteristic time scale t c = H/(2U g ).

The viscosities and densities considered here correspond approximately to a light oil or hydrocarbon (fluid 1) and water (fluid 2). These configurations are interesting in so far as they provide interface and turbulent characteristics similar to those of jet atomization with simpler and more easily controlled initial conditions, and are as a result more convenient for numerical simulations. These liquid-liquid two-phase flows make it possible to select easily a Reynolds number range smaller, or even equivalent, to those encountered in real atomization processes while exhibiting equivalent dynamics. The characteristics of the fluids and related dimensionless numbers are given in Tables 1, 2 and 3. Case 1 corresponds to a configuration with moderate fragmentation. The box used in case 2 has a larger size and thus yields a larger Reynolds number (by a factor of 30 compared to case 1). The Weber number is also higher in case 2, (by a factor of 10 compared to case 1) but is increased less than the Reynolds number since surface tension was also artificially increased. This was done specifically to avoid having an exceedingly large Weber number which would lead to excessive fragmentation of the interface, making it even harder than it is now to capture most of the small-scale structures and droplets.

The small scales of the flow may be characterized by the Kolmogorov scale η K = Re -3/4 H and the Hinze scale η H = We -3/5 H. Both scales are given in Table 3. While the Kolmogorov scale is very small, and prevents any attempt at DNS with the grids used in the paper (up to 2048 grid points in each direction), the Hinze scale is larger, and as η H /H > 1/100, it should be easily possible to resolve droplets at the Hinze scale. We note, however, that droplet sizes much smaller than the Hinze scale have been observed in numerical simulations of atomization [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF] and that one can estimate the diameter of the Taylor-Culick rims involved in the atomization process to be d T = We -1 H η H (see [START_REF] Marcotte | Density contrast matters for drop fragmentation thresholds at low Ohnesorge number[END_REF] for a detailed derivation of d T and its significance). Notice that for both reference simulations below, the cell size ∆x is such that ∆x < d T /4.

Quantities of interest

Several quantities characterizing the evolution of the flow field are of primary interest. Among these physically relevant quantities we select the following.

• The potential energy which helps to measure the advancement of the relaxation from a high energy state with the light fluid located below the heavy fluid to the equilibrium state with the light fluid on top, and provides a characterization of the density stratification.

• The kinetic energy which characterizes the agitation of the system and allows to monitor the amount of dissipation.

• The total mechanical energy minus the surface energy which is the sum of the previous two. As the Weber numbers are large, the surface energy is negligible and the above sum characterizes well the dissipation of energy in the system. It can also be used to characterize how far the system is from the possible universal high Re and We state discussed above.

• The time evolution of the volume integral of the enstrophy in both fluids. This quantity is defined as E r,n = 1 2 Ω C n ω 2 dV , with ω = ∇ × u denoting the vorticity.

Model and numerical methods

As is now well established, incompressible two-phase flows involving fluid-fluid interfaces and Newtonian fluids can be modeled by a single set of incompressible Navier-Stokes equations with phase-specific densities and viscosities and capillary forces, together with the transport equation of the phase function C. The resulting model takes implicitly into account the mass and momentum jump relations at the interface [START_REF] Delhaye | Jump conditions and entropy sources in two-phase systems. Local instant formulation[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF], whereas the continuity of the fluid-fluid and fluid-solid interfaces is taken into account by the C equation. The entire set of equations reads

∇ • u = 0 , (6)
ρ ∂u ∂t + (u • ∇)u = -∇p + ρg + ∇ • µ ∇u + (∇u) T + F st , (7)
∂C ∂t + u • ∇C = 0 , ( 8
)
where p is the pressure, F st is the interfacial force per unit volume and ρ and µ are the local density and viscosity of the two-phase medium, respectively. Capillary effects are inserted in the source term F st in the form F st = σκn i δ i , where σ denotes the surface tension, κ is the local mean curvature of the interface, n i is the unit vector normal to the interface and δ i is the interface delta function [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF].

The above one-fluid system is the classical model for multiphase incompressible flows with sharp interfaces and effects of uniform surface tension. Localizing the interface requires solving a transport equation for the characteristic function C. This can be performed using the volume fraction (also noted C in VOF methods) or using a continuous level-set function φ such that φ = 0 at the interface and φ > 0 (resp. < 0) in fluid 2 (resp. 1). In this case, the characteristic function is obtained as C = H(φ) where H is the Heaviside function [START_REF] Osher | Level set methods: An overview and some recent results[END_REF]. Note that the above equations are distinct from those involved in various models of multiphase turbulence, that may involve Reynolds stresses, terms accounting for unresolved eddies, turbulent viscosities and diffusivities, and non-sharp interfaces.

The Thétis code

Thétis is a CFD code developed in the TREFLE department of the I2M laboratory. It solves the one-fluid Navier-Stokes equations discretized with implicit finite-volumes on an irregular staggered Cartesian grid.

A second-order centered scheme is used to approximate the spatial derivatives while a second-order Euler or Gear scheme is used for time integration [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. All terms are written at time (n + 1)∆t, except the inertial term which is expressed in the following semi-implicit manner:

u n+1 • ∇u n+1 ≈ 2u n -u n-1 • ∇u n+1 . (9)
It has been shown that this approximation allows to reach second-order convergence in time [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. The coupling between velocity and pressure is ensured by using an implicit algebraic adaptive augmented Lagrangian method [START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF]. The augmented Lagrangian methods used in this work are independent of the chosen discretization and could for instance be implemented in a finite-element framework [START_REF] Bertrand | A three-dimensional fictitious domain method for incompressible fluid flow problems[END_REF]. In two dimensions, the standard augmented Lagrangian approach [START_REF] Fortin | Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites[END_REF] can be used to deal with two-phase flows since direct solvers [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF] are efficient in this case. However, as soon as three-dimensional problems are considered, the linear system resulting from the discretization of the augmented Lagrangian terms has to be treated with a BiCG-Stab II solver, preconditioned by a Modified and Incomplete LU method [START_REF] Van Der | BiCGSTAB: a fast and smoothly converging varian of BiCG for the solution of non-symmetric linear systems[END_REF].

As for the interface tracking and advection of C, two different volume of fluid (VOF) methods have been implemented in Thétis [53] [54]. They are evaluated here. The above numerical methods and the one-fluid model have been validated in previous works, e.g. [START_REF] Randrianarivelo | Numerical modelling of the solid particle motion using a new penalty method[END_REF] and [START_REF] Vincent | Local penalty methods for flows interacting with moving solids at high Reynolds numbers[END_REF].

The Gerris Flow Solver code

The Gerris Flow Solver (GFS) is an open source code implementing finite volume solvers on an octree adaptive grid together with a piecewise linear VOF interface-tracking method [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The variables are located on the octree grid in a collocated manner. Time-advancement is achieved through a Bell-Collela-Glaz scheme for advection and a Crank-Nicolson algorithm for viscous stresses. Incompressibility is satisfied at the end of each time step through a projection method, with a second MAC projection of the face velocities. The resulting code is second-order accurate in both time and space for single-phase flows. The simulations reported in this paper use a VOF method of the piecewise linear type, in which the interface segments are reconstructed using the mixed-Youngs-centered (MYC) approximation [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry[END_REF]. While not the most accurate for very fine grids, the MYC approximation is easily implemented when using only information from the nearest-neighbour cells, an important advantage when domain-decomposition on octrees is used. Advection of C is performed using the Lagrangian-Explicit or "CIAM" method first published by Li [START_REF] Li | Calcul d interface affine par morceaux (piecewise linear interface calculation)[END_REF] and discussed in [START_REF] Scardovelli | Interface reconstruction with least-square fit and split Lagrangian-Eulerian advection[END_REF][START_REF] Aulisa | A geometrical area-preserving volume of fluid advection method[END_REF][START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. While not volume-conserving to machine accuracy, it has very good volume and mass conservation properties. Surface tension is a vexing question in multiphase flow simulations. As density and viscosity ratios become very large or very small, the simulations become increasingly difficult with standard methods [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. The problem has been considerably improved in Gerris as a result of the use of Height-Function methods [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] and a so-called balanced-force algorithm [START_REF] Francois | A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework[END_REF][START_REF] Renardy | PROST: A parabolic reconstruction of surface tension for the volumeof-fluid method[END_REF]. For a review of surface-tension methods including a discussion of the differences between those used in Gerris and in other codes, see [START_REF] Popinet | Numerical models of surface tension[END_REF]. Gerris has been used with success in two-dimensional [START_REF] Fuster | Instability regimes in the primary breakup region of planar coflowing sheets[END_REF] and three-dimensional [START_REF] Fuster | Simulation of primary atomization with an octree adaptive mesh refinement and VOF method[END_REF][START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF] atomization and droplet impact studies.

The Jadim code

JADIM is a versatile code developed for a number of years at IMFT. In JADIM, the momentum equations are discretized on a staggered orthogonal grid using a finite-volume approach. Spatial discretization is performed using second-order centered differences. Time-advancement is achieved through a third-order Runge-Kutta algorithm for advection/source terms and a Crank-Nicolson algorithm for viscous stresses.

Incompressibility is satisfied at the end of each time step through a projection method. More details may be found in [START_REF] Calmet | High-Schmidt number mass transfer through turbulent gas-liquid interfaces[END_REF]. The resulting code is second-order accurate in both time and space for single-phase flows.

In two-and three-phase configurations, a VOF method with no interface reconstruction is used. The advection equation for C is solved using a Zalesak flux-corrected scheme [START_REF] Zalesak | Fully multidimensionnal flux-corrected transport algorithms for fluids[END_REF] split into successive onedimensional steps [START_REF] Bonometti | An interface-capturing method for incompressible two-phase flows. validation and application to bubble dynamics[END_REF] in order to avoid the well-known tendency of the multidimensional version of this scheme to distort iso-C surfaces. Owing to the splitting procedure, the overall transport scheme is not rigorously conservative. Therefore a local mass error control which improves upon the global control technique described in [START_REF] Bonometti | An interface-capturing method for incompressible two-phase flows. validation and application to bubble dynamics[END_REF] is employed here. The corresponding strategy is based on the detection of drop or bubbles (through an index function and a topological monitoring of each bubble or drop) and on an iterative solution which modifies the volume fraction within the transition regions (i.e. those in which 0 < C < 1) in order to keep the volume of each bubble or drop constant. Since no interface reconstruction step is involved, smearing of interfaces frequently occurs in high-shear regions. The original strategy described in [START_REF] Bonometti | An interface-capturing method for incompressible two-phase flows. validation and application to bubble dynamics[END_REF] is employed to keep this smearing within reasonable bounds. In the most difficult cases (e.g. break-up and coalescence), this approach is supplemented by an antidiffusion technique in which the local volume is redistributed in the direction normal to the interface so as to eliminate spurious values of the volume fraction outside interfacial regions.

The capillary force is represented using a modified version of the Continuum Surface Force (CSF) model [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. This modification consists in defining the interfacial region as that in which the capillary force takes non-zero values rather than that in which the volume fraction takes intermediate values. By doing so, it improves the control of the thickness of the transition region between the two fluids. The above strategy for the transport of C makes the approach selected in JADIM intermediate between VOF and Level Set techniques: it is based on the transport of the local volume fraction of one of the fluids but this transport is not achieved in a strictly conservative manner and no explicit reconstruction of the interface is carried out, the position of the interface being merely identified a posteriori with the C = 0.5 iso-surface.

The Archer code

Archer is a CFD code developed at CORIA laboratory, mainly devoted to multiphase flows. It was previously applied to compute atomization [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF], vaporization and mixing [START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF], and other complex interfacial flows [START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF]. The Level Set (LS) method [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] is used for tracking interfaces and shapes. It is based on a continuous distance function φ defined as the signed distance between any point of the domain and the interface. Similar to the volume fraction C, φ obeys a pure advection equation.

To avoid singularities in the φ field, the fifth-order conservative WENO [START_REF] Shu | Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws[END_REF] scheme is applied to discretize convective terms. When the LS advection is carried out, high velocity gradients can cause wide spreading or stretching of φ which then no longer remains a distance function. A redistancing algorithm [START_REF] Sussman | Improved level set method for incompressible two-phase flows[END_REF] is thus applied at every time step to restore the distance property of φ, i.e. |∇φ| = 1. Advancement of the φ-equation and the redistancing algorithm can induce mass loss in under-resolved regions [START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF]. This is the main drawback of LS methods. To improve mass conservation, the Coupled Level-Set Volume of Fluid (CLSVOF) method [START_REF] Sussman | A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows[END_REF] is used. The main idea of this approach is to benefit from the advantages of each tracking strategy: minimize the mass loss using VOF and keep a fine description of interface properties with the smooth LS function. (A study of the performance of various VOF variants was recently published [START_REF] Aniszewski | Volume of fluid (VOF) type advection methods in two-phase flow: A comparative study[END_REF] using Archer in the CLSVOF context; the phase inversion problem is also discussed therein.) The coupling between LS and VOF is maintained using a correction scheme based on a geometric reconstruction of both VOF and LS interfaces [START_REF] Aniszewski | A new approach to sub-grid surface tension for LES of two-phase flows[END_REF]. The LS method itself is coupled with a projection method for the incompressible Navier-Stokes equations, where the density and the viscosity depend on the sign of the LS function (with appropriate interpolations used in interfacial cells). To finalize the description of the two-phase flow, jump conditions across the interface are taken into account using the Ghost Fluid (GF) approach. In this approach, ghost cells are defined on each side of the interface [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] [77] and appropriate discretization schemes are applied to the jump of each quantity. As defined above, the interface is characterized through the distance function, and jump conditions are extrapolated on some nodes on each side of the interface. Following the jump conditions, the discontinued functions are extended continuously and then derivatives are estimated. Heaviside functions are designed according to the distance function in order to provide a characteristic function for the two-phase medium. They follow the functions proposed in formula ( 6) of [START_REF] Towers | Finite difference methods for approximating Heaviside functions[END_REF] . Discrete delta functions are built in a similar way and are regularized using a properly set cosine function. These techniques have been presented and validated in previous works, e.g. [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF].

The DyJeAT code

DyJeAT( Dynamics of Jet ATomisation) is an in-house computational fluid dynamics library developped at ONERA to deal with atomization processes in aircraft or launcher engines. It is based on the onefluid formulation for incompressible two-phase flows with specific interface capturing features. Classical projection methods are used to ensure the incompressibility constraint [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires[END_REF]. The spatial discretization is based on staggered uniform Cartesian grids for the velocity components, all others quantities such as density, pressure and level-set function φ being cell-centered. Advection terms in the momentum equations are approximated in a conservative way with 5 th -order accurate WENO schemes [START_REF] Shu | Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws[END_REF]. This particular choice was motivated by the robustness and low numerical dissipation of these schemes observed in DNS studies [START_REF] Trontin | Direct numerical simulation of a freely decaying turbulent interfacial flow[END_REF]. Viscous terms are discretized with second-order centered scheme. Time integration is performed with a second-order accurate Adams-Bahsforth scheme. The Poisson equation for the pressure is solved by a fast multigrid preconditioned conjugate gradient method [START_REF] Tatebe | Mgcg method, a robust and highly parallel iterative method[END_REF]. A level-set method [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] is used to capture the interface, implicitly defined as the zero level of the smooth function φ. Similar to the momentum equation, the level-set equation is solved with a 5 th -order WENO scheme for the spatial discretization and a 3 rd -order TVD Runge-Kutta scheme for time advancement. During its transport by the fluid velocity, the level-set function no longer remains a distance function. Since redistancing is mandatory for curvature computation, and as is standard practice, the level-set implicit function is repeatedly reinitialized to the signed distance. Redistancing in the level-set scheme, including the associated resolution issues and mass conservation properties, is discussed in detail in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles. Application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF] (see in particular figures 3.8 and 3.9 therein). Jump conditions at the interface for both pressure and viscous terms are handled with a Ghost Fluid approach [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson equation on irregular domains[END_REF]. All the codes used in the present work are based on finite volumes but with different variants of the methods. Obviously, each code can lead to a different result on a given grid, due to variations in the methods used. Table 4 is given to clarify the comparison of the model and numerical methods used. Schematically, the codes cover the broad range of numerical methods that can be found in the literature for simulating multiphase flows using the sharp interface approximation described above. This approach is sometimes called direct or detailed, to contrast it with other approaches that depart from the sharp interface approximation. Four of the codes use staggered grids while one (Gerris) is based on collocated unknowns.

Pros and cons of the different codes

Staggered grids have several advantages: a more accurate pressure solution, and the avoidance of spurious "red-black" velocity oscillations. However, Gerris has demonstrated high accuracy on test cases [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF] and ten years of experience with Gerris show almost no occurrence of spurious oscillations. Four codes are using a projection technique achieving a velocity-pressure coupling that contains a time splitting error, while one code (Thetis) uses an exact augmented Lagrangian approach. This technique is exact if the residual of the iterative solver is zero at machine accuracy. In terms of interface tracking, Gerris and Thetis are using a geometrical VOF-PLIC technique while Jadim employs a FCT scheme to directly approximate the advection of the volume fraction without reconstructing explicitly the interface geometry (this is why the corresponding approach may be thought of as Implicit VOF). DyJeAT uses the Level Set method whereas Archer couples this technique with VOF-PLIC in order to improve mass conservation and also to better evaluate advective effects in the momentum transport. Concerning capillary effects, the Ghost Fluid approach based on jump relations at the interface and the CSF formulation with height functions [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] are the most accurate. Thetis and Jadim use less sophisticated capillary force approximations that are known to generate larger spurious currents than the other methods. To finish with numerical methods, all projection steps except for Gerris are achieved with an iterative BiCGStab II solver preconditioned with a multigrid algorithm, while Thetis and the coupled augmented Lagrangian method use BiCGStab II and an incomplete LU preconditioner because the linear system is not symmetric. The major drawback of the Thetis solver is that its parallelization has a good speed-up until 1000 processors but can hardly handle in its present form a larger number of processors due to the ILU preconditioner. In contrast, the multigrid preconditioner used in the projection approach is nicely extendable in parallel computations until 100, 000 processors. In Gerris, the projection step is achieved thanks to an in-code multigrid Poisson solver.

Detailed numerical simulations using DyJeAT

We used the DyJeAT code to obtain results at very high resolution. The grids we used are still too coarse to reach the Kolmogorov scale but can be considered suitable to provide "detailed numerical simulations" offering a view of the large-scale vortical structures. The meaning of such detailed simulations for interfacial structures may be understood from a study of the convergence of the flow properties. As we shall see, some of the characteristics of the flow linked to its multiphase character, such as the PDF of droplet sizes, also seem to converge. We can thus consider that the detailed numerical simulations also offer a view of the large-scale properties of the interfacial structure. As a preliminary result relevant to both cases, we show in Figure 3 the evolution of the total mechanical energy

E m = E p,1 + E k,1 + E p,2 + E k,2 ,
obtained on a 512 3 grid in both cases which, as we shall see, are well converged. The initial energy is

E m (0) = E i p .
As t → ∞ one expects the system to converge to a state in which the light fluid occupies exactly the cuboid Ω up and the kinetic energy vanishes, so we expect

E m → E f m = E f p .
We rescale E m (t) for both cases as

E m (t) = E m (t) -E i p E f p -E i p .
The plot shows a remarkable similarity in both cases, especially at the end of the process where the two curves nearly superpose. Energy in case 1 is somewhat lower, presumably because a fraction of it was transferred to surface energy. For the same reason, one observes more fluctuations in the curve corresponding to case 1. The similarity between the two cases is an indication of the possibility of the Reand We-independent fully developed turbulent regime discussed in the introduction.

Case 1: a phase inversion problem with few fragmentation events

In this section we present the numerical results pertaining to case 1 obtained on 128 3 , 256 3 and 512 3 grids with DyJeAT. The characteristics of this case are summarized in Tables 123. Dimensionless data are considered according to the definitions indicated in Table 5.

A representation of a sequence of interfacial shapes is shown in Figure 4. It is seen that the interface forms a thin sheet with expanding holes that lead to the formation of ligaments and then droplets. The expanding holes are surrounded by a characteristic Taylor-Culick rim [START_REF] Taylor | The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets[END_REF][START_REF] Culick | Comments on a ruptured soap film[END_REF]. The potential and kinetic energies are mostly dependent on the large scales of the flow and are thus "macroscopic" quantities which, by definition, are related to the large scales and should for that reason be well resolved on moderately fine grids. A simple verification that this is indeed the case is achieved by plotting the total mechanical energy as done in Figure 5 for case 1. It is seen that the energy evolves in a nearly monotonic manner, and that it asymptotes to the expected value E f p at long times. The corresponding results for the kinetic and potential energies are reported in Figures 6 and7, respectively. There is still some difference between the 256 3 grid and the 512 3 grid for fluid 2 around time t * = 2.5, although not in fluid 1. We note that, since the kinetic energy is scaled by (1/16)ρ n U 2 g H 3 for fluid n (see Table 5), the maximum kinetic energy of fluid 2 in the scaled variables is 8(U K /U g ) 2 (ρ 1 /ρ 2 ) = 11.1375. This is close to the maximum of the curve for fluid 2 in Figure 6. The fact that the kinetic energy is close to its upper bound also means that U K is close to the L 2 -norm velocity, and that our Reynolds and Weber number estimates are consistent with the observed L 2 -norm velocity.

Parameter

Value Units The grid convergence results for the enstrophy are illustrated in Figure 8. Although the enstrophy seems to converge in fluid 1, its peak increases with the grid resolution in fluid 2 and no convergence is reached even on the 512 3 grid. Therefore, although primary moments such as potential and kinetic energies or volume ratio of the light fluid suggest that a real DNS was achieved, the analysis of higher-order moments such as enstrophy invalidates this hope. The interface shape and isovalues of vorticity magnitude at t * = 3 (which corresponds to the peak within fluid 2 in Figure 8) are reported in Figure 9. Finer vortical structures are clearly captured by the finer grid inside fluid 2 (near the bottom corners for instance). As pointed out in Table 3, the 512 3 grid does not reach the Kolmogorov scale η K and thus does not resolve the finest vortical structures. In addition, as we are considering two-phase flows with jumps in the physical properties, especially viscosity, properly capturing the interfacial vortical layers requires that a sufficient number of grid points be located around the interfaces. Therefore, we may suspect that the way the local viscosity is estimated as a function of the local volume fraction plays a role in the local vorticity magnitude. Indeed, an arithmetic average between µ 1 and µ 2 is generally used but there are good reasons to rather favor an harmonic average [START_REF] Ritz | A numerical continuous model for the hydrodynamics of fluid particle systems[END_REF][START_REF] Pianet | Simulating compressible gas bubbles with a smooth volume tracking 1-fluid method[END_REF][START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF]. To avoid this possible influence, we re-computed case 1, still with DyJeAT, while considering the same viscosity in both fluids. The corresponding results are discussed in Appendix B. They show that grid convergence of enstrophy is still not achieved when the viscosity is set to 0.1P a.s (the corresponding Reynolds number is 137) although the enstrophy excursions are less violent. A similar study was performed by some of us [START_REF] Sayadi | A convergence study of the one-fluid formulation in a phase inversion application at moderate Reynolds and Weber numbers[END_REF] with the conclusion that sheet breakup may be responsible for spurious large enstrophy even at reduced Re and La. The obvious conclusion of these additional computations is that the differences found among the various enstrophy evolutions are not due to the effects of an inaccurate numerical evaluation of the viscosity jump but to other effects requiring further study. We also note that in the atomization simulations performed by some of us [START_REF] Ling | A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability[END_REF], enstrophy as measured statistically is not diverging upon grid refinement. Hence the results do not correspond to a true Direct Numerical Simulation: convergence is achieved on quantities dominated by the large-scale motions, such as the kinetic and potential energies, but it is not on enstrophy for which the main contribution is from the small scales.

t * = t/t c t 0.643 - E f p,1 the potential energy in fluid 1 for t → ∞ 0.1035 J E f p,2 the potential energy in fluid 2 for t → ∞ 0.3755 J E * k,1 = E k,1 (1/16)ρ 1 U 2 g H 3 E k,1 0.000341 - E * k,2 = E k,2 (1/16)ρ 2 U 2 g H 3 E k,2 0.000378 - H 3 /

Droplet sizes

We now briefly analyze the evolution of the droplets size observed in results for case 1 obtained with DyJeAT. We saw in Figure 4 that the interface forms a thin sheet with expanding holes that lead to the formation of ligaments and then droplets. This hole formation mechanism was already observed in simulations of atomizing jets [START_REF] Shinjo | Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[END_REF] or mixing layers [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup[END_REF] and is critical in the atomization or fragmentation process [START_REF] Lohse | Double threshold behavior for breakup of liquid sheets[END_REF][START_REF] Villermaux | Fragmentation versus cohesion[END_REF]. It is likely that the control of the hole formation process is essential for obtaining a convergent droplet size distribution upon grid refinement [START_REF] Chirco | Manifold death: the implementation of controlled topological changes in thin sheets by the signature method[END_REF]. Figure 10 shows how the distribution of the largest drops observed with DyJeAT varies with the grid resolution and dimensionless time. At later times, no droplets are found. The number of droplets decreases in time as droplets reach the top of the box and merge with the bulk of fluid 1. Overall, at a given t * , the number of droplets also decreases as the resolution increases. Similarly, for a given resolution, the average size of the droplets decreases as time proceeds.

Case 2: a phase inversion problem involving many fragmentation events

In this section we present the numerical results pertaining to case 2. The characteristics of this case are also provided in Tables 123. Dimensionless data are considered according to the definitions indicated in Table 6. In order to examine grid convergence, the numerical simulations were run on five grids i.e. 128 3 , 256 3 , 512 3 , 1024 3 and 2048 3 . (It is worth noting that over 30 million hours of CPU time were necessary to perform the simulation on the finest 2048 3 grid.)

A plot of the total mechanical energy is shown in Figure 11. As in case 1, the energy decreases monotonically. It is seen that the convergence is irregular, as the 1024 3 simulation is further away from the reference 2048 3 simulation than the 512 3 one. This effect will be seen in more detail in the separate plots for the kinetic and potential energies shown below. At long times, the computed energy is clearly above the expected asymptotic value E f p . The results for the potential and kinetic energies are displayed in Figures 12 and13, respectively. The data obtained on the three finest grids superpose until t * = 4.5 approximately. We characterized the rate of convergence by comparing the potential energies at t * = 2.63. The corresponding conclusions are reported in Table 7. The order of convergence is near p = 2 for the low-resolution grids and near p = 1 for the highresolution ones. It is likely that at low resolution, the error on the flow characteristics away from interfaces dominates, while at high resolution the error near the interfaces is dominant. Since discontinuities in the pressure gradients and the derivatives of the velocity field are present near the interfaces, the discretization schemes lose an order of accuracy there and the decrease of the order of convergence to p = 1 is expected.

After that time, the comparison of the three finest grids shows that a well-characterized convergence no longer exists: the predicted kinetic energy varies irregularly with the grid size, the difference between the 2048 3 and 1024 3 grids being often larger than that between the 1024 3 and the 512 3 grids. This behavior starts at t * = 1.75 and becomes pronounced beyond t * = 5. This can only indicate one thing: despite the kinetic energy resulting mostly from large-scale motions, the small scales influence strongly the evolution of the kinetic energy. These small scales may for instance be involved in coalescence events, which themselves display the perforation of thin liquid sheets. Whether the sheets perforate and are destroyed or not has a strong influence on the large scales. Because of the observed convergence behavior of DyJeAT, its 2048 kinetic energy can be taken as a reference, with certainty until t * = 1.75, and as a reasonable estimate until t * = 5.

E * k,1 = E k,1 (1/16)ρ 1 U 2 g H 3 E k,1 3.41 - E * k,2 = E k,2 (1/16)ρ 2 U 2 g H 3 E k,2 3.78 -
An oscillating behavior is still present in the kinetic energy evolution, even though it is less regular than in case 1. In all simulations, the kinetic energy decays as t * -2 in fluid 2, as already found with case 1. In fluid 1, the exponential Stokes decay law (see Appendix A) is clearly more difficult to obtain, as large amplitude variations are observed. Compared to case 1, the turbulent intensity in fluid 2 is larger, especially in the bottom part of the box, and thus provides a strong forcing to the large blob of light fluid that stands on top of it. This is why the oscillatory motion of fluid 1 is more complex and its decay does not strictly follow the purely viscous Stokes law.

The time histories of enstrophy are plotted in Figure 14 enstrophy in fluid 2 (compared to fluid 1) was not observed in case 1. This newly found delay is a consequence of the higher Reynolds number: vortical layers develop quite quickly in the more viscous fluid 1, but a significantly longer time is required for shear regions to develop in fluid 2, owing to the prevalence of inertial effects.

Again, the enstrophy magnitude increases with the grid resolution. Compared to case 1, this tendency is reinforced by the large number of break-up events which result in a large population of droplets of fluid 1 that modulates the motion of fluid 2. Whatever the grid resolution, including 2048 3 , convergence is not achieved and the finer the grid, the larger the enstrophy magnitude. As Figure 14 shows, the difference between the peak magnitudes obtained with two successive resolutions increases as the grid is further refined. This is an indication that the highest resolution considered here is still far from that required to achieve a true DNS. However, it is observed that the width of the time interval in which the enstrophy differs is decreasing with the refinement of the grid, testifying that grid convergence may be reached on a finer mesh.

A view of the interfacial structure is provided in Figure 15. Many lenticular (or "flattened" ) droplets are seen, an effect of the large rising speed that results in inertial effects frequently stronger than surface tension effects. This relative weakness of capillary effects occurs at fairly large scales, with a threshold that seems to be close to 0.05H. This length is close to the Hinze scale η H = 0.03H given in Table 3, which was to be expected. Smaller droplets are seen to remain nearly spherical. The interface rendered in Figure 15 is also seen to be much smoother than in Figure 2, indicating that most of the droplets are well-resolved. We will return to the issue of the resolution of the smallest scales in Section 4.4.

Figure 15: Topology of interfaces for case 2 at t * = 3.5 (near the enstrophy peak) in fluid 2 using DyJeAT on a 1024 3 grid. The interface is much smoother than in Figure 2, indicating that most of the droplets are well-resolved.

Droplet size distributions in case 2

We performed a study of the droplet size distributions in case 2 at various times. The corresponding results are summarized in Figure 16.

The droplet volumes V are measured in the simulations and a corresponding equivalent radius is obtained for each drop as r e = (3V /(4π)) 1/3 . The equivalent radius may belong in any of a number of adjacent bins B i defined as the interval B i = (r i -∆ i /2, r i + ∆ i /2) where r i and ∆ i are the bin center and width, respectively. We consider regularly spaced bins with ∆ i = 1/256, r i = (i -1/2)/256 and B i = ((i -1)/256, i/256). Then the number of droplets whose equivalent radius lies in B i is noted N i . N i is approximately proportional to an asymptotic probability f such that f (r i ) = N i /∆ i ; however this approximation is inappropriate for the first bin because ∆ 1 /r 1 is far from small. The function f (r i ) defines the Probability Distribution Function (PDF) of the droplet sizes. In Figures 16 and17, we plot N i versus r i at various times or various resolutions.

The observation of these distributions leads to several remarks. First, as time progresses (see Figure 16), the number of droplets increases, but the number of large drops increases faster than the number of small droplets. Beyond a radius of 0.02 m, the distribution roughly follows an exponential law with an average droplet size r that increases in time.

In Figure 17, a convergence study of the droplet size distribution is performed at time t * = 4. First, it is seen that the 256 3 grid result does not match the results obtained on finer grids. Instead, droplet populations obtained with that grid are systematically less numerous than with the other grids, indicating that the total mass of droplets obtained with that grid is smaller. This can be easily explained by the fact that poorly resolved droplets "evaporate" in a Level Set method such as that used in DyJeAT. Second, the 512 3 , 1024 3 and 2048 3 grid results almost superpose in a region of variable size, going from a minimum radius r m (N ) to a maximum radius r M (N ), with r m and r M both depending on the grid resolution N . To be specific, the minimum radius of agreement decreases with increasing N until approximately r m (2048) 0.005 m which indicates a range of convergence extending to very small droplet radii for the 2048 3 grid. The upper boundary r M (N ) of the range of convergence is harder to determine, because large radii are affected by statistical noise. However, it is seen that few additional very large droplets (r > 0.07 m) are identified as the resolution increases from 256 3 to 2048 3 . Moreover, grids 1024 3 and 2048 3 are in agreement over a wide range with r M (2048) 0.05 m. Finally, we note that the behavior at very small radii differs between the 2048 3 grid and the other grids: while the latter have a cutoff (a fall in droplet counts N 1 and/or N 2 ) at small sizes, the 2048 3 grid does not. This is due to the fact that in the 2048 3 case, the fall in droplet count can only be seen if instead of the peculiar bin B 1 defined above, one subdivides the interval (0, 1/256) into several smaller bins. We have actually done that and did observe the cutoff even for the 2048 3 case. Similar effects of the grid refinement, including cutoffs or drops in frequency, have been reported for atomization of liquid jets, especially with Level Set methods [START_REF] Herrmann | On simulating primary atomization using the refined level set grid method[END_REF] and less markedly for VOF approaches [20]. This indicates that 512 3 grids are a minimum requirement to reach an approximately accurate droplet count at any scale, and that 256 3 simulations are poorly resolved in the whole range of scales. On the other hand, the absence of convergence of the enstrophy does not imply a non-convergence of the droplet sizes in the intermediate range r m (N ) < r < r M (N ). Another interesting fact is that the large scales are influenced by how the small scales are computed: more large droplets are seen with finer grids, probably because a higher grid resolution implies less breakup of thin layers or filaments, thus reducing the rate at which large structures break into smaller ones. In order to better understand the PDF, we plot the computed frequencies for the two most refined grids as a graph of ln(r i N i ) versus ln r i in Figure 18. The first bin B 1 was excluded due to its particular character, and only results for the most refined grids 1024 3 and 2048 3 are shown. Using these coordinates has the advantage that the log-normal PDF appears as a parabola. This specific PDF is defined as

f (r) = A r exp - (ln r -ln μ) 2 2σ 2 , ( 10
)
where A is a normalization constant, ln μ is the logarithmic average and σ is the logarithmic standard deviation. A parabolic curve corresponding to such a log-normal distribution is shown on Figure 18 to aid the interpretation. This parabola is built with μ = 10 -3 m, σ = 1.7 and A = 3. While the value of μ is very uncertain, it is clear that extrapolating the trend of the distribution leads to expect a larger value of rf (r) at smaller r. These larger values of rf (r) could be revealed in simulations with yet more refined grids than the current ones. This strongly reinforces the expectation that many smaller droplets would be seen in a true DNS. However, these droplets would contribute little to the interfacial area. Indeed the interfacial area scales as

Σ S = ∞ 0 4πr 2 f (r) dr , (11)
which for a log-normal f (r) or a simpler f (r) ∼ 1/r converges at the r = 0 bound. Note that we took the assumption of spherical droplets to estimate Σ S . The latter assumption is likely to be a good approximation for the smallest droplets. Thus, the absence of such droplets in the integral does not affect much the interfacial area.

A very simple Pareto distribution of the form f (r) ∼ A/r would be a horizontal line in the variables of Figure 18. That approximation is in fact a good fit in the intermediate range of droplet sizes with A 1, and is reminiscent of distributions seen in other contexts [START_REF] Balachandar | Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines[END_REF][START_REF] Pairetti | On shear layer atomization within closed channels: Numerical simulations of a cough-replicating experiment[END_REF]. As a third alternative, an exponential distribution f (r) = B exp(-r/r m ) is also plotted on Figure 18 with r m = 0.02 m. The exponential distribution provides a better fit for the largest r. With this distribution, the interfacial area integral converges for ∆x r m = 0.02 m and thus for 256 3 grids.

Finally, we note that none of the distributions yields a very good fit throughout the entire range of scales. It is possible that the PDF is in fact bimodal, with different "modes" corresponding to different physical mechanisms being superimposed at various scales.

Results of the benchmark using all codes at moderate resolution

The teams involved in the present benchmark ran their own code each. Throughout, 256 3 grids were used.

Case 1 results for all codes

We first compare the performance of all codes in terms of mass conservation in Figure 19. DyJeAT is also run at 256 3 resolution for this purpose. In the case of DyJeAT, the issue of mass conservation is quite specific, since the mass decrease or increase is computed at each time step before the level-set function is shifted to bring the mass deviation back to zero. In doing so, mass may be transferred in a non-local manner. Thus the interface is rearranged locally, making for instance small droplets disappear. What is plotted in Figure 19 in the case of DyJeAT is the cumulative absolute value of the deviation before the level-set shift. All codes have fair mass conservation properties with the worse volume error always smaller than 1.6% of the volume of light fluid. This is a nevertheless quite significant error that should be reflected in energy conservation, since conservation of the total mechanical energy is closely related to volume and mass conservation. Here, since almost all of fluid 1 goes to the top, the final potential energy of the system, equal to the final total mechanical energy minus the surface energy, is determined to a very good approximation by the final mass of fluid 1 in the top region. Thus, any error in the final volume (hence the mass) of fluid 1 in the top region should result in an error on the final mechanical energy.

Figure 20 shows that all codes predict a final mechanical energy close to the reference value. More precisely, normalizing differences by the total energy variation from t = 0 to t → ∞, it is found that the departure from the reference value is less than 1.4% for all codes, which is consistent with the error found on mass conservation. More specifically, it is seen in Figure 20 that Archer and Jadim agree with the reference within the thickness of the line, while Thetis and Gerris perform somewhat worse, with Thetis even below the theoretical asymptotic value. This ranking of the codes may be partially inferred from Figure 19 which confirms the good mass conservation behavior of Archer and Jadim but does not show a significantly worse behavior for Thetis. Indeed there are other aspects that can affect the final potential energy: droplets of light fluid that remain attached to the bottom wall by capillary forces increase the final potential energy, as do droplets of heavy fluid attached to the top. This effect is counterbalanced by mass conservation: an excess (resp. deficit) of light fluid at the end of the simulation, near the top, decreases (resp. increases) potential energy. It is seen that the mass of light fluid increases with Thetis, which agrees with the undershoot seen in the corresponding energy. Finally, DyJeAT with a 256 3 resolution is the only code that is slightly above the theoretical asymptotic value. This is probably related to the specific treatment of mass conservation described above. for the two versions of Thetis, and to a lesser degree the Jadim code. The behavior of the kinetic energy is more complex as differences between the behavior in fluids 1 and 2 are noticed. However, at large times the kinetic energy predicted by the Archer code for both fluids is the closest to the reference solution. This may be due to the non-converged character of the DyJeAT reference solution for t * > 2.5, which would make the agreement with Archer easier since both codes are at least partially based on the Level Set approach. In other words, Archer being the code most similar to DyJeAT is the most likely to agree with the reference. Figure 21 indicates that phase separation is achieved after a dimensionless time t * ≈ 15. The dimensionless frequency of the waves observed at the surface of the light fluid is about 0.6. Again, it is observed that all codes predict evolutions close to each other and to the reference solution.

To better view the various predictions of kinetic and potential energies by the different codes in another perspective, we plot in Figure 23 the deviation of the energies of all codes with respect to the DyJeAT reference. This deviation is shown for the kinetic energy, the potential energy, and the sum of these two components (total mechanical energy minus the surface energy). Interestingly, for some of the codes (especially Thetis) the differences are oscillatory while for others there is a less coherent drift. Note also that the energy deviations shown in Figure 23 are of O(10 -3 ) and thus much smaller than the total potential energy which is of O(0.1) (This value of the potential energy in Joules can be deduced from the characteristic energies used for the normalization defined in Table 5.) The code-to-code comparison provided by Figures 23 indicates some phase shifts and various damping rates at large times. In this limit, Thetis is the furthest from the DyJeAT reference, followed by Gerris, while Jadim is very close to that reference. Archer has fluctuations that make the energy at large times slightly unsteady. The codes exhibit the largest difference with the reference at time t * = 2.5, close to the kinetic energy peak (see Figure 22). These differences, while all of the same order of magnitude are clearly the largest for Thetis followed by Archer, then Jadim and Gerris, and are the smallest for DyJeAT 256 3 . A possible explanation for these code-to-code differences could be the quality of mass conservation in each of them, since it directly influences the amount of fluid in the top part. Another explanation would be the accuracy of energy conservation since it drives the observed oscillations. A third explanation is the accuracy of the computation of the surface tension force, since it can remove or add energy proportional to (i) the interfacial area and (ii) the "spurious currents" generated by small imbalances with the other terms in the momentum equation (see [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]).

Figures 24 and25 display the instantaneous values of the enstrophy in both fluids. Although all codes were shown to agree well on other macroscopic quantities, the magnitude of the enstrophy peak differs dramatically for each of them. The location of this peak is found to occur at t * ≈ 2 (resp. 3) for all codes in fluid 1 (resp. fluid 2). However, differences in magnitude up to 80% are observed, depending on the grid and the code. Despite the large discrepancies noticed between all codes, it is interesting to observe that they all provide a t * -3 temporal decay law which may be shown to correspond to a turbulent scaling

Case 2 results for all codes

We first compare all codes in terms of mass conservation in Figure 26. DyJeAT is also run at 256 3 resolution for this purpose. The same remark as above applies to the technique involved for mass conservation in DyJeAT. The codes have very dispersed mass conservation properties. As above, conservation of the total mechanical energy approximately deduced from the analysis of Figure 27 indicates how each code conserves volume, since the final potential energy is related to the final volume of fluid 1 in the top region. It is seen in Figure 27 that Archer and Gerris agree with the reference within the thickness of the line, while Thetis and Jadim perform worse, with Thetis again below the theoretical asymptotic value. Unlike case 1, this ranking of the codes can also be deduced from the mass conservation error plotted in Figure 26, where Gerris and Archer perform much better than the other codes. In this case, both Gerris and DyJeAT 256 3 end up above the DyJeAT reference value at late times. The reason why mass conservation deteriorates in Jadim beyond t * = 5 is that the local mass error control strategy described in Section 3.3 was not designed to handle topological changes, and presumably fails during coalescence and breakup events. Such events being much more frequent in case 2 than in case 1, departures from exact mass conservation are much larger in the former case.

Similar to the observations reported for case 1, the results provided by the five codes are in approximate agreement as far as the time evolution of the potential and kinetic energies is concerned, even if larger discrepancies are observed compared to case 1. These evolutions are displayed in Figures 28 and29. The behaviors of the codes are diverse. For example, the Gerris code performs well for the kinetic energy of fluid 2. However, it predicts a significantly different oscillatory dynamics for the kinetic energy of fluid 1 beyond t * = 1.5. In fact, for fluid 1, all codes are close to the reference, except Gerris, until time t * = 2.5.

Then they all diverge from the reference. For fluid 2, Gerris and Thetis behave better than the other codes within the time interval 2 < t * < 3.

Similar to case 1, to view the different predictions of kinetic and potential energies by the different codes in another perspective, we plot the deviation of the energies predicted by all codes with respect to the DyJeAT reference in Figure 30. This deviation is shown for the kinetic energy, the potential energy, and the sum of these two components (total mechanical energy minus the surface energy). For some of the codes (especially Gerris), the sum of the kinetic and potential energies remains close to the reference at early times, with both components deviating in an oscillatory manner and in phase opposition. In the case of Gerris, the plot also includes the deviation of the surface energy from the "reference" DyJeAT case. This deviation mainly serves to show that the surface energy fluctuations are of the same order as the deviations in kinetic and potential energies. It does not seem that the performance of the codes in terms of prediction of the kinetic and potential energies is related to their performance in the accounting for surface tension forces, since, at least in the case of Gerris, the deviations in these various energies are not correlated.

For some other codes, the deviation from the reference is much less oscillatory, indicating that the error is dispersed in a non-oscillatory manner. Finally, in some cases the total mechanical energy (minus the surface energy) is not conserved, indicating a possible dissipative error or a mass conservation error. In the long time limit, Thetis and to a lesser extent Jadim have energies that deviate from the reference, while Gerris and Archer are very close to it. This behavior is consistent with the mass conservation of these codes. DyJeAT 256 3 is close to the reference (obtained with DyJeAT itself). Unlike case 1, there is no clear maximum of the error around t * = 2 or 3 for all the codes, although maxima are found around t * = 2 for some of them and either the total, potential or kinetic energies. These "maximum errors" are of the same order of magnitude, between 2.5 and 5 Joules. Note also that these energy deviations need to be compared with the total kinetic energy which is of the order of 40 Joules (see Table 6 and Figure 29). The time histories of enstrophy are plotted in Figures 31, 32 and 14. Again, all codes find the peak value nearly at the same time, namely t * ≈ 2.5 in fluid 1 and t * ≈ 4 in fluid 2, although the code-to-code differences are significantly larger than in case 1. Again, the enstrophy magnitude increases with grid resolution and is code-dependent.

The uncertainty surrounding the 256 3 grid results for the PDF of droplets sizes leads us not to discuss the results provided by the other codes for this PDF, since these codes were not run at higher resolution. When droplets are counted in these low-resolution simulations, an exceedingly large number of small droplets resulting from numerically dominated breakup, -so called debris, flotsam or jetsam or wisp droplets -is found, precluding a meaningful analysis in the small droplet range. This type of debris droplet is absent or rare in pure Level Set simulations.

Summary and future work

A computational benchmark based on the phase inversion of two immiscible fluids of different densities confined in a closed cubic box has been set up. The purpose of this benchmark was to check the capabilities and limitations of current codes and grid resolutions. The first observation made with the reference code, DyJeAT, is that the total mechanical energy converges easily with grid resolution and can be rescaled so that it exhibits a very similar behavior in both cases. This effect was unexpected and it would be interesting to see if it applies to a wider range of cases. The second observation, put simply, is that the various codes, all run at much lower resolution than the reference, reproduce easily the results of the reference for quantities tied to the large scales, such as kinetic and potential energies, but fail to reproduce quantities related to the small scales, such as the PDF of droplet sizes or the enstrophy integral.

The first problem we considered, case 1, is characterized by an integral scale Reynolds number Re = 1.37 10 4 . Computational results have been analyzed by considering the time evolution of several volumeaveraged indicators, namely the potential and kinetic energies in each fluid, the relative volume of light fluid in the top part of the box and the enstrophy in each fluid. The results reveal that all codes provide close evolutions for the first three quantities, indicating that reliable predictions are obtained when quantities essentially determined by large-scale motions are considered. The long-term decay rates of the kinetic energies have been found to be in good agreement with the viscous Stokes law (in fluid 1) and the turbulent decay law (in fluid 2) discussed in Appendix A, respectively. In addition, all codes agree on the sloshing frequency of the light fluid. The situation has been found to be much more problematic when enstrophy is considered: all codes provide markedly different values of the enstrophy maximum (although they agree on the time at which this maximum occurs) and none of them converges when the grid is refined.

A second phase inversion test case, involving a higher Reynolds number, Re= 4.33 10 5 , was subsequently considered. The same general conclusions apply to this configuration. All codes provide converging predictions as long as only kinetic and potential energies are considered, whereas large discrepancies are observed on enstrophy. Here again, the true DNS conditions are not satisfied and much finer grids should be considered to expect convergence on quantities such as enstrophy and dissipation which are governed by small-scale processes, especially multiple break-up events. It has to be noticed that despite the discrepancies observed on enstrophy, approximate convergence is obtained for PDFs of droplet sizes as soon as a 1024 3 grid is used. Clearly, case 2 is an implicit LES rather than a DNS. As a conclusion, this benchmark illustrates the fact that implicit LES simulations on fine grids allow to reach convergence on first-order moments such as kinetic or potential energies while they fail to converge as soon as small-scale statistics are investigated. Explicit LES models may be implemented in this case [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF][START_REF] Larocque | Parametric study of LES subgrid terms in a turbulent phase separation flow[END_REF][START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF][START_REF] Aniszewski | A new approach to sub-grid surface tension for LES of two-phase flows[END_REF][START_REF] Vincent | A priori filtering and LES modeling of turbulent two-phase flows -application to phase separation[END_REF][START_REF] Saeedipour | Large eddy simulation of turbulent interfacial flows using approximate deconvolution model[END_REF], or much finer meshes should be employed to achieve a true DNS.

Future work of some of the partners will first be devoted to the consideration of cases with smaller Re and We in order to perform true DNS with convergence of the small scales. Multiscale Eulerian VOF or Level Set approaches coupled to a Lagrangian description of small droplets may also be introduced by some of them with the goal of reaching convergence on all physical quantities on "reasonable" grids. Last, adaptive mesh refinement (AMR) techniques will be considered to concentrate numerical efforts in flow regions where the vorticity magnitude is expected to approach its maxima.

Acknowledgements

This work was part of the MODEMI ANR project (Modélisation et Simulation Multi-échelle des Interfaces, ANR-11-MONU-0011) devoted to the multiscale modeling of two-phase flows. We thus acknowledge the support of the Agence Nationale de la Recherche (ANR) for the I2M, D'ALEMBERT, CORIA and IMFT laboratories concerning the Multiscale Modeling of Interfaces. The authors wish to thank the Midi-Pyrénées and Aquitaine Regional Councils for the financial support dedicated to a PhD thesis at ONERA and IMFT and a 256-processor cluster investment, located in the TREFLE team of the I2M laboratory. The D'ALEMBERT team benefited from the ERC grant TRUFLOW. We are grateful for access to the computational facilities of the French CINES (National computing center for higher education) and TGCC, notably Irene-Rome, granted by GENCI under project numbers x20142b6115, x20152b6115, x20162b6115, A0012b06115 and A0092B07760. The project also benefited from the PRACE grant TRUFLOW number 2020225418 for a large amount of CPU hours on TGCC Irene-Rome in 2021. This work was also partly granted access to the HPC resources of CALMIP under the allocation 2013-P0633. The CORIA team would like to express its gratitude to the CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, www.criann.fr) for providing CPU resources. We thank the technical and administrative teams of these supercomputer centers and agencies for their kind and efficient help.

A Appendix A: Scaling laws for the kinetic energy decay

After the acceleration stage induced by buoyancy forces, during which most of the lighter fluid goes to the top part of the box, the phase inversion problem is characterized in a second stage by a wavy behavior which makes it look quite similar to a sloshing flow progressively damped by shearing and viscous effects.

The classical analysis of Stokes regarding the viscous damping of gravity waves (see e.g. [START_REF] Lamb | Hydrodynamics[END_REF], pp. 623-624), makes it possible to predict the time evolution of the kinetic energy in a weakly viscous flow driven by a surface wave. First, the time evolution of the mechanical energy, which is the sum of the potential and kinetic energies, is known to result from the internal viscous dissipation, so that

dE m dt = - 1 2 µ Ω ∇u + (∇u) T 2 dV . (12)
If we assume the shear layer at the free surface to have negligible effect owing to the moderate velocity gradients expected in this region, the flow can be considered irrotational. Then the velocity potential, Φ, and rate of change of the mechanical energy, E m , read Φ = Φ 0 e kz cos(kx -ωt) ,

dE m dt = -8µk 4 Ω Φ 2 dV , (13)
where k and ω are the wave number and radian frequency, respectively, and the overbar denotes the time average value. In the framework of linear wave theory, the potential and kinetic energies are equal. However the wave amplitude is not small in the present phase inversion problem. Nevertheless, as soon as most light fluid stands in the top part of the box, the potential energy stays almost constant, making the time variation of E m mainly governed by that of the kinetic energy. Therefore, still in the linear approximation, we can approximately write

E m = 1 2 ρ Ω u 2 dV ≈ 1 2 ρ Ω ∇Φ • ∇Φ dV = ρk 2 Ω Φ 2 dV . (15)
Now, two markedly different flow situations can be met in the phase inversion problem, namely a "gentle" configuration in which break-up and coalescence events of the light fluid scarcely occur and a "violent" configuration in which such events are much more numerous. Combining ( 14) and ( 15), the time evolution of the kinetic energy in the first regime is found to obey

E m = 1 2 ρ Ω u 2 dV = Ke -8νk 2 t = Ke -8ν ω 4 g 2 t , (16)
where K is a constant, ν = µ/ρ, and use has been made of the dispersion relation ω 2 = gk.

Let us now consider the "violent" configuration. In this case, averaging throughout the whole volume Ω, ( 13)-( 14) simply yields

u 2 = ∇Φ • ∇Φ = k 2 Φ 2 , ( 17
)
dE m dt = -8µk 2 u 2 , (18)
where • denotes the mean value over Ω. In other words, one has

d u 2 dt = -16νk 2 u 2 . ( 19
)
In this configuration, the flow is expected to be turbulent and the molecular viscosity is no longer relevant to estimate the damping rate of the flow. As a crude surrogate, we can introduce an effective turbulent viscosity ν t scaling as l u, where l is a characteristic length of the large-scale flow, which can be taken as the box size. Replacing ν by ν t in [START_REF] Vincent | Direct numerical simulation of phase separation and a priori two-phase LES filtering[END_REF], we then obtain

d u 2 dt = -16lk 2 u 3 . ( 20
)
Assuming the volume-averaged velocity to follow a power law, i.e. u to be of the form u 0 t -n , yields n = 1, from which we infer

E m ∼ t -2 . ( 21
)
This prediction is to be compared with the exponential decay predicted by [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] in the "gentle" regime. The result ( 21) is reminiscent of the approximate decay law of the kinetic energy in decaying homogeneous isotropic turbulence (e.g. [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]). This is expected since, according to Taylor's estimate, the dissipation rate is known to scale as u 3 0 /l 0 , where u 0 and l 0 stand for the large-scale velocity and length scales, respectively. Therefore, when l 0 is constant, ( 21) is immediately recovered and the dissipation rate is predicted to decay as t -3 .

B Appendix B: influence of viscosity on the flow dynamics

In order to better understand the origin of code-and grid-dependencies of the volume-averaged enstrophy in cases 1 and 2, we performed an extra series of computations with DyJeAT with the same viscosity in both fluids. In this way, any possible influence of the averaging procedure selected to compute the local viscosity as a function of the volume fraction and of the numerical treatment of the viscosity jump in the interfacial grid cells is removed. We make use of the physical parameters of case 1, except for viscosity which is set to 0.1P a.s (case 1a), 0.01P a.s (case 1b) and 0.001P a.s (case 1c), respectively. The convergence study was performed on three different grids, namely 128 3 , 256 3 and 512 3 in case 1a, whereas four grids ranging from 128 3 to 1024 3 were considered in cases 1b and 1c. The evolution of the volume-averaged enstrophy is reported in Figures 33,34 and 35 for cases 1a, 1b and 1c, respectively. In case 1a, the Reynolds number is 130 and the enstrophy while still not converging as the grid is further refined, has less violent excursions. In this configuration, the enstrophy decay at large time follows an exponential law, a behavior typical of the viscous decay of gravity waves (see Appendix A). In case 1b, grid convergence is again not achieved in both fluids, the Reynolds number in fluid 2 being 1300. At large time, the enstrophy obeys a t * -3 decay law typical of turbulent conditions (see Appendix A). In case 1c, enstrophy convergence is not achieved in either fluid. The Reynolds number is 1.37 10 4 and it is observed that, even though the results concerning energy and relative volume in the upper part of the box converge, as shown with case 1, the grid is not thin enough to capture the small-scale eddies. Again, the decay law at large times is found to correspond to turbulent conditions.

Figure 1 :

 1 Figure 1: Initial setup of the phase inversion problem in a closed box. The light liquid or oil is fluid 1 and the more dense, less viscous liquid is fluid 2.

  ρ

Figure 3 :

 3 Figure 3: Evolution of the rescaled (see text) total mechanical energy for both cases using DyJeAT.

Figure 4 :

 4 Figure 4: A sequence of interface shapes for case 1, using DyJeAT at 512 3 resolution. The dimensionless times are (from left to right) t * = 2, t * = 2.25 and t * = 2.5. One can see a large sheet pierced by two expanding holes (arrows).

Figure 5 :

 5 Figure 5: Evolution of the dimensional total mechanical energy (without the surface energy), E m , as a function of dimensionless time for case 1. The horizontal line is the expected value E f p .

8 2 Table 5 :

 825 the final volume of fluid 1 in the top part of the box 0.000125 m 3 Maximum of enstrophy in fluid 1 (DyJeAT on a 512 3 grid) 0.0733 m 3 .s -2 Maximum of enstrophy in fluid 2 (DyJeAT on a 512 3 grid) 1.3759 m 3 .s -Parameters used to define the dimensionless variables in case 1.

Figure 6 :

 6 Figure 6: Grid convergence of kinetic energies for case 1.

Figure 7 :

 7 Figure 7: Grid convergence of potential energies for case 1.

Figure 8 :

 8 Figure 8: Grid convergence of enstrophy for case 1.

Figure 9 :

 9 Figure 9: Vorticity structures and interface shape for case 1: iso-60 vorticity magnitude (blue) and iso C = 0.5 interface (orange surface) obtained with DyJeAT on the 128 3 grid (left) and 256 3 grid (right).

Figure 10 :

 10 Figure 10: A representation of the sizes of the largest droplets (assuming H = 0.1 m) found in case 1, as a function of time and grid resolution.

H 3 / 8 3

 383 the final volume of fluid 1 in the top part of the box 0.125 m Maximum of enstrophy in fluid 1 (DyJeAT on a 2048 3 grid) 103 m 3 .s -2 Maximum of enstrophy in fluid 2 (DyJeAT on a 2048 3 grid) 1049 m 3 .s -2

Figure 12 :

 12 Figure 12: Grid convergence of potential energies for case 2.

Figure 13 :

 13 Figure 13: Grid convergence on kinetic energies for case 2 in log-linear (left) and log-log (right) coordinates.

Figure 14 :

 14 Figure 14: Convergence of enstrophy for case 2.

Figure 17 :

 17 Figure 17: Droplet radius frequency N i obtained with DyJeAT at t * = 4. The droplet radii are in meters, assuming the box size to be H = 1 m.

Figure 18 :

 18 Figure 18: Droplet radius frequencies (PDF) in logarithmic coordinates, with rf (r) versus r. A log-normal distribution and an exponential one are shown to guide the eye.

Figure 19 :

 19 Figure 19: Mass conservation for case 1. The dimensional volume deviation of fluid 1 (light fluid) is plotted as a function of the dimensionless time.

Figure 20 :Figure 21 :

 2021 Figure 20: Evolution of the dimensional total mechanical energy in case 1 (without the surface energy), E m , for all codes as a function of the dimensionless time. The horizontal line is the expected value E f p .

Figure 22 :Figure 23 :

 2223 Figure 22: Kinetic energies for case 1 in log-linear (left) and log-log (right) coordinates.

Figure 24 :

 24 Figure 24: Enstrophy for case 1.

Figure 25 :

 25 Figure 25: Enstrophy for case 1 in log-log coordinates.

Figure 26 :Figure 27 :

 2627 Figure 26: Mass conservation for case 2.

Figure 28 :

 28 Figure 28: Potential energies for case 2.

Figure 29 :

 29 Figure 29: Kinetic energies for fluid 1 in log-linear coordinates (left) and fluid 2 in log-log coordinates (right) (case 2).

Figure 30 :

 30 Figure 30: Deviation of the energies from the reference in case 2. Results shown for Gerris also include the deviation of the surface energy from the "reference" DyJeAT case. DyJeAT 256 3 is compared to the DyJeAT 2048 3 reference.

Figure 31 :

 31 Figure 31: Enstrophy for case 2.

Figure 32 :

 32 Figure 32: Enstrophy for case 2 in log-log coordinates.

Figure 33 :

 33 Figure 33: Grid convergence of enstrophy for case 1a in fluid 1 (left) and fluid 2 (right) -linear-linear (top) and log-linear (bottom) coordinates.

Figure 34 :

 34 Figure 34: Grid convergence of enstrophy for case 1b in fluid 1 (left) and fluid 2 (right) -linear-linear (top) and log-log (bottom) coordinates.

Figure 35 :

 35 Figure 35: Grid convergence of enstrophy for case 1c in fluid 1 (left) and fluid 2 (right) -linear-linear (top) and log-log (bottom) coordinates.

  1 gH 4 /32 and E i p,2 = 15ρ 2 gH 4 /32 while for t → ∞, E f p,1 = 15ρ 1 gH 4 /128 and E f p,2 = 49ρ 2 gH 4 /128. An upper bound on the kinetic energy is

Table 1 :

 1 Fluid and physical properties common to both cases.

	1	ρ 2	µ 1	µ 2	g
	(kg.m -3 ) (kg.m -3 ) (P a.s) (P a.s) (m.s -2 )
	900	1000	0.1	0.001	9.81
		Case	σ	H
			(kg.s -2 ) (m)
		1	0.045	0.1
		2	0.45	1

Table 2 :

 2 Fluid and physical properties that depend on the considered case.

Table 3 :

 3 Dimensionless numbers for the considered cases. η K is the Kolmogorov scale at which inertia balances viscous stresses in a cascade, η H is the Hinze scale at which inertia balances capillary effects in a cascade, and d T is the diameter of a Taylor-Culick rim.

	Case	Re	We La=Re 2 /We	η K /H	η H /H d T /H
	1	1.37 10 4 41.6	4.5 10 6	7.9 10 -4 0.107	0.02
	2	4.33 10 5 416	4.5 10 8	5.9 10 -5	0.03	0.002

Table 4 :

 4 Summary of the numerical methods used in the various codes.

	Code	Interface	Navier-Stokes	Capillary	Grids	Solvers
		tracking	velocity-pressure	forces
			coupling
	Archer	CLSVOF	Projection	Ghost	Cartesian BiCGStab
			method	fluid	staggered	Multigrid
	DyJeAT	Level set	Projection	Ghost	Cartesian BiCGStab
			method	fluid	staggered	Multigrid
	Gerris	Geometrical	Projection	CSF with	Collocated BiCGStab
		VOF	method	height functions and AMR	Multigrid
	Jadim	Implicit	Projection	CSF	Cartesian BiCGStab
		VOF	method		staggered	Multigrid
	Thetis	Geometrical	Augmented	CSF with	Cartesian BiCGStab
		VOF	Lagrangian	smooth VOF	Staggered	ILU

  Figure 11: Evolution of the dimensional total mechanical energy (without the surface energy), E m , as a function of dimensionless time. The horizontal line is the expected value E f p .

		4880
		4860			DyJeat -2048
					DyJeat -1024
					DyJeat -512
					DyJeat -256
	Total energy	4820 4840			DyJeat -128 E pot, inf
		4800
		4780
		0	2	4	6	8	10
				Dimensionless time
		Parameter				Value	Units
	t * = t/t c p,1 the potential energy in fluid 1 for t → ∞ E f	t 2.03 1035	-J
	E f p,2 the potential energy in fluid 2 for t → ∞	3755	J
							3

Table 6 :

 6 Parameters used to define the dimensionless variables in case 2.

Table 7 :

 7 . The delay observed in the development of Case 2. Order of convergence p of the dimensionless potential energy of fluid 1 at t

	n	cells	|E * p (n + 1) -E * p (n)|	p
	1 128 3	0.0281	1.6674
	2 256 3	0.0088	2.4120
	3 512 3	0.0017	0.9675
	4 1024 3	0.0008	-
	5 2048 3	-	-

* = 2.63.

  Figure 16: Droplet radius frequency N i at various dimensionless times obtained with DyJeAT at maximum resolution. The droplet radii are in meters, assuming the box size to be H = 1 m.

						t*=1.84 t*=1.98 t*=2.48 t*=3.21 t*=4
		100
	Droplet count N(r)	10
		1
		0	0.02	0.04	0.06	0.08	0.1
				Droplet radius r (meters)
						256 grid 512 grid 1024 grid 2048 grid
		100
	Droplet count N(r)	10
		1
		0	0.02	0.04	0.06	0.08	0.1
				Droplet radius r (meters)