Aude Bernard-Champmartin

Florian De Vuyst
email: devuyst@cmla.ens-cachan.fr

Florian De

A low diffusive Lagrange-remap scheme for the simulation of violent air-water free-surface flows

Keywords: numerical method, multiphase flow, air-water flow, free boundary, interface capturing, compressible fluid, finite volume, Lagrange-remap solver, solver for parallel computing, wave breaking, sloshing, impact problem, wave formation and wave breaking;, wall wave impacts, local pressure peaks and pressure loadings;

] in the context of compressible multifluid flows using a five-equation model. In this paper, we explore a simplified variant approach for gas-liquid three-equation models. The Eulerian numerical scheme has two ingredients: a robust remapped Lagrange solver for the solution of the volume-averaged equations, and a low diffusive compressive scheme for the advection of the gas mass fraction. Numerical experiments show the performance of the computational approach on various flow reference problems: dam break, sloshing of a tank filled with water, water-water impact and finally a case of Rayleigh-Taylor instability. One of the advantage of the present interface capturing solver is its natural implementation on parallel processors or computers.

Introduction

Simulation of free surface flows knows an increasing interest as an essential predictive tool for innovative Engineering designs into many fields of applications, and a complementary analysis tool compared to physical experiments. This includes for instance the safety study of water dams, tsunamis, the extraction of offshore petroleum, the sizing of Liquified Natural Gas (LNG) carriers, processes of phase separation, waste water treatment, flocculation processes, bio-Engineering, medical applications, etc. The evolution of the interfaces between phases and the consecutive complex dynamics need to be simulated for the understanding of the flows and the process optimization in the industrial case. For gas-liquid applications involving fast dynamics, beyond quantities of interest, numerical models have to be able to capture most of the flow features that drive the dynamics, e.g.:

• formation of air pockets;

• ejection, fragmentation of liquid droplets;

• Archimedes buoyancy effect with rising of bubbles and fall of droplets;

• effects of gas compressibility inducing a gas-to-liquid response by a pressure wave, etc.

In this paper, we consider immiscible gas-liquid two-phase flow problems. The strong ratio of mass density between gas and liquid (typically 1:1000) is known to be a source of numerical stiffness and numerical instability. Therefore robust computational approaches supporting high density ratio have to be considered. Among the family of conservative Finite Volume methods (FVM), the Lagrange-Remapped solvers [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to godunov's method[END_REF]; Woodward and [START_REF] Colella | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]; [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF]; [START_REF] Bailey | Lagrange-Remap methods for the Euler equations for single and multi gas flows[END_REF]; [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]; [START_REF] Arber | A Staggered Grid, Lagrangian-Eulerian remap Code for 3-D MHD Simulations[END_REF]...] provide both robustness and stability with achievement of mathematical properties of positiveness and entropy compatibility.

Lagrange-remap numerical schemes (also referred to as Euler-Lagrange schemes) are a particular family of Eulerian FVM where, at each time step, the equations are solved according to a Lagrangian evolution with a mesh that is convected by the flow itself, then the "Lagrange solutions" are remapped on the initial mesh into a conservative way by estimating both fluxed mass and momentum. Because of the Lagrangian step of these methods, code coupling or coupling of different physics is made easier against conventional FV methods. Moreover, the Lagrangian description is very practical for multi-material flows of multi-phase flows because we are a natural control of the fluxed quantities material-by-material. Lagrange and Lagrange-remap solvers still know strong developments today with major contributions as e.g. energy-preserving compatible schemes for staggered methods Caramana et al. (1998b), collocated variables and cell-centered entropy-satisfying schemes, see [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF]; [START_REF] Maire | A cell-centered lagrangian scheme for twodimensional compressible flow problems[END_REF]; [START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF].

In this paper, we rather consider a simpler staggered Lagrange-Remap solver with a directionby-direction remapping. More precisely, the 2D multidimensional Lagrange step (operator L ∆t xy ) completely solves the fluid equations while an operator splitting alternating direction (AD) for the projection (operators R x and R y respectively) is used to interpolate on the reference Cartesian Eulerian grid into a conservative manner, involving convective flux balances. Symmetrized operator splitting can be used to ensure second order accuracy. Both linear and nonlinear numerical stability are ensured by the use of standard pseudo-viscosity (viscous pressure) terms, detailed into an Appendix at the end of this paper. Actually, we use this simple Euler solver because of its simplicity of code implementation and because it can be vectorized/parallelized into a natural manner. Moreover, the aim of this paper is not about the Hydrodynamics solver: the article mostly focuses on numerical antidiffusive methods for interface capturing, as part of a global multifluid Hydrodynamics solver. At the present time, the antidiffusive approaches assumes a direction-by-direction remapping which leads to a simpler derivation of the andiffusive fluxes.

The issue of an interface tracking/capturing algorithm providing expected properties like robustness, accuracy, conservation of volume and mass while not being too much computationally intensive is still the object of today's active Research. Pure Lagrangian approaches like Smoothed Particle Hydrodynamics (SPH) methods naturally captures the moving interfaces because each macro-particle moves with the flow. Each particle is also attached to a given material with its own equation of state (EOS). For liquid-gas flows, we have liquid particles and gas particles and the interface of nothing else but the discrete interface separating liquid particles from gas particles. In the last decade, we have seen in the literature major contributions of improvement in the SPH world with improved accuracy, stability, and ability to tackle multiphase flow problems with high density ratios for violent flow applications, see for example [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF]; [START_REF] Marrone | δ-sph model for simulating violent impact flows[END_REF]; [START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF]. Despite these improvements, SPH still know some issues. Moreover parallelization techniques for SPH are quite technical and require a great expertise (as dealing with strong density ratio between phases). That's why we rather choose a more classical computational approach. Notice also that there are also mesh-based Lagrangian methods [START_REF] Maire | A cell-centered lagrangian scheme for twodimensional compressible flow problems[END_REF]; [START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]). But for flows with large deformations, cells may become degenerate, and both regularization or remapping procedures are needed.

For Eulerian methods, the family of interface tracking methods try to reconstruct the free boundary according to some tracking procedure (level set methods [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]] for example). Interface reconstruction methods try to recontruct a moving interface according to some incomplete information: volume-of-fluid (VOF) methods [START_REF] Noh | SLIC (Simple Line Interface Calculation)[END_REF]; [START_REF] Youngs | Time-depend multi-material flow with large fluid distortion[END_REF]] or MOF (Moment-of-Fluid) methods [START_REF] Dyadechko | Moment-of-fluid interface reconstruction[END_REF]). The family of interface capturing methods involves at the continuous level the transport equation of an indicator function that distinguishes the location of the different materials:

∂ t z + u • ∇z = 0, z ∈ {0, 1}
(where u is the velocity and z the color function or a phase indicator), expressing that the interface property between the two phases is advected with the local fluid velocity. From the numerical point of view, finite volume cells may host different materials. The so-called "mixed cells" need an additional closure but on the other hand the mass fraction c g ∈ [0, 1] of the gas fluid may be used as the numerical indicator function:

∂ t c g + u • ∇c g = 0, c g ∈ [0, 1].
It can be set up in conservative form in order to conserve the gas mass:

∂ t (ρc g ) + ∇ • (ρc g u) = 0
with ρ representing the mean cell "mixture" density, but of course it induces a "diffuse" interface, which has to be kept as less diffusive as possible. Let us emphasize that this diffuse feature is only a pure numerical artefact. This kind of methods uses more or less sophistication levels including high-order schemes, compressive flux limiters, artificial compression stages, local adaptive mesh refinement (AMR) [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF]], a posteriori methods in which an anti-diffusion phase is added after the projection (or the advection) of the quantities, use of pseudo-velocities to correct the truncation error of the numerical scheme [START_REF] So | Anti-diffusion method for interface steepening in two-phase incompressible flow[END_REF]; [START_REF] Hill | A multidimensional positive definite remapping algorithm for arbitrary Lagrangian Eulerian methods[END_REF][START_REF] Navaro | Aéroacoustique numérique d'un écoulement tourbillonnaire[END_REF]...], Eulerian methods with Lagrangian tracking of the interface (VFFC-ENIP [START_REF] Lagoutière | A totally Eulerian finite volume solver for multi-material fluid flows: Enhanced natural interface positioning (ENIP)[END_REF]]), etc. These methods have been proved to be very efficient. But the price to pay is the relative important implementation and computational effort. For parallel computer architectures, the parallel implementation may be tricky or require a tedious work. To summarize what can be said is that each method in the literature shares both advantages and drawbacks. The combination of the three properties conservation-accuracy-robustness is actually highly constrained. The parallel computing aspect may also add strong constraints with heavy programming efforts. Regarding interface capturing methods with a color function, there is two points of view: either a threshold on z (say z = 1 2 ) which discriminates the fluid zones as level sets do, or one considers a possibly smoothed colored function involving a smoothed transition between both fluids. In this case, some "regularization" closure has to be defined into the model. The "uncertainty" related to a z belonging to on the open interval (0, 1) has to be expressed, for example by a volume-averaged mixture closure. This is discussed in the next section.

In this work, we have decided to investigate the use of a recent antidiffusive advection scheme for stepwise solutions, initially proposed by Després and Lagoutière. The idea is to combine both upwinding and downwinding discretizations for the gradient operators. The upwinding process is known to provide strong stability in L p norm under a standard CFL condition while pure downwinding is unconditionally unstable because of its over-compressive nature. Després and Lagoutière then proposed a combination of both upwinding and downwinding "at the limit of stability", providing the most compressive solver while ensuring stability with a local discrete maximum principle. This computational approach has been considered in Kokh and Lagoutière [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF]] and more recently in [START_REF] Billaud Friess | A second order anti-diffusive lagrange-remap scheme for two-component flows[END_REF][START_REF] Billaud Friess | An anti-diffusive lagrange-remap scheme for multi-material compressible flows with an arbitrary number of components[END_REF]] for multiphase flow problems with a "five-equation" model that includes the transport of a color function. In fact the color variable z acts for the interface location whereas another variable y, a gas mass fraction is necessary for conservation purposes. Thus we present here an adaptation of the Després-Lagoutière advection scheme to the case of a simpler "three-equation" gas-liquid volume averaged model. In our case, the gas mass fraction c g also acts as the color variable and is andiffused. We rather reformulate the mass conservation equations using a gas volume fraction α which is of course a function of c g for given phase densities. We use isentropic pressure law equations of state per phase, the liquid phase being seen as a (weakly) compressible fluid. To close the system, we here assume a simple pressure equilibrium in each mixed cell. Methods based on the instantaneous pressure equilibrium usually produce overshoots of some material quantity at the interfaces because the instantaneous pressure equilibrium does not correspond to the behavior of two pure material cells facing each other (see for example [START_REF] Kamm | A pressure relaxation closure model for one-dimensional, twomaterial lagrangian hydrodynamics based on riemann problem[END_REF]; [START_REF] Bachmann | Random sampling remap for compressible two-phase flows[END_REF] for a discussion about this issue). But there are more realistic closure mechanisms like pressure relaxation mechanisms [START_REF] Kamm | A pressure relaxation closure model for one-dimensional, twomaterial lagrangian hydrodynamics based on riemann problem[END_REF] that improve the behavior at the interface. Once again, this paper does not really focus on this topic. We mostly want to evaluate the both advantages, limitations and drawbacks of such an antidiffusive approach, in particular in the case of stiff gas-liquid multiphase flows.

Another motivation that justifies our choice of that type of interface capturing method based on the solution of a transport equation of an Eulerian Cartesian grid is that we believe they are very promising for manycore coprocessor parallel computing, like Intel Xeon PHI or Graphics Processing Units (GPU). While standard PDE discretization methods (mesh, array data structures, sparse matrices, memory) have been designed and optimized a few decades ago on the assumption of a sustainable CPU-global memory model, parallel many-core processing architecture is completely different such that data structures organized by grids of blocks are particularly suitable for finite difference/volume methods on Cartesian grids. At the present time, our computations are performed sequentially on a standard CPU for accuracy assessment of the numerical method and its improvement. Ongoing works are dedicated to the GPU programming of such a method. Let us finally remark that the last update of the NVIDIA CUDA Software Development Kit includes an Application Programming Interface (API) for Adaptive Mesh Refinement (AMR). AMR methods of course could be used with the low-diffusive interface capturing technique for a more accurate level of interface resolution.

Three-equation two-fluid model and equations of state

There are numerous models for gas-liquid two-phase flows subject to gravity forces. In our targeted field of "violent flow" applications with strong topology changes involving wave breaking and air pockets, the Saint-Venant shallow-water equations involving a variable height of the interface are not relevant. We must resort to a multi-fluid system of equations. Let us consider here a phaseseparated model composed of inviscid fluids, separated by a free boundary Γ(t). We use a unique fluid velocity u leading to a momentum conservation equation. The resulting system is supposed to be hyperbolic according to the choices of pressure laws and mixture closure. To simplify, we do not here consider surface tension effects. Let us consider regular smooth solutions with smooth interface boundaries. At instant t, for each spatial point x, one has to consider the continuity equation for the phase k being present:

∂ t ρ k + ∇ • (ρ k u) = 0
where ρ k is the mass density of the fluid k, k = g or ℓ (g (resp. ℓ) stands for the gas (resp. for the liquid)) and u is the fluid velocity. The momentum balance equation reads, for the phase k being present at point

x ∂ t (ρ k u) + ∇ • (ρ k u ⊗ u) + ∇p = ρ k g
considering a pressure p = p(ρ k ) and a constant gravity field g. Denoting by z = z(x, t) the indicator function that returns 0 if the current fluid is the liquid and 1 if it is the gaseous phase, one can define a general mass density ρ, ρ = zρ g + (1 -z)ρ ℓ .

Thus the mass conservation equations write under a condensed form for any phase

∂ t ρ + ∇ • (ρu) = 0.
For the same reason, we have the global momentum equation

∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ρg.
The indicator function z is simply convected by the flow field because its Lagrangian derivative is zero:

D t z = 0 (D t is the Lagrangian derivative, i.e. D t = ∂ t + u • ∇), leading (if preferred) to a conservation law ∂ t (ρz) + ∇ • (ρzu) = 0.
From the numerical point of view, using a fixed mesh that does not move with the flow, we will necessary have some finite volume cells hosting both phases, with respective volume fractions α and (1 -α) for gas and liquid. Averaging operations should also involve phasic velocities u g and u ℓ . But, assuming that gas and liquid phases are sufficiently separated, we consider a unique velocity u.

Introducing the gas mass fraction c g per cell, we will then have c g ∈ [0, 1] numerically. The mass conservation for the gas phase reads

∂ t (c g ρ) + ∇ • (c g ρu) = 0.
The mass density now is seen as a mixture-like density according to the volume fraction variable α ∈ [0, 1]:

ρ = αρ g + (1 -α)ρ ℓ . (1)
We expect the mass conservation of the mixture, that is

∂ t ρ + ∇ • (ρu) = 0. ( 2
)
In particular, we get a transport equation for the gas mass fraction

D t c g = ∂ t c g + u • ∇c g = 0.
Actually, the gas mass fraction c g can play the role of a "smoothed" indicator function. The mass conservation of the different phase can also be written using volume fraction variables, for example for the gas phase ∂ t (αρ g ) + ∇ • (αρ g u) = 0.

(3)

We then have the relationship between α and c g

c g = αρ g ρ = αρ g αρ g + (1 -α)ρ ℓ or α = c g ρ ℓ c g ρ ℓ + (1 -c g )ρ g . (4)
We still have a momentum balance equation considering the mixture momentum on the volume. That implies the definition of pressures in this volume. As a simple closure model, we here consider pressure equilibrium between phases p = p g (ρ g ) = p ℓ (ρ ℓ ).

(5)

As discussed in the Introduction, more realistic closures involving for example pressure relaxation [START_REF] Kamm | A pressure relaxation closure model for one-dimensional, twomaterial lagrangian hydrodynamics based on riemann problem[END_REF]) may be used, but this is not at the aim of this paper.

Equations of state, speeds of sound

We discuss here the choice of the equation of state (EOS) for each phase. We assume an isentropic evolution for each fluid, with a phasic pressure p k as a function of the phasic mass density, i.e., p k = p k (ρ k ). To get hyperbolicity, we will assume

∂p k ∂ρ k = c 2 s,k > 0
where c k denotes the speed of sound into the fluid k. For the gas, we decide to use an isentropic perfect gas law :

p g (ρ g ) = p 0 ρ g ρ 0 g γg (6)
with constant reference pressure p 0 , density ρ 0 g and specific heat ratio γ g ∈ (1, 3]. The gas speed of sound c s,g is given by c s,g = γ g p g ρ g .

An isentropic approximate Tait equation of state is used for the liquid phase [START_REF] Bailey | Lagrange-Remap methods for the Euler equations for single and multi gas flows[END_REF]]:

p ℓ (ρ ℓ ) = p 0 + p 0 K ρ ℓ ρ 0 ℓ γ ℓ -1 , (7)
with reference pressure p 0 and density ρ 0 ℓ , γ ℓ > 1 and where K =

ρ 0 ℓ c2 s,ℓ
γ l p 0 is the bulk modulus (c s,ℓ is a reference liquid speed of sound). With K ≫ 1, this EOS assumes that the liquid is quasiincompressible since an high change of pressure entertains very slow variations on the density 1,2 . For the liquid speed of sound, we have

c 2 s,ℓ = c2 s,ℓ ρ ℓ ρ 0 ℓ γ l -1 , (8)
that we can rewrite as :

c 2 s,ℓ = c2 s,ℓ ρ 0 ℓ ρ ℓ + γ ℓ (p -p 0 ) ρ ℓ
1 For instance, with a speed of sound cs,ℓ = 350 m.s -1 , ρ 0 ℓ = 1000 kg.m -3 and γ ℓ = 7: from a pressure of reference p0 = 10 5 Pa, a division by 2 of this pressure yields to a variation of the density of the order of 0.004%.

2 The treatment of the liquid as a compressible fluid avoids to invert a Poisson system and gives us a real pressure, the drawback is the high celerity of sound of the water which restrains the time step of the numerical scheme.

showing the stiffened gas-like behaviour. The validity domain of this EOS is expressed by c 2 s,ℓ > 0, so we will assume that the pressure will always satisfies the inequality

p > max 0, p 0 - ρ 0 ℓ c2 s,ℓ γ ℓ .
Let us add a comment on the whole hyperbolicity of the system. Of course, under pure phase conditions, the existence of speeds of sound for the EOS guarantees the property. Let us consider now mixed cell conditions that arise into the diffuse interface (α ∈ (0, 1)) and consider the 1D case for simplicity. Setting up the system in variables (ρ, c g , u) (c g is the gas mass fraction), we have for smooth solutions

∂ t ρ + u∂ x ρ + ρ∂ x u = 0, ∂ t c g + u∂ x c g = 0, ∂ t u + u∂ x u + 1 ρ ∂ x p = g.
Considering p as a function of ρ and c g , i.e. p = p(ρ, c g ), the characteristic polynomial is given by

P (λ) = (u -λ) (u -λ) 2 -∂p ∂ρ .
Under pressure equilibrium conditions, we respectively have

α = ρc g ρ g (p) , 1 -α = (1 -c g )ρ ρ ℓ (p) . Using the identity α + (1 -α) = 1, we get c g ρ g (p) + (1 -c g ) ρ ℓ (p) = 1 ρ .
By derivating w.r.t ρ this expression, we get

∂p ∂ρ - c g ρ ′ g (p) (ρ g (p)) 2 - (1 -c g )ρ ′ ℓ (p) (ρ ℓ (p)) 2 = - 1 ρ 2 . As ρ ′ g (p) = 1 c 2 s,g and ρ ′ ℓ (p) = 1 c 2 s,ℓ
, then we derive the mixture speed of sound c m,s of the system:

c 2 m,s = ∂p ∂ρ = 1 c g ρ ρg(p) 2 1 c 2 s,g + (1 -c g ) ρ ρ ℓ (p) 2 1 c 2 s,ℓ .
The system is thus hyperbolic, with real eigenvalues u -c m,s , u and u + c m,s .

Summary

To summarize, there are two equivalent ways to represent the system, either using the volume fraction α:

∂ t αρ g + ∇ • (αρ g u) = 0, ( 9
) ∂ t (1 -α)ρ ℓ + ∇ • (1 -α)ρ ℓ u = 0, ( 10
)
∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ρg, (11)
or using the mass fraction c g :

∂ t ρ + ∇ • (ρu) = 0, ( 12
) ∂ t (ρc g ) + ∇ • (ρc g u) = 0, ( 13
)
∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ρg, (14)
adding a pressure equilibrium closure p = p g (ρ g ) = p ℓ (ρ ℓ ) for mixed cells. This system is hyperbolic.

EOS parameters

In our targeted applications, we consider air-water flows. For simplicity, the air phase will be modeled using a perfect gas law, with γ g = 1.4. For the liquid water phase, we will use γ ℓ = 7. We will work near under atmospheric conditions, and thus we will practically use the following parameters: ρ 0 g = 1.28 kg.m -3 , ρ 0 ℓ = 1000 kg.m -3 , p 0 = 10 5 Pa, c ℓ = 1500 m.s -1 . We can consider during the simulation a constant celerity of sound since c s,ℓ ∼ cs,ℓ thanks to (8), and due to the weakly compressible caracter of the water. Moreover, in order to lower the stiffness between fluid velocities and sound speed (low Mach number conditions), we can artificially lower the liquid speed of sound, for example cℓ = 350 m.s -1 (keeping the weakly compressible character of the water, cf [START_REF] Monaghan | Simulating free surface flows with sph[END_REF]]).

Numerical scheme

In this section we describe a staggered Lagrange+Remap scheme (also referred to as remapped Lagrange scheme). The solver here is a variant version of the BBC scheme discussed in Woodward and Collela [Woodward and Colella (1984)] also discussed in [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF][START_REF] De Vuyst | Staggered lagrange-remap schemes in conservation form[END_REF]]. This Eulerian scheme uses a two-dimensional Cartesian staggered grid where the velocity variables are defined at the edge midpoints while all the other quantities are defined at cell centers (Fig. 1). In the following we write the constant spatial steps ∆x and ∆y such that the coordinates of the center of the cell are x i = i∆x and y j = j∆y (for simplicity we write x i,j = (x i , y j )). Note that the volume of the Eulerian cell is constant, given by V i,j = ∆x∆y.

Cell i,j x ∆ y ∆ 1/2, i j u + 1/2, i j u - , 1/2 i j v - , 1/2 i j v +
Figure 1: Staggered mesh of the BBC scheme: the first component of the velocity u is defined at the vertical edge midpoints and the second v at the horizontal edge midpoints. All the other physical quantities (like α and the pressure p) are defined at the center of the cell.

Lagrangian step

This step allows us to solve the Euler equations written in Lagrangian form:

ρ D t 1 ρ -∇ • u = 0, ( 15
)
ρ D t u + ∇p = g, (16)
D t c g = 0, ( 17
)
where D t = ∂ t + u • ∇ is the Lagrangian (particle) derivative. In that scheme, the derivatives in space are centered. To achieve numerical stability, a pseudo-viscosity pressure q is added. Artifical viscosity only acts on compression areas (∇ • u < 0) allowing us to keep second-order accuracy into expansion zones [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF], [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], Caramana et al. (1998a)] (in AppendixB the exact form of the pseudo-viscosity used is described).

To get second order of accuracy in time we use a multi-step scheme with the help of two intermediate time steps t n+1/4 = t n + ∆t 4 and t n+1/2 = t n + ∆t 2 before the final time step t n+1 = t n + ∆t. To summarize, the Lagrangian step consists of:

• A first step of prediction of an intermediate velocity at the time t n+1/4 :

u n+1/4,L i+1/2,j = u n i+1/2,j - ∆t 4 ∆y m n i+1/2,j (p + q) n i+1,j -(p + q) n i,j , (18)
v n+1/4,L i,j+1/2 = v n i,j+1/2 - ∆t 4 ∆x m n i,j+1/2 (p + q) n i,j+1 -(p + q) n i,j + ∆t 4 g, (19)
where the total masses at the edges are defined by: m n i+1/2,j =

m n i,j +m n i+1,j 2 and m n i,j+1/2 = m n i,j +m n i,j+1 2 (m n i,j = m n g i,j + m n ℓ i,j
, the sum of the mass of gas and the mass of liquid into the finite volume cell (i, j)).

• A second step of prediction of the quantities at the time t n+1/2 using the new velocities at t n+1/4 :

V n+1/2,L i,j = V n i,j + ∆t 2 ∆y u n+1/4,L i+1/2,j -u n+1/4,L i-1/2,j + ∆t 2 ∆x v n+1/4,L i,j+1/2 -v n+1/4,L i,j-1/2 ; ( 20
)
from the knowledge of the partial mass of each fluid m n g i,j , m n ℓ i,j and the updated volume

V n+1/2,L i,j
, the pressure equilibrium assumption requires the solution on an algebraic problem set up in variable α:

p n+1/2,L i,j = p g m n g i,j α n+1/2,L i,j V n+1/2,L i,j if m n ℓ i,j = 0 ⇒ α n+1/2,L i,j = 1, p n+1/2,L i,j = p ℓ m n ℓ i,j (1-α n+1/2,L i,j )V n+1/2,L i,j if m n g i,j = 0 ⇒ α n+1/2,L i,j = 0, p n+1/2,L i,j = p g m n g i,j α n+1/2,L i,j V n+1/2,L i,j = p ℓ m n ℓ i,j (1-α n+1/2,L i,j )V n+1/2,L i,j
otherwise.

(21)

Here, we remind that ρ g = mg αV and ρ ℓ = m ℓ (1-α)V and the conservation of the partial masses of each phase comes from (17) as c g = mg mg+m ℓ . This algebraic problem is numerically solved by a Picard fixed-point algorithm. In AppendixA we give details on the numerical solution which is key for the overall performance of the method.

To summarize, thanks to the pressure equilibrium assumption, we get the values of

α n+1/2,L i,j , p n+1/2,L i,j , ρ n+1/2,L g i,j
and ρ n+1/2,L ℓ i,j at the new time step (n + 1/2). This Lagrangian step keeps both gas and liquid masses constant. Thanks to the new pressure p n+1/2,L i,j

, the velocities at time step t n+1/2 are then updated:

u n+1/2,L i+1/2,j = u n i+1/2,j -∆t 2 ∆y m n i+1/2,j (p + q) n+1/2,L i+1,j -(p + q) n+1/2,L i,j , (22)
v n+1/2,L i,j+1/2 = v n i,j+1/2 -∆t 2 ∆x m n i,j+1/2 (p + q) n+1/2,L i,j+1 -(p + q) n+1/2,L i,j + ∆t 2 g. (23)
• A third step that enables us to get the variables on the distorted cells (of volumes V n+1,L i,j

(24)), using the predicted quantities at t n+1/2 :

V n+1,L i,j = V n i,j + ∆t∆y u n+1/2,L i+1/2,j -u n+1/2,L i-1/2,j + ∆t∆x v n+1/2,L i,j+1/2 -v n+1/2,L i,j-1/2 . ( 24
)
The new values of α n+1,L i,j

, p n+1,L i,j , ρ n+1,L g i,j
and ρ n+1,L ℓ i,j

are then obtained thanks to the solution of the pressure equilibrium assumption (Eq. ( 21)) in which we replace the volume

V n+1/2,L i,j
and α n+1/2,L i,j

by their updated values V n+1,L i,j and α n+1,L i,j

, the masses remaining constant.

Finally, the velocities at time t n+1,L are given by the extrapolations procedure:

u n+1,L i+1/2,j = 2u n+1/2,L i+1/2,j -u n i+1/2,j , (25)
v n+1,L i,j+1/2 = 2v n+1/2,L i,j+1/2 -v n i,j+1/2 . ( 26
)
Expressions ( 25)-( 26) are condensed form of conservative momentum evolution.

Projection step

At the end of the Lagrangian step, we get updated values on a moved mesh. The so-called remap step is required to get the quantities on the initial cartesian grid. The conservative projection step is performed at a given order accuracy (current work only deals with a first order accuracy), however we take a particular attention to the projection related to the volume fractions in order not produce too much artificially mixture by a numerical diffusion.

The remap step is performed in two sub-steps following an alternating direction fractional step procedure:

1. Projection along the x direction. The vertical left and right edges are turned back on their initial positions, giving a cell of volume:

V n+1, * i,j = V n+1,L i,j -∆t∆y(u n+1/2,L i+1/2,j -u n+1/2,L i-1/2,j ), (27)
2. Projection along the y direction. We thus recover the initial cartesian grid:

V n+1 i,j = V n+1, * i,j -∆t∆x v n+1/2,L i,j+1/2 -v n+1/2,L i,j-1/2 (28) = ∆x∆y. ( 29
)
In each subset, the projection of the partial masses and the velocities are performed without equilibrium in pressure requirement. Consistance and stability properties are asked to the remap numerical scheme by imposing conditions to find a suitable value for α at the interface, appearing in the following formulas ( 30)-( 31) for the first step (resp. ( 33)-( 34) for the second) in order not to create extremum. After the remap, a step of equilibrium in pressure is performed varying the partial volume fraction of each phase in order to get the volume fraction α at equilibrium. The conservative projection of masses are given by:

• Projection step along x: (the quantities are obtained on the intermediate volume V n+1, * ij defined by ( 27)):

m n+1, * g i,j = m n g i,j -∆t∆y u n+1/2,L i+1/2,j αρ g n+1,L i+1/2,j -u n+1/2,L i-1/2,j αρ g n+1,L i-1/2,j , (30)
m n+1, * ℓ i,j = m n ℓ i,j -∆t∆y u n+1/2,L i+1/2,j (1 -α)ρ ℓ n+1,L i+1/2,j -u n+1/2,L i-1/2,j (1 -α)ρ ℓ n+1,L i-1/2,j , (31)
with the quantities αρ g n+1,L i±1/2,j for the gas (resp. (1 -α)ρ ℓ n+1,L i±1/2,j for the liquid) corresponding to an interpolation of the quantity αρ g (resp.

(1 -α)ρ ℓ ) at the interface (left or right) i ± 1/2, j of the cell i, j. As the volume and the partial masses of each phase have changed during the projection step, the solution of the pressure equilibrium relation is needed to get (among others) the value of α n+1, * i,j

that we need for the next step:

p n+1, * i,j = p g m n+1, * g i,j α n+1, * i,j V n+1, * i,j if m n+1, * ℓ i,j = 0 ⇒ α n+1, * i,j = 1, p n+1, * i,j = p ℓ m n+1, * ℓ i,j (1 -α n+1, * i,j )V n+1, * i,j if m n+1, * g i,j = 0 ⇒ α n+1, * i,j = 0, p n+1, * i,j = p g m n+1, * g i,j α n+1, * i,j V n+1, * i,j = p ℓ m n+1, * ℓ i,j (1 -α n+1, * i,j )V n+1, * i,j otherwise, ( 32
)
where m n+1, * g i,j , α n+1, * i,j

, m n+1, * ℓ/g i,j and V n+1, * i,j

are the intermediate projected values. We still refer to the AppendixA for details on the algebraic solution.

• Projection step along y: at the end of this step, the quantities are known on the volume V n+1 ij = ∆x∆y, and thus on the cartesian grid: for the liquid) which approximate the quantities αρ g (resp.

m n+1 g i,j = m n+1, * g i,j -∆t∆x v n+1/2,L i,j+1/2 αρ g n+1, * i,j+1/2 -v n+1/2,L i,j-1/2 αρ g n+1, * i,j-1/2 , (33)
m n+1 ℓ i,j = m n+1, * ℓ i,j -∆t∆x v n+1/2,L i,j+1/2 (1 -α)ρ ℓ n+1, * i,j+1/2 -v n+1/2,L i,j-1/2 (1 -α)ρ ℓ n+1, * i,j-1/2 , ( 34
(1 -α)ρ ℓ ) at the interface (top or bottom) i, j ± 1/2 of the cell i, j.

Once again, after the projection step, we need to solve the pressure equilibrium assumption:

p n+1 i,j = p g m n+1 g i,j α n+1 i,j V n+1 i,j if m n+1 ℓ i,j = 0 ⇒ α n+1 i,j = 1, p n+1 i,j = p ℓ m n+1 ℓ i,j (1 -α n+1 i,j )V n+1 i,j if m n+1 g i,j = 0 ⇒ α n+1 i,j = 0, p n+1 i,j = p g m n+1 g i,j α n+1 i,j V n+1 i,j = p ℓ m n+1 ℓ i,j (1 -α n+1 i,j )V n+1 i,j otherwise, ( 35
)
and we get the final values at time t n+1 on the Eulerian cartesian grid: α n+1 i,j , p n+1 i,j , ρ n+1 g i,j and ρ n+1 ℓ i,j .

In the formulas of projection of the masses ( 30)-( 31) and ( 33)-( 34), for first order accuracy, the most obvious stable choice to define the interpolated values of α, ρ g and ρ l at the edges is an upwind strategy:

αρ g n+1,L/ * i+1/2,j,up = α n+1,L/ * i,j ρ n+1,L/ * g i,j if u n+1/2,L i+1/2,j > 0 (36) = α n+1,L/ * i+1,j ρ n+1,L/ * g i+1,j if u n+1/2,L i+1/2,j ≤ 0 (37)
and

αρ g n+1,L/ * i,j+1/2,up = α n+1,L/ * i,j ρ n+1,L/ * g i,j if v n+1/2,L i,j+1/2 > 0 (38) = α n+1,L/ * i,j+1 ρ n+1,L/ * g i,j+1 if v n+1/2,L i,j+1/2 ≤ 0, ( 39
)
and the same for the quantities related to the liquid phase. Such a choice is actually irrelevant to simulate immiscible fluids since the projection on volume fraction is too much diffusive and does not enable to follow a thin interface (cf. Fig. 2a). Let us emphasize that second order projection is still not sufficient to keep a numerical thin interface (cf. Fig. 2b). With the first order upwind projection (Fig. 2a), the interface between the liquid (in blue) and the air (in red) is too much diffusive to be followed with accuracy (e.g. droplets and air pockets are hidden because of the diffusion). With the second order upwind projection (Fig. 2b), the interface is less diffuse and we begin to see an air pocket which were previously concealed by the diffusion. But the diffusion is still present and we do not observe formation of droplets.

For details on the initialization, please refer to the Test case II of the section 5 dedicated to the numerical experiments.

We did not detail the projection formulas for the momentum. This can be found in [START_REF] De Vuyst | Staggered lagrange-remap schemes in conservation form[END_REF]].

In the next section, we explain how we improved this remap step in order to keep a numerical thin interface between gas and liquid.

A low-diffusive procedure to choose the fluxes at the edges in the remap numerical scheme

As initially discussed into the introduction, our strategy here is to design a computational approach which is quite simple to implement but can also be implemented in parallel in a rather straightforward way. Our construction follows ideas initially proposed by [START_REF] Lagoutière | A totally Eulerian finite volume solver for multi-material fluid flows: Enhanced natural interface positioning (ENIP)[END_REF]; Després and [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF]] to build low-diffusive advection schemes.

The principle is to take advantage of both upwind scheme (known to be strongly stable but artificially diffusive) and downwind scheme (low-diffusive but unstable) in order to obtain a stable scheme with minimal numerical diffusion. This family of anti-diffusive schemes has been extended to the case of multiphase flows [START_REF] Billaud Friess | A second order anti-diffusive lagrange-remap scheme for two-component flows[END_REF]; [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF]]. We here adapt the methodology to the present "three-equation" two-phase system. First, we recall the conservation form of the mass conservation (Eq. 2) in terms of c g

∂ t (ρc g ) + ∇ • (ρc g u) = 0. ( 40
)
It is equivalent to the mass conservation of the gas phase (3):

∂ t (αρ g ) + ∇ • (αρ g u) = 0.
We immediately deduce the nonconservative transport equation for the gas mass fraction (thanks to equations ( 1) and ( 40):

D t c g = ∂ t c g + u • ∇c g = 0. ( 41
)
This linear transport equation implies a local maximum principle on the variable c g . It is easy to check that α and c g (4) are linked according to the dual relations:

α(c g ; ρ g , ρ ℓ ) = c g ρ ℓ c g ρ ℓ + (1 -c g )ρ g , c g (α; ρ g , ρ ℓ ) = αρ g αρ g + (1 -α)ρ ℓ . (42)
From the numerical point of view, the maximum principle checked by c g into (41) should also be verified. The idea is to get the "best" interface value of c g in order to ensure the local maximum principle while being as compressive as possible. But, because of the dependency of c g on α, ρ g and ρ ℓ (second formula in ( 41)), we need to express the maximum principle of c g in terms of interface values of α, ρ g and ρ ℓ :

c g i+1/2,j = α n+1,L i+1/2,j (ρ g ) n+1,L i+1/2,j α n+1,L i+1/2,j (ρ g ) n+1,L i+1/2,j + (1 -α n+1,L i+1/2,j )(ρ ℓ ) n+1,L i+1/2,j . (43)
Following once again ideas from [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF]], we decide to simply set the upwind values for ρ g and ρ ℓ . The remaining degree of freedom is the value of α n+1,L i+1/2,j which has to be optimized, subject to consistency and stability requirements (on c g ). For simplicity, we only describe the first step of the x-projection, considering variables at stage (n+1, L) in order to get intermediate projected variables at (n + 1, * ). From now on, the sought value of α is denoted α LD i±1/2,j (LD stands for low diffusive). The process is:

1. first, to define a trust interval I as the intersection of intervals in which α LD i+1/2,j must be in order to check consistency and stability properties for the projection. The intersection is ensured not to be empty since the diffusive value α n+1,L i+1/2,j,up is known to belong to all of these intervals; 2. then to take for α LD i+1/2,j the nearest value from the downwind one α n+1,L i+1/2,j,down (most compressive choice) while remaining in the trust interval 3 .

3 We recall that for a generic variable z, the upwind and downwind values correspond to:

z n+1,L i+1/2,j,up = z n+1,L i,j if u n+1/2,L i+1/2,j > 0 = z n+1,L i+1,j if u n+1/2,L i+1/2,j ≤ 0,
The definition of the trust interval I, in which we will choose the value α LD i+1/2,j at the edges (i + 1/2, j), is subject to some consistency and stability requirements.

Consistency requirement. The value at the edge c g i+1/2,j must be between the values at the left and right cells of the edges, which means

c n+1,L g i+1/2,j ∈ [a n+1,L i+1/2,j , A n+1,L i+1/2,j ], (44)
where a n+1,L i+1/2,j = min(c n+1,L g i,j

, c n+1,L g i+1,j ), A n+1,L i+1/2,j = max(c n+1,L g i,j

, c n+1,L g i+1,j ), ( 45)

with at the edges, c g i+1/2,j expressed thanks to its definition (43). In order to define the first condition for the trust interval, we have to find a sufficient condition α LD i+1/2 ∈ [d n+1,L i+1/2,j , D n+1,L i+1/2,j ] to check the condition (44) 4 . We refer to [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF] for details, we get the same value for [d n+1,L i+1/2,j , D n+1,L i+1/2,j ] using the definition of c g (43) and of ρ (1):

α n+1,L i+1/2,j ∈ I 1 := [d n+1,L i+1/2,j , D n+1,L i+1/2,j ] ⇒ (44), (46)
with the lower bound d n+1,L i+1/2,j :

d n+1,L i+1/2,j = ρ n+1,L ℓ i+1/2,j,up a n+1,L i+1/2,j ρ n+1,L g i+1/2,j,up (1 -a n+1,L i+1/2,j ) + ρ n+1,L ℓ i+1/2,j,up a n+1,L i+1/2,j (47)
and the upper bound D n+1,L i+1/2,j :

D n+1,L i+1/2,j = ρ n+1,L ℓ i+1/2,j,up A n+1,L i+1/2,j ρ n+1,L g i+1/2,j,up (1 -A n+1,L i+1/2,j ) + ρ n+1,L ℓ i+1/2,j,up A n+1,L i+1/2,j . (48)
Note that the interval [d n+1,L i+1/2,j , D n+1,L i+1/2,j ] ⊂ [0, 1] is well defined since ρ g > 0 and ρ l > 0 in all the cells if the pressure p > 0 thanks to the definition of the EOSs ( 6)-( 7) and thus the denominator is strictly positive (recording that a n+1,L i+1/2 and A n+1,L i+1/2 all belong to [0, 1], see footnote 4). Thus we have defined the first interval I 1 (46) to define the final trust interval I. It can be easily proved that the upwind value α n+1,L i+1/2,j,up belongs to I 1 . Note that this interval I 1 does not differ from the original referred article (since only the definition of the mass fraction c g (43) is used).

Stability requirement. We now define the second interval I 2 used for the definition of the final trust interval. This interval depends on the sign of the velocity at the edges u n+1/2,L i+1/2,j which determines in which cell the stability condition is calculated. We detail the calculus of the interval in the positive case (that we denote I + 2 ), the other case can be found in AppendixC and is denoted

I - 2 . and z n+1,L i+1/2,j,down = z n+1,L i+1,j if u n+1/2,L i+1/2,j > 0 = z n+1,L i,j if u n+1/2,L i+1/2,j ≤ 0.
4 Note that a n+1,L i+1/2,j and A n+1,L i+1/2,j belong to the interval [0, 1].

The projected value c n+1, * g i,j

must respect the following stability condition (in the Leroux-Harten meaning): if the velocity at the edge i + 1/2 is positive (u n+1/2,L i+1/2,j > 0) and if u n+1/2,L i-1/2,j > 0, to ensure the stability of c g in cell i, j5 , c n+1, * g i,j

must check:

a n+1,L i-1/2,j ≤ c n+1, * g i,j ≤ A n+1,L i-1/2,j , (49)
with a n+1,L i-1/2,j and A n+1,L i-1/2,j defined by ( 45). This new condition (49) allows us to have the positivity of the partial masses during the projection step since a n+1,L i-1/2,j ∈ [0, 1] and A n+1,L i-1/2,j ∈ [0, 1]. The projection of the masses per phase can be equivalently rewritten in terms of ρc g and ρ thanks to the definition of c g (4): ( 30)-( 31) are equivalent to:

(c g ρ) n+1, * ij = V n+1,L ij V n+1, * ij (c g ρ) n+1,L ij -∆t∆y V n+1, * ij u n+1/2,L i+1/2,j (c g ) n+1,L i+1/2,j ρ n+1,L i+1/2,j -u n+1/2,L i-1/2,j (c g ) n+1,L i-1/2,j ρ n+1,L i-1/2,j , (50)
(ρ) n+1, * ij = V n+1,L ij V n+1, * ij (ρ) n+1,L ij -∆t∆y V n+1, * ij u n+1/2,L i+1/2,j ρ n+1,L i+1/2,j -u n+1/2,L i-1/2,j ρ n+1,L i-1/2,j , (51)
with the intermediate V n+1, * i,j defined by ( 27):

V n+1, * i,j = ∆x∆y + ∆t∆x v n+1/2,L i,j+1/2 -v n+1/2,L i,j-1/2 . ( 52
)
We now characterize the interval I 2 (I + 2 in the case of positive velocities) by the following successive results:

Proposition 1. If u n+1/2,L i+1/2,j > 0 and if u n+1/2,L
i-1/2,j > 0, the stability condition is calculated on the cell (i, j) and if we take

α n+1,L i+1/2,j ∈ I + 2 := [b i+1/2,j , B i+1/2,j ], with b n+1,L i+1/2,j = α n+1,L i,j + ρ n+1,L i,j ρ n+1,L g i+1/2,j,up (1 -A n+1,L i-1/2,j ) + A n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up >0 A n+1,L i-1/2,j -c n+1,L g i,j u n+1/2,L i-1/2,j u n+1/2,L i+1/2,j - V n+1, * i,j ∆t∆y u n+1/2,L i+1/2,j ≤0 under the condition (55) , (53)
and

B n+1,L i+1/2,j = α n+1,L i,j + ρ n+1,L i,j ρ n+1,L g i+1/2,j,up (1 -a n+1,L i-1/2,j ) + a n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up >0 a n+1,L i-1/2,j -c n+1,L g i,j u n+1/2,L i-1/2,j u n+1/2,L i+1/2,j - V n+1, * i,j ∆t∆y u n+1/2,L i+1/2,j ≥0 under the condition (55) , (54)
we ensure c n+1, * g i,j

to stay in the interval defined by ( 49) and thus we guarantee the stability of the remap scheme. Moreover, under the CFL-like cell strain limitation condition

V n+1, * i,j -∆t∆y u n+1/2,L i-1/2,j ≥ 0 (55)
the upwind value α n+1,L i,j belongs to both I + 2 and I 1 .

Proof. We first consider the inequality:

a n+1,L i-1/2,j ≤ c n+1, * g i,j (56)
by multiplying it by ρ n+1, * i,j

and by using ( 50)-( 51), we get the condition:

c n+1,L g i+1/2,j ρ n+1,L i+1/2,j ≤ V n+1,L i,j ∆t∆y ρ n+1,L i,j u n+1/2,L i+1/2,j c n+1,L g i,j -a n+1,L i-1/2,j + a n+1,L i-1/2,j ρ n+1,L i+1/2,j . ( 57
)
This condition is a sufficient one: indeed the definition of a n+1,L i-1/2 (45) involves that the quantity

ρ n+1,L i-1/2,j u n+1/2,L i-1/2,j u n+1/2,L i+1/2,j (c n+1,L g i-1/2,j -t n+1,L i-1/2,j
) ≥ 0, which should appear in the right side of (57), can be erased to give a sufficient condition (57) independent of the edge i -1/2 to check the stability condition (56). Thus, replacing c g ρ by αρ g and using the definition of ρ (1) we rewrite the sufficient condition (57) to check (56) in terms of conditions on the value of α at the edges:

α n+1,L i+1/2,j ρ n+1,L g i+1/2,j,up (1 -a n+1,L i-1/2,j ) + a n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up :=Q 2 ≤ V n+1,L i,j ∆t∆y ρ n+1,L i,j u n+1/2,L i+1/2,j c n+1,L g i,j -a n+1,L i-1/2,j + ρ n+1,L ℓ i+1/2,j,up a n+1,L i-1/2,j :=Q 1 . ( 58
)
We thus prove that the upwind value α n+1,L i+1/2,j,up = α n+1,L i,j

(since u n+1/2,L i+1/2,j > 0 in that case) is lower than the quantity Q 1 /Q 2 , by displaying α n+1,L i+1/2,j,up in the term Q 1 , in order to construct an nonempty trust interval by intersection of the two intervals I 1 and I 2 . This should be done by rewriting term Q 1 as follows:

Lemma 4.1. The term Q 1 defined in (58) can be rewritten equivalently as:

Q 1 = α n+1,L i,j ρ n+1,L g i+1/2,j,up (1 -a n+1,L i-1/2,j ) + a n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up +ρ n+1,L i,j c n+1,L g i,j -a n+1,L i-1/2,j ≥0 V n+1, * i,j -∆t∆y u n+1/2,L i-1/2,j ∆t∆y u n+1/2,L i+1/2,j . ( 59
)
Proof. First, we can rewrite the term Q 1 of (58) as:

Q 1 = ρ n+1,L i,j c n+1,L g i,j -a n+1,L i-1/2,j + a n+1,L i-1/2,j ρ n+1,L ℓ i,j +ρ n+1,L i,j c n+1,L g i,j -a n+1,L i-1/2,j V n+1,L i,j ∆t∆y u n+1/2,L i+1/2,j -1 . ( 60
)
Thanks to the definition of c g (4) and ρ (1), and since the upwind values on the edges i + 1/2 correspond to the values of the variables in the cell i, j (u n+1/2,L i+1/2,j > 0), the following relation is checked:

ρ n+1,L i,j c n+1,L g i,j -a n+1,L i-1/2,j + a n+1,L i-1/2,j ρ n+1,L ℓ i,j = α n+1,L i,j ρ n+1,L g i+1/2,j,up (1 -a n+1,L i-1/2,j ) + a n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up . (61)
Then, using ( 61) in ( 60) and the definition of the intermediate volume V n+1, * i,j

(27), we can rewrite Q 1 in the form (59) of Lemma 4.1.

Let's go back to the proof of Proposition 1. Besides, as ρ n+1,L g/ℓ i+1/2,j,up

> 0 if p > 0, then Q 2 = ρ n+1,L
g i+1/2,j,up (1 -a n+1,L i-1/2,j ) + a n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up > 0 and a sufficient condition to respect the condition of stability ( 56) is to consider

α n+1,L i+1/2,j ≤ B n+1,L i+1/2,j , (62)
with B n+1,L i+1/2,j defined by ( 54). Finally, if the CFL condition ( 55) is respected, we see immediately that α n+1,L,upw i+1/2,j ≤ B n+1,L i+1/2,j . The treatment of the other inequality is similar (we still have ρ n+1,L g i+1/2,j,up (1 -

A n+1,L i-1/2,j ) + A n+1,L i-1/2,j ρ n+1,L ℓ i+1/2,j,up > 0 since ρ n+1,L g/ℓ i+1/2,j,up
> 0 if the pressure p > 0) and we get the definition (53) for the lower boundary b n+1,L i+1/2,j of the interval I + 2 . Thus, if the condition ( 55) is checked, then

α n+1,L i+1/2,j,up = α n+1,L i,j ∈ [b n+1,L i+1/2,j , B n+1,L i+1/2,j ] := I + 2 (63)
and a new interval I + 2 is defined and is added to the trust interval for the choice of the low diffusive value α LD i+1/2,j . Finally, we summarize the results obtained in that section, that is the construction of a nonempty trust interval depending on the sign of the edges velocities to find a low-diffusive value for the flux α LD i+1/2,j used in the formula of projections ( 30)-( 31).

Theorem 4.2. Under the condition to be respected by the time step (in which s = sign(u

n+1/2,L i+1/2,j )) V n+1, * i+1/2,j,upw -s∆t∆y u i+1/2-s,j ≥ 0, ( 64
)
when u n+1,L i+1/2,j u n+1,L i+s/2,j > 0 (i.e. when the velocities at the edges of the cell where the stability condition is calculated are of the same sign), the value of α LD i+1/2,j can be taken in the following trust interval I:

I = I 1 consistency for cg ∩ I s 2 stability for cg := [ω n+1,L i+1/2,j , Ω n+1,L i+1/2,j ] ∈ [0, 1], ( 65
)
which is nonempty since the upwind value α n+1,L i+1/2,j,up ∈ I, where the interval I 1 are defined by (46) and I s 2 by (C.1)-(C.2). Moreover, taking α LD i+1/2,j ∈ I ensures to respect maximum principle on c g and especially to keep the positivity of the masses of each phases during the projection6 .

Remark. Note that the trust interval is only defined in two cases: if

u n+1/2,L i+1/2,j > 0 and u n+1/2,L i-1/2,j > 0 or if u n+1/2,L i+1/2,j < 0 and u n+1/2,L
i+3/2,j < 0. In the other cases, the procedure described above takes the upwind choice without the need of a trust interval.

Remark. Contrary to [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF], in our case the value of α after each step of projection naturally respects the maximum principle (α stays in [0, 1] during the solution of the pressure equilibrium algebraic problem, cf. AppendixA.2) without having to add extra conditions. Besides, the volume fraction α (obtained in terms of a balance of partial volumes) depends on the other quantities of the system ρ, c g and u, whereas the color function of [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF] (which is advected) only depends on the velocities.

Remark. For the second step of projection, the procedure to define the trust interval is exactly the same by replacing the quantities at time n + 1, L used in the previous formulas or the lowdiffusive value by their updated values at time n + 1, * ; by substituting the quantity ∆yu n+1/2,L i±1/2,j by ∆xv n+1/2,L i,j±1/2 and the intermediate volume

V n+1, * i,j by V n+1 i,j = ∆x∆y.
For the choice of the low-diffusive value at the edges α LD i+1/2,j thanks to the knowledge of the trust interval I, we refer to [START_REF] Kokh | An numerical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model[END_REF] using the trust interval defined above (or see the remainder in the AppendixD).

Time step restriction and CFL condition

The implemented low-diffusive procedure does not restrict the time step although it needs to respect the condition (64). Indeed, this will be checked if we impose that the vertical edges of the cells cannot move of more than ∆x 4 and the horizontal edges of the cells of more than ∆y 4 . This is in particular valid if we impose the following restriction on the time step:

∆t max(|u|, |v|) ≤ min(∆x, ∆y) 4 . ( 66
)
Note that this condition also ensures that all the intermediate volumes constructed during either the Lagrangian phase V n+1,L i,j

(24) or during the projection phase V n+1, * i,j

(27) are positives. We remind the classical CFL condition for the Euler equation:

∆t max(|u|, |v|, |c|) < min(∆x, ∆y) 2 , ( 67
)
where c is the speed of sound of the flow that we choose constant (taken at the celerity of sound of the water c 0 ℓ in all our simulations). In our case, due to the quasi-incompressible nature of the fluid, the Mach number has to be always less than 0.1 (in time and space): [START_REF] Monaghan | Simulating free surface flows with sph[END_REF] for instance). In particular, thanks to the last relation, the numerical speed of sound checks:

M = |u| c < 0.1 (cf.
c > 2 max(|u|, |v|), (68)
and thus the classical CFL condition (67) completed with the condition (68) deduced from the low Mach hypothesis is efficient to ensure the condition (66). Thus the condition (64) (or (55) in the case of positive edge velocities) which is necessary to define a nonempty trust interval for the lowdiffusive choice of α LD i+1/2,j is checked: the low diffusive procedure does not conduct to a reduction of the time step.

Numerical experiments

The resulting two-dimensional numerical code (ODYSSEY) has been tested on various dam break test cases, selected since they allow us for comparisons with other numerical codes and, when available, with real water experiments realized in tank. Only wall slip boundary conditions are currently implemented in the ODYSSEY code. Extension to more complex boundary conditions will soon be realized.

Test case I. We use the parameters from [START_REF] Cruchaga | Collapse of a liquid column: Numerical simulation and experimental validation[END_REF] enabling a comparison with both real experiments and numerical results (obtained with a Finite Element method coupled with an interface capturing method for the Navier-Stokes incompressible equations). We consider a numerical box (0.44 m × 0.42 m) filled up with gas except for a column of water at the left bottom corner of width 0.144 m and twice higher. We assume that the gate which retains the water instantaneously disappears at the initial time. The parameters for the EOS are: γ g = 1.4, γ ℓ = 7, c 0 = 350 m.s -1 , p 0 = 10 5 Pa, ρ 0 g = 1 kg.m -3 , ρ 0 ℓ = 1000 kg.m -3 ; the number of cells in each direction is N x = N y = 300. The results are presented in Fig. 3 and Fig. 4 for different times. Note that the experiments are realized on three-dimensional tanks which can lead to transverse effects that cannot be rendered by a two-dimensional numerical code. On Fig. 3, we observe that except at the initial time where we see that the rise of the gate has a little impact on the experiment (not taken into account by the codes), the collapse of the column with our code gives results in very good agreements with those of Cruchaga et al.. Up to the time t = 0.6 s, the numerical results are indeed very similar. In the next snapshots, discrepancies are mainly due to the droplets which fall from the upper wall (in particular at the left bottom of the box), which is not the case in the simulation of Cruchaga et al. since they have taken open boundary conditions for the upper wall. On Fig. 4: at these longer times (especially the three last shots), our code reproduces better the experiment than the numerical code of Cruchaga et al. (we just point out that we use a more refined mesh (300x300) than the one used in Cruchaga et al. (2007) (100x75)). The evolution of the water in time is well reproduced by our code, the agitation we still observe is due to the fact that our model does not take into account viscosity.

Test case II. In this numerical experiment, we add an obstacle in the middle of the box. This conducts to the formation of a long wave which impacts the right wall. Using sizes from Greaves (2006) (Fig. 5) we can compare with both experiments [START_REF] Koshizuka | A particle method for incompressible viscous flow with fluid fragmentation[END_REF]] and numerical simulations [START_REF] Greaves | Simulation of viscous water column collapse using adapting hierarchical grids[END_REF]] at different times. The EOS parameters are identical to the previous test case. The number of points of discretization are N x = 600 and N y = 600. Comparative results of the volume fraction at different times are shown at Fig. 6 and Fig. 7 7 . On Fig. 6, small deviations are observed at the beginning between numerical and physical experiments originate from the finite time necessary for the removal of the gate. Otherwise, it is essentially in good agreement with referred results. On Fig. 7, the results are in agreement with the experiments and the other simulation. In particular, in the last snapshots (t = 0.645 s), our code developed the same patterns near the obstacle and a gas pocket on the top right corner, similarly to the experiment. Then, in Fig. 8, we compare results obtained with our 2D code with the first order projection with (case B) and without (case A) the low-diffusive procedure. We plot the volume fraction of gas at three different times. The first order low-diffusive projection allows us to keep a thin interface between water and air. We get a detailed interface and physical phenomenons including pockets of gas, ejection of droplets, etc.

Besides, due to the high density ratio (1/1000) between both liquid and gas phases, gas mass fraction c g is a better indicator of numerical diffusion. Indeed, if the volume fraction of gas α is α = 1 -10 -8 , the mass fraction c g is c g ∼ 1 -10 -5 (taken ρ g ∼ 1 and ρ l ∼ 1000). In Fig. 9, we present this quantity at different time steps in the whole box of simulation (without the cutoff at h = 0.6 m). After the formation of a wrapped motion, we observe some numerical diffusion due to the ejection of droplets initiated from the edge of the wave. Outside of this zone, the diffusion of the interface of the free surface flow is restricted to one cell. At the final time (t = 0.645 s), these small amounts of liquid begin to slightly disturb the free surface on the bottom right due to the fall of droplets due to gravity. But this diffusion appears only in terms of mass and not in term of volume of water which remains negligible. In short, droplets (which are initially ejected from the wave edge) generate a small volume fraction of water in mixed cells. Since our code considers an unique mean velocity u, we then keep all created mixed cells within the numerical simulation. Snapshots of the conservative quantities (α g ρ g ) et (α ℓ ρ ℓ ) [kg.m -3 ] in log 10 scale are displayed on Figures 10 and11.

In particular, one can observe a "numerical spray" developing through the domain, mostly due to the appearance and fragmentation of fine structures (filaments, high-frequency instabilities). This is a shortcoming of the method but that doesn't avoid us to track the free boundary with sufficient accuracy since the fluid interface is kept sharp. In Fig. 12 we superimpose both gas volume fraction and velocity field, thanks to the representation of this vector field, we better understand the motion of the free surface flow. Finally in Fig. 13 the pressure field is plotted, in particular at t = 0.258 s and t = 0.387 s, we see a pressure peak due to the obstacle at the middle of the box.

Test case III. We here consider the sloshing of a liquid under a horizontal excitation, i.e. the tank translates horizontally (surge acceleration). The tank position moves according to the periodic motion x(t) = A cos(2πt/T ) where A = 0.032 m is the maximum amplitude of the excitation and T the period (cf. Fig. 14). In [START_REF] Akyildiz | Sloshing in a three-dimensional rectangular tank: Numerical simulation and experimental validation[END_REF], it is analyzed that the quantity of energy is expected. The first natural frequency of the liquid contained in a tank could be obtained thanks to the formula (cf. e.g. [START_REF] Akyildiz | Sloshing in a three-dimensional rectangular tank: Numerical simulation and experimental validation[END_REF]):

ω fluid = g π L tanh( π L h w ), ( 69
)
where L is the length of the tank, h w the initial height of the water in the tank. In the following, we take a forced period of T = 1.3 s (i.e. ω forced = 4.83 rad/s) and a height of water of 0.6 m. The relation ( 69) gives a natural frequency for the fluid of ω fluid = 3.77 rad/s. For this test case, both experimental results performed by [START_REF] Faltinsen | Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth[END_REF] and numerical one are available [START_REF] Shao | An improved SPH method for modeling liquid sloshing dynamics[END_REF]. In Fig. 15 we present the gas volume fraction at different times: due to the horizontal excitation of the tank, a wave initiates and displaces from left to right and right to left with an amplitude which increase up to t ≃ 3.5 s before decreasing and having a situation more stable at t = 6.54 s. After, the amplitude of the oscillation increases again (cf. Fig. 16).

In the next figure 16, we follow the time evolution of the free surface elevation of the water at a distance of 0.05 m of the left of the box and we superpose to our results the scanned experimental results of [START_REF] Faltinsen | Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth[END_REF]. We obtained both the correct period for the oscillation and the correct wrapping of the signal. The slight overestimation of the amplitude of the height of the water can be explained by the fact that we have neglected the viscosity in our model. Without viscosity, the amplitude of the wave is not slowed down and the results are quite overestimated compared to the ones obtained in Fig. 5 a) of [START_REF] Shao | An improved SPH method for modeling liquid sloshing dynamics[END_REF] with their SPH code with viscosity. Anyway, we recover the expected nonlinear characteristics: the upwind sloshing amplitude is larger than the downwind sloshing amplitude. By performing a fit of the curve given by our code (in plain blue line in the Fig. 16), we recover the curve as a superposition of the forced and natural frequencies (Fig. 17).

Test case IV. Here, we study the free fall of a disk of water (surrounded by air) and its impact with some water at rest (Fig. 18). We show the evolution of the volume gas fraction at different time steps (cf. Figs. 19 and 20). On Fig. 19, the smashing of the sphere conducts to the formation of two jets on either side of the initial disk. Then, on Fig. 20, we see that the gap generated by the fall of the disk is progressively filled up due to the falling of the jet of water at right and left of the box due to the gravity. This motion conducts at the end of the presented snapshots to the formation of a new bouncing jet at the middle of the box due to the meeting of the two lateral ones. The expected symmetry is well-conserved during quite a long time. What can be observed is that the computational quality and accuracy are almost as good as pure Lagrangian methods (SPH, particle-based), at least for the characteristic time of strong dynamics. A) is compared to the (also first order but) low-diffusive procedure (case B). It is observed that numerical diffusion is mainly suppressed and complex interface details can be tracked.

Test case V. As a complementary test case, we show the long-time development of air-water Rayleigh-Taylor instabilities generated by gravity. At initial time, the liquid is put over the gas with a sine-shaped boundary between the two fluids (the boundary is given by Γ(x) = Ly/2 + Ly/20 cos(2πx/Lx), cf. Fig. 21). Wall boundary conditions are enforce in all the box's sides. In Figs. 22 and 23, we observe the gas volume fraction α. Because the surface tension is not taken into account, instabilities at all wavelengths should emerge and grow. From the numerical point of view, only wavelengths greater than the mesh size can be captured by the code. We are aware that this computation is mesh-dependent but its goal is to show that we are able to capture free-boundary air-water flows with strong topological changes with a small amount of numerical diffusion. At the first times of the simulation (t = 0.07 s and t = 0.09 s), we observe a lot of secondary instabilities. This is due to the discretization of the initial cosine shape of the perturbation of the interface (they are clearly mesh-dependent). These instabilities should be reguralized by taking into account viscous or surface tension effects. Moreover, we emphasize that we do not observe the classical mushroom shape during our simulations. Indeed, due to the instabilities develppment's characteristic time τ ∼ 0.1 s in our simulation (τ = 1 Agk , where A = ρ l -ρg ρ l +ρg ∼ 1 is the Atwood number, g = 9.81m.s -2 the gravity and k = 2π λ , where λ = L x = 0.5 m is the wavelength), we only are at 3τ when the pike of water reaches the bottom of the box. We are thus closer to a linear regime than to a non linear one during the simulation. On the contrary, on the right and left sides, we reach faster the non linear regime since the wavelengths of the secondary instabilities are smaller and thus the characteristic time too (see for instance at t = 0.30 s). We point out that we are also able to reproduce Kelvin-Helmholtz instabilities due to the friction between the two fluids arising when the heavy fluid sink towards the bottom. This is at the source of the vortex that we observe for instance at t = 0.30 s or t = 0.53 s on the boundary of the main pike.

Discussion and conclusions

In this paper we have presented a numerical solver for immiscible fluid two-phase flow problems with numerical solution and experiments for air-water flows. The hydrocode solver is based on a remapped Lagrange discretization which appears to have several advantages in this context. After a pure conservative finite volume reinterpretation, the strategy for free boundary capturing is to design anti-diffusive phasic mass fluxes. We adapt the antidiffusive approach by Després-Lagoutière and Lagoutière-Kokh to our system of equations. We derive accuracy-stability trust regions to select the best interface gas fractions ("best" means most compressive but stable). Roughly speaking, what can be observed is that the computational quality and accuracy are almost as good as pure Lagrangian methods (SPH, particle-based) as soon as details of the moving interface are not of the order of the mesh step. In case of fine structures like filamentation, sprays, high-frequency instabilities, the numerical method inherently creates artificial phase mixing and numerical diffusion of the interface. A way to improve accuracy would be to use adaptive mesh refinement AMR strategies in regions of strong gas fraction gradients, but this has not be done yet. This works is a milestone toward a more physical air-water simulation code including fluid viscosity and free boundary surface tension and more "multiphase effects". The use of two phasic velocities for example would allow for fluid sliding at interfaces, but also for phase reseparation (by buoyancy). Viscosity and surface tension should be helpful and good for interface capturing accuracy because they have interface regularizing effects. Another milestone is the GPU code parallelization. The antidiffusive interface capturing scheme has been held up because we are confident on its natural and straightforward parallelization without particular specific interface treatment. Figure 16: Test case III. Tank partially filled up with water submitted to a horizontal excitation, with h w = 0.6 m and a period for the tank of T = 1.3 s. We show the free surface elevation, i.e. the height of the water relatively to the initial height h w . We superimpose our results (plain blue line) over the scanned experimental curve of [START_REF] Faltinsen | Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth[END_REF] (dotted black line). Our code reproduces with good accuracy the frequency of the oscillations of the height of the water and the envelop of these oscillations is well conserved. We notice small defects in the amplitude of the relative height of the water compared to the experiment. In particular, in the previous formula (A.18), α n+1,L i±1/2,j = I(α) |V i±1/2,j represents the interpolation of α on the volume V i±1/2,j defined by the displacement of the edge (i±1/2, j) during the time step ∆t at the velocity u n+1/2,L i±1/2,j . In Eq. ( A.18), α n+1,L i±1/2,j is given by the low-diffusive value α LD i±1/2,j calcultated in Sec. 4. Then, starting from this initial point, the Newton algorithm modifies the partial volume of each phase up to reach the equilibrium in pressure. We ensure the guess to be in the interval [0, 1] by defining α guess,remap = I(α) n+1, * i,j

)). For the second step of projection, the formula to initiate the Newton algorithm is similar: the interpolation of α on the Eulerian volume V n+1 i,j is given by: on the volumes V i,j±1/2 defined by the displacement of the vertical edges i, j ± 1/2 at the velocities v n+1/2,L i,j±1/2 during the time step ∆t. It corresponds to the low-diffusive value α LD i±1/2,j calculated during the second step of remap to get the projected masses thanks to Eqs. ( 33)-(34).

I(α) n+1 i,j = V n+1

AppendixB. Comments on the artificial viscosity

During the Lagrange phase (cf. Sec. 3.1), we use an artificial viscosity q (the so-called pseudoviscosity) in order to stabilize the staggered scheme which is centered in space. The pseudoviscosity q is a viscous pressure. It is designed to only act into compression zones (div(u) < 0) and shock waves. On the contrary, for smooth expansion zones the pseudo-viscosity is set to zero to keep second order accuracy. This conducts to change the momentum balance by means of the pressure gradient but also the mean density by means of the divergence of the velocity. The pseudoviscosity is a combination of a linear term (q lin ≈ div(u)) which acts on the linear stability and a quadratic one q quad ≈ |div(u)|div(u) provides nonlinear (large-amplitude shock) stability. We refer to [START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF]; [START_REF] Wilkins | Use of artificial viscosity in multidimensional fluid dynamic calculations[END_REF]; Caramana et al. (1998a); [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]; [START_REF] Robinson | ALEGRA: An Arbitrary Lagrangian-Eulerian multimaterial, multiphysics code[END_REF]; [START_REF] Sprague | The numerical treatment of simple hydrodynamic shocks using the von Neumann-Richtmyer method[END_REF]] for more details on this subject.

As we deal with the two-dimensional problems, we have to take into account a 2D approach of the velocity divergence and in order to check the previous requirements, we choose the following form of the pseudo-viscosity in the code: the linear term is expressed as q lin = -a 1 ρcs √ ∆y∆x (∆u∆y + ∆v∆x), if div(u) < 0 (B.1) = 0 elsewhere, where ∆u represents the jump in velocity along the x direction: ∆u i,j = u n i+1/2,j -u n i-1/2,j , and ∆v is the jump along the y direction: ∆v i,j = v n i,j+1/2 -v n i,j-1/2 , naturally defined in the center of a cell (i, j) thanks to the staggered grid (as the real pressure term p). And then, we get: (divu) ij = ∆u i,j ∆x + ∆v i,j ∆y . The pseudo-viscosity coefficient a 1 > 0, a 1 = O(1) is a constant to be defined depends on hydrodynamics quantities. For the quadratic term, it has to be on the form of q quad ≈ -a 2 ρ V D t V |D t V | (cf. [START_REF] Sprague | The numerical treatment of simple hydrodynamic shocks using the von Neumann-Richtmyer method[END_REF]) and as the total volume V = m ρ follows the equation D t V = V ∇ • u in the Lagrangian phase cf. Eq. ( 15), we choose the following form: q quad = -a 2 ρ ∆x∆y |∆u∆y + ∆v∆x|(∆u∆y + ∆v∆x) if div(u) < 0 (B.2) = 0 elsewhere, which has the dimension of a pressure (a 2 = O(1)). The values of the two a 1 and a 2 are put to 0.15 and 0.1 respectively in all the simulations.

Figure 2 :

 2 Figure2: Effect of upwind projection without any specific treatment of the volume fraction, due to numerical diffusion. Example of collapse of a column of water with an obstacle in the middle of the box. We represent the numerical volume fraction α at a certain time obtained with two different upwind projections. With the first order upwind projection (Fig.2a), the interface between the liquid (in blue) and the air (in red) is too much diffusive to be followed with accuracy (e.g. droplets and air pockets are hidden because of the diffusion). With the second order upwind projection (Fig.2b), the interface is less diffuse and we begin to see an air pocket which were previously concealed by the diffusion. But the diffusion is still present and we do not observe formation of droplets. For details on the initialization, please refer to the Test case II of the section 5 dedicated to the numerical experiments.

Figure 3 :

 3 Figure 3: Test case I. Collapse of a column of water, evolution of the volume fraction α at different times. The experimental (resp. numerical) results of Cruchaga et al. (2007) are at left (resp. in the middle) and the results with our code ODYSSEY at right are compared at different time step.

Figure 4 :Figure 5 :

 45 Figure 4: Test case I. Same as Fig. 3 at longer time.

Figure 6 :

 6 Figure 6: Test case II. Collapse of a water column with an obstacle. At left, the results obtained with the code of Greaves, in the middle, the results of Koshizuka et al. and at right the results with our code ODYSSEY (we observe the volume fraction α).

Figure 7 :

 7 Figure 7: Test case II. Collapse of a water column with an obstacle at longer time.

Figure 8 :

 8 Figure 8: Test case II. Accuracy of interface capturing: we present the gas volume fraction α. A first order projection (caseA) is compared to the (also first order but) low-diffusive procedure (case B). It is observed that numerical diffusion is mainly suppressed and complex interface details can be tracked.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Test case II. Collapse of a water column with an obstacle. Snapshots of the gas mass fraction c g (4). We represent the whole numerical box (without the cut at h = 0.6 m).

Figure 12 :

 12 Figure 12: Test case II. Dam Break test case with an obstacle with the ODYSSEY code. We superimpose both gas volume fraction and velocity field.

Figure 13 :Figure 14 :

 1314 Figure 13: Dam Break test case with an obstacle with the ODYSSEY code. We show the pressure field.

Figure 15 :

 15 Figure15: Test case III. Tank partially filled up with water submitted to a horizontal excitation, with h w = 0.6 m and a period for the tank of T = 1.3 s. We present the volume fraction of the gas α (always gas in red and water in blue) at different times. A wave generated by the horizontal motion of the tank moves on both sides of the tank with a varying amplitude.

Figure 17 :Figure 18 :

 1718 Figure17: Test case III. Fit of the curve representing the time evolution of the free surface elevation of the water given by our code with a function f (t): f (t) = A sin(ω 1 t + φ 1 ) + B sin(ω 2 t + φ 2 ) + off (curve given by our code in dotted red line and fitted curve in plain blue line). We get for the two frequencies: ω 1 = 3.74 rad/s and ω 2 = 4.83 rad/s (and φ 1 = 1.55 rad/s, φ 2 = 1.59 rad/s, A = -0.08, B = 0.07, off = 0.01). The first one is very closed to the natural frequency of the fluid (ω fluid = 3.77 rad/s) and the second to the induced frequency of the box (ω forced = 2π T = 4.83 rad/s).

Figure 19 :Figure 20 :

 1920 Figure 19: Test case IV. Free fall of a disk of water and impact with water at rest. Evolution of the volume fraction of gas α at different time steps.

Figure 21 :

 21 Figure 21: Test case V. Configuration: Lx = 0.5 m, Ly = 1. m, N x = 300, N y = 600. The border between air and liquid is disturbed at a height of Ly/2 by a sine-shaped signal: Ly/2 + Ly/20 cos(2πx/Lx).

Figure 22 :

 22 Figure 22: Test case V. Rayleigh-Taylor instabilities. Heavy (water) fluid above the light one (gas) with a sine-shaped perturbation of the initial interface, evolution of the gas volume fraction α.

Figure 23 :

 23 Figure 23: Test case V. Rayleigh-Taylor instabilities. Development of a main needle followed by secondary fine needles and topology changes. The overall symmetry of the computational is globally fulfilled.

  In the formula(A.19), the interface values α n+1, * i,j±1/2 = I(α) |V i,j±1/2 represent the interpolations of the quantity α n+1, *

	, * V n+1 i,j i,j	α n+1, * i,j	-	∆t ∆x	v	n+1/2,L i,j+1/2 α n+1, * i,j+1/2 -v	n+1/2,L i,j-1/2 α n+1, * i,j-1/2 .	(A.19)
	i,j

When u n+1/2,L i+1/2,j < 0, and if u n+1/2,L i+3/2,j < 0, the stability condition must be calculated in the cell i + 1 since the information comes from the right cell.

Note that the arbitrary choice to calculate maximum principle conditions on the quantities related to the gas induces the same properties for the liquid due to the relations cg + c ℓ = 1, and thus min(cg i,j , cg i-1,j ) = 1max(c ℓ i,j , c ℓ i-1,j ) and max(cg i,j , cg i-1,j ) = 1 -min(c ℓ i,j , c ℓ i-1,j ).

The tank used by[START_REF] Koshizuka | A particle method for incompressible viscous flow with fluid fragmentation[END_REF] is opened at a height of 0.6 m whereas the numerical simulations are performed on a box of height of 4a = 1 m.

We recall that for a generic variable z, z i+1/2,j,up corresponds to the value zi,j if u n+1/2,L i+1/2,j > 0 and otherwise to the value zi+1,j.

Acknowledgements

This work is partly supported by the joint Lab LRC MESO between CMLA and the French Nuclear Agency CEA, DAM DIF. We would like to thank J.-M. Ghidaglia for fruitful advices and suggestions and L. Desvillettes for his support. A.B.C. acknowledges partial support of CBDif-Fr, Collective behaviour & diffusion : mathematical models and simulations ANR-08-BLAN-0333-01.

AppendixA. Numerical solution of the pressure equilibrium equation

AppendixA.1. Solution of pressure equilibrium equations and equations of states

Local pressure equilibrium conditions usually lead to a local algebraic system to solve, and one has to take care of the existence and uniqueness of a solution of the system, but also on the numerical solution using a fixed point Newton-type method (see the interesting work from [START_REF] Flatten | On solutions to equilibrium problems for systems of stiffened gases[END_REF] on the case of a mixture of stiffened gases). In this Appendix, we discuss the case of pressure equilibrium using an isentropic perfect gas law for gas and a modified Tait equation for water.

From the conservative variables W g = αρ g and W ℓ = (1 -α)ρ ℓ , we have to compute both pressure p and gas volume fraction α. The pressure equilibrium assumption (barotropic closure) leads to a scalar algebraic equation to solve p = p g (ρ g ) = p ℓ (ρ ℓ ), i.e.

for α ∈ (0, 1). As an example, let us consider an isentropic perfect gas law for the gas

and the modified Tait equation for the liquid

where K = ρ 0 ℓ c 2 ℓ p 0 γ ℓ , c ℓ being the speed of sound in the liquid. The modified Tait equation is a rather good, local approximation of the compressibility of the water near references conditions (ρ 0 ℓ , p 0 ). It is assumed that flow conditions will not exceed the validity domain of this EOS, i.e. p ℓ (ρ ℓ ) > 0 or equivalently

Under "atmospheric" conditions, one can use the following numerical parameters: p 0 = 10 5 P a, ρ 0 g = 1.28 kg m -3 , ρ 0 ℓ = 1000 kg m -3 , c ℓ = 1500m s -1 , γ g = 1.4, and γ ℓ = 7. For these values, K = 3214.3 and the validity domain of the Tait equation (A.4) is

(showing the low compressibility of the water). With these EOS, the pressure equilibrium equation (A.1) writes ρ g ρ 0

Introducing the dimensionless variables µ g = Wg ρ 0 g and µ ℓ = W ℓ ρ 0 ℓ , we get the algebraic equation to solve

AppendixA.2. Uniqueness of the root

For mixed cell, i.e. µ g , µ l > 0 fixed, we search the zero of the following function Φ:

By derivating, we get:

(we are looking for values of α ∈]0, 1] since we consider mixed cells). As φ(α, µ g , µ ℓ ) -→

Moreover this solution ensures to get positive pressures in the mixed cells. Practically we use volume fraction thresholds by defining two critical values α g c = 10 -8 and α ℓ c = 10 -11 : when α * < α g c , then α * is forced to 0 and when 1-α * < α ℓ c , then α * is forced to 1 to avoid to deal with too small numbers (α * such that φ(α * , µ g , µ ℓ ) = 0). Note that these cuts conduct to neglect mass fraction of gas less than c g = αρg ρ < c cutof f g ≈ 10 -11 and mass fraction of the liquid less than

(for ρ g ≈ 1 and ρ ℓ ≈ 1000). Notice that the cutoff is chosen not symmetrically due to the high density ratio between the two phases (ρ ℓ /ρ g ∼ 1000). In such a way, this cutoff leads to the same minimal mass of each material in the cell.

AppendixA.3. Efficient Picard fixed point method strategy

We can rewrite (A.6) in a more appropriate form

. In order to solve f (α) = 0 numerically, we need a convergent fixed point algorithm. Usually, a Newton method is used because of its quadratic convergence rate. Unfortunately, we experienced a poor convergence rate on the equation (A.8). The reason is that α → f (α) has a derivative at the root which is (numerically) close to zero and thus the convergence rate is quasi-linear. On figure A.24, we show a typical profile of the function α → f (α). One can observe a "flat" region near the root, making the Newton method slowly convergent. Our strategy is to correctly initialize the Newton method with some "good" initial guess by means of a predictor step.

Initial guess strategy for the Newton algorithm

Because of the stiffness of the pressure equilibrium equation, a convenient initial guess for Newton's algorithm is needed to ensure both fast convergence and for the iterates to stay into the admissible interval [0, 1]. Below, we describe two different ways (depending on the step of the numerical scheme considered, i.e. Lagrange or remap) in order to find a systematic initial guess for Newton's algorithm.

AppendixA.4.1. Lagrangian step

During the Lagrangian step, we look for an equation checked by the volume gas fraction α which includes the pressure equilibrium assumption. This has been already discussed in papers dealing with relaxation strategies (cf. [START_REF] Colella | Multifluid algorithms for Eulerian finite difference methods[END_REF]; [START_REF] Robinson | ALEGRA: An Arbitrary Lagrangian-Eulerian multimaterial, multiphysics code[END_REF] for instance). In the following, we will denote α g the volume fraction of gas and α ℓ the one of the liquid, with α g + α ℓ = 1. From the EOSs of the two phases ( 6)-( 7), we obtain:

(A.9)

where p is the pressure of the mixture of gas and liquid at the equilibrium: p = p g (ρ g ) = p ℓ (ρ ℓ ). Moreover, from the equations of conservation of masses of each phase ( 9)-( 10), we have:

Using (A.9)-(A.10) in (A.11), we get:

Since α g + α ℓ = 1, we find the non-conservative equation check by the pressure p in the mixed cells:

depending on the coefficients of each EOS and on the volume gas fraction α. And thus using (A.12), we obtain the equation followed by the volume gas fraction α including the pressure equibirum assumption:

Thanks to (A.15), we get (using that the total volume checked D t V = V ∇ • u thanks to Eq. 15):

and we finally take for the initial value in the Newton algorithm during the Lagrangian phase:

AppendixC. Definition of the interval I 2 to ensure stability conditions for c g

In the section 4, we have treated the case when the velocities u n+1/2,L i±1/2,j > are positive and the stability condition for the remap scheme on c g is calculated on the cell i, j. Without detailing the other case when u n+1/2,L i+1/2,j < 0, (and if the velocity of the other vertical edge is negative u n+1/2,L i+3/2,j < 0), by defining s = sign(u n+1/2,L i+1/2,j ) we can write the interval

] in a generic form . The lower bound is expressed as:

and the upper bound as: Thus, under the generic restriction condition on the time step:

i+1/2,j < 0 and u n+1/2,L i+3/2,j < 0. We remind here that in the case of two velocities of different sign in the cell where the stability is calculated, an upwind choice is performed for α at the edges, which always ensure the stability requirement (see the Remark after Theorem 4.2).

AppendixD. Reminder of the choice of the low-diffusive value at the edges α LD i+1/2,j

In the section 4, we have define the trust interval in which α LD i+1/2,j has to be taken. We just recall here the procedure of low-diffusion that can be found in Kokh and Lagoutière (2010) using our trust interval I (65) to choose α LD i+1/2,j as close to the downwind value as possible in order to limit the diffusion:

• If u n+1/2,L i+1/2,j > the objective is to take a value for α at the edges as close as possible to the downwind one. Stability of α and thus positivity of the masses of each cell must be ensured in the cell i + 1, j in that case:

if ω n+1,L i+1/2,j < α n+1,L i+1,j < Ω n+1,L i+1/2,j , α LD i+1/2,j = n+1,L i+1 if α n+1,L i+1

≥ Ω n+1,L i+1/2,j , α LD i+1/2,j = Ω n+1,L i+1/2,j (D.1) 2. if u n+1/2,L i-1/2,j ≤ 0, we take the upwind value α LD i+1/2,j = α n+1,L i,j .

• if u n+1/2,L i+1/2,j < 0, the objective is to take a value for α at the edges as close as possible to the downwind one. Stability of α and thus positivity of the masses of each cell must be ensured in the cell i, j:

≤ ω n+1,L i+1/2,j , α LD i+1/2,j = ω n+1,L i+1/2,j if ω n+1,L i+1/2,j < α n+1,L i,j

< Ω n+1,L i+1/2,j , α LD i+1/2,j = α n+1,L i,j if α n+1,L i,j

≥ Ω n+1,L i+1/2,j , α LD i+1/2,j = Ω n+1,L i+1/2,j (D.2) 2. if u n+1/2,L i+3/2,j ≥ 0, we take the upwind value: α LD i+1/2,j = α n+1,L i+1,j .