Florian Desmons
email: florian.desmons@u-bordeaux.fr

Mathieu Coquerelle
email: mathieu.coquerelle@bordeaux-inp.fr

A generalized high-order momentum preserving (HOMP) method in the one-fluid model for incompressible two phase flows with high density ratio

Keywords: Navier-Stokes, two phase flows, numerical method, consistent transport, momentum, high-order method Highlights

A generalized high-order momentum preserving (HOMP) method in the one-fluid model for incompressible two phase flows with high density ratio

Florian Desmons, Mathieu Coquerelle

Introduction

Numerous applications in fluid dynamics involve the interaction of two incompressible (or almost incompressible) phases with a high density ratio around 10 3 and viscosity ratio around 10 2 , such as quotidian air and water. Water wave breaking, droplets impacts, bubbles trajectories, thin films formation, tank sloshing, etc. are common phenomena that are widely studied with CFD codes. May it be at micrometer or oceanic scale, inertia usually plays a dominant role, particularly for flows at high Froude, Reynolds or Weber numbers. The importance of the precise simulation of the mechanical interactions between liquid and gases prompts for well-suited models and accurate numerical methods.

The physical transition between two immiscible phases takes place at the nanometer scale, relying out of scope of the Navier-Stokes (NS) equations. Hence, for CFD simulation, the common hypothesis is to consider it as infinitely small. The NS equations can be solved separately in each phase with appropriate boundary conditions -namely, at the interface, no transfer of mass, the pressure jump due to surface tension and the viscous stress continuity -, given that the surface between them is precisely tracked. A very common alternative, in a Eulerian framework is usually adopted: herein the interface is implicitly captured by a volume fraction and the NS equations are solved in the whole domain. Under this framework, two principal families of methods exist: the Ghost Fluid Method (GFM) and the One-Fluid Method (OFM, also known as the whole-domain formulation), originating principally from Kataoka [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF]. The former is based on the injection of the boundary conditions into the discrete schemes (i.e. for incompressibility, the convection and the viscous terms); it requires a precise location of the interface and appropriate numerical methods, usually limited to second order. On the other hand, the latter relies on a regularized (i.e. smooth) transition of the density and viscosity at the mesh level (e.g. the density climbs from 1 to 1000 • kg • m -3 in approximately 6 cells around the interface). In that case, classical numerical schemes to compute the momentum equation, and in particular the critical convection term, have been widely used. While a first order scheme is very diffusive and will quickly damp momentum, it is noticeably very stable; conversely high-order schemes are much less diffusive, better capture small vortices but might induce numerical instabilities. However, as the mesh is refined and the gradient of density is going to zero, deriving the momentum across the interface, if no particular care is taken, will eventually lead to severe numerical errors and thus unstable simulations.

Additionally, the OFM is based on an underlying interface method that is used to ensure a controlled (i.e. fixed) interface thickness. Several methods exist in the literature and are available in CFD codes; we focus here on the principle ones. For that purpose, Volume of Fluid (VOF) has been introduced by Hirst and Nichols [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF], enhanced with the geometric approach Piecewise Linear Interface Calculation (PLIC) by Rider and Kithe [START_REF] Rider | Reconstructing Volume Tracking[END_REF]. One disadvantage of that method the fact that it is only first order accurate. More recently, the Moment Of Fluid (MOF) approach has been introduced by Dyadechko and Shashkov [START_REF] Dyadechko | Moment-of-fluid interface reconstruction[END_REF] for its increased accuracy (up to second order) and further developed in [START_REF] Jemison | Filament capturing with the Multimaterial Moment-of-Fluid method[END_REF][START_REF] Milcent | Moment-of-fluid analytic reconstruction on 3D rectangular hexahedrons[END_REF] for sub-cell filaments and analytical optimization. Finally, Level Set Methods (LSM), originating from Sethian and Osher [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], are widely used for their inherent smooth property making them particularly well fitted for surface tension computation (see Coquerelle et al. [START_REF] Coquerelle | A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces[END_REF]). LSM however is known to suffer from mass loss and deformation: several techniques are employed to counteract those problems, among others re-initialization algorithms [START_REF] Russo | A Remark on Computing Distance Functions[END_REF][START_REF] Hartmann | The constrained reinitialization equation for level set methods[END_REF]. A detailed review of LSM has been conducted by Gibou et al. in [START_REF] Gibou | A review of level-set methods and some recent applications[END_REF].

The coupling of the continuity and the momentum equations with the interface representation leads to a dilemma where, when solved at the discrete level, momentum can be lost and even worse transferred from one phase to another. For instance, because of numerical errors -e.g. approximations, truncature terms, decoupling, linearization, inconsistency, etc. -, when the momentum of the heaviest phase (i.e. the liquid) is transferred towards the lighter one (i.e. the gas), spurious high velocities arise in the latter; conversely when the lighter fluid transfers momentum towards the heaviest one, even though the velocity might rise in a small manner, the total energy is increased drastically. While these phenomena were not particularly detected for years, these numerical errors have been observed to be increasing as the mesh is refined as the growth of computational power allows the simulation up to billions of cells, thus bringing a new challenge to tackle. In consequence, capturing accurately and with stability the dynamic of two-phase flows at large and small scales simultaneously has become a real challenge which requires the development of adequate and robust methods.

Recently, multiple works have addressed with this problem. In a VOF framework, following the precursor works of Rudman [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF], then Bussman et al. [START_REF] Bussmann | Modeling High Density Ratio Incompressible Interfacial Flows[END_REF], Fuster et Popinet [START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF] have recently proposed a method that is adapted for general compressible flows but that is limited to first order. Owkes et Desjardin [START_REF] Owkes | A mass and momentum conserving unsplit semi-Lagrangian framework for simulating multiphase flows[END_REF] on their side have developed a semi-Lagrangian scheme adapted coupled with VOF. In a MOF framework, Jemison et al. [START_REF] Jemison | Compressible, multiphase semi-implicit method with moment of fluid interface representation[END_REF] have provided the CISL (Cell Integrated Semi-Lagrangian) numerical scheme for compressible flows. Both methods are based on a PLIC reconstruction of the interface. In the LSM framework, several works [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Olsson | A conservative level set method for two phase flow[END_REF] have been conducted in order to synchronize the momentum and mass flux permitting to diminish the problems inherent to the spreading of the LSM interface approach. More recent articles [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF][START_REF] Manik | A generic algorithm for three-dimensional multiphase flows on unstructured meshes[END_REF][START_REF] Kumar | A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids[END_REF][START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF][START_REF] Ghods | A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods[END_REF] have focused on a more generic formulation, i.e. designed independently of the underlying interface transport, based on the advection of an auxiliary smooth density field. Manik et al. [START_REF] Manik | A generic algorithm for three-dimensional multiphase flows on unstructured meshes[END_REF] have addressed the problem for unstructured meshes, based on an iterative approach for matching fluxes. The latest work of Nangia et al. p. 2 [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF], in an approach similar to ours, takes advantage of a fully coupled computation to compute consistent fluxes, with a third order spatial scheme and an underlying LSM. All those articles point the finger on the need for a consistent mass and momentum advection approach when the density ratio is larger than 100.

Most approaches are critically linked to the interface transport/advection method and can be delicate to implement in details and/or extend to high-order schemes and 3D. At the price of less precise conservation of momentum, more generic approaches ensure the synchronized and consistent preservation of both mass and momentum.

Our approach follows that principle based on a consistent solving of the continuity and momentum equations in a smooth one-fluid Eulerian framework. It is independent on the interface method and is compatible with conventional high-order schemes such as WENO and Runge-Kutta, thus being easily adaptable to existing CFD codes in 2D and 3D. WENO schemes are particularly suited for capturing a fast varying quantity, such as density and momentum near the interface, and in consequence our approach relies on using a unique scheme for the advection equation, with no particular treatment done in the interfacial region. This strategy ensures to converge towards the expected continuous solution from a discrete point of view, at the price of a non-exact conservation of momentum in each phase, hence the named momentum preserving (MP) principle. We demonstrate over various examples that the proposed method effectively reduces transfers of momentum due to numerical errors.

We first present the context of the incompressible immiscible two-phase equations in the one-fluid model in section 2.1. Then, in section 2.2 we detail and analyze the origins of erroneous transfers of momentum between phases from the discrete point of view. The proposed method, named HOMP for High-Order Momentum Preserving, is developed and discussed in section 3.1. The underlying numerical schemes used for solving the Navier-Stokes equations and the interface transport are presented in section 3.2. Results of the proposed method are detailed and studied in section 4 wherein we have proposed verification and validation cases as well as the application to more complex two phase flows with gravity and surface tension forces. The HOMP method shows a very satisfactory reduction of erroneous momentum transfers on all interface methods, even when coupled with high-order advection schemes such as the fifth-order WENO.

Context, difficulties and objectives

In this section we firstly present the governing equations used to model two-phase flows withing the one-fluid model and secondly we present the principal numerical schemes that are used throughout this work.

Navier-Stokes equations for two-phase flows within the one-fluid model

The Navier-Stokes equations with variable density and viscosity can be written in a conservative form as:

∂ρ ∂t +∇ • (ρu) = 0 (1a) ∂ρu ∂t + ∇ • (ρu ⊗ u) = ∇ • (2µD) -∇p + f (1b)
with u = (u(x, t), v(x, t), w(x, t)) the fluid velocity, D the deformation tensor defined as D = (∇u + ∇u T )/2, p the pressure field and f external forces (such as gravity, surface tension, etc.). ρ ≡ ρ(x, t) is the fluid density and µ ≡ µ(x, t) the dynamic viscosity. For the sake of clarity, (x, t) will be dropped in the rest of the document. For an incompressible two immiscible phases flow, the physical properties ρ and µ are constant in each phase. The widely used and so-called one-fluid model [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF] consists in that the system of equations ( 1) can be solved -with the common dynamic boundary conditions at the interface -within a single fluid framework where ρ and µ are constant for each phase. The regularity (i.e. the smoothness) of both quantities is a key factor of the good convergence and stability of the approach. Moreover, when the density ratio becomes high (as it is the case for common liquid/gas applications where it is of the order 10 3 ) numerical errors increase as they are function of |∇ρ|.

In order to model the incompressible property of both phases, a characteristic function χ can be defined such that χ(x) = 1 (resp. χ(x) = 0) in phase 1 (resp. phase 2); ρ and µ are thus obtained by interpolation through the characteristic function:

p. 3 ρ = ρ 1 χ + ρ 2 (1 -χ) (2a) µ = µ 1 χ + µ 2 (1 -χ) (2b)
with ρ 1 and µ 1 (respectively ρ 2 and µ 2 ) the properties of the phase 1 (resp. phase 2). We present in section 2.2 how χ evolves with time.

Still today, as it is commonly used in the literature on the base of historical frameworks [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF][START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Scardovelli | Direct Numerical Simulation of Free-Surface and Interfacial Flow[END_REF][START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF], two-phase flows are solved with the approximation of the One-Fluid Model (OFM) that is based on the decomposition of the inertial term mixed with the continuity equation, stating that:

∂ρu ∂t + ∇ • (ρu ⊗ u) = ρ ∂u ∂t + (u • ∇) u + u ∂ρ ∂t +∇ • (ρu) = ρ ∂u ∂t + ∇ • (u ⊗ u) -(∇ • u) u + u ∂ρ ∂t +∇ • (ρu) = ρ ∂u ∂t + ∇ • (u ⊗ u) -(∇ • u) u where the identity (u • ∇) u = ∇ • (u ⊗ u) -(∇ • u) u has
been used in the second line in order to write the inertial term in a conservative form and wherein the continuity equation 1a cancels the second right-hand side term to obtain the last formula. Under the assumption of incompressibility, i.e. ∇ • u = 0, the last term also cancels and equations 1 can finally be rewritten as:

ρ ∂u ∂t + ∇ • (u ⊗ u) = ∇ • (2µD) -∇p + f ∇ • u = 0 (3)
which is algebraically equivalent to 1a in its continuous form and only in that case. The mass conservation equation was simplified as ∇ • u = 0 because ∂ρ ∂t + ∇ • (ρu) = ∂ρ ∂t + u • ∇ρ + ρ∇ • u and, in each phase, ρ is constant. These equations are easier to discretize and very commonly used in many CFD codes.

Origin of the problem

The use of equation 3 has been found to generate an important transfer of momentum through the interface due to inconsistency between continuous and discrete framework (we refer the reader particularly to [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF] for detailed discussion). By the erroneous discretization through the splitting of momentum and continuity equations, momentum fluxes appear at the interface between the two phases. Those transfers can be dreadful as, concerning air and water interaction, they impact the velocity field in a phase that is as 1000 times denser than the other. Hence, the rapid variation of density in the vicinity of the interface has to be taken care adequately both for the continuity and the momentum equations.

Trivial illustration

To depict the problem, consider the 1D problem where two Lagrangian particles are located on each side of the interface between the two phases. The first one, P 1 , in air (of density ρ 1 = 1 kg/m 3 ) moves at the velocity u 1 (t) towards the other one, P 2 , that resides in water (of density ρ 2 = 1000 kg/m 3 ) and travels at velocity u 2 (t) .

At time t+δt, discretization errors induced by the method (such as the inconsistency between the momentum/continuity equations and the discrete phase representation) or the numerical schemes (i.e. principally diffusion), will impact the movement of particle P 2 because of P 1 velocity. In other terms, the velocity

u 2 (t + δt) is impacted by u 1 in such a way that u 2 (t + δt) = u 2 (t + δt) + (u 1 , δt)
where u 2 would be the ideal solution when phases only interact through boundary conditions and (u 1 , δt) the numerical error due to the one-fluid approach. Hence, the momentum (ρu) 2 of particle P 2 is increased by ρ 2 . As the density difference between the two phases is quite high (ρ 2 /ρ 1 = 1000) even a small error in the velocity field can produce a very important error in the momentum and thus, more dreadfully, in the overall energy of the system. We call this phenomenon a discrete transfer of momentum.

In practice, we have observed the transfer of air vortices inside the water phase (and reciprocally) during the simulation of, for instance, ocean waves, as illustrated later on in figure 14. The discrete transfer of momentum was inducing non-physical under-water vortices which eventually destabilized the simulations. As the errors are rising with smaller meshes, we believe that the discrete transfer of momentum is proportional to |∇ρ| ∼ O δx -1 . p. 4

Temporal error

From an algebraic point of view, the problem arises from the fact that, for practical matters, we separate the resolution of equations 1, i.e. we do not solve them simultaneously and, as they are coupled, thus violate continuity. In consequence, when one integrates the momentum equation 1b, both velocity and density has to be known. For the sake of simplicity, we drop the viscous and external force terms such that the equation becomes:

∂ρu ∂t + ∇ • (ρu ⊗ u) = -∇p. (4)
If we integrate eq. 4 over time between t n = t and t n+1 = t + δt, we get:

t n+1 t n ∂ρu ∂t + ∇ • (ρu ⊗ u) dt = t n+1 t n -∇p dt (ρu) n+1 -(ρu) n + t n+1 t n ∇ • (ρu ⊗ u) dt = t n+1 t n -∇p dt (5)
where (ρu) n = ρ n u n is the known momentum at time t n and (ρu)

n+1
= ρ n+1 u n+1 the unknown momentum at time t n+1 to be computed. The manner one approximates the remaining integrals leads to the various well known temporal integration schemes. Under that form, the only temporal errors are due to the integration of the inertial and the pressure gradient terms. We will show in section 3.1 how, in practice, we discretize this equation.

Applying a similar approach to the continuity equation 1a leads to:

ρ n+1 -ρ n + t n+1 t n ∇ • (ρu) dt = 0 (6)
The difficulty comes from that we have two unknowns, ρ n+1 and u n+1 , and two equations containing integrals of time varying terms ρ (t) and u (t). In consequence, in order to find a solution, we have to make assumptions, usually using Taylor expansion over time. This is where the one-fluid method described earlier in its simplified form (refer to coupled equations 3) is limited to first order in time as, to solve equation 6 we need to exactly know the velocity field u (t) and for the equation 5 to exactly know the density field ρ (t) (as well as u (t)). As this cannot be done simultaneously, the continuity assumption that lead to the simplified equation is violated.

In practice, it is usually assumed that the density field ρ (t) does not vary much during time and one can take the first order approximation ρ (t) ρ n when solving equation 5. Hence, the term (ρu) n+1 is approximated by ρ n u n+1 and the inertial term by ρ n t n+1 t n ∇ • (u ⊗ u) dt . That assumption permits simplifying the problem by dividing the momentum by ρ n . As with common immiscible two-phase flows methods the density field is constant in almost the whole domain but in the vicinity of the interface, the approximation has no influence in the bulk and only where ρ varies in space and time.

However, this first order approximation limits the overall methods and particularly introduce errors near the interface where a lot of important physical interactions take place. Moreover, the one-fluid model relies on the fact that the interface thickness -i.e. the volume in which ρ (x) varies smoothly between the two phases, as exhibited by equation 2 -is as small as possible. If, ideally, this thickness -noted -should be zero, for spatial continuity matters, it has a finite support. It is commonly defined as = k • δx where δx is the local mesh size and k is a constant, usually between 3 and 6.

In consequence, if the interface advances in the normal direction with a CFL condition -i.e. the interface is moving at the rate of around one cell per time step -, the variation of density during time in that cell is non-negligible and the approximation that ρ n+1 ρ n will lead to a lot of errors. And this problem is particularly important as the ratio between the density of the two fluids is high, as for example for air and water for which it is around 10 3 . If one managed to obtain a good match between ρ * and ρ n+1 , the integration error would be much reduced. It is the case with the interface based momentum conservation techniques such as [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Jemison | Compressible, multiphase semi-implicit method with moment of fluid interface representation[END_REF][START_REF] Popinet | Numerical Models of Surface Tension[END_REF] that are, however, usually limited to low-order.

Γ n Γ n+1 Γ *
Figure 1: Illustration of the impact of the first order approximation for the density field ρ through the position of the interface Γ over time. The superscripts refer to the temporal positions of Γ, where Γ n (resp. Γ n+1 ) is the interface at time t n (resp. t n+1 ); Γ * is a good prediction (approximation) of the position of the interface at time t n+1 . The control volume is shown in gray. The first figure presents the overall shape of the fluid's interface, the second is zoomed near the gray control volume (a mesh cell for example). The three sub-figures on the right represent the density for each approximation: the opaque color stands for the heaviest fluid, the translucent color for the lighter one. The difference between the density fields at time t n+1 and the prediction through Γ * is much smaller than with the density field at time t n , thus illustrating the need for a good prediction.

Spatial error

Additionally, the spatial discretization of the equations brings more errors in the momentum equation. In this section we give a glimpse of the origin of those errors; the numerical details are beyond the scope of this article as they are very much linked to the various frameworks and advection schemes part of each CFD code.

Mean momentum definition.

First, the definition of the mean momentum inside a control volume (i.e. a mesh cell) Ω is defined as ρu = |Ω| -1 Ω ρ (x) u (x) dx where the over-line denotes the mean value. The density and the velocity can be split, inside Ω, into a mean value and a variation such as: ρ (x) = ρ + ρ (x) and u (x) = u + u (x). We can rewrite the mean momentum as:

ρu = ρ u + ρ u
because the means of the variation of the density and velocity are null, while the mean of the product of those variations is not. Most of the time the momentum is written as: ρu ≡ ρ u which is only a second order approximation. When the variation of ρ and u in a cell is high, as it is expected in the one-fluid model in the cells containing the interface, the approximation can become relatively weak. It's important to note that ρ usually scales with ∆x -1 as the interface numerical thickness usually scales with ∆x. On the other hand, u scales with the derivative of the velocity at the interface which is discontinuous with variable viscosity fluids as the jump condition µ ∂un ∂n = 0 holds. In practice, as the viscosity field is smoothed, this condition also scales as ∆x -1 . Hence, under this framework, the product ρ u should scale as ∆x -2 and should not be neglected with mesh refinement.

However, when the viscosity ratio is relatively small, and the interface regularized over several cells, the velocity variation remains quite negligible compared to the density variation. The question of the convergence of those approximations is beyond the scope of this article. Nonetheless, it is worth to note that an interface thickness diminishing at a rate smaller that ∆x 1 (e.g. scaling as ∆x 1/2 ) would ensure a proper convergence of the mean momentum. In our simulations, we have found the approximation ρu ρ × u to give sufficiently good results.

p. 6

Splitting and advection schemes. A second source of spatial error inherent of the standard non preserving approach comes from the splitting of the incompressible inertial term ∇ • (u ⊗ (ρu)) (first form) into two parts: ρ∇ • (u ⊗ u) + (∇ • (uρ)) u (second form, usually used in practice, where we have used the fact that ∇ • u = 0). In the latter, the left term is the density times the convection term and the second one is the density advection term (the same as in the continuity equation) times the velocity. While this decomposition holds on the algebraic part, it suffers from important approximations when solving at the discrete level. For the sake of simplicity, we base our reasoning on a collocated grid in 1D, using staggered grids and higher dimension would lead to similar conclusions.

In a control volume Ω i , we can easily compute the fluxes at the faces with traditional schemes. We apply the finite volume approach to the momentum flux in the x direction for the first form as:

|Ω i | -1 Ωi ∂ ∂x (ρuu) dx = (ρuu) i+1/2 -(ρuu) i-1/2 ∆x
where (ψ) i-1/2 , for example, denotes the value of the ψ field at position i -1/2 (e.g. the west face of the control volume). The second form yields to:

|Ω i | -1 Ωi ρu ∂ ∂x (u) + u ∂ ∂x (ρu) dx = (ρu) i (u) i+1/2 -(u) i-1/2 ∆x + u i (ρu) i+1/2 -(ρu) i-1/2
∆x where ρ (resp. u) holds for the mean value of the density (resp. velocity field) in Ω i . Needless to say that the approximation of the integral of the product has an impact on the approximation, as explained in the previous paragraph. But furthermore, the two forms do not equate numerically. If we assume a positive velocity and a first order upwind scheme for simplicity, the flux at the i + 1/2 (resp. i -1/2) face is approximated by (ψ) i+1/2 ψ i (resp. (ψ) i-1/2 ψ i-1 ), i.e. from the control volume before the face, for ψ being uu or uρ. In consequence, the two forms differ when developing the fluxes:

ρ i u i u i -ρ i-1 u i-1 u i-1 ∆x = ρ i u i u i -u i-1 ∆x + u i ρ i u i -ρ i-1 u i-1 ∆x .
The truncature terms coming from the approximations are not the same and the study of their impact on the solution may not be trivial. When knowing the origin of the flux, one could rearrange the equations to balance the inequity in order to reduce the inconsistency.

Furthermore, the problem is more substantial when using non-linear schemes (such as WENO) as the interpolation ρu of the momentum would never be equal to the product ρ u of the interpolations ρ of the density and u of the velocity.

Prescribed criteria

In conclusion, the simplifications usually used in the splitting approach induce supplementary error terms that can be avoided by the use of the initial conservative form of the coupled equations. Indeed, completing the latest article [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF], two criteria are necessary to ensure the coherent and consistent computation of mass and momentum:

1. the use of a synchronized temporal integration, 2. the use of consistent spatial conservative schemes both for the advective part of the momentum and mass equations.

Coupling with an interface method

It is clear that solving equations 1b and 1a with a common advection scheme would be sufficient to ensure those two criteria. For example, a high-order WENO explicit scheme could be used to advect the density field and the momentum precisely. However, in practice, doing so would introduce over time an undesirable smearing or sharpening of the interface region between the two phases, thus violating the immiscible fluids physical assumption. Hence, the coupling with an interface transport method is undoubtedly necessary. The literature proposes a wide range of representations and associated numerical methods such as Volume-of-Fluid (VOF), Moment-of-Fluid (MOF) and Level Set Method (LSM). By construction, those representations p. 7 ensure the immiscibility at the discrete level while, for some of them, are not fully conservative for the mass. Their detailed description is beyond the scope of this document.

Even though they are very commonly used, their carelessness use violates the second criterion as the solution of pure mass advection (i.e. the continuity equation giving the expected ρ * ) and interface transport (i.e. the effective density ρ n+1 for the next time step) do not match. The use of consistent high-order conservative schemes to accurately capture momentum fluxes is much more delicate and when coupled with various interface transport methods such as LSM, VOF and MOF. Associated particular methods have been developed for the momentum conservation. However, they are inherently linked to the interface representation and need specific algorithms that match exactly the associated mass flux. Our approach is parallel to theirs as it is generic and works accurately whatever the representation. Furthermore, it is not restricted to low-order schemes.

Proposed approach

High-Order Momentum Preserving (HOMP)

In this section we present the HOMP method proposed in order to preserve the momentum, i.e. reduce discrete momentum transfers. Our approach is prolonging the initial work of Bussman et al [START_REF] Bussmann | Modeling High Density Ratio Incompressible Interfacial Flows[END_REF] and relates to the recent article [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF]. The generalized MP approach is developed to be independent of interface representation, adaptable to many high-order advection scheme and it can be easily implemented in 2D and 3D as a modification of conventional CFD codes using the one-fluid approach. The proposed method takes advantage of high-order temporal and spatial schemes, Runge-Kutta with WENO 5, 3 (as explained in section 3.1.5), without necessitating specific treatment near the interface, and shows satisfactory results with VOF-PLIC, LSM and MOF. We develop below the different steps of the algorithm; the differences between the non preserving method, referred thereafter as standard method, and the MP method are shown in the table 2.

Benefits

We propose a discrete framework for incompressible two phase flows, named High-Order Momentum Preserving (HOMP), that: 1. takes advantage of existing high-order advection schemes, both spatially (WENO) and temporally (RK); 2. is independent on the underlying interface representation; 3. is almost effortlessly adaptable to an existing Navier-Stokes 2D and 3D code by the use of existing schemes; 4. drastically reduces the momentum transfer across the interface.

Algorithm

Step 1: Phase advection (interface transport) First of all, we solve the phase advection equation for χ n+1 thanks to the chosen interface method (VOF, LSM, MOF) and deduce from it the (smooth) final physical properties ρ n+1 and µ n+1 of the fluid at the time step t n+1 . For that purpose, as it is necessary for high-order temporal integration of the interface, we extrapolate the velocity at time t n+1 such as: u n+1 = 2u nu n-1 (see 3.2.3 for more details).

Step 2: Continuity equation (prediction) We solve the Navier-Stokes equations 1 in a split manner and start by the continuity equation. The objective of this first step is to compute a prediction of the density field ρ * ρ(t+δt), where δt = t n+1 -t n . One could directly solve the equation as-is but, as for incompressible two phase flows, ρ (resp. µ) is defined as a linear interpolation of ρ 1 and ρ 2 (resp. µ 1 and µ 2 , refer to equations ( 2)), we propose to solve the conservative advection equation on the characteristic function χ from which we deduce the local physical properties ρ * and µ * (that will appear in the momentum equation):

∂χ ∂t + u • ∇χ ≡ ∂χ ∂t + ∇ • (uχ) = 0 p. 8
written from the indentity ∇ • (uχ) = u • ∇χ + χ∇ • u and using ∇ • u = 0 for incompressible flows. When integrated over t n , t n+1 , the conservative advection equation leads to:

χ * -χ n + t n+1 t n ∇ • (uχ) dt = 0. (7)
This equation can be solved with high-order temporal and spatial schemes as it is traditionally done, for example, for the energy equation. In practice, the χ function is defined as the volume fraction of one phase (sometimes noted C in the literature), which is deduced from the underlying interface representation. We will discuss these two points in section 3.1.4. ρ * and µ * are predicted density and viscosity fields at the time step t n+1 , but they do not replace the final ones (noted ρ n+1 and µ n+1 ) that will be computed through the chosen interface transport method, may it be VOF, MOF, LSM, etc. When dealing with a staggered grid method, the cell based ρ * and µ * fields have to be interpolated/reconstructed onto the faces as it is commonly done before solving the momentum equation (see discussions in 2.2 and remarks in 3.1.6).

For example, equation 7 can be solved by a first order explicit scheme, such that:

χ * = χ n -δt ∇ • (u n χ n )
Higher order schemes will, of course, be used as shown in the tests of the method in section 4.

Step 3: Momentum equation (prediction) The predicted momentum (ρu) * is obtained in a split manner [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or threedimensional cavity flows[END_REF], after integrating equation 1b over the temporal interval, by solving:

(ρu) * -(ρu) n + t n+1 t n ∇ • (u ⊗ (ρu)) dt = t n+1 t n ∇ • (2µD) + f dt
where the pressure gradient term has been dropped and will be reintroduced in step 4. By taking the definition of (ρu) n ≡ ρ n u n (which is known at that point) and (ρu) * ≡ ρ * u * where ρ * has been calculated in step 1, we finally get the equation:

ρ * u * -ρ n u n + t n+1 t n ∇ • (u ⊗ (ρu)) dt = t n+1 t n ∇ • (2µD) + f dt (8)
that can be solved with various advection schemes.

The key aspect of our method is that the inertia term ∇ • (u ⊗ (ρu)) has to be integrated in the exact same numerical manner as for the continuity equation (step 1). The viscous term can be integrated in various manners. We will discuss these two points in section 3.1.4.

To illustrate further discussion in the document regarding the inertial term, we will note the advection equation part of 8 as:

ρ * u * = ρ n u n - t n+1 t n ∇ • (u ⊗ (ρu)) dt (9)
We refer the reader to section 3.2 for detailed descriptions of how to integrate the various terms of the equation.

Step 4: Pressure/Poisson equation The pressure gradient term that has been dropped in previous step is used to force the incompressibility of the flow. Thus, we can write:

(ρu) n+1 -(ρu) * = - t n+1
t n ∇p dt which can be implicitly approximated at first order by:

(ρu) n+1 -(ρu) * -δt ∇p n+1 p. 9
Now, a priori, nothing is known about (ρu) n+1 that could help us to solve this equation as-is. However, we assume that, as discretization steps (δx and δt) go to zero, ρ * → ρ n+1 : hence, by taking the approximation (ρu) n+1 ≡ ρ n+1 u n+1 ρ * u n+1 , we modify the previous equation as:

ρ * u n+1 -ρ * u * = -δt ∇p n+1
which finally yields to:

u n+1 -u * = - δt ρ * ∇p n+1
By incompressibility, the property ∇ • u n+1 = 0 is assumed, resulting in the well known Poisson equation:

∇ • u * = ∇ • δt ρ * ∇p n+1 (10)
The solution p n+1 of this equation in used in the next step to ensure incompressibility.

Step 5: Velocity correction The correction of the predicted velocity field u * permits to compute the incompressible velocity field u n+1 :

u n+1 = u * - δt ρ * ∇p n+1 (11)

Summary

Figure 2 presents the differences between the non preserving and the momentum preserving method. The first step is added for density prediction ρ * and can be reduced to the solving of a classic advection equation. After that, there are only very few modifications inside the other steps: basically, the replacement of ρ n by the predicted ρ * . The key improvement of the proposed HOMP method is the effective use of high-order spatial and temporal advection schemes for both density prediction and momentum inertial term in steps 2 and 3, in the whole domain and without specific treatment near the interface. Therefore, as we will see in the results, the algorithm permits to drastically reduce discrete momentum transfer while necessitating only few simple modifications inside the original code, independently of the interface representation method, and for a low computational cost.

Discussion on the need for a consistent numerical scheme for advection equations

The key idea behind the method is the adjoined use of a single numerical scheme for the solution of the advection terms in 7 and 8. When solving the advection term on ρ and ρu, it is mandatory to have very close numerical approaches as already discussed. The reason behind this necessity is that, when solving the momentum equation, whether the integration scheme is implicit or explicit, we divide the momentum ρ * u * by ρ * to obtain the velocity u * . If there is a large diffusion and/or anti-diffusion that do not match in the resolution of the advection of the density and the momentum, the mismatch could induce instabilities.

A simplified 1D case.. Let's assume that the velocity is constant: u = c and ρ is variable at the interface. Then, the first order explicit continuity equations reads:

ρ * = ρ n -δt c
∂ρ n ∂x and, with a similar discretization, the advection part of the momentum is treated as:

ρ * u * = ρ n u n -δt c 2 ∂ρ n ∂x = ρ n -δt c ∂ρ n ∂x c
assuming that u n = c. Hence, by replacing ρ * in the second equation by the continuity solution, we obtain:

ρ n -δt c ∂ρ n ∂x u * = ρ n -δt c ∂ρ n ∂x c p. 10

Momentum Preserving method

Step 1: Interface transport + Continuity

χ n+1 ⇒ ρ n+1 , µ n+1 (see section 3.2.3)
Step 2: Continuity equation (prediction)

χ * ⇒ ρ * , µ * (see. eq 7)
Step 3: Momentum equation (prediction

) (ρu) * (ρu) * -(ρu) n + t n+1 t n ∇ • (u ⊗ (ρu)) dt = t n+1 t n ∇ • (2µD) + f dt (see. eq 8)
Step 4: Pressure/Poisson equation p n+1

∇ • u * = ∇ • δt ρ * ∇p n+1
(see. eq 10)

Step 5: Velocity correction

u n+1 = u * -δt ρ * ∇p n+1
(see. eq 11)

Non preserving method

Step 1: Interface transport + Continuity

χ n+1 ⇒ ρ n+1 , µ n+1

Empty

Step 2: Velocity equation (prediction

) u * u * -u n + t n+1 t n ∇ • (u ⊗ u) dt = 1 ρ n t n+1 t n ∇ • (2µD) + f dt Step 3: Pressure/Poisson equation p n+1 ∇ • u * = ∇ • δt ρ n ∇p n+1
Step 4: Velocity correction which gives the expected result u * = u n = c if and only if the derivative operator ∂•/∂x (∇ in multi-dimension) is the same for both equations. Said in a different way: if a numerical error is made on the advection of the density (left part), the same numerical error should be done on the advection of the momentum (right part) to counteract the variation. This imposes the condition of consistency in the numerical schemes when computing the two advection terms.

u n+1 = u * -δt ρ n ∇p n+1

The importance of using WENO schemes

Within the one-fluid model, the sharp transition of the volume fraction (and hence the density and viscosity) from one phase to the other leads to a challenging problem. While the immiscible two-phase flow dynamic is based on an infinitely small interfacial region, the one-fluid principle is on the other side based on a differentiable density field, ensuring straightforward numerical discretization of both continuity and momentum equations. Hence, as discussed in the previous section, the finite support of the interfacial region, i.e. where χ varies from 0 to 1, has an important impact on spatial discretization errors. The sharper the transition, the more it is bound to be spread by numerical diffusion over neighbouring cells and hence induce momentum transfer.

In consequence, dedicated numerical schemes has to be used in order to ensure that under advection, ρ remains as close as possible to the ideal analytical solution. For that matter, WENO schemes, first introduced in [START_REF] Shu | Efficient implementation of essentially non-oscillatory shockcapturing schemes[END_REF] for hyperbolic conservative equations, have been widely used, particularly for compressible flows with shocks, and improved over the last decades in numerous articles. They have been designed to reduce numerical diffusion while avoiding most of oscillatory errors. Numerical experiments show that they are capable of capturing accurately the advection of a step function; hence, they are a genuine, reliable choice for resolving the advection of ρ and ρu.

It is important to note that, if the interface thickness reduces linearly with the spatial step, the high-order p. 11 convergence of WENO schemes is not achieved per se. However, numerical errors are still much lower than, for example, a first order upwind scheme. As we will demonstrate in the results section, within the momentum preserving framework, using a WENO 5, 3 scheme instead of a first order upwind and third order WENO for computing fluxes at the interface permits reducing numerical damping, capturing small flow structures and converging faster towards the expected result. While WENO schemes do not ensure boundedness of advected quantities, they are still very robust and limit oscillations. Indeed, adding limiters or developping taylor-made spatial schemes, as in [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF] for instance, will surely benefit to the overall accuracy and stability of the method. This topic is the subject of ongoing works.

3.1.6. Five important remarks Remark 1. Even though the proposed method results in the reduction of transfer of momentum between the two phases, the last step of the algorithm replaces the ideal ρ * (here we use ideal in the sens that, as-is, if one could use the solution ρ * for the phase transport, there would be no additional error) by ρ n+1 which is a result of the phase advection. This part is usually done with level set, VOF or any other traditional method that permits to keep a sufficiently sharp and controlled transition between the two phases. In practice, ρ n+1 is the solution of the translation of the interface position, transferred into a Eulerian volume fraction (which is inherent to VOF method or obtained through a cut cell algorithm or, more simply, a Heaviside function for level set). In the end, it is not the momentum ρ * u n+1 that is the solution at the next time step but rather ρ n+1 u n+1 . Hence, a variation/transfer of momentum, proportional to ∆ (ρu) u n+1 ∆ρ = u n+1 ρ n+1ρ * , will nevertheless be introduced by the method, as it is the case with the standard approach when no specific momentum conservation algorithm is used. However, this spurious momentum is much less than the one without the proposed method which removes the error inherent of the convection process. [START_REF] Popinet | Numerical Models of Surface Tension[END_REF][START_REF] Ghods | A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods[END_REF] alleviate this problem by designing a method adapted to the interface representation; mostly being limited to first or second order precision. On the other hand, our method is designed to be independent of that interface representation while permitting to use high-order advection schemes as we will show in the results section 4.

Remark 2. The continuity equation can be efficiently implemented as far from the interface the density field ρ (equivalently the characteristic function χ) is constant (because of the incompressible hypothesis). Therefore, the advection term ∇ • (uρ) is numerically only non-null near the interface, where the stencil of the advection scheme crosses variations of ρ. For order schemes, the stencil is larger and thus requires more computational efforts. Nonetheless, the relative number of cells close to the interface where we need to compute the term ∇ • (uρ) compared to the total cells in the volume where we need to compute ∇ • (u ⊗ (ρu)) reduced when the mesh is refined. Hence, the additional cost of the HOMP method is one dimension less than the cost of the mandatory inertial term.

Remark 3. Most modern finite volume codes are based on staggered grids where the density is defined at the center of cells while the velocity components at the associated faces. For computational efficiency, unless otherwise stated, we have made the choice to solve the continuity equation on the cells and the momentum on the faces. This requires the interpolation of the density field ρ * from the cells onto the faces control volumes after the prediction (step 1 of the algorithm). Even though this approach introduces additional numerical errors, we haven't found this approximation to cause any problem during the development of the method and the simulation of fluids with density ratio up to 10 3 ; a second order centered interpolation was found to be sufficient. We also have implemented the advection of a face-based density field in order to link more closely the solution to the momentum equation. This discretization strategy has the main drawback to be more computationally costly as it requires the solution of two more advection equations in 3D (1 for each component). Nevertheless, the increased cost can be controlled following the previous remark.

Remark 4. There is a conjoint necessity for the use of the same temporal and spatial schemes to integrate the continuity (step 1) and the inertial term (step 2) in order to ensure consistency. Alternatively to the proposed method, one could think of using the solution ρ n+1 obtained after the interface transport equation as the predicted ρ * . However, that approach would violate the principle of having a paired numerical method for the continuity equation and the convection term. We have found it to give very poor and unstable results, particularly when using an interface transport that is not strongly mass conserving.

Remark 5. The reader can note that, in the splitting made for the NS equations and the interface advection step, we have chosen to start by the advection of the interface (step 1 of the algorithms) instead than at the end. While the study of that matter is beyond the scope of this article, by experience, we have found p. 12 that, when using the standard method, the simulations were much more stable when doing so. Regarding the proposed method, both choice give similar results.

Numerical schemes and discretization

The Navier-Stokes equations are solved on a staggered grid using the time splitting correction method [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or threedimensional cavity flows[END_REF] to account for the velocity-pressure coupling. The overall temporal approach is thus first order. However, as inertia plays a key-role in our momentum preserving method, high-order schemes are important to solve the advection equation, as detailed below. A validation with various numerical schemes is performed in section 4.

Linear systems of momentum and pressure equation are solved with the generalized minimal residual method (GMRES) initially preconditioned by a left Jacobi method. We have used the HYPRE library which provides precise results with very good efficiency and scalability when using multi processors parallel computation.

The proposed algorithm is independent of the interface representation, hence the verification and validation tests are performed using three different methods: VOF-PLIC, MOF and LSM. The three of them and keys of their implementation are described subsequently below.

All of the presented numerical methods and algorithms have been developed inside the Open Source Notus CFD code 1 which has been, among others, validated for two-phase flows. The HOMP method is available for testing and reproducing results starting from Notus v0.5.0.

Conservative advection equations

Any kind of schemes can be used to integrate the two advection equations 7 and 9, may they be implicit or explicit. In practice, as we want to reduce the discrete momentum transfer as much as possible, we will use high-order temporal and spatial schemes such as a 2 nd order Non-Strong Stability Preserving explicit Runge-Kutta method (RK-NSSP 3, 2) and a 5 th order WENO scheme. We refer the reader to the article [START_REF] Wang | Linear Instability of the Fifth-Order WENO Method[END_REF] for a detail description of those schemes. We also have tested our method with a 3 rd order temporal integration scheme (NSSP 5, 3) but this hasn't shown a remarkable difference in the results' accuracy while necessitating more computational time. Nevertheless, there is no limit regarding the temporal integration scheme -which can be conveniently increased when the time step becomes large -and any type of explicit scheme could be used as long as it is consistent for both mass and momentum equations.

Viscous term

The viscous term

t n+1 t n ∇ • (2µD
) dt in the momentum equation can be approximated in an explicit or implicit form. As the time step restriction can become very problematic for small δx, the use of an implicit scheme is usually advised. Hence, we approximate the integral by

t n+1 t n ∇ • (2µD) dt ∇ • 2µ * D n+1
where the deformation tensor D is added to the linear system and µ * is the predicted viscosity field obtained through χ * .

Interface transport and characteristic function / volume fraction

Volume-of-Fluid method (VOF), is a two-fluid framework directly based on the volume fraction function χ. We use the geometric Piecewise-Linear Interface Calculation scheme (PLIC) [START_REF] Rider | Reconstructing Volume Tracking[END_REF] coupled with a first order in time advection scheme. Moment-of-Fluid method (MOF) is an extended VOF method which incorporates material centroid in addition to material volume fraction for interface reconstruction, originating from [START_REF] Dyadechko | Moment-of-fluid interface reconstruction[END_REF].

In order to reduce the computational cost of MOF, we use the analytical reconstruction proposed by Milcent et al [START_REF] Milcent | Moment-of-fluid analytic reconstruction on 3D rectangular hexahedrons[END_REF]. Level-Set (LS), introduced by Osher et Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] is a two-phase representation based on a signed function φ, usually a distance function. The function is transported by the fluid through an advection equation; it requires applying frequently a re-initialization procedure. The χ function is then computed through a Heaviside function (described below). Hence, the LSM method appears to be the closest approach to the continuity equation as it is inherently based on the advection of a scalar field related to the density field via linear interpolation. The advection equation is solved using a WENO 5, 3 scheme and a 2 nd order 1 Notus CFD code: http://www.notus-cfd.org is developed in the I2M Laboratory p. [START_REF] Gibou | A review of level-set methods and some recent applications[END_REF] Runge-Kutta NSSP 3, 2 temporal integration method; the frequent re-initialization is done by solving the classical Eikonal equation in the essence of [START_REF] Russo | A Remark on Computing Distance Functions[END_REF] in a band around the surface.

As the interface transport algorithms are beyond the scope of this article, and because the HOMP method is independent of them, we refer the reader to Notus CFD documentation for more information about the implementation details of the associated algorithms.

These methods are regularized near the interface in order to avoid high gradients (namely ∇ρ) and hence damp eventual numerical instabilities due to the one-fluid inner representation. In consequence, the transitional region between the two phases (the number of cells over which ρ ∈]ρ 1 , ρ 2 [) is -unless stated otherwise -made of 6 cells. To do so with VOF and MOF methods, we have smoothed the volume fraction with an inverse distance weighting method, using 3 consecutive iterations (i.e. spreading the interface in the 2 × 3 surrounding cells). Regarding LSM, unless stated otherwise, we have simply used a regularized Heaviside function to translate the level set function φ into a smooth volume fraction field χ (x) = H (φ (x)) where = 3δx and the Heaviside function is defined as:

H (φ) =          0 if φ > 1 if φ < - 1 2 1 + φ + sin(π φ ) π otherwise (12)

Pressure gradient

In practice, in order to obtain more precision on the pressure term, the pressure gradient in step 3 is split in two terms as proposed by Goda [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or threedimensional cavity flows[END_REF], such that: ψ = p n+1p n . Hence, an explicit term -∇p n is added to the right-hand side of 8. After that, the solution of the Poisson equation 10 will instead be the pressure increment ψ, which is then used to correct the velocity field 11. Finally, the pressure field is incremented with p n+1 = p n + ψ.

Gravity term

We have chosen to implement a 2 nd order Crank-Nicolson scheme for approximating the gravitational term in the momentum equation such as:

f = ρ n + ρ * 2 g ( 13
)
where g is the gravity vector. This scheme permits to obtain a very good precision on the gravity force, particularly important on the standing wave test case 4.2. It has for principal advantages of being much more stable than the explicit scheme (i.e. using ρ n alone) while being much less diffusive than the implicit scheme (i.e. using ρ * alone).

Surface tension term

The surface tension term is introduced in the momentum equation following Brackbill et al. [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] in a Continuum Surface Force (CSF) manner such as:

f = γκ∇χ ( 14
)
where γ is the surface tension coefficient and κ the interface curvature computed via the particular interface representation.

Results

In order to demonstrate the effectiveness of the HOMP method, we present here the results obtained on various test cases, from the more numerical one towards more physical applications. In section 4.1 we show the convergence of the method and its ability to preserve momentum on the advection of a heavy droplet. We believe that the gravity and capillary waves test cases are very much representative of the momentum transfer problems as they focus on the dynamic of one wave length. Even though non-linear interactions between gas and liquid phases are very important in daily applications, the study of separate wave lengths can, on a p. 14 simplified basis similar to linear theory, help to exhibit numerical errors and their dynamic implications. In consequence, we have studied into more details the associated test cases in sections 4.2 and 4.3 in which we show how the HOMP method permits to drastically reduce the transfer of momentum on standing gravity and capillary waves, in 2D and 3D. Finally, in sections 4.4 and 4.5 we apply the method to the simulation of more complex dynamic cases: the splashing drop and the plunging of a wave.

Ad hoc mixed method. In order to assess the importance of using a high-order non linear scheme such as WENO 5, 3 in the whole domain, even near the interface, and to compare various numerical schemes, we have designed an ad hoc strategy that permits chosing a different numerical scheme near the interface and in the phases bulk. Desjardin et al. [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF] have proposed a similar strategy and used a first order scheme in the interfacial region. In our article, we reference this method as ad hoc as it is made for comparison purposes and does not represent the prescribed numerical strategy, i.e. we will see in the following results that they confort the relevance of using a unique scheme in the whole domain.

This ad hoc method, refered as mixed Momentum Preserving thereafter, is used when computing the fluxes of mass and momentum in the advection equation where ρ varies around the control volume. First, we define an ensemble of cells E Γ surrounding the interface Γ, i.e. a band, as all cells adjacent to a variation of density. In algorithmic terms, in 2D, a cell Ω i,j ∈ E Γ if there exists a direct neighbouring cell Ω i ,j for which ρ i,j = ρ i ,j . This ensemble is computed after each update of the density ρ n+1 in the first step of the MP algorithm. When solving the advection term for the density and momentum predictions, we can then chose a particular scheme depending on the proximity of the interface, i.e. for the x direction, at a face i + 1/2 if either Ω i,j or Ω i+1,j is part of E Γ , ρ and ρu are reconstructed with a scheme Θ Γ , otherwise a scheme Θ b (for bulk) is used. The couple (Θ Γ , Θ b ) represents the mixed numerical schemes. Unless stated otherwise, in this section, the mixed scheme is not applied, i.e. Θ Γ ≡ Θ b simply denoted by HOMP and, when the mixed approach is used, we set a fifth order WENO scheme in the bulks, i.e. Θ b ≡ WENO 5, 3, as it is quite usualy used for monophasic flows.

Post processing. The momentum is measured in the whole domain and in each phase through the computation of the integral:

M i = Ωi ρ i u • dx = Ω χ i ρ i u • dx
for the i th phase, where χ i is the associated characteristic function. From the discrete point of view, as we are using a staggered grid, the characteristic function is interpolated onto the faces and the integral is approximated as the sum over each face control volume. Similarly, the total kinetic energy in a domain Ω is computed as:

E k Ω = Ω 1 2 ρu 2 • dx.
We define the density-vorticity as the product ρω where ω = ∇ × u is the vorticity vector (i.e. reducing to a scalar in 2D). It is derived from the momentum and, in the hypothesis of separated phases where ρ is constant per phase, is equal to momentum rotational. This is an important measure as it can exhibit spurious vortices that arise in each phase, as we will see in the following sections. The weighting of the vorticity by ρ is relevant as even a low intensity vortex convey a lot of energy in a dense fluid such as water.

Convection of a high density droplet

Description

The spurious transfer of momentum due to model or numerical approximations is directly related to the density ratio of the two fluids. Increasing this ratio will increase the phenomena. Thus, the numerical ideal test case of the convection of a high density droplet, with a ratio of density ρ l /ρ g = 10 6 , can be used to exhibit and detect spurious transfers of momentum. This tailor-made case has been studied in several references such as [START_REF] Ghods | A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods[END_REF][START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF][START_REF] Mirjalili | Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[END_REF][START_REF] Kumar | A novel consistent and well-balanced algorithm for simulations of multiphase flows on unstructured grids[END_REF][START_REF] Manik | A generic algorithm for three-dimensional multiphase flows on unstructured meshes[END_REF][START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF][START_REF] Bussmann | Modeling High Density Ratio Incompressible Interfacial Flows[END_REF]. It is not based on a physical experiment and should be taken with care as it is very sensitive to initialization and numerical methods (see below for explanations). However, in order to compare to the literature, we study here the impact of the use of the HOMP method.

The droplet is initialized with a uniform velocity and the gas at rest. Ideally, as the droplet is very dense, its shape should remain the same and its total momentum should be conserved. If erroneous momentum is transferred into the droplet, we expect it to deform and, conversely, the gas is expected to be highly perturbed by this exchange. It is important to note that, in the continuity of the remarks in Nangia et al. [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF], within the incompressible framework and as the liquid density is finite (though very high), we should anyhow (i.e. even without any numerical error in the momentum equation) expect a deformation of the drop. This is due to four facts: p. 15

• the initial velocity field imposed in the whole domain is more and more discontinuous as the mesh is refined.

• as well, the flow should be incompressible; hence, there is no trivial analytical solution for instantiating the test.

• finally, as when studying the flow around a sphere, the pressure in the gas will be high in front of the droplet and low in the wake. That will undeniably lead to a compression of the liquid phase upstream and its elongation downstream. Moreover, the periodic translation of the droplet in the gas medium induces non-linear fluid dynamic in the latter, making the appearance of vortices inevitable and their interaction with the liquid hardly predictable.

• even though the density ratio is important, the dynamic of the liquid is not the one of a rigid body.

In consequence, the velocity in the droplet will not be uniform unless a more strict initialization (such as with penalization) would be applied at the discrete level, which is beyond the scope of the article. However, as the density ratio is very high, the variation of velocity should not be too important, while still enough to deform a bit the liquid domain. Without viscous effects, small structures are expected to arise when refining the mesh.

As the two incompressible fluids are non-viscous, and in the absence of gravity and surface tension, the total kinetic energy should be conserved. We will see that it is not the case with the standard method. The problem is even more pregnant with such high density ratio ρ l /ρ g as, for example, if a small erroneous momentum is transferred from the gas to the liquid, the resulting momentum in the liquid is increased by ρ l ρg . Moreover, with the standard approach, the result even diverges as the mesh is refined. On the other hand, the HOMP method permits conserving well the momentum by reducing momentum transfers across the interface.

In conclusion, we will first study the simulations in 2D with a constant interface thickness, then with a reducing thickness. We will finally present the results obtained with the 3D simulation of a spherical droplet to show the capacity of the method to well resolve a 3D flow.

Physical parameters, initial and boundary conditions.

A droplet (spherical in 3D and cylindrical in 2D) of radius R = 0.2 m holds into a domain Ω = [0, L] 2 of size L = 1 m, located initially at the center. The boundary conditions are periodic in the x direction and slip on top. As the problem is purely symmetric, in practice, we have only computed the top half of the domain. The heaviest fluid's density is set to ρ 1 = 10 6 kg/m 3 and the lighter one to ρ 2 = 1 kg/m 3 . In order to focus on the inertial term and so that the viscous term doesn't affect the solution, both fluids are inviscid.

The initialization of the velocity field is consequential as it drives the whole case. At t = 0, u 0 = u (t = 0) is set to u 0 (x) = χ 0 (x) • (1, 0) T where χ 0 (x) is the initial characteristic function (i.e. the volume fraction) of the droplet, given by the underlying interface representation (see below for the corresponding results). This initialization is consistent with the finite volume framework and is necessary to regularize the velocity field which would otherwise be discontinuous. Additionally, the velocity field is projected onto the divergence free space of solutions in order to use an initial condition that is physically meaningful (otherwise, it would be imposed only at the end of the first time step). Hence, the velocity inside the droplet is not exactly uniform and would in any case deform it. We however expect this deformation to be as small as possible. The result will clearly show the inadequacy of the standard method to capture the expected dynamic.

Numerical methods. Concerning the computation of ρ * (equation 7) and the inertial term (equation 9), we have studied the impact of three numerical schemes: 1 st order, 3 rd order WENO 3, 2 and 5 th order WENO 5, 3 coupled with 2 nd order (NSSP3, 2) Runge-Kutta integration. As, with such a high density ratio, the truncature terms in the continuity equation will bring a lot of errors, we have decided to use the resolution of that equation based on the faces of the mesh. This helps by removing the approximations that are made when using the interpolation from cells to faces after the advection step and enforcing the coherence with the momentum equation. This appeared to be less necessary for smaller density ratios (see next test cases).

We have compared the use of three different underlying interface representations: LSM, VOF and MOF. As prescribed before, we have made sure that the interface thickness is of smoothly varying through several cells (see below for precision), thanks to a standard smoothing method described in paragraph 3.2.3.

p. 16

Spatial discretization. We have applied a regular and uniform grid spacing for the discretization of the domain with N 2 cells (i.e. N × N 2 cells with symmetry) such that δx = L N . We denote by S the scale of the simulation such that N = 64 • S. Integration time step. The time step has been set to δt = 10/(31.25 • N ), which follows a CFL-like condition. With the imposed initialization, this corresponds to a Courant number σ = δt|u| max δx around σ 0.3 (in practice, this would rise particularly in the presence of vortices in the gas, without disturbing the results). This is 10 times larger than the step used in [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF]. With the second order temporal scheme for the density prediction and the inertial term that we have used, we have found the difference to be relatively small while permitting to save a lot of computational time. It also demonstrates the robustness of the HOMP method when using large time steps close to the explicit CFL restriction.

Results

Constant interface thickness. First, in order to elude the problem of a discontinuous velocity field when refining the mesh, and thus reduce the inevitable deformations of the droplet, we first study the convergence of the methods with a constant interface thickness (i.e. the physical width over most of the interface Dirac mass is held). It is imposed through the smoothing process applied to the sharp volume fraction computed by the interface algorithms. In the particular case of LSM, in order to match as precisely as possible to VOF and MOF, we have used the sharp volume fraction computation as proposed by Min and Gibou [START_REF] Min | Robust second-order accurate discretizations of the multi-dimensional Heaviside and Dirac delta functions[END_REF]. The smoothing procedure is applied 3 • S 2 times thus ensuring an interface thickness of around ∼ 0.1 • m, whatever the scale is. Hence, for the finest mesh N = 256, there are approximately 24 cells for the transition of density from gas to liquid. In practice, it is impossible to define and control very precisely the interface thickness, but the proposed approach permits to obtain a sufficient approximation for our purposes.

Results of the simulations with various not-mixed (i.e. Θ Γ ≡ Θ b ) schemes and interface representations are shown in figure 3. Globally, we can see that, as expected, the shape of the droplet does not converge to a perfect disk: as the horizontal velocity increases in the wake of the droplet (at the line of symmetry), the liquid is pushed toward the right, inducing a small pinch. We clearly see in the left sub-figure of 3 that the first order scheme is much more diffusive and brings an important deformation of the liquid phase, even stretching the interface in the opposite direction. However, thanks to the HOMP method, the droplet is much more compact because less momentum has been transferred in the gaseous phase. Regarding that, the blue contours captures important momentum in the lighter phase -note that, at those locations not far from the interface, because of the smoothing of the density field, ρ > 1, and thus the x momentum ρu is impacted. It is worth noting that momentum is also impacted by lower-order schemes in the bulk of the phases; we will study the impact of using WENO 5, 3 scheme in the bulk and different schemes near the interface with the ad hoc mixed method in the next test cases. Figure 3: Convection of a high density droplet with constant interface thickness: comparison of the numerical convergence of the shape and x momentum for various numerical and interface methods. Time is t = 1 s. For each simulation, the top half part is the result using standard method, the bottom half one with HOMP method. Purple filled color is for momentum up to 10 6 (which is the momentum inside the drop at initialization), white is for ρu = 0. The contours highlight particular values: black for the interface position, blue (in the gas phase) for ρu = -10 and red (inside the droplet) for ρu = 10 6 .

p. 17

Reducing interface thickness. As the one-fluid model aims to capture immiscibility and sharp transition between the phases when the interface thickness goes to zero, we study the convergence of our method compared to the standard one. For that purpose, we used the same smoothing procedure as previously, but here we have set the number of iterations to be 3 • S. This way, is still very much comparable from an interface method to another but however, it scales as the square root of the mesh size, in other terms:

∼ O δx 1/2 . While in [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF] the authors have used a very sharp interface thickness (of 2 • δx in LSM), as many authors have noticed, we believe that the interface thickness has to scale slower than linearly with the mesh in order to obtain converging and stable results.

In this configuration, as the gradient of density scales at a slower rate than the mesh refinement, we expect the transfers of momentum with the standard method to be much higher than in the previous paragraph. As δx → 0, both the velocity and the density fields are getting sharper; hence, the momentum being the product of both, will see the high-order truncature terms in the advection equations increase at doubled rate. Again, as there is no trivial analytical/physical solution to this test case, we do not expect a regular and smooth solution regarding the interface and the flow around the droplet.

Moreover, high-order numerical schemes, showing very small numerical diffusion, can eventually induce more error in the simulations as the mesh is refined. It is important to note that it is an expected result on that particular test case that is not correlated with the HOMP method but is due to the inherent stiffness of the case. For that matter, lower order schemes might be preferred but, for general and more physically based experiences as shown in the next sections, the smoothness of the incompressible flow permits to properly use high-order schemes.

The results are shown in figure 4. We first see that, as the mesh is refined, the dynamic of the fluids is more complex and brings more perturbations to the shape of the droplet. As explained in the description, this behavior is expected as the velocity field is getting close to a discontinuity. In the absence of forces that would keep the spherical shape such as surface tension, the interface eventually deforms. As it has been observed in the literature, the standard method (top part of each simulation) leads to unstable and destructive dynamic that eventually make the code diverge. This is mostly due to erroneous transfers of momentum from the gaseous phase into the liquid phase; a phenomenon that is increased as the density gradient goes to infinity.

On the other hand, all the results using the HOMP method (bottom part of each simulation) demonstrate a much better handling of the underlying physics. Even when using the non-linear fifth order WENO scheme that is bound to induce more vortices in the inviscid flow, the simulations stay stable for all the interface methods for N = 64 and N = 128. For the finest mesh (N = 256), the initial discontinuity of both velocity and density induces a lot of high-order variations in the flow: we can see that the standard method produces much more deformation of the interface while the HOMP method manages to keep the droplet more compact.

We also see that our method introduces much less momentum extrema in both phases, particularly in the gas, demonstrating a reduced exchange of momentum across the interface. The reader could note that, as the LSM is smoother by essence and its advection closer to the continuity equation, the interface is less deformed as we refined the mesh compared to VOF and MOF.

While the use of the first order scheme gives very satisfactory results, even with coarse meshes, the WENO schemes bring high-order terms in the flow that are only lightly damped; thus resulting in more perturbations of the interface. Our results differ from the one exposed in [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF] where the interface seems smoothed and closer to a sphere. We explain this difference by three facts: first, the result of the simulation is very sensitive to the initial condition of the velocity field, secondly, their LSM interface thickness is much smaller thus permitting to enclose the momentum more sharply inside the droplet, and finally, the authors have adapted a particular advection scheme with limiters that enhances the preservation of momentum while we have made the choice to use more common schemes. Developing high-order advection schemes that are adapted to the coupled continuity/momentum equations is an important factor. Nevertheless, the results obtain with HOMP is much better than without it. We will see the accuracy of our method in the next sections with more common density ratios of 10 3 and more physically relevant cases.

The plots of fig. 5 show how the total kinetic energy E k Ω in the whole domain accompanied by the liquid phase volume. First, we observe a clear link between kinetic energy and mass conservation; it is straightforward that a gain in mass is accompanied by a gain in kinetic energy, and conversely. We also clearly see that the standard method (dashed lines) conserves much less the kinetic energy. While numerical diffusion can explain the important decrease with the first order scheme, we see that even higher order p. 18 schemes, particularly with coarse meshes, show a lot of reduction of E k . This can be explained by the transfer of momentum from the more dense phase towards the less dense one, as it is illustrated by the increasing momentum in fig. 4. For the finest meshes (N = 512) the simulations become unstable with time as more small structures arise and spurious momentum.

On the other hand, we see that the HOMP method enables much better the conservation of the kinetic energy with all three numerical schemes. The fifth order WENO shows the best results while permitting capturing fine structures in the flow. The increase in E k with the first order scheme as well as with the low-order VOF-PLIC can be explained by the important diffusivity of these methods that leading to small transfer of momentum from the gaseous phase towards the liquid one, coupled with the important gain in the liquid phase volume, increasing more importantly the energy. The reader can note that those quantitative results are very similar for LSM and MOF which both benefit from high-order transport equation solving. The latter ensures a better mass conservation and hence for the energy. The HOMP method keeps the simulations very stable even with very fine meshes. We will demonstrate this important property later on in more complex applications such as the breaking wave 4.4 and the drop impact 4.5. For each simulation, the top half part is the result using standard method, the bottom half one with HOMP method. Purple filled color is for momentum up to 10 6 (which is the momentum inside the drop at initialization), white is for ρu = 0. The contours highlight particular values: black for the interface position, blue (in the gas phase) for ρu = -10 and red (inside the droplet) for ρu = 10 6 .

3D simulations. We have experimented the same test case as above in 3D where we have used symmetric conditions for the bottom y = 0 and the back z = 0 planes. The interface thickness as in the previous paragraph, is linearly following the mesh size, i.e. ∼ 3δx. The time step is the same as in the 2D case.

We have used a WENO 3, 2 -NSSP3, 2 RK scheme for the advection terms. Figure 6a illustrates the shape of the interface using LSM (reinitialization at every time step) and VOF for that mesh, without (top half parts) and with the use of HOMP method (bottom half parts). We see that the method grants a reduced deformation of the interface and much less transfer of momentum into the gaseous phase. In 3D, the liquid is sucked in the wake by the low pressure zone (which is more important than in 2D), as commonly observed when studying rain drop free fall. Quantitatively, we can see with the isocontours in fig. 6b a reduced transfer of x momentum inside the air phase through time with the HOMP method. Moreover, in the continuity of the 2D results, the total kinetic energy E k Ω , shown in fig. 6b, is also much better conserved. The decrease of kinetic energy with the standard method is mostly due to erroneous transfer of momentum across the interface from the liquid to the gaseous phase where the density is 10 6 times lower.

Standing gravitational waves 4.2.1. Description

This verification case was designed in order to show the impact and accuracy of the proposed HOMP method on a very simple, yet classical, problem. A sinusoidal standing gravitational wave is the sum of two progressive waves of same amplitude and length. They exhibit a trivial periodic dynamic that captures well p. 19 any discretization and numerical error. Under linear assumption, for high Reynolds number and neglected viscous effects, a standing wave should keep its sinusoidal shape, its amplitude and its frequency of oscillation.

Still, according to linear theory, we expect the fluid's velocity to be maximum at the interface; this clearly justifies the elementary interest of the test case as the numerical errors should be very much focused at that place, hence highlighting the enhancements brought by the HOMP method.

When written in the vorticity form, the two-phase inviscid flow equation of momentum can be simplified to:

Dw Dt = (w • ∇) u + ∇ρ × ∇p ρ 3
where w = ω/ρ, D • /Dt the particle derivative and (w • ∇) u is the vorticity stretching term, canceling in 2D. Hence, when the density variation across the interface is not orthogonal to the pressure gradient, vorticity appears. Regarding the standing wave, as |∇ρ × ∇p| = 0, particularly around the poles of the sinusoidal, we expect to have a regular and natural increase in the vorticity around the interface that will eventually create a thin boundary layer around the interface. At the crest and trough, the vorticity has to cancel as ∇p and ∇ρ are orthogonal and also because the flow is symmetric. The initial velocity field is set to 0 (which corresponds to a standing wave at its maximum potential energy and null kinetic energy). The initial pressure field is also set to 0, it will be naturally set to hydrostatic after p. 20 the first iteration of the algorithm. The wave maximum vertical velocity is approximately U max 0.8 m/s and its is T 1.3 s.

Physical parameters, initial and boundary conditions.

Forces in the momentum equation.

As we want to focus on inertial effects, and as we want to take the same assumption as for linear theory (Re = ρwaterUmaxλ µwater 8 • 10 6 ), we have canceled the viscous term in the momentum equation. On the other hand, as gravitational effects are dominant (F r II = U 2 max gH/2 6.5•10 -3 ), we have chosen a a 2 nd order Crank-Nicolson scheme for approximating the gravitational term in the momentum equation (as presented in section 3.2).

Numerical methods. Concerning the computation of ρ * 7) and the inertial term (equation 9), we have studied the impact of three numerical schemes: 1 st order, 3 rd order WENO 3, 2 and 5 th order WENO 5, 3 coupled with a 2 nd order (NSSP3, 2) Runge-Kutta integration. Finally, we have compared the use of three different interface representations: LSM, VOF and MOF. As prescribed before, we have made sure that the interface thickness is of smoothly varying through 6 cells, thanks to appropriate regularization and smoothing.

Spatial discretization. We have applied a regular and uniform grid for the discretization of the domain with N × 2N cells such that δx = L N . We define the number of cells per wave length as N λ = λ δx = N . For ease of convergence study, we denote by S the scale of the simulation such that N = 64 • S. Integration time step. An important point concerns the time step: when we refine the grid, as δx → 0, the interface thickness of width 6δx is shrinking. Hence, a CFL-like condition has to be applied in order to keep the number of cells traveled by the interface per second constant. It permits an appropriate study of numerical spatial and temporal convergence at once. We have thus chosen to make the time step decrease linearly with the mesh:δt = 0.05/S seconds. Relatively to the maximum velocity of the wave reached when the interface is flat, u max 0.8 m/s, this corresponds to a CFL number of approximately σ = δt δx • u max

0.05•64 L • 0.8 1 4 .

Results

The validation of the method and its accuracy is demonstrated through the comparison of the various numerical methods and convergence tests. By symmetry, the x component of the momentum integral is null (in practice: very close to machine precision, unless stated otherwise). We expect the y component to oscillate between a maximum (positive) and a minimum (negative) value in a periodic trend, when the interface is flat, while being null when the wave is fully deployed. If energy is conserved, the interface should remain exactly the same on a regular period T of oscillation. Moreover, there should be no transfer of momentum between the two phases and hence one phase should not disturb the other's.

p. 21

The reader has to note that, as air has a very much smaller density than water, the air phase will contain much more vortices than the water phase. Regarding the total momentum, there is a factor of roughly 400 between the expected maximum values of M 1,y and M 2,y . Also, as spurious vortices are pushed inside both phases and evolve with time, leading to a non-null momentum temporal mean; this is particularly observable in the air momentum plots and the visualization of vortices.

Due to numerical errors, there exists an effective and significant transfer of momentum which will be exhibited thereafter, which might not be apparent if only looking at the wave height evolution. As we will show, the HOMP method manages to stabilize numerical simulation by drastically reducing the numerical transfer of momentum near the interface and the intensity of spurious vortices.

Advection schemes spatial convergence. First, we now compare the impact of the use various numerical schemes for the computation of the advection equations in the momentum and the density prediction. For that purpose, we study the spatial convergence for the first order, the WENO 3, 2 and the WENO 5, 3 schemes, with and without the ad hoc mixed method. The interface representation used is the LSM with a regularization parameter of = 3δx. The temporal integration for the advection terms is done through a NSSP3, 2 scheme. Figure 7 shows the total momentum temporal variation in each phase after t = 90 s (i.e. after around 70 periods) for various methods and meshes, while figure 10a shows the result with the mixed (O1, WENO5) and the most precise (WENO5, WENO5) for the whole simulation. Therein, we see that the first order scheme has significant impact on the damping rate of the gravity wave, particularly for coarse meshes. Table 1 summarizes the damping coefficients for both schemes, showing an absolute ratio of 1/3 for the finest mesh in favour of the HOMP method and a greater convergence rate.

We recall that spatial convergence is studied under constant CFL number (i.e. δt scales linearly with δx); for clarity, we have only showed scales S = {1, 4, 16}. It is important to note that the total momentumwhich will be studied first -is a criterion that has to be completed with the effective momentum distribution as well as the vortices inside each phase -which will be studied afterwards. Also, for all methods, at the finest mesh, the wave's amplitude are of the same order of magnitude and all exhibit (see figure 8) a persistent sinusoidal shape. Despite that, we show below that the standard method fails to reproduce a correct dynamic and even diverges while the HOMP method exhibits very good converging properties.

We can note in sub-figure 7b that the total air momentum seems to converge at the finest mesh around a value of 1 kg • m -2 • s -1 for all methods. However, compared to that value, the standard method coupled with the first order scheme exhibits a significant total momentum in the air phase for the coarse meshes. This is due to an excessive diffusion of momentum that leads to transfers between the two phases. While being reduced when refining the grid, this increase is very much unwelcome; it will eventually lead to numerical instabilities in the simulation as energy is added to the system by numerical errors, and thus to undesired erroneous dynamic behaviors. On the other hand, the use of the proposed MP method with the first order scheme shows values of momentum that is much smaller, proving a reduced amount of transfer. When using high-order schemes, the total momentum in the air phase is much more acceptable and the MP method is always lower than the standard method, as much as 2 times smaller with WENO 3, 2 for the coarsest mesh. On the other hand, we see in sub-figure 7c that the total water momentum is not properly converging when using the standard method while the MP method exhibits (in 7d) a very good convergence by the bottom (i.e. with increasing value towards convergence when the mesh is refined). The reader can note that the most precise WENO 5, 3 gives qualitatively the same result for both methods, which is one expected property of WENO schemes, made to behave well in presence of discontinuities (hence near the interface). The mixed ad hoc schemes give globaly satisfactory results. However, we can see in the water phase that even though the (O1, WENO5) coupling reduces damping compared to (O1, O1); still, using third or fifth order WENO at the interface is significantly less dissipative, highlighting a higher capacity of these latter schemes to preserve momentum.

Furthermore, when comparing the momentum in figure 8 and vortices distribution in figure 9, we clearly see that the standard method leads to substantial spurious currents that perturb both the bulk of air and water phases. Even worse, the underwater vortices get more dreadful as the mesh is refined, exhibiting increasing numerical errors. On the other hand, the MP method permits attaining converging results as the mesh size becomes smaller. Moreover, the use of high-order schemes such as the WENO 5, 3, even across the interface, brings more precision in the results by capturing fine structures as expected, which is a very p. 22 important criterion for the accurate simulation for complex simulations (as presented in the next sections).

Those results all show that whatever the spatial scheme used for the advection terms, one can expect a very good reduction of momentum transfer between phases when using the proposed approach. Hence, they prove that the MP method and its high-order version is very well suited for coupling with any classical advection scheme without the need of particular treatment near the interface, and can be plugged into an existing code without too much effort. (a) Air phase total momentum using the standard method. (c) Water phase total momentum using the standard method. (d) Water phase total momentum using the MP method.

Figure 7: Standing gravity wave: spatial convergence of the total y momentum as a function of time for LS+NSSP3, 2 using various advection schemes, with and without MP method, after around 70 wave periods. With the standard method, the total momentum is less well conserved. At the contrary, the HOMP method well captures the momentum in each phase and shows very good numerical convergence. The use of the first order mixed method reduces the overall damping of the wave.

Temporal convergence. We have studied the impact of the time step on numerical errors. For that purpose, we have used the LS+NSSP3, 2+WENO5, 3 with HOMP method and fixed the mesh with S = 4. The different time steps have been set to δ T S = 0.05 respectively. The simulation with T S = 1 has p. 23 already been presented in the previous spatial convergence paragraph and serves as reference. Figure 10b shows, as expected, a very good temporal convergence; for very small time steps, the residual error is due to spatial discretization. The exponential functions drawn in the plot are defined as M 0 e -(t-t0)µ damping where M 0 is the maximum momentum value for the first period at time t 0 and µ damping is the associated measured damping coefficient, as summarized in table 1. We see that, even after an important number of time steps, for long time, the temporal error is very small and the results given at T S = 2 (hence σ = 1 8 ) gives a sufficiently good result compared to smaller time steps. Hence, it seems that it is not necessary to choose a CFL number below 1/10 while a CFL number around 1/4 still gives very acceptable results, making the method very reliable even with big time steps. p. 24 1). Mixed scheme (O1, WENO5) on top and (WENO5, WENO5) HOMP on bottom. Over time, the former damps the gravity wave noticeably more than the later. µ damping 1.7 • 10 -2 5.5 • 10 -3 1.9 • 10 -3 8.2 • 10 -4 2.9 • 10 -4

Order -1.6 1.5 1.2 1.5

Table 1: Standing capillary wave: parameters for the envelop exponential decay (in the form M 0 e -(t-t 0 )µ damping ) of measured total y momentum, as drawn in figure 10a. Results for the mixed (O1, WENO5) and (WENO5, WENO5) HOMP scheme shows a valuable gain for the latter.

Interface methods. We now propose to validate the reliability and flexibility of the HOMP method by comparing the results obtained with LSM, VOF-PLIC and MOF interface transport methods as presented in section 3.2.3. The two later are so-called sharp methods as they rely on an explicit geometric discretization of the interface, hence localizing the volume fraction transition region in only one cell. The interface geometry is advected by the underlying fluid and permits, then, reconstructing the volume fraction at the next time step. MOF is known to be more precise than VOF-PLIC, while being more computationally intensive. A smoothing procedure is applied to VOF-PLIC and MOF volume fraction as to relate closely to the regularized Heaviside function associated to the LS. These three interface methods are conveniently available in Notus and usable without any more effort than changing the case numerical parameters. Hence, every other parameters are set identical to the previous simulations with LS. We show in figure 11 the convergence results obtained with VOF-PLIC, MOF and LSM (for reference, the same as in previous paragraph) using the HOMP method coupled with the WENO 5, 3 scheme. While VOF and MOF both methods show quite similar results in early times (i.e. after a few periods), we observe for longer time that MOF dissipates much less than VOF-PLIC. This concurs to the fact that MOF has better precision than the latter method. At the opposite, the use of VOF-PLIC exhibits more damping in the momentum for long time which reflects a characteristic of a low-order method. Regarding MOF, we observe a small increase in momentum over time which might be explained by the important difference between ρ n+1 and ρ * induced by the volume fraction prediction based on an underlying smooth representation and the sharp advection of the interface, a behavior that is typical to high-order methods. This anti-diffusive behavior can also explain the better result of MOF against LSM for the coarsest mesh.

It is worth to note that, within the LSM framework, small jumps in the momentum might be observed at regular intervals. This artifact is due to the re-initialization procedure that is applied to the LSM every 10 time steps. Without the re-initialization, the LSM equation is very close to the continuity equation, i.e. ρ n+1 ρ * , particularly if they share the same numerical schemes, as we did using WENO 5, 3 scheme. This makes the LSM particularly convenient for two-phase flows. However, during time, the LSM is contracted or spread (i.e. caused by variations of ∇φ) by the underlying fluid, so is the volume fraction/density field. When the LSM is reinitialized at the end of a time step, the associated regularized volume fraction/density ρ n+1 (computed through the Heaviside function) is instantaneously and harshly modified (sharpened or spread). It thus departs much more from the predicted ρ * than if no re-initialization had been done. Indeed, one could apply the re-initialization procedure more often in order to spread the numerical error due to the correction more frequently; nevertheless, this would not remove all accumulated error as the deformation is inherently induced by the underlying fluid and the smooth interface assumption. p. 26 Figure 11: Standing gravity wave: spatial convergence of the total y momentum as a function of time for VOF-PLIC, MOF and LS, using NSSP3, 2+WENO5, 3 with HOMP. As expected, the lower order VOF-PLIC is more diffusive while MOF and LSM give similar results.

Standing capillary waves 4.3.1. Description

This verification case was designed in order to confirm the results shown on the previous case for interfacial flows at very small scales, i.e. at the capillary scale. Hence, it is very similar to the previous one, while we focus here on surface tension forces rather than gravitational forces that have thus been removed. In order to reflect the eventual viscous effects at fine mesh scales and to test the method in a general framework, we have set viscous forces. Unlike gravity which acts as a volumetric force, in the one-fluid model, surface tension is singular at the interface and thus interacts differently with the flow. Its discretization is described here-after.

Physical parameters, initial and boundary conditions. The domain has been down-scaled by a factor of 10 Numerical methods and surface tension. Concerning the computation of ρ * and the inertial term, we have used the most precise WENO 5, 3 scheme coupled with 2 nd order (NSSP3, 2) Runge-Kutta integration. The interface method is LSM with a regularization parameter of 3 • δx. The surface tension forces is modeled in the momentum equation in a CSF manner. Following the conclusions of the study done by Denner et al. [START_REF] Denner | On the convolution of fluid properties and surface force for interface capturing methods[END_REF], we find that it is important, as well as having smooth density and viscosity fields, to use a sharp representation of the characteristic function χ for computing surface tension forces such that the volumetric external force added to the momentum equation is:

f = γκ∇χ sharp p. 27
where χ sharp is computed with the same procedure as described for the heavy droplet case, thanks to the sharp method described in [START_REF] Min | Robust second-order accurate discretizations of the multi-dimensional Heaviside and Dirac delta functions[END_REF], and κ is the curvature computed with 4 th order precision in the normal direction as prescribed in [START_REF] Coquerelle | A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces[END_REF].

Spatial discretization and integration time step. Capillary waves celerity, by linear theory, when neglecting viscous effects, scale as λ -1/2 ; it can be approximated for deep waves by: c γ ρ 2π λ . We those expect velocity magnitudes of the order of v = O (0.2) m • s -1 . In practice, the velocity in the water phase is less than the celerity, and we have observed that |v| max ∼ 0.16 m • s -1 , even in the air phase. Additionally, surface tension methods are subjected to time step constraint more restrictive than advection CFL. It scales as δx 3/2 rather than δx 1 . However, we still would like to use the larger time step as possible for evident computational efficiency. Large time step induce more temporal error, and we expect the HOMP method to perform, as proven on the gravity wave, much better than the standard approach.

We have applied a regular and uniform grid for the discretization of the domain with N × 2N cells such that δx = L N . As previously, for ease of convergence study, we denote by S the scale of the simulation such that N = N λ = 32 • S. For all the following simulation, we have chosen δt = 10 -4 /S such that at the maximum, the CFL number is σ max ∼ 1/20.

Results

Figure 12 shows the resulting density-vorticity (i.e. ρω) fields for varying spatial scales. Qualitatively, we clearly see that, with the standard method the quantity of erroneous underwater vorticity is much higher than when using the HOMP method. Under the crest, for the former, we observe a much larger zone (related as pocket in figure 12) of high spurious vorticity under the crest (a similar pocket can be observed under the trough). On the other hand, for the latter, the pocket is less spread and the maximum value herein is also reduced by approximately half, as specified in table 2. Furthermore, under the trough, we also observe a much better convergence, with HOMP, towards the finer mesh where the vorticity is confined closer to the interface. It is important to note that, as the mesh is refined, small spurious currents and vortices induced by the momentum transfers -as observed in the gravity wave case -are damped by viscous effects, non-negligible at that spatial scale. This explains the reduction of maximum vorticity and convergence of the simulations, even with the standard method; at finer meshes (S > 4), as shown in the right-most results in figure 12, these spurious vorticity pockets eventually disappear for all methods.

Plunging breaking wave 4.4.1. Description

To study the breaking wave, the standard method is still commonly used with VOF-PLIC [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Lubin | Numerical simulations of three-dimensional plunging breaking waves: generation and evolution of aerated vortex filaments[END_REF] or LSM [START_REF] Iafrati | Modeling of ocean-atmosphere interaction phenomena during the breaking of modulated wave trains[END_REF]. We will investigate the use of the proposed method on this validation case to show the influence of using the proposed method for simulating a complex multiphase flow on fine mesh. The Bond number of the wave simulated is Bo = ∆ρg γk = 34 005 and the wave steepness is: Ak = 0.45. Deike et al. [START_REF] Deike | Capillary effects on wave breaking[END_REF] have created a diagram of the regime of breaking wave in deep water ( d λ ≥ 0.5) in function of the Bond number and the wave steepness. For this set of parameters, the breaker type is plunging. The wave starts to break, a plunging jet is generated which will impact the surface and turn into a roller. We will focus the results on the first steps of the breaking event, until the jet generation. During the jet generation, a high shear arises above it. As a high shear is generated, this case is a good demonstration of the limit to use the standard method with a very fine mesh and the influence of using a high order scheme such as WENO 5, 3 for a coarse mesh.

Physical parameters, initial and boundary conditions. The periodic sinusoidal wave is initialized regarding the first order stokes wave (airy) equations:

η 0 (x) = A • cos(kx) u 0 (x, y) = α cosh(k(y+d))
sinh(kd)

• cos(kx)

sinh(k(y+d)) sinh(kd)

• sin(kx) with η 0 the initial wave surface height, k = 2π λ the wave-number, λ the wave length, u 0 = (u 0 , v 0 ) the initial velocity inside the water, d the constant wave depth, and α = A (gk + γ ρwater k 3 )tanh(kd). The wave is initialized in 2D box of size L × H = 1 m × 1 m. The boundary conditions are periodic in the x direction and wall (no-slip) on top and bottom. The wall on top has been set such that there is no influence on the breaking event. The wave is initialized as a sinusoidal function of amplitude A = 0.072 m, wave length λ = 1.0 m and wave depth d = 0.5 m . The surface tension coefficient is taken as γ = 73 • 10 -3 N • m -1 and the constant gravity acceleration is taken as

g = 9 • 81 m • s -2 .
Spatial and temporal discretization. The simulation is performed with a WENO 5, 3 scheme coupled with NSSP 3, 2. An adaptive time step procedure is used by requiring the CFL number to be σ = 0. The smoothing procedure (presented in section 3.2.3) is done 3 times; hence, the effective interface thickness between air and water is reduced linearly with the spatial step.

Numerical method for surface tension. The simulations were carried out using the height function method to compute the curvature within the VOF-PLIC framework. The implemented algorithm is based on a hybrid method: it consists in using the initial method developed by Popinet [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] for well represented interface and, otherwise, a local mesh rotation is performed as presented in [START_REF] Owkes | A mesh-decoupled height function method for computing interface curvature[END_REF]. p. 29

Results

When a plunging breaker wave starts its breaking, a jet is generated in the front side of the wave. During this generation, a high shear is located at the top of the jet. Figure 13 shows a close view of the breaking event. At t = 0.612 s., the jet is starting to plunge, falling under gravity, before impacting the free surface and capturing a roller. The validation of the method and its accuracy is demonstrated through the comparison of the breaking event for different meshes. The simulations performed with the proposed method (HOMP), using WENO 5, 3, demonstrate convergence as the mesh is refined. On the other hand, the standard method, using WENO 5, 3, generates non-physical instabilities where the highest shear is located. These instabilities are increasing with time until the simulation blows up around t = 0.62 s, and are even worse when refining. These spurious currents are generated by transfers of momentum across the interface, as was illustrated previously on the case of the standing wave. Figure 14 shows the comparison of the mixed MP (O1, WENO5) and HOMP (WENO5, WENO5) for each mesh at the same time. A reference solution is shown in a black dashed line, corresponding to the surface profile using HOMP on a very fine mesh. The behavior of the wave surface is becoming noticeably different as soon as the wave starts to break. The jet is highly influenced by the numerical scheme used near the interface. Using the mixed MP (O1, WENO5) modifies the dynamic of the jet during the breaking process: it is much thicker than the reference solution. This is effect is due to more important diffusion of the first order numerical scheme used near the interface, as discussed in previous sections. By contrast, using a WENO 5, 3 scheme near the interface improves the jet profile: even on the coarsest mesh, there is a good adequation with the reference solution, exhibiting a thinner width. However, as expected, both methods converge to the reference solution as shown in figure 15. These results show the influence of the underlying numerical scheme used at the interface, poiting out the relevance of using (WENO5, WENO5) for capturing the overall dynamics, more especially for coarse meshes. p. 30

Drop impact 4.5.1. Description

Finally, we study the capacity of the method to capture correctly the dynamic of the complex case of the impact of a drop over a flat surface of water. This problem has been widely studied in the literature, and we refer the reader to Liow [START_REF] Liow | Splash formation by spherical drops[END_REF], Morton et al. [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] and Cole's thesis [START_REF] Cole | The Splashing Morphology of Liquid-liquid Impacts[END_REF] for detailed experiments at increasing Froude (Fr ≡ U 2 gD ), Weber (We = ρ l U 2 D σ ) and Reynolds (Re = ρ l U D µ l ) numbers, where D is the droplet diameter, U its mean velocity before impact and subscript l stands for the liquid's property.

In order to validate our simulations, we have reproduced Morton et al. validation case B (corresponding to Fr = 220, We = 248 and Re = 7250) that shows a thick jet formation and Cole's experiment number III (corresponding to Fr = 67, We = 62 and Re = 3082), referred to a case of "formation of a primary vortex ring / pre-entrapment jetting transition". This case exhibits very interesting phenomena: the formation of a small mushroom-like cavity at early stage, followed by an intermediate cavity with sharp edges that eventually collapses towards a pinch. Afterwards, an air bubble is captured below the surface while a thick jet rises. The use of a dye permits to highlight the formation of a vortex ring that pushes the bubble underwater. Those phenomena are illustrated in the close up view of our simulation in fig. 18c.

It is important to note that the various phenomena that arise after the impact of the droplet on the pool are complex and very sensitive to initial conditions (as well as uncertainty in the experimental measures). Surface tension and the capturing of very small structures (of the order of 10 -5 m) that trigger some consequent phenomena is a key factor for accurate simulations. For these reasons, this application is still today a challenging problem for numerical computation and its detailed analysis beyond the scope of this article.

p. 31

We show in this section the capacity of the HOMP method to avoid spurious momentum and vorticity that permits to obtain high-quality results.

Physical parameters, initial and boundary conditions. Using the same physical domain size as in the experiments isn't necessary as the deep pool hypothesis is largely valid with a smaller domain. Also, the impact of the lateral walls can be neglected in the early stages of the experiment. Also, the dynamic of the impact and its consequence is mostly axisymmetric. Hence, we have taken a box of size L × H = 9D × 18D with axisymmetric condition, no-slip on the bottom, slip on the right and a Neumann condition on the top boundary. This permits to let air freely flow as in the experiment while confining the liquid inside the computational domain. The liquid pool is initialized at the vertical center of the domain, i.e. at y = 0, hence the depth of the pool is 9D. An illustration of the setup is given in figure 16.

Morton's case B (referred by superscript B) and Cole's case III (referred with superscript III) physical parameters are summarized in table 3. These choices lead to slightly different numerical values for Fr, We and Re but this hasn't showed to give noticeable difference in the dynamic of the impact while corresponding more closely to the experimental fluid property measures. Gravity is set to g

= 9.81 • m • s -2 .
While the initial condition of the fluid's velocity inside and near the droplet is a complicated problem, after several tests, we have found that it wasn't necessary to position it too high above the pool, hence we have set y 0 D = 2D. Furthermore, we have initially imposed the velocity U inside the droplet and, in order to alleviate any discontinuity matter, have linearly reduced the velocity in the surrounding air such that:

v y =      U if r ≤ R U 1 -r-R 0.4R else if r ≤ 1.4R 0 otherwise
where v y is the vertical component of the velocity and r is the distance to the droplet's center.

Case B from [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] Spatial and temporal discretization. As the dynamic of the flow is less important far from the impacting zone and the cavity, in order to save computational time, the physical domain has been discretized with a grid p. 32 with a uniform spacing centered around the zone of interest of size 3D × 6D and an exponentially increasing spacing farther away, as depicted in figure 16. In the uniform zone, the spatial discretization is such that δ r = δ y = 6D N where N (resp. 2N ) is the total number of cells in the r (y) direction. In order to be able to capture accurately the very small physical phenomena that particularly arise near the interface (e.g. the mushroom-like cavity, the pinch, the entrapped bubble and secondary droplets) we have chosen N = 1024 (i.e. the mesh consists of 1024 × 2048 = 2 097 152 cells). This corresponds to a number of cells per diameter of δ r /D 170, and hence gives δ B r = 1.7 • 10 -5 m for case B and δ III r = 1.541 • 10 -5 m for case III. The choice of the time step is very restrictive at those small spatial scales as capillary waves can travel very fast. In order to avoid any instabilities, following the prescriptions in [START_REF] Denner | Numerical time-step restrictions as a result of capillary waves[END_REF], we have imposed step constraint in the form:

δt ≤ σ δ • c γ + u •
where σ is the desired Courant number, • denotes the most restrictive spatial direction (i.e. here r or y), c γ = 2π (ρ l +ρg)δ• is the fastest capillary wave captured by the mesh and u • is the flow velocity. This corresponds to the classic CFL condition where the considered velocity is the maximum velocity of a fluid particle at the interface where the smallest wave can eventually arise. However, at the scales of the order of 10 -5 m, viscous forces start to play a non-negligible role and diffuse capillary waves that are quickly damped. Even though we have found that using time steps as large as δ B t = 1.27 • 10 -6 • s and δ III t = 2.3 • 10 -6 • s have given very satisfactory results, for the sake of rigor, we have preferred to follow the CFL condition with σ = 1 with the prescribed constraint.

Numerical methods. For capturing accurately the surface tension effects that are a key and triggering point of most of the expected phenomena, we have chosen to use the LSM and a WENO 5, 3 + NSSP 3, 2 scheme for integrating the advection equations. The regularization parameter associated to the LSM interface thickness has been set to = 1.5 • δr and we have opted for the same numerical method for computing surface tension forces as for the capillary standing wave case 4.3.

Results

Morton et al. case B.

Figure 17 summarize the results obtained for the simulation of [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] case B. In fig. 17a, without the use of the HOMP method, the very early stages of the formation of the cavity are very much different: shortly after the impact (around t = 1 ms) the thin film is less well captured, the crown has less energy and thus grows less high while spreading more. The reader can also remark a non-negligible difference in the vorticity field at the top of the droplet on the first image that is later on transferred at the center of the cavity; the HOMP method captures an important layer of negative vorticity that is due to the important inertia inside the drop. Those differences lately impact the overall dynamic of the simulation that make the standard method to fail to capture accurately the process.

Using the HOMP method, the resulting cavity and jet formation are qualitatively compared to the experiments from [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] in fig. 17b: we observe a very good match of the simulation regarding the shape of the cavity. Here, the dimensionless time is obtained as T = (t simt impact ) U/D and the dimensionless lengths as L = l/D (where l is either the cavity depth or the jet height). The measured cavity depth plotted in fig. 17c is following the same curves as the literature's results. However, starting from T = 15, a small capillary wave has converged towards the center of the cavity while not disappearing: this causes the noteworthy formation of a small teat (illustrated in fig. 17d) for a very short period of time, where the surface profile induces a very high curvature captured by our simulation that cannot be captured with coarser meshes. That phenomenon makes the measure of the water level go back upwards then downwards in fig. 17c; measuring the surrounding water plateau depth recovers the same trend as reference results.

Regarding the jet, in our simulation we have observed a higher maximum height of around 5.8 before it releases a first droplet at T = 41. The release time is in good agreement with the simulations from [START_REF] Ray | Regimes during liquid drop impact on a liquid pool[END_REF] while it is later than in the experimental results. Herein, the reader can note a possible mistake in the reported experimental and simulation times of [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] as the last picture (see 17b top row, last column) T = 41 contradicts in the experimental measures (see 17c, filled disks) for the discontinuity around T = 37. Nevertheless, the diameter of the secondary droplet has been measured to D = 2.6 mm, hence a ratio of 0.897 regarding the initial one, which is of the same order as reported by previous works. Because of existence of several capillary waves along the jet, shortly after the secondary droplet, we observe the release of a third one.

p. 33

The difference in the simulations can be explained by the use of our more momentum conserving method coupled with high-order schemes for the interface representation, inertial and capillary forces. Also, our mesh is around 3.4 times finer than the one used by Ray et al, which permits capturing smaller structures.

Cole case III. The numerical results shown in fig. 18 compare very well to the pictures taken by Cole in [START_REF] Cole | The Splashing Morphology of Liquid-liquid Impacts[END_REF]. Similar disparity has been observed for the standard method as in the previous case, justifying again the necessity to use the HOMP method with WENO 5, 3.

A dye has been added inside the droplet in order to trace the initial liquid and exhibit the vortices. In the last picture of 18b, at t = 18 ms, we see that the entrapped bubble is pushed towards the bulk while the jet is going upwards. In our simulation, the captured bubble has a diameter of approximately 0.5 mm (∼ 0.3D) which is twice smaller than in the experiment. This difference is not surprising as this bubble entrapment phenomenon is, as explained previously and in the literature, very sensitive to initial conditions and error measurements, as well as numerical errors. When the pinch has collapsed, a very strong air current is flowing out of the cavity towards the top at a velocity of 16 m • s -1 . Two small droplets of diameter 0.2 mm and 0.3 mm (see close up view of fig. 18c) are also rapidly ejected in the early stage of the jet formation, around t ∼ 15 ms at a velocity of approximately 7.3 m • s -1 and 1.65 m • s -1 . This very small phenomena has also been observed in [START_REF] Cole | The Splashing Morphology of Liquid-liquid Impacts[END_REF].

Again, the advantage of the HOMP method with WENO 5, 3 is clearly illustrated in fig. 18a (wherein we used a constant time step for the sake of comparison). We see that the standard method fails to capture the mushroom-like initial cavity (left images, at t = 3.27 ms) and the sharp shape of the cavity (at t = 10.7 ms) that leads to the entrapment of the bubble (t = 13 ms then t = 17 ms). The mixed method using (O1, WENO5) also fails to capture important key features of the dynamic and the following early capture of a bubble. We believe that this is due to a more important diffusion of the momentum at the interface before the impact, when inertial forces are largely dominant. Around t = 13 ms, the bottom of the cavity collapses and capture correctly the large bubble but, later on, as the surrounding vortex has not formed appropriately earlier, the jet grows faster while the expected underlying vortex is not present.

Conclusion

We proposed a discrete framework for the preservation of momentum in incompressible two phase flows with high density ratio that enables the use of accurate and robust advection schemes, up to high-order WENO 5, 3 coupled with second order Runge-Kutta temporal integration. The HOMP method can thus be almost effortlessly integrated into an existing code in 2D and 3D and essentially requires the solving of a supplementary advection equation for the mass prediction, consistently with the inertial term of the momentum equation, as now acknowledged in the domain [START_REF] Nangia | A robust incompressible Navier-Stokes solver for high density ratio multiphase flows[END_REF]. This strategy greatly reduces spurious exchanges of momentum across the interface, problem which the standard one-fluid model is suffering, and yields to converging solutions with mesh refinement, unraveling a convergence problem that has been observed in the literature.

Furthermore, the MP approach is independent on the underlying interface representation, following early works of [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF][START_REF] Ghods | A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods[END_REF]; we have demonstrated that it works convincingly when coupled with VOF-PLIC, MOF and LSM, and particularly benefits from the proposed High-Order MP improvements as exposed in the various results, without necessitating any particular algorithm as opposed to [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF][START_REF] Fuster | An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF][START_REF] Owkes | A mass and momentum conserving unsplit semi-Lagrangian framework for simulating multiphase flows[END_REF][START_REF] Ghods | A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods[END_REF]. We successfully tested the method on selected numerical validation cases, including gravitational and capillary waves, that showed very satisfactory results even with relatively coarse meshes and large time steps. Additionally, we successfully confronted our approach to complex applications at large scales -with the breaking wave -as well as small scales -with the droplet impact on a liquid pool. Our simulations show good agreement to literature results and the expected dynamic of the flows, where the non preserving method failed and low-order schemes are more dissipative, and this without the necessity of a specific, complex or costly numerical strategy other than the underlying advection scheme. Moreover, the use of high-order advection schemes even near the interface eases the overall consistency and accuracy of the two phase Navier-Stokes solver. p. 34 We believe that any kind of mass / momentum advection coupled scheme can work well within our method as long as they are consistent with each other. Also, other interface representations such as front-tracking could also benefit from the HOMP method. We are considering working on specifically designed advection schemes that would take full advantage of the known characteristic function profile in order to increase the accuracy at the interface, where the density field becomes steeper as the mesh is refined. The three left most plots are using the standard method at scales S = {1, 4, 16} while the three right most plots are the results using the MP method with the same meshes. Using a finer mesh permits to reduce momentum diffusion; however, for the standard method, spurious momentum appears underwater while the MP method converges well for all simulations. The results using the mixed method are not shown as being very close, regarding the momentum, to the non mixed version.

CRediT

p. 39 with three different advection schemes, including the ad hoc mixed method. For all methods, scales are S = {1, 4, 16}. For the standard method, we observe intense underwater vortices that are even more apparent when refining the mesh, while the HOMP method converges well for all simulations. Results with the ad hoc mixed method show less diffusion because of the use of WENO 5, 3 in the bulk. However, the MP(O1, WENO5) scheme shows more damping at lower resolution while the HOMP(WENO3, WENO3) and particularly HOMP(WENO5, WENO5) schemes conserve a more dense vorticity near the interface.

p. 40 (b) Comparison of the observed free surface profiles at various time. Top row is extracted from [START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF] while the bottom row is for our simulation. Be aware that the point of view of the photos can hide the shape of the cavity, particularly for the crown; also, the photographs have different zoom levels. For the simulations, only the last image has been zoomed out in order to capture the full jet with the almost released secondary droplet. p. 42

Figure 2 :

 2 Figure 2: Algorithm comparison between Momentum Preserving (left) and standard method (right).

  (a) First order upwind scheme. (b) WENO 3, 2 scheme. (c) WENO 5, 3 scheme.

  (a) First order upwind scheme. (b) WENO 3, 2 scheme. (c) WENO 5, 3 scheme.

Figure 4 :

 4 Figure4: Convection of a high density droplet with reducing interface thickness: comparison of the numerical convergence of the shape and x momentum for various numerical and interface methods. Time is t = 1 s. For each simulation, the top half part is the result using standard method, the bottom half one with HOMP method. Purple filled color is for momentum up to 10 6 (which is the momentum inside the drop at initialization), white is for ρu = 0. The contours highlight particular values: black for the interface position, blue (in the gas phase) for ρu = -10 and red (inside the droplet) for ρu = 10 6 .

  = 64 standard, N = 128 standard, N = 256 standard, N = 512 HOMP, N = 64 HOMP, N = 128 HOMP, N = 256 HOMP, N = 512 (c) MOF.

Figure 5 :

 5 Figure 5: Convection of a high density droplet with reducing interface thickness: comparison of the numerical convergence of the total kinetic energy E k Ω against time for the three studied advection (not mixed, i.e. Θ Γ ≡ Θ b ) schemes and interface representations, and water phase volume associated to the WENO 5, 3 scheme results (right most column).

  A 2D horizontal wave is initialized in a rectangular box of size L × H = 10 m × 20 m. The boundary conditions are periodic in the x direction and wall (no-slip) on top and bottom. The wave is initialized as a sinusoidal function of amplitude A = 0.05 m and wave length λ = L = 10 m centered in the domain. The heavier fluid (water) is of density ρ 1 = 1000 kg/m 3 and the lighter one (air) is of density ρ 2 = 1 kg/m 3 . Gravity is set to g = 9.81 m/s 2 towards the negative y direction.

4

 4 (a) 3D interface representation in blue and x momentum isocontours on the plane x, y. Time is t = 1 s and N = 128. Top row is for the standard method and bottom row using HOMP; left column is for LSM and right column for VOF. The 2D contours highlight particular values for the momentum: ρu = -0.1 (a hundred times less than in 2D) and red (inside the droplet) for ρu = 10 6 . Total kinetic energy standard, N = 64 standard, N = 128 standard, N = 256 HOMP, N = 64 HOMP, N = 128 HOMP, N = 256 (b) Convergence of the total kinetic energy E k Ω as a function of time for the LSM method. The initial value has been plotted with a light gray dotted line.

Figure 6 :

 6 Figure 6: Convection of a high density 3D droplet: comparison of the numerical convergence of the shape and x momentum integral, using a WENO 3, 2 + NSSP 3, 2 RK scheme.

  Air phase total momentum using the MP method.

  Spatial convergence in the interval t ∈ [0, 100] s. The associated envelopes for periodic minimum and maximum values of are shown as colored decreasing exponential functions (whose parameters are given in table

  Close up in the interval t ∈ [7.2, 7.4] s for the temporal convergence at scale 4 with decreasing Courant number.

Figure 10 :

 10 Figure 10: Standing gravity wave: total y momentum as a function of time for LS+NSSP3, 2+WENO5, 3 with HOMP method. The method clearly converges towards a solution capped by spatial discretization errors.

  p

  VOF-PLIC, MOF and LSM (already presented) results after around 70 periods. Comparing the three interface methods at scale S = 16.

  -3 : a 2D horizontal wave is initialized in a box of size L × H = 10 mm × 20 mm. The boundary conditions are periodic in the x direction and wall (no-slip) on top and bottom. The wave is initialized as a sinusoidal function of amplitude A = 0.2 mm and wave length λ = L = 10 mm. The relative amplitude has been increased compared to the previous case in order to match more closely to capillary scale applications such as shown in the next sections. The surface tension coefficient is taken as γ = 73 • 10 -3 N • m. Density and viscosity are set to ρ l = 10 3 • kg • m -3 and µ l = 10 -3 • P a • s (resp. ρ g = 1 • kg • m -3 and µ g = 2 • 10 -5 • P a • s) for the liquid phase (resp. gaseous phase).

Figure 12 :

 12 Figure 12: Standing capillary wave: spatial convergence of density-vorticity ρω, using LSM+NSSP3, 2+WENO5, 3 with the standard method (top row) and with HOMP method (bottom row), at time t = 0.2 s. Meshes are, from left to right, obtained with: S = {1, 2, 4, 16}. Spurious vorticity pockets (located by white circles) are particularly observable at coarser meshes, notably under the crest; as expected they are much smaller and disappear faster with the HOMP method. The color map goes from -10 4 • kg • m -3 • s -1 (blue) to +10 4 • kg • m -3 • s -1 (red).

  3. We have applied a uniform grid discretization of the domain with 512 × 512 = 262 144 cells for the coarse mesh, 1024 × 1024 = 1 048 576 cells for the medium mesh and 2048 × 2048 = 4 194 304 cells for the fine mesh. The spatial step are respectively δx coarse = 1 512 = 1.95 • 10 -3 m, δx medium = 1 1024 = 0.98 • 10 -3 m and δx f ine = 1 2048 = 0.49 • 10 -3 m. A reference solution is given on a very fine mesh with δx ref = 1 4096 = 0.24 • 10 -3 m.

Figure 13 :

 13 Figure13: Plunging breaking wave: comparison of the breaking dynamics with the standard and proposed HOMP method using VOF-PLIC interface representation for every meshes (coarse, medium and fine). The interface thickness is reduced linearly with the spatial step.

Figure 14 :Figure 15 :

 1415 Figure 14: Plunging breaking wave: comparison of the surface profil with mixed MP (O1, WENO5) (red line) and with HOMP (WENO5, WENO5) (blue line) using VOF-PLIC interface representation with WENO 5, 3 for coarse (left), medium (center) and fine (right) meshes. The black dashed line corresponds to the reference solution calulated with HOMP.

Florian Desmons:

  Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing -Original Draft, Visualization. Mathieu Coquerelle: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing -Original Draft, Visualization, Supervision.

  (a) First order scheme. (b) WENO 3, 2 scheme. (c) WENO 5, 3 scheme.

Figure 8 :

 8 Figure 8: Standing gravity wave: visualization of y momentum for LS+NSSP3, 2 at time t = 48 s (around 47 periods) with three different advection schemes.The three left most plots are using the standard method at scales S = {1, 4, 16} while the three right most plots are the results using the MP method with the same meshes. Using a finer mesh permits to reduce momentum diffusion; however, for the standard method, spurious momentum appears underwater while the MP method converges well for all simulations. The results using the mixed method are not shown as being very close, regarding the momentum, to the non mixed version.

  (a) First order: standard method (left), MP(O1, O1) (center), MP(O1, WENO5) (right). (b) WENO 3, 2: standard method (left), HOMP(WENO3, WENO3) (center), MP(WENO3, WENO5) (right).

Figure 9 :

 9 Figure 9: Standing gravity wave: visualization of density of vorticity ρω for LS+NSSP3, 2 at time t = 48 s (around 47 periods)with three different advection schemes, including the ad hoc mixed method. For all methods, scales are S = {1, 4, 16}. For the standard method, we observe intense underwater vortices that are even more apparent when refining the mesh, while the HOMP method converges well for all simulations. Results with the ad hoc mixed method show less diffusion because of the use of WENO 5, 3 in the bulk. However, the MP(O1, WENO5) scheme shows more damping at lower resolution while the HOMP(WENO3, WENO3) and particularly HOMP(WENO5, WENO5) schemes conserve a more dense vorticity near the interface.

  Comparison of density-vorticity without (left side of each snapshot) and with HOMP (right side of each snapshot). The interface is drawn with a thick black line and the color scale goes from -10 -5 kg • m -2 • s -1 (blue) to 0 (white) to+10 -5 kg • m -2 • s -1 (red). From left to right: snapshots at t = {1.02, 2.04, 5.1, 10.6} ms.

  Cavity depth and jet temporal evolution compared to[START_REF] Ray | Regimes during liquid drop impact on a liquid pool[END_REF][START_REF] Morton | An investigation of the flow regimes resulting from splashing drops[END_REF]; all values are adimentionalized. Our results are in filled squares, Ray et al. in open circles, Morton et al. simulation in open triangles and their experiment in filled disks. (d) Cavity retraction phenomenon: 2D free surface profile (thick black line) with density-vorticity coloring (same scale as in fig. 17a) overlapping the semi transparent 3D revolution.A small teat appears at the center of the cavity.

Figure 17 :

 17 Figure 17: Drop impact: Morton et al. case B simulation: qualitative comparison and quantitative results.

  (a) Comparison of density-vorticity. Top line: with standard method (left side of each snapshot) and with HOMP (WENO5, WENO5) (right side of each snapshot). Bottom line: with mixed MP (O1, WENO5) (left side of each snapshot) and with HOMP (WENO5, WENO5) (right side of each snapshot, same as above). The interface is drawn with a thick black line and the color scale goes from -10 -5 kg • m -2 • s -1 (blue) to 0 (white) to+10 -5 kg • m -2 • s -1 (red). From left to right: snapshots at approximately t = {3.27, 10.7, 13, 18} ms.(b) Free surface profile from our simulation (bottom rows, turquoise color) compared to pictures extracted from[START_REF] Cole | The Splashing Morphology of Liquid-liquid Impacts[END_REF] (top rows, gray-scale) at the same instants. A dye is used to trace the droplet's liquid. Focus on the free surface profile and droplet dye from our simulation at time t = 15.18 ms and illustration of the different phenomena following the impact. Here, fast (with velocity superior to 5 m/s) secondary droplets are ejected by the jet in formation after the cavity retraction while the entrapped bubble is traveling slowly downards. The dye exhibits the presence of relatively important underwater vortices.

Figure 18 :

 18 Figure 18: Drop impact: Cole case III simulation: qualitative results.

  . 25µ damping 2.8 • 10 -2 1.1 • 10 -2 5.3 • 10 -3 2.5 • 10 -3 1.1 • 10 -3

	Scale (S)		1	2	4	8	16
	(O1, WENO5)	M 0	319	334	354	367.9	375.1
	(WENO5, WENO5)	Order M 0	-328.8	1.3 338.7	1.1 356.6	1.1 369	1.2 375.6

Table 2 :

 2 4 • kg • m -3 • s -1 (red). Standing capillary wave: approximate maximum density-vorticity values (ρω) max at the center of the spurious vorticity pocket under the crest (see figure12for its localization) at time t = 0.2 s, with WENO 5, 3 standard method and WENO 5, 3 HOMP. Units are 10 4 • kg • m -3 • s -1 .

	p. 28

Table 3 :

 3 D is the drop diameter in m, U its mean impacting velocity in m • s -1 , ρ l (resp. ρg) the liquid's (resp. gas) density in kg • m -3 , µ l (resp. µg) the liquid's (resp. gas) viscosity in P a • s and γ the surface tension in N • m -1 .

		9D
	9D
	D 2D	6D → N cells
		2N cells
	9D	3D → N 2 cells
		N cells
	Figure 16: Physical domain and the underlying mesh (with N = 16 for illustration). The uniform part of the mesh has been highlighted in light gray.

Drop impact: parameters (left table) and setup (right figure) of the simulations.

We present in figure1an illustration of the difference between a first order and a higher order temporal approximation of the density field ρ (through the interface position). We there assumed, for simplicity, that p. 5 the density field is sharp, i.e. = 0 and thus the density is constant on each side of the interface Γ. If one assumes that the velocity is constant in this region over space and time, i.e. u = u, the difference of computed momentum ∆ (ρu) is thus linearly proportional to the difference in computed density ∆ρ =ρ n+1ρ * .Γ n Γ n + 1 Γ * Γ n Γ n+1 Γ * Γ n Γ n+1 Γ * Γ n Γ n+1Γ *

Acknowledgements

Graphical in this article were rendered with PGFPlots [34] and simulation visualizations were created using the VisIt visualizing and analysis tool [3]. VisIt is supported by the Department of Energy with funding from the Advanced Simulation and Computing Program and the Scientific Discovery through Advanced Computing Program.

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux -SysNum (ANR-10-IDEX-03-02) and was granted access to the HPC resources of the Mésocentre de Calcul Intensif Aquitain (MCIA) and TGCC under the allocations 2018-A0082A06104 and 2020-A0082A06104 made by GENCI.