Ziqiang Zou
email: zou@limsi.fr

Nicolas Grenier

Samuel Kokh

Christian Tenaud

Edouard Audit

Compressible solver for two

Keywords: low Mach correction, two-phase flow, sharp interface, truncation error analysis, compressible approach, Level Set

Introduction

In the present work, we are interested in flows that involve two phases both represented by two compressible materials separated by an infinitely thin interface. The interface is considered sharp and treated as a discontinuity [START_REF] Hu | On the hllc riemann solver for interface interaction in compressible multi-fluid flow[END_REF]. The Ghost Fluid Method (GFM) [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid method)[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with Level Set/Ghost Fluid method[END_REF][START_REF] Caiden | A numerical method for two-phase flow consisting of separate compressible and incompressible regions[END_REF][START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF][START_REF] Liu | The ghost fluid method for compressible gaswater simulation[END_REF][START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF][START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] is a very popular family of methods that propose discretization strategies for simulating such flows. They rely on a Level Set function to describe the position of the interface [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF][START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. A common approach in the GFM consists of coupling single-fluid numerical schemes like Godunov-type methods across the interface in order to accurately approximate the fluid motion in the vicinity of the interface.

Unfortunately, in the case of single-fluid flows it is now well-established that standard Godunov-type methods may dramatically lose their accuracy when the material velocity is small compared to the sound velocity of the medium and when the computation domain is discretized over a quadrangle (resp. hexahedral) grid in 2D (resp. 3D). This low-velocity regime is referred to as the low Mach regime. This flaw has been extensively studied in the literature and several possible fixes are available (see for example [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Murrone | Behavior of upwind scheme in the low mach number limit: Iii. preconditioned dissipation for a five equation two phase model[END_REF][START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF][START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF][START_REF] Boniface | Rescaling of the roe scheme in low mach-number flow regions[END_REF][START_REF] Xie | An accurate and robust HLLC-type Riemann solver for the compressible Euler system at various Mach numbers[END_REF][START_REF] Pelanti | Wave Structure Similarity of the HLLC and Roe Riemann Solvers: Application to Low Mach Number Preconditioning[END_REF]). This question has also been addressed concerning two-phase flows for diffuse interface models [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustictransport splitting scheme for two-phase flows[END_REF] where the interface is considered as a numerically diffused zone that can reduce the density and sound speed ratios. Comparing with the sharp interface method, the challenge associated with large discontinuities across the interface could be reduced as these discontinuities are smoothed by the mixing layer.

In the case that we consider here, the properties of the medium may experience a severe jump across the interface. Therefore one cannot assume that the orders of magnitude associated with the fluid parameters are the same for both phases. Consequently, it is necessary to define a Mach regime that pertains to each fluid separately. In this context, the classic low Mach fixes are no longer sufficient to propose an accurate discretization of the flow equations when both Mach numbers of each fluid are low.

In the present work, we propose a new low Mach regime analysis that is adapted to compressible flows with sharp interfaces including surface tension effects. Moreover, we also propose a new low Mach correction of the GFM presented in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] that accounts for the jump of characteristic magnitude with respect to density and sound velocity that preserves its accuracy in the low Mach regime. The truncation error is not only uniform with respect to the Mach number as many low Mach corrections for two-phase flows [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF][START_REF] Peluchon | A robust implicit-explicit acoustictransport splitting scheme for two-phase flows[END_REF][START_REF] Pelanti | Low Mach number preconditioning techniques for Roe-type and HLLC-type methods for a two-phase compressible flow model[END_REF], but also uniform with respect to the density ratio and sound speed ratio between both phases.

Before going any further, let us mention that the study of the low Mach regime that we propose in the present work borrows elements from the analysis of the Asymptotic Preserving properties (AP) [START_REF] Cordier | An asymptotic-preserving all-speed scheme for the euler and navier-stokes equations[END_REF][START_REF] Degond | All speed scheme for the low mach number limit of the isentropic euler equations[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equations[END_REF][START_REF] Bispen | Asymptotic preserving imex finite volume schemes for low mach number euler equations with gravitation[END_REF][START_REF] Liu | An asymptotic preserving scheme for the two-dimensional shallow water equations with coriolis forces[END_REF][START_REF] Yan | A successive penalty-based asymptotic-preserving scheme for kinetic equations[END_REF]. However, we would like to emphasize that we do not claim that our method is a genuine AP numerical scheme. Indeed, proposing an AP method would first require to exhibit a clear limit system when the Mach number vanishes in both phases while accounting for discrepancies in the fluid properties and Mach numbers. Such task is beyond the scope of this work: we choose here to focus on studying the behavior of the error with respect to the Mach number and characteristic quantities like the density and the sound velocity of each fluid.

The paper is structured as follows: the overall system is described in Section 2. We briefly recall the numerical scheme proposed by [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] for this model in Section 3 and define the low Mach regime in Section 4 for our interface model. The low Mach behavior of the numerical scheme [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] far from the interface is studied in Section 5, suggesting that the low Mach fixes proposed in [START_REF] Peluchon | A robust implicit-explicit acoustictransport splitting scheme for two-phase flows[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] will enhance the behavior of the scheme far from the interface. In Section 6 we investigate these low Mach schemes in the vicinity of the interface and we expose truncation errors that are driven by the jump of magnitude with respect to both the density and the sound velocity of the material that dramatically affects the accuracy of the low Mach schemes. A new low Mach scheme is then proposed in Section 7 to address this issue. This new numerical method is tested in Section 8 against several 2D tests that show significant improvement in the simulation results over the previous low Mach fixes.

Governing equations for compressible two-phase flows

We consider a space volume D ⊂ R d occupied by two compressible fluids i = 1, 2 that are both equipped with a Stiffened Gas Equation of State (EOS) of the form

(ρ, e) → p EOS i (ρ, e) = (γ i -1)ρe -π ∞ i γ i , (1)
where ρ and e denote respectively the density and the specific internal energy of the material. The parameter γ i is the heat capacity ratio of fluid i, π ∞ i is a constant pressure representing the molecular attraction between molecules of fluid i. The sound speed c i is then defined for the fluid i by

c 2 i = (∂p EOS i /∂ρ) + (p EOS i /ρ 2 )(∂p EOS i /∂e) = γ i p i + π ∞ i ρ .
We suppose that both materials are separated by a sharp interface Γ(t) that is defined at the instant t > 0 by the points x ∈ R d where a function (x, t) → ϕ(x, t) vanishes, i.e. Γ(t) = {x ∈ D | ϕ(x, t) = 0}. The function ϕ is called the Level Set function and we shall suppose that in the vicinity of Γ(t) the Level Set function ϕ coincides with the signed distance to Γ(t) [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF][START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. The domain where the fluid i = 1, 2 lies will be noted D i (t) = {x ∈ D | (-1) i ϕ(x, t) < 0}. In this stage, we neglect the heat transfer phenomenon. We shall suppose that each fluid i = 1, 2 verifies the compressible Navier-Stokes in D i (t) that is to say:

t > 0, x ∈ D i (t)            ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ∇ • S, ∂ t (ρE) + ∇ • ((ρE + p)u) = ∇ • (Su), (2a) (2b) (2c)
where u is the fluid velocity and p = p EOS i (ρ, e) is the fluid pressure. For the sake of simplicity, we shall assume in the following that both fluid 1 and 2 share the same viscosity µ 1 = µ 2 = µ so that the viscous stress tensor S is defined by

S = µ ( ∇u + (∇u) T ) -2 3 µ(∇ • u)I in D i , i = 1, 2.
The total energy per unit of mass is noted E = e + 1 2 |u| 2 . If n is the unit normal vector to the interface oriented from fluid 1 to fluid 2, we suppose that the motion equations for the fluid 1 and the fluid 2 are coupled across the interface thanks to the following three hypotheses:

(H1) the fluid velocity normal to the interface (u • n) is continuous across the interface, (H2) no mass transfer occurs between the two fluids, (H3) surface tension effects are accounted for by supposing that there is a pressure jump across the interface that verifies the Laplace relation [START_REF] Landau | courseof theoretical physics[END_REF]:

lim x→y x∈D 1 p -lim x→y x∈D 2 p = σκ(y, t), for y ∈ Γ(t), (3)
where σ is the surface tension coefficient and κ is the curvature at point y of the interface. The hypotheses (H1) and (H2) imply that ϕ verifies an advection equation with velocity u :

∂ t ϕ + u • ∇ϕ = 0, ∀x ∈ D, t > 0. (4)
Finally let us recall that the unit normal vector n and the curvature κ at the interface can be computed as follows for x ∈ D

n = ∇ϕ |∇ϕ| , κ = ∇ • ( ∇ϕ |∇ϕ| ) . (5)

Numerical solver for compressible two-phase flows

In this work, both fluids 1 and 2 are considered compressible. To adapt the global system to "all-regime Lagrange-Projection-like strategy "proposed in [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF], acoustic, transport and viscous terms are split during a small time interval.

Acoustic/transport/viscous diffusion operator splitting

We will now exhibit three subsystems the original model is split into. The first one is the acoustic system that is set by accounting for pressure-related terms and surface tension effects. We consider the following acoustic subsystem defined on the whole domain D:

t > 0, x ∈ D                                ρ∂ t (1/ρ) -∇ • (u) = 0, ρ∂ t u + ∇p = ( σκ |∇ϕ| ) δ(ϕ)∇ϕ, ρ∂ t E + ∇ • (pu) = ( σκ |∇ϕ| ) δ(ϕ)u • ∇ϕ, ∂ t ϕ = 0, (6a) (6b) (6c) (6d)
where δ(ϕ) is the Dirac measure whose support is Γ(t). System (6) allows to encompass altogether both bulk fluid motion and the jump relation (3) thanks to the introduction of the Dirac measure that handles the surface tension effects.

Then we consider a subsystem that involves the transport terms for both fluids and the Level Set function:

t > 0, x ∈ D i                    ∂ t ρ + u • ∇ρ = 0, ∂ t (ρu) + (u • ∇)ρu = 0, ∂ t (ρE) + u • ∇(ρE) = 0, ∂ t ϕ + u • ∇ϕ = 0, (7a) (7b) (7c) (7d)
and a last subsystem for viscous effects:

t > 0, x ∈ D i                    ∂ t ρ = 0, ∂ t (ρu) = ∇ • S, ∂ t (ρE) = ∇ • (Su), ∂ t ϕ = 0. (8a) (8b) (8c) (8d)
The overall numerical solver proposed in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] is a three-step strategy: given a fluid state (ρ, ρu, ρE, ϕ) n j at time t n = n∆t, the updated value (ρ, ρu, ρE, ϕ) n+1 j at instant t n+1 is built as follows:

• Update the fluid state (ρ, ρu, ρE, ϕ) n j to the value (ρ, ρu, ρE, ϕ) n+ j by approximating the solution of the acoustic subsystem [START_REF] Liu | The ghost fluid method for compressible gaswater simulation[END_REF].

• Update the fluid state (ρ, ρu, ρE, ϕ) n+ j to the value (ρ, ρu, ρE, ϕ) n+1- j by approximating the solution of the transport subsystem [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF].

• Update the fluid state (ρ, ρu, ρE, ϕ) n+1- j to the value (ρ, ρu, ρE, ϕ) n+1 j by approximating the solution of the viscous subsystem (8).

Approximate Riemann solution for the acoustic subsystem

The structure of the acoustic system (6) shares similarities with the acoustic system studied in [START_REF] Chalons | A large time-step and wellbalanced Lagrange-projection type scheme for the shallow water equations[END_REF]. This allows [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF] to propose an approximate Riemann solver W RP that is consistent in the integral sense with the full system (6) including the surface tension terms.

Let W T = (ρ, u, v, w, E, ϕ, π) where π acts as an approximated pressure value for the system, we briefly recall hereafter the definition of W RP . We have

W RP (ζ; W l , W r ) =                    W l , if ζ < -a l , W * l , if -a l < ζ < 0, W * r , if 0 < ζ < a r , W r , if a r < ζ, (9)
for intermediate states W * l and W * r defined by

π * l = (a l p r + a r p l ) a l + a r - a l a r (u r -u l ) a l + a r + σκa l a r + a l [H (ϕ r ) -H (ϕ l )], (10a)
π * r = (a l p r + a r p l ) a l + a r - a l a r (u r -u l ) a l -a r - σκa r a l + a r [H (ϕ r ) -H (ϕ l )], (10b)
1 ρ * l = 1 ρ l + a r (u r -u l ) + p l -π r -σκ[H(ϕ r ) -H(ϕ l )] a l (a l + a r ) , ( 10c
)
1 ρ * r = 1 ρ r + a l (u r -u l ) + p r -p l + σκ[H(ϕ r ) -H(ϕ l )] a r (a l + a r ) , ( 10d
)
u * = u * l = u * r = a l u l + a r u r a l + a r - p r -p l + σκ[H(ϕ r ) -H(ϕ l )] a l + a r , (10e)
v * l = v l , v * r = v r , (10f)
w * l = w l , w * r = w r , (10g)
E * l = E l - π * l u * -p l u l a l , (10h)
E * r = E r + π * r u * -p r u r a r , ( 10i
)
where H is the Heaviside function and κ is the local value of interface curvature. The parameters a l and a r are user-chosen constants that need to be sufficiently large approximations for ρ l c l and ρ r c r so as to satisfy stability constraints for the numerical scheme. The Suliciu-type relaxation approximation adapted to the vacuum that allows to derive W RP is recalled in Appendix A. Indeed, according to Bouchut [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws: and Well-Balanced schemes for Sources[END_REF], if one notes {b} + = max(b, 0), it is possible to show that defining a l and a r as follows

if π r -π l + σκ [ H(ϕ r ) -H(ϕ l ) ] ≥ 0, set                a l = ρ l c l + γ l +1 2 ρ l { π r -π l + σκ[H(ϕ r ) -H(ϕ l )] ρ r c r + u l -u r } + a r = ρ r c l + γ r +1 2 ρ r { π l -π r -σκ[H(ϕ r ) -H(ϕ l )] a l + u l -u r } + if p r -p l + σκ [ H(ϕ r ) -H(ϕ l ) ] ≤ 0 set                a r = ρ r c l + γ r +1 2 ρ r { p l -p r -σκ[H(ϕ r ) -H(ϕ l )] c l ρ l + u l -u r } + a l = ρ l c l + γ l +1 2 ρ l { p r -p l + σκ[H(ϕ r ) -H(ϕ l )] a r + u l -u r } +
ensures positivity for the density. In the present work we will consider two alternate choices: a first simple choice that allows to recover the classic Suliciu relaxation solver [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. i. rate-type fluids[END_REF][START_REF] Chalons | Navier-stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws: and Well-Balanced schemes for Sources[END_REF] by setting

a l = a r = K max(ρ l c l , ρ r c r ), (11)
where K > 1 is a constant and a second heuristic choice obtained by defining

a l = K ′ ρ l c l , a r = K ′ ρ r c r , K ′ = 1.1. (12)

Finite volume discretization

We now recall the finite volume approximation for each of the acoustic, transport and viscous effect systems. We suppose that the computational domain is discretized over a Cartesian grid. If Ω j and Ω k are two cells of the mesh, we note ∂Ω jk = Ω k ∩ Ω j . We define the neighbourhood of a cell Ω j thanks to the set of indices N( j) = {k | ∂Ω jk ∅} and we note n jk the unit normal vector to ∂Ω jk oriented from Ω j to Ω k . The center of the cell Ω j will be denoted by x j and for k ∈ N( j) the center of ∂Ω jk will be noted x jk = x k j as depicted in figure 1. With a slight abuse of notation, for a function x → b, when k ∈ N( j) we will note

b(x jk ) = lim x→x jk x∈Ω j b(x), b(x k j ) = lim x→x jk x∈Ω k b(x).
We define N ± ( j) = {k ∈ N( j) | ± ϕ n k ϕ n j > 0} and we also note

∆x jk = x j -x k = ∆x k j , so that ∆x jk /2 = x j -x jk = x k j -x k and ∆x k j /2 = x k -x k j = x k j -x j .
x k j x j cell j

x k cell k Figure 1: Mesh for a multi-dimensional problem with an intersection point x jk at the cell face.

Approximation of the acoustic subsystem

By adapting the lines of Chalons et al. [START_REF] Chalons | A large time-step and wellbalanced Lagrange-projection type scheme for the shallow water equations[END_REF] thanks to the rotational invariance of acoustic system, we use the approximate Riemann solver described in Section 3.2 to build a Finite Volume approximation of the acoustic system. It reads

L j ρ n+ j = ρ n j , L j (ρu) n+ j = (ρu) n j - ∆t |Ω j | ∑ k∈N( j) |∂Ω jk |π * jk n jk , L j (ρE) n+ j = (ρE) n j - ∆t |Ω j | ∑ k∈N( j) |∂Ω jk |π * jk u * jk , ϕ n+ j = ϕ n j , L j = 1 + ∆t |Ω j |         ∑ k∈N( j) |∂Ω jk |u * jk         , (13a) (13b) (13c) (13d) (13e) (13f)
where

                   u * jk = n jk • (a j u n j + a k u n k ) a j + a k - p n k -p n j + σκ n jk [H(ϕ n k ) -H(ϕ n j )] (a j + a k ) , π * jk = a k p n j + a j [p n k + σκ n jk (H(ϕ n k ) -H(ϕ n j )] a j + a k - a j a k a j + a k n jk • (u n k -u n j ). (14a) (14b)
The flux in [START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF] gives an explicit update for the acoustic system. The term κ jk is an approximation of the interface curvature in the vicinity of ∂Ω jk . Using [START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF] we note that κ can be evaluated by:

κ = ∇ • (∇ϕ) |∇ϕ| = 2ϕ x ϕ y ϕ xy -ϕ 2 x ϕ yy -ϕ 2 y ϕ xx (ϕ 2 x + ϕ 2 y ) 3/2 . ( 15
)
Relation ( 15) is used to compute a value κ j for any cell Ω j in the vicinity of the interface by approximating first and second derivatives of ϕ thanks to centered second order discretizations. Then κ jk is defined as the following average value

κ jk = κ j |ϕ k | + κ k |ϕ j | |ϕ j | + |ϕ k | . ( 16
)
The estimate κ jk of the curvature in the vicinity of ∂Ω jk allows to account for surface tension effects thanks to the flux formula (14b). A Courant-Friedrichs-Lewy (CFL) constraint on the time step choice for (13)-( 14) pertaining to the acoustic effects reads ∆t acoustic max j, k∈N( j)

( |∂Ω jk | |Ω j | max ( a j ρ j , a k ρ k )) ≤ 1 2 . ( 17
)
Following the study of [START_REF] Elmore | Physics of waves[END_REF], a capillary wave with a wave number K travels at velocity c cw given by

c cw = √ σK ρ 1 + ρ 2 . ( 18
)
If we note Π the Archimedes' constant, the largest wave number that can be described over our mesh is K = Π max j, k∈N( j) |∂Ω jk |/|Ω j |. Following [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] we impose the stability for the approximation of the capillary waves by requiring that:

c cw ∆t capillary |∂Ω jk | |Ω j | ≤ 1 2 . ( 19
)
By injecting ( 18) into [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF], we impose a restriction on the time step accounting for surface tension effects:

∆t capillary max j, k∈N( j) √ 4Πσ (ρ j + ρ k ) ( |∂Ω jk | |Ω j | ) 3 ≤ 1. (20)
3.3.2. Approximation of the transport subsystem for the fluid variables As a precise interface description is highly requested for two-phase flow simulations, the Level Set advection will be treated independently with a high-order scheme and will be introduced later. We write the transport subsystem with conservative variables which causes conservative terms to appear:

∂ t b + ∇ • (bu) -b∇ • u = 0, b ∈ {ρ, ρu, ρv, ρw, ρE}.
This system is approximated with an upwind scheme using the values of the normal velocities u * jk computed at each cell edge ∂Ω jk during the acoustic step. The discretized equations then read

b n+1- j = b n+ j - ∆t |Ω j | ∑ k∈N( j) |∂Ω jk |u * jk b n+ jk + b n+ j ∆t |Ω j | ∑ k∈N( j) |∂Ω jk |u * jk , (21)
with b n+ jk =              b n+ j if u * jk > 0, for k ∈ N( j), b n+ k if u * jk ≤ 0, for k ∈ N + ( j), b n+ k,ghost if u * jk ≤ 0, for k ∈ N -( j), ( 22
)
where b n+ k,ghost is a ghost value that associated with the so-called ghost-cells Ω k when k ∈ N -( j) as represented in figure 2. Let us briefly recall how the ghost values are computed (see [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] for a detailed presentation). For the sake of simplicity we shall suppose here that the computational domain is discretized over a uniform Cartesian grid of space step ∆x. We consider b ∈ {ρ, ρE} and we note (b n ) = (n • ∇b) the differential of b along n. The values of b are propagated from the regions where ϕ j < 0 to the ϕ j > 0 regions using the following steps: • compute an approximate value (b n ) j = (n • ∇b) j for (b n ) in the cell Ω j for j such that ϕ j < -∆x thanks to second-order central differences formulas,

• extrapolate the value of b n in the region ϕ > -∆x by solving the evolution equation

∂ t b n + H(ϕ + ∆x)n • ∇b n = 0, (23)
thanks to a WENO5 scheme as presented in [START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF],

• propagate the value of b in the ghost region ϕ > 0 by solving the PDE

∂ t b + H(ϕ)(n • ∇b -b n ) = 0, (24)
with a method of lines [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF].

Stability is ensured under the following classic CFL that only involves u * jk :

∆t transport max j, k∈N( j) ( |∂Ω jk | |Ω j | max(|u j |, |u * jk |) ) ≤ 1 2 . ( 25
)
Let us note that in the low Mach regime, ( 25) is less restrictive than the acoustic CFL (17).

Approximation of the transport for the Level Set function

As detailed in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF][START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF], the discrete values of the Level Set function ϕ are updated by means of a high-order one-step numerical method. This method provides 5 th -order accuracy with respect to both time and space approximation of (7d). Stability for this numerical scheme is also granted under the CFL condition [START_REF] Pelanti | Low Mach number preconditioning techniques for Roe-type and HLLC-type methods for a two-phase compressible flow model[END_REF]. This advection step requires to be completed by a classic redistancing stage. Following [START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF], we achieve this task by performing a numerical resolution of a Hamilton-Jacobi equation with a WENO5-RK3 scheme.

Cell state update

As the interface is advected by the fluid velocity, the value of ϕ evolves with time and may change its sign. For a cell that changes its sign, its real state is replaced by its ghost state in that cell.

Viscous subsystem discretization

The viscous subsystem is approximated by using a classic centered secondorder finite-volume method.

The time step constraint for the viscous step reads

∆t viscosity max j, k∈N( j) ( µ ρ j |∂Ω jk | 2 |Ω j | 2 ) ≤ 1 2 . ( 26
)
In the next sections, we will examine the behavior of the present scheme and propose a low Mach fix to limit the truncation error in this regime.

Definition of the low Mach regime for the two-phase flows with interface

By following the work of [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF][START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF][START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], the behavior of this system of equation (2) will be investigated with respect to the Mach regime at continuous level. We will perform an asymptotic analysis in the regime M i ≪ 1 that will provide a set of conditions that must be fulfilled by the flow parameters during the motion and also at the initial conditions [START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF][START_REF] Cordier | An asymptotic-preserving all-speed scheme for the euler and navier-stokes equations[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equations[END_REF].

Definition of non-dimensional variables

In order to study the behavior of numerical schemes in the low Mach regime, we first need to characterize the flow with respect to the Mach number. To that end, we introduce characteristic values for length, time and velocity that will be respectively noted x, û, t = x/û. This allows to define the following nondimensional quantities and nabla operator that are common to both fluids

t = t/ t, x = x/ x, κ = κ x, ũ = u/û, ∇ = ( ∂ ∂ x1 , . . . , ∂ ∂ xd ) T = x ∇.
The non-dimensional domain occupied by the fluid i and the interface will be respectively noted Di and Γ. We now define a set of non-dimensional variables for each fluid i = 1, 2 by setting

⌊ ρ⌋ i = ρ/ ρi , ⌊ p⌋ i = p/ pi , ⌊c⌋ i = c/ĉ i , ⌊ẽ⌋ i = e/ê i , ⌊ σ⌋ i = σ/( pi x), ⌊ S⌋ i = 1 Re i S/( ρi û2 ) (27)
where ρi , ĉi , pi = ρi ĉ2 i , êi = pi / ρi , σi and Re i denote respectively a characteristic density, sound velocity, pressure, internal energy, surface tension coefficient and Reynolds number associated with the fluid i = 1, 2. In the following, we will enclose an expression between the signs ⌊•⌋ i in order to indicate that it has to be considered as a non-dimensional expression with respect to the characteristic quantities of the fluid i. For example, ⌊ ρ⌋ 1 is the dimensionless density variable with respect to the characteristic density ρ1 of fluid 1. We note M i the Mach number associated with the fluid i that is defined by

M i = û/ĉ i .
The dimensionless total energy is then defined as:

⌊ Ẽ⌋ i = E/ê i = ⌊ẽ⌋ i + M 2 i . (28)
Furthermore, we introduce the following non dimensional parameters α i ∈ R + and β i ∈ R + that characterize the magnitude of characteristic densities and sound velocities with respect to the fluid 1

α i = ρi / ρ1 , β i = ĉi /ĉ 1 .
And finally, we also note

α m = min(α 2 , 1/α 2 ) ≤ 1, β m = min(β 2 , 1/β 2 ) ≤ 1. (29)

Dimensionless evolution equation and low Mach regime

For the sake of readability, when there is no ambiguity, we shall replace ⌊•⌋ i by • to refer to a non-dimensional expression with respect to the fluid i. With this simplified notation, the governing equation ( 2) for the bulk flow in D i can be written in a dimensionless form:

t > 0, x ∈ Di (t)                          ∂t ρ + ∇ • ( ρ ũ) = 0, ∂t( ρ ũ) + ∇ • ( ρ ũ ⊗ ũ) + 1 M 2 i ∇ p = 1 Re i ∇ • S, ∂t( ρ Ẽ) + ∇ • (( ρ Ẽ + p) ũ) = M 2 i Re i ∇ • ( S ũ). (30a) (30b) (30c)
In this work, we aim at deriving a numerical method that preserves its accuracy in the low Mach regime. In order to clearly distinguish between the improvement of our method and possible damping of spurious effects due to the viscosity, we only consider flows that involve small viscosity effects. More precisely we shall suppose that O(Re) = 1/M i . In this case, one can see that the discretization of the dissipative terms is uniform with respect to the Mach number and does not affect the performance of the scheme. Therefore for the sake of brevity, we shall only consider the analysis of the convective and capillary terms.

The non-dimensional form of the interface evolution ( 4) is

∂tϕ + ũ • ∇ϕ = 0, t > 0, x ∈ D. (31)
Concerning the boundary conditions across the interface Γ: (H1) implies that ( ũ • ñ) is continuous across Γ and the relation (3) yields a non-dimensional jump relation that reads

             lim x→ỹ x∈ D1 p             1 -              lim x→ỹ x∈ D2 p             1 = ⌊ σ⌋ 1 κ(ỹ, t), for ỹ ∈ Γ( t). (32)
We want to note that as the capillary effects are treated as the jump condition, the non-dimensional jump relation is always valid whether the reference is fluid 1 or 2. Indeed, we have also:

             lim x→ỹ x∈ D1 p             2 -              lim x→ỹ x∈ D2 p             2 = ⌊ σ⌋ 2 κ(ỹ, t), for ỹ ∈ Γ( t). (33)
Before going any further, let us underline that the overall two-phase system formed by [START_REF] Liu | An asymptotic preserving scheme for the two-dimensional shallow water equations with coriolis forces[END_REF], [START_REF] Yan | A successive penalty-based asymptotic-preserving scheme for kinetic equations[END_REF] and the boundary conditions on Γ depends on both M 1 and M 2 .

In this context, we will now characterize what will be referred to in this paper as the low Mach regime for our two-phase interface model. We suppose that M 1 ≪ 1 and M 2 ≪ 1, and for a dimensionless variable Ãi of fluid i, we consider an asymptotic expansion with respect to M i as follows:

Ãi (x, t) = M 0 i Ã(0) i (x, t)+ M 1 i Ã(1) i (x, t)+ M 2 i Ã(2) i (x, t)+• • • , x ∈ Di ( t), t > 0. ( 34
)
We now examine the behavior of Ã(k) i in the limit regime M i → 0: introducing expression (34) in [START_REF] Liu | An asymptotic preserving scheme for the two-dimensional shallow water equations with coriolis forces[END_REF] and collecting the terms with equal power of M i , we obtain:

• Order M -2 i : ∇ p(0) = 0, x ∈ Di ( t), t > 0, (35)
• Order M -1 i : ∇ p(1) = 0, x ∈ Di ( t), t > 0, (36)
• Order M 0 i :

t > 0, x ∈ Di (t)                    ∂t ρ(0) + ∇ • ( ρ(0) ũ(0) ) = 0, ∂t( ρ(0) ũ(0) ) + ∇ • ( ρ(0) ũ(0) ⊗ ũ(0) ) + ∇ p(2) = 1 Re i ∇ • S(0) , ∂t( ρ(0) Ẽ(0) ) + ∇ • (( ρ(0) Ẽ(0) + p(0) ) ũ(0) ) = 0, (37a) (37b) (37c)
where p(2) is the dynamic pressure related to the fluid motion, and p(0) is the thermodynamic pressure. Thanks to [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic euler and navier-stokes equations[END_REF], we can see that Ẽ(0) i = ẽ(0) i , so that p(0) can be expressed as a function of ρ(0) and Ẽ(0) thanks to the EOS (1):

p(0) = p EOS i ( ρ(0) , Ẽ(0) ), in Di ( t), t > 0. (38)
Then, using [START_REF] Elmore | Physics of waves[END_REF] and (37c) we obtain

∂t p(0) + ũ • ∇ p(0) + ρc 2 ( ρ(0) , ẽ(0) ) ∇ • ũ(0) = 0, in Di ( t), for t > 0. (39)
We can now further characterize the limit regime M i → 0: equations ( 35) and ( 36) that pertain respectively to terms of order M -2 i and M -1 i imply that the pressure is homogeneous with respect to the space variable up to a fluctuation of magnitude of M 2 i :

∇ p = O(M 2 i ), x ∈ Di ( t), t > 0. ( 40
)
We now add two supplementary hypotheses. First, we assume that the boundary conditions are chosen such that the space-homogeneous p(0) is also independent of t. Consequently, [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] implies that

∇ • ũ(0) = 0 in Di ( t), for t > 0. (41)
Let us emphasize that if such hypothesis is classic in the case of single-material flow, it is much stronger in our case due to the presence of the two-phase interface. Second, we suppose that the above property is also true across the interface, that is to say:

⌊ ∇ • ũ(0) ⌋ 1,2 = 0, on Γ. (42)
Before going any further, let us sum up what shall be referred to as a flow in the low Mach regime for our two-phase interface model : such flow shall verify

                               ũ is regular (C 1 function) in D, p and ρ are regular (C 1 function) in D i (t), ∇ • ũ = O(M i ), t > 0, x ∈ D i (t), ∇ p = O(M 2 i ), t > 0, x ∈ D i (t), ⌊ ∇ • ũ⌋ 1,2 = O(M i ), t > 0, x ∈ Γ(t). (43a) (43b) (43c) (43d) (43e)
An accurate scheme for low Mach two-phase flows should be able to preserve the low Mach regime as previously defined. In the following sections, we will examine the behavior for initial data that verifies a discrete equivalent of (43).

Low Mach behavior of the numerical scheme

In this section we will present a set of conditions that the discretization strategy of Section 3 needs to verify in order to accurately approximate the two-phase flow presented in Section 2 when the flow matches the low Mach definition [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF]. Let us mention that we employ an uniform Cartesian mesh, which implies ∆x jk = ∆x.

We first recast the first two subsystem ( 6) and [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF] we considered in Section 3 using non-dimensional form in each domain Di . With non-conservative dimensionless variables ( 1 ρ , ũ, Ẽ, ϕ), the acoustic system (6) for bulk flows reads

t > 0, x ∈ Di                          ρ∂t(1/ ρ) -∇ • ũ = 0, ρ∂t ũ + 1 M 2 i ∇ p = 0, ρ∂t Ẽ + ∇ • ( p ũ) = 0, ∂tϕ = 0. (44a) (44b) (44c) (44d)
As the Mach number M i is explicitly involved in [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF], it is also interesting to consider the evolution of ρ(0) , ũ(0) and Ẽ(0) through the acoustic phase. We get

t > 0, x ∈ Di                      ∂t ρ(0) = 0, ρ(0) ∂t ũ(0) + ∇ p(2) = 0, ∂t Ẽ(0) = 0, ∂tϕ = 0. (45a) (45b) (45c) (45d)
The transport subsystem [START_REF] Fechter | A sharp interface method for compressible liquid-vapor flow with phase transition and surface tension[END_REF] in non-dimensional form becomes

t > 0, x ∈ Di ∂t b + ∇ • ( b ũ) -b( ∇ • ũ) = 0, b ∈ { ρ, ρũ, ρṽ, ρ w, ϕ}. (46)
In order to examine the behavior of the numerical scheme, we will follow the lines of [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF] by considering smooth solutions of (30) that fulfill the low Mach regime requirements [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF]. We note (x, t) → b any rescaled fluid parameter associated with this solution. We commit an abuse of notation by setting b(x j , t) = b j . Then we shall inject such solutions into the rescaled numerical scheme and study the resulting truncation error.

Low Mach behavior of the acoustic step in the bulk

We start by examining the acoustic step by considering a cell

Ω j ⊂ D i . Let us note a i = max {r | Ω r ⊂D i }
(ρ r c r ), for the sake of simplicity we suppose that the artificial acoustic impedance are defined using [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] with the uniform choice a j = a i for all cell Ω j ⊂ D i . Then, the non-dimensional form of ( 13) is

                                     ρn j        1 ρn+ j - 1 ρn j        -∆ t ∑ k∈N( j) |∂ Ω jk | |Ω j | ũ * jk = 0, ρn j ( ũn+ j -ũn j ) + ∆ t M 2 i ∑ k∈N( j) |∂ Ω jk | |Ω j | π * jk n jk = 0, ρn j ( Ẽn+ j -Ẽn j ) + ∆ t ∑ k∈N( j) |∂ Ω jk | |Ω j | π * jk ũ * jk = 0, (47a) (47b) (47c)
where

                 ũ * jk = n jk • ( ũn j + ũn k ) 2 + 1 M i πn j -πn k 2ã i , π * jk = πn j + πn k 2 -M i ãi ( ũn k -ũn j ) 2 n jk . (48a) (48b)
The low Mach hypothesis [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF] implies that there exists three smooth functions A i , B i and C i of magnitude O(M 0 i ) such that we have tn > 0, for any xk and x j in Di

       ũ(x k , tn ) • n = ũ(x j , tn ) • n + A i (x j , xk , tn , n)|x k -x j |, p( xk , tn ) = p(x j , tn ) + M 2 i B i (x j , xk , tn )|x k -x j |. (49a) (49b)
Relations ( 48) and (49) yield that

                       ũ * jk = ũ(x jk , tn ) • n jk - 1 2ã i M i B i (x j , xk , tn )|x k -x j | + O(M i ∆ x) + O(∆ x2 ), π * jk = p(x jk , tn ) - ãi 2 M i A i (x j , xk , tn , n)|x k -x j | + O(M 2 i ∆ x2 ), π * jk ũ * jk = p(x jk , tn ) ũ(x jk , tn ) • n jk -M i C i (x j , xk , tn , n)|x k -x j | + O(∆ x2 ), (50a) (50b) (50c)
so that by injecting (50) into (47), we get

                                         ρn j ρn+ j -1 - ∆ t ∆ x          ∑ k∈N j ũn jk • n jk + O(M i ∆ x) + O(∆ x2 )          = 0, ( ũn+ j -ũn j ) + ∆ t ρ n j ∆ xM 2 i          ∑ k∈N j pn jk + O(M i ∆ x) + O(M 2 i ∆ x) + O(M 2 i ∆ x2 )          = 0, ( Ẽn+ j -Ẽn j ) + ∆ t ρn ∆ x         ∑ k∈N( j) pn jk ũn jk • n jk + O(M i ∆ x) + O(∆ x2 )         = 0. (51a) (51b) (51c)
Classic computations yield

∑ k∈N( j) ũ(x jk , tn ) • n jk = ∇ • ũ(x j , tn )∆ x + O(∆ x2 ), ∑ k∈N( j) p(x jk , tn ) • n jk = ∇ p( x j , tn )∆ x + O(M 2 i ∆ x2 ), ∑ k∈N( j) p(x jk , tn ) ũ(x jk , tn ) • n jk = ∇ • ( p(x jk , tn ) ũ(x jk , tn ) ) ∆ x + O(∆ x2 ). ( 52a
) (52b) (52c)
With (52) at tn , we get:

                           ρn+ j = ρn j / ( 1 + ∆ t ∇ • ũ + O(∆ tM i ∆ x) + O(∆ t∆ x2 )
) ,

ũn+ j = ũn j - ∆ t M 2 i ρ j ∇ p + O( ∆ t∆ x M i ) + O(∆ x2 ∆ t), Ẽn+ j = Ẽn j - ∆ t ρ j ∇ • ( p ũ) + O(M i ∆ x∆ t) + O(∆ x2 ∆ t). ( 53a
) (53b) (53c)
In the low Mach regime, the rescaled discretization of the acoustic system is consistent with:

t > 0, x ∈ Di                          ∂t ρ + ρ ∇ ũ = O(∆ xM i ) + O(∆ t), ∂t ũ + 1 M 2 i ρ ∇ p = O(∆ xM 0 i ) + O ( ∆ x M i ) + O(∆ t), ∂t Ẽ + 1 ρ ∇ • ( p ũ) = O(∆ xM i ) + O(∆ t). (54a) (54b)
(54c)

Low Mach behavior of the transport step in the bulk

The transport step in non-dimensional form for a cell Ω j ⊂ D i reads

Ω j ∩ Di ∅, bn+1- j = bn+ j - ∆ t | Ω j | ∑ k∈N j |∂ Ω jk |ũ * jk bn+ jk + bn+ j ∆ t | Ω j | ∑ k∈N j |∂ Ω jk |ũ * jk , (55)
where bn+ jk are non-dimensional form of ( 22) with respect to the fluid i. By injecting the low Mach solution into (55), standard computations show that:

bn+1- j = bn+ j - ∆ t ∆ x ∑ k∈N j ( ũn jk • n jk + O(M i ∆ x) + O(∆ x2
)

) bn+ jk + bn+ j ∆ t ∆ x          ∑ k∈N j ũn jk • n jk + O(M i ∆ x) + O(∆ x2 )          , (56) hence bn+1
- j = bn+ j -∆ t ∇ • ( ũb ) + ∆ t b ∇ • ũ + O(∆ x∆ t) + O( xM i ∆ t). (57)
In the low Mach regime, the rescaled transport system is consistent with:

t > 0, x ∈ Di , ∂t b+ ∇•( ũb )-b( ∇• ũ) = O(M i 0 ∆ x)+O(∆ t)+O(M i ∆ x), b ∈ { ρ, ρũ
, ρṽ, ρ w, ϕ}.

(58)

Low Mach regime accuracy analysis and fix for the bulk flow

Results presented in Sections 5.1 and 5.2 are coherent with the literature: in the acoustic step, a truncation error term of magnitude O(∆ x/M i ) that appears in (54b) suggests that important loss of accuracy may occur when M i ≪ 1, while in the transport step, the truncation error remains uniform with respect to M i .

Many works over the past years have proposed low Mach corrections for single-fluid flows [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Li | An all-speed Roe-type scheme and its asymptotic analysis of low Mach number behaviour[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Bruel | A low Mach correction able to deal with low Mach acoustics[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes[END_REF]. The fix usually amounts to use a centered pressure evaluation at the cell interfaces in the low Mach regime. In the case of a two-phase flow, let us recall two modifications that have been proposed in the literature.

Noting (M i ) jk a local evaluation of the Mach number M i at the interface ∂Ω jk and setting θ jk = min((M 1 ) jk , (M 2 ) jk , 1), Peluchon et al. [START_REF] Peluchon | A robust implicit-explicit acoustictransport splitting scheme for two-phase flows[END_REF] proposed to replace π * by altering terms related to velocity jump as follows:

π * ,θ,AW jk = a k π j + a j (π k + [ H(ϕ k ) -H(ϕ j )
] σκ)

a j + a k -θ jk a j a k a j + a k n jk (u k -u j ). (59)
Chalons et al. [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] studied a slightly different low Mach correction with a centered pressure that can be expressed as:

π * ,θ,CP jk = (1 -θ jk ) π j + π k + [ H(ϕ k ) -H(ϕ j ) ] σκ 2 + θ jk π * jk . (60)
Let us underline that the two-phase models studied in both Peluchon et al. [START_REF] Peluchon | A robust implicit-explicit acoustictransport splitting scheme for two-phase flows[END_REF] and Chalons et al. [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase flows on unstructured meshes[END_REF] are different from our framework as they involve potential mixture regions where both fluids can simultaneously be present. Applying either (59) or (60) in our case succeeds in controlling the error term in the momentum update of the acoustic step in pure fluid i region. Indeed, either correction will enable a new truncation error estimate

t > 0, x ∈ Di ρ∂t ũ + 1 M 2 i ∇ p = O(∆ xM 0 i ) + O ( θ∆ x M i ) + O(∆ t), (61)
where θ is an upper bound for all θ jk . Proof. For the bulk flow of fluid i, both low Mach corrections can take the simple form:

π * ,θ jk = πn j + πn k 2 -M i θ jk ãi ( ũn k -ũn j ) 2 n jk . (62)
There exists a function D i of magnitude 1 with respect to M i such that:

π * ,θ jk = p(x jk , tn ) -

D i ãi 2 M i A i (x j , xk , tn , n)|x k -x j | + O(M 2 i ∆ x2 ). (63)
Following the similar analysis given in (50)-(54), then we can get:

ũn+ j = ũn j - ∆ t M 2 i ρ j ∇ p + O(θ jk ∆ t∆ x M i ) + O(∆ x2 ∆ t). (64)

Asymptotic behavior of existing low Mach schemes across the interface

We have seen in Section 5 that in the low Mach regime, our discretization strategy may suffer a loss of accuracy in each pure fluid region D i of the computational domain. This drawback is a classic pathology in the case of single fluid flows that can be improved by modifying the pressure discretization. We shall now investigate in the present section the behavior of the numerical scheme across the two-phase interface in the low Mach regime. We will see that the approximation of several flow parameters will be impacted in the low Mach regime due to the abrupt jump of the medium properties across the interface. Once again, we consider a smooth solution of (30) that fulfill the low Mach regime requirements [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF]. We consider a cell Ω j such that ϕ n j > 0 that lies in the vicinity of Γ so that N -( j)

∅. Let V be a neighborhood of x j so that V ∩ D 1 ∅ and V ∩ D 2 ∅. Thanks to the smoothness of the velocity we know that ⌊ũ⌋ i and its derivatives are of magnitude O(M 0 i ) within the entire region V for i = 1, 2. Thus there exist a smooth function A such that the variation of the normal velocity at the boundary Ω j reads

                               ⌊ ũ(x j , tn ) • n jk ⌋ 1 = ⌊ ũ(x jk , tn ) • n jk + A(x jk , x j , tn , n k j ) |∆x jk | 2 + O(∆ x2 ) ⌋ 1 for k ∈ N( j), ⌊ ũ(x k , tn ) • n jk ⌋ 1 = ⌊ ũ(x k j , tn ) • n jk + A(x k j , xk , tn , n jk ) |∆x k j | 2 + O(∆ x2 ) ⌋ 1 for k ∈ N + ( j), ⌊ ũ(x k , tn ) • n jk ⌋ 2 = ⌊ ũ(x k j , tn ) • n jk + A(x k j , xk , tn , n jk ) |∆x k j | 2 + O(∆ x2 ) ⌋ 2 for k ∈ N -( j). ( 65a
) (65b) (65c)
Similarly according to (43d), ⌊ ∇ p⌋ i is of magnitude O(M 2 i ) within V ∩ D i so that there exist two smooth functions B 1 and B 2 respectively of magnitude

O(M 0 1 ) and O(M 0 2 ) such that                                ⌊ p(x j , tn ) ⌋ 1 = ⌊ p(x jk , tn ) + M 2 1 B 1 (x jk , x j , tn ) |∆x jk | 2 + O(M 2 1 ∆ x2 ) ⌋ 1 for k ∈ N( j), ⌊ p(x k , tn ) ⌋ 1 = ⌊ p(x k j , tn ) + M 2 1 B 1 (x jk , xk , tn ) |∆x k j | 2 + O(M 2 1 ∆ x2 )) ⌋ 1 for k ∈ N + ( j), ⌊ p(x k , tn ) ⌋ 2 = ⌊ p(x k j , tn ) + M 2 2 B 2 (x jk , xk , tn ) |∆x k j | 2 + O(M 2 2 ∆ x2 )) ⌋ 2 for k ∈ N -( j). ( 66a
) (66b) (66c)
The jump condition across the interface can be expressed as:

⌊ p(x jk , t) -p(x k j , t) ⌋ 1 = ⌊[ H(ϕ j ) -H(ϕ k ) ] σκ(x jk , t) ⌋ 1 = ⌊ σκ(x jk , t) ⌋ 1 when k ∈ N -( j). ( 67
) It is important to note that ⌊ p(x k , t) ⌋ 1 (resp. ⌊ã k ⌋ 1 ) when k ∈ N -( j) is not of magnitude O(1) as p(x k ) (resp. a k ) is evaluated in D 2 but
it is rescaled using the characteristic values of the fluid 1. Nevertheless when k ∈ N -( j), one can express ⌊p(x k )⌋ 1 and ⌊ã k ⌋ 1 using terms of magnitude O(1) as follows:

⌊ p(x k , t) ⌋ 2 = ⌊ p(x k , t) ⌋ 1 /(α 2 β 2 2 ), (68)
and

⌊ ãk ⌋ 1 = ⌊ ãk ⌋ 2 (α 2 β 2 ). (69)
Thus (66c) can now be expressed using characteristic values of the fluid 1, that is to say

⌊ p(x k , t) ⌋ 2 = ⌊ p(x k , t) ⌋ 1 /(α 2 β 2 2 ) = ⌊ p(x k j , t) + M 2 1 B 2 (x jk , xk , t) |∆x k j | 2 α 2 + O(M 2 1 ∆ x2 )α 2 ⌋ 1 . (70)
Following the same standard lines as in Section 5 for the bulk flow, we can now reinject the low Mach solution into the expression of the numerical scheme in order to express the truncation error. We examine the cell Ω j and consider that it is mostly occupied by the fluid 1 as we supposed that ϕ j > 0. We will study the behavior of two elements of the numerical scheme: the discretized pressure at ∂Ω jk and the update of ρ(0) j when k ∈ N -( j) ∅.

⌊ π * ,θ,CP jk ⌋ 1 =            πn j + πn k + [ H(ϕ n k ) -H(ϕ n j ) ] σκ n jk 2            1 +θ jk            π * jk - πn j + πn k + [ H(ϕ n k ) -H(ϕ n j ) ] σκ n jk 2            1 .

(71)

By ( 65), (66), ( 67), ( 69), ( 70) and ( 60), tedious calculations show that

⌊ π * ,θ,CP jk, j ⌋ 1 = ⌊ p(x jk, j , tn ) ⌋ 1 + O(M 2 1 ∆ x2 ) + θ jk M 1 ⌊ ã j ⌋ 1 ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 A(x j , xk , tn , n jk )|x k -x j | -θ jk          ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 B 1 (x j , x jk , tn )|x jk -x j | - ⌊ ã j ⌋ 1 α 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 B 2 (x jk , xk , tn )|x k -x jk |          M 2 1 - (1 -θ jk ) 2 ( B 1 (x j , x jk , tn )|x jk -x j | -B 2 (x k , x jk , tn )α 2 |x k -x jk | ) M 2 1 . ( 72
)
Then we can see that in the third line of (72), when α 2 becomes large, it may produce important errors for the estimation of the pressure gradient in Ω j . Indeed, let us consider

⌊ π * ,θ,CP jk -π j ⌋ 1 ∆ x = ⌊ p( x jk ) -p( x j ) ⌋ 1 ∆ x + O(α 2 M 2 1 ) + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 )
.

Using (66) we get ⌊ π * ,θ,CP jk -π j ⌋ 1 ∆ x = O(α 2 M 2 1 ) + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 ). ( 73
)
Therefore, when α 2 = O(1/M 1 ), this suggests that the discrete gradient estimation in the vicinity of the interface may grow to reach a O(1/M 1 ) magnitude which violates the low Mach hypothesis (43d). The analysis given in Appendix B shows that the AW type correction could involve a truncation error of magnitude O(α 2 ) to the discrete pressure gradient in the momentum update, which also contravenes the low Mach hypothesis (43d).

Low Mach flow across the interface: evolution of the density near the interface during the acoustic step

We now evaluate the discretized normal velocity

⌊ ũ * jk ⌋ 1
for k ∈ N -( j): by combining (65), (66), ( 67), ( 69), ( 70) and (14a) we obtain

⌊ ũ * jk ⌋ 1 = ⌊ ũ(x jk , tn ) ⌋ 1 • n jk +          ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 - 1 2          A(x j , xk , tn , n jk )|∆x k j | + O(∆ xM 1 ) + O(∆ x2 ) - 1 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 β 2 α 2 ( B 1 (x j , x jk , tn ) |∆x k j | 2 + B 2 (x k , x jk , tn )α 2 |∆x k j | 2 ) M 1 . (74)
We can now inject (74) and (50a) into (47a), we get

ρn+ j = ρn j /          1 + ∆ t ∇ • ũ + ∑ k∈N -( j)          ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 - 1 2          + O(∆ tM i ∆ x) + O(∆ t∆ x)          . ( 75
)
Compared with (53a), one can then see that due to the discrepancy of characteristic quantities that occurs when α 2 ≫ 1 and β 2 ≫ 1, relation (75) suggests that the leading order of density ρ(0) can no longer remain constant up to a term of magnitude O(1). This error is associated with the estimation of the velocity divergence:

1 ∆ x ∑ k∈N( j) ⌊ ũ * ,θ,N jk ⌋ 1 = ∇ • ũ + O(1) + O(M 1 ). ( 76
)
Without detailed analysis, in the acoustic step update the CP type or AW type correction gives the same truncation error on density and energy update because the same intermediate velocity u * is employed. The CP type intermediate pressure π * ,θ,CP or AW type intermediate pressure π * ,θ,AW could involve a truncation error of magnitude of O(α 2 ) to the momentum update.

Truncation error for the transport update in the vicinity of the interface

We now turn to the transport step and by using (74) we first remark that the normal velocity estimation verifies

⌊ ũ * jk ⌋ 1 = ⌊ ũ(x jk , tn ) ⌋ 1 • n jk + ⌊ O(∆ xM 0 1 α 0 m β 0 m ) + O(∆ xM i ) ⌋ 1 . (77)
Now, for b ∈ {ρ, ρu, ρv, ρw, ρE} and k ∈ N -( j), following [START_REF] Pelanti | Wave Structure Similarity of the HLLC and Roe Riemann Solvers: Application to Low Mach Number Preconditioning[END_REF] we know that b n+ jk is obtained either by the upwind choice within the same fluid either by a linear Ghost Fluid extrapolation. This means that in either case

bn+ jk = b(x jk , tn ) + O(∆ xM 0 1 α 0 m β 0 m ). (78)
Therefore, in the neighborhood of Γ, we obtain a similar truncation error as in the bulk fluid case (58), if b ∈ { ρ, ρũ, ρṽ, ρ w, ϕ}, we have bn+1-

j = bn+ j - ∆ t ∆ x ∑ k∈N j ( ũn jk • n jk + O(M i ∆ x) + O(∆ x) ) bn+ jk + bn+ j ∆ t ∆ x          ∑ k∈N j ũn jk • n jk + O(M i ∆ x) + O(∆ x)          , ( 79
) hence bn+1- j = bn+ j -∆ t ∇ • ( ũb ) + ∆ t b ∇ • ũ + O(∆ x∆ t) + O(M i ∆ t). (80)
Consequently, (80) suggests that the transport step does not perturb the low Mach regime, the resealed transport system is also consistent with (58).

New low Mach correction

We have seen in Section 6 that within cells neighboring the interface Γ, even the classic low Mach corrected schemes failed to preserve the low Mach regime defined by [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF]. The cause of the error seems to originate from the jump of the characteristic values that are associated with the fluid across Γ. Indeed, if one considers the expression of the dimensional pressure in D i , we get from the asymptotic expansion with respect to M i that

p = pi p = p(0) ρi ĉ2 i + p(1) ρi ĉi û + p(2) ρi û2 + • • • in D i . ( 81
)
In the low Mach regime defined by [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF] we know that the first two terms in (81) are constant in D i so that the variations of p are dominated by the dynamic pressure with a magnitude ρi û2 . Across the interface Γ, the dynamic pressure undergoes a jump that is characterized by

ρ2 û2 ρ1 û2 = α 2 . ( 82
)
This corroborates the fact that the error term of magnitude O 73) is issued from the discretized centered pressure term. For similar reasons, relation (76) shows that the the presence of α 2 and β 2 in (74) does not allow to recover an accurate estimate of ∇ • ( ũ(0) ) = 0 in the acoustic step. This results in variations of ρ(0) of magnitude O(1) that generate spurious oscillations in the thermodynamic pressure. In order to cure these problems, we propose to use a new numerical discretization for the pressure and the normal velocity terms at the cell interfaces. Our new numerical scheme should account for the variation of characteristic values associated with the dynamic pressure and also provide a reliable discretization of ∇ • ( ũ(0) ) in the vicinity of Γ. We introduce the following discretization of the pressure and the normal velocity at the interface ∂Ω jk of a cell

( 1 α m ) in (
Ω j for k ∈ N( j) u * ,θ,N jk = (1 -θ jk )n jk • u n j + u n k 2 + θ jk n jk • a i u n j + a k u n k a i + a j + π n j -π n k a j + a k , (83a)
π * ,θ,N jk = (1 -θ jk ) ρ n k π n j + ρ n j [π n k + σκ n jk (H(ϕ n k ) -H(ϕ n j )] ρ n j + ρ n k + θ jk π * jk . (83b)
We can now apply the lines presented in Section 6. Supposing that the flow verifies the low Mach hypotheses [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF] and considering a cell

Ω j that is crossed by Γ, if k ∈ N -( j) we have that ⌊ ũ * ,θ,N jk ⌋ 1 = ⌊ ũ(x jk , t) ⌋ 1 • n jk + θ jk          ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 - 1 2          A(x j , xk , t, n jk )|∆x k j | + O(∆ xM 1 ) + O(∆ x2 ) - 1 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 β 2 α 2 ( B 1 (x j , x jk , t) |∆x k j | 2 + B 2 (x k , x jk , t)α 2 |∆x k j | 2 ) M 1 , (84)
and

⌊ π * ,θ,N jk, j ⌋ 1 = ⌊ p(x jk, j , t) ⌋ 1 + O(M 2 ∆ x2 ) + θ jk M 1 ⌊ ã j ⌋ 1 ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 A(x j , xk , t, n jk )|∆x k j | -θ jk          ⌊ã k ⌋ 2 α 2 β 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 B 1 (x j , x jk , t) |∆x jk | 2 - ⌊ ã j ⌋ 1 α 2 ⌊ ã j ⌋ 1 + ⌊ã k ⌋ 2 α 2 β 2 B 2 (x jk , xk , t) |∆x k j | 2          M 2 1 - α 2 (1 -θ jk ) 1 + α 2 ( B 1 (x j , x jk , t) |∆x k j | 2 -B 2 (x k , x jk , t) |∆x k j | 2 ) M 2 1 . (85)
Relation (85) gives a discrete pressure gradient:

⌊ π * ,θ,N jk -π j ⌋ 1 ∆ x = ⌊ p( x jk ) -p( x j ) ⌋ 1 ∆ x + O(M 2 1 ) + O(θ jk M 1 ) + O(θ jk α 2 M 2 1 ). (86)
The truncation error is now uniform with respect to the density ratio, sound speed ratio and Mach number. Injecting the estimates (84) into (47) yields a discrete velocity divergence to the density evolution for the acoustic step:

1 ∆ x ∑ k∈N( j) ⌊ ũ * ,θ,N jk ⌋ 1 = ∇ • ũ + O(M 1 ) + O(M 2 1 ). (87)
Compared with (76), the error that disturbs the constant thermodynamic pressure profile is eliminated. As a conclusion, (86) and (87) provides a truncation error that is uniform respect to M 1 , α m and β m to the discrete gradient estimation and density evolution, the low Mach hypothesis ( 43) is generally maintained without truncation error of large magnitude such as O(

1 M 1 ) and O( 1 α m ).
As for the transport step, the dimensionless corrected intermediate velocity ũθ,N jk (85) can be written into the same form as (50a). Therefore, we can get the same conclusion about the transport step: the discrete transport step with the new low Mach correction is capable of satisfying preserving the low Mach regime defined by [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF].

An important feature of the new low Mach corrected solver for the acoustic step is that it can be associated with an approximate Riemann solver

(ζ; W l , W r ) → W N RP (ζ; W l , W r ), defined by W N RP (ζ; W l , W r ) =                    W l , if ζ < -a l , W * ,N,θ l , if -a l < ζ < 0, W * ,N,θ r , if 0 < ζ < a r , W r , if a r < ζ, (88)
with intermediate states W * ,N,θ l and W N,θ, * r defined by

u * ,θ,N = u * ,θ,N l = u * ,θ,N r = (1 -θ) u n l + u n r 2 + θ jk a l u n l + a r u n r a l + a r + π n l -π n r -σκ n lr [H(ϕ n r ) -H(ϕ n l )] a l + a r , (89a)
π * ,θ,N l = (1 -θ) ρ n r π n l + ρ n l [π n r + σκ n lr [H(ϕ n r ) -H(ϕ n l )] ρ n l + ρ n r + θπ * l , (89b)
π * ,θ,N r = (1 -θ) ρ n l π n r + ρ n r [π n l -σκ n lr [H(ϕ n r ) -H(ϕ n l )] ρ n l + ρ n r + θπ * r , (89c)
1 ρ * ,θ,N l = 1 ρ l + a l +a r 2 (u r -u l ) + π l -π r + σκ [ H(ϕ l ) -H(ϕ r ) ] a l (a l + a r ) , (89d)
1 ρ * ,θ,N r = 1 ρ r + a l +a r 2 (u r -u l ) + π r -π l -σκ [ H(ϕ l ) -H(ϕ r ) ] a r (a l + a r ) , (89e)
v * ,θ,N l = v l , v * ,θ,N r = v r , (89f)
w * ,θ,N l = w l , w * ,θ,N r = w r , (89g)
E * ,θ,N l = E l - π * ,θ,N l u * ,θ,N -π l u l a l , (89h)
E * ,θ,N r = E r + π * ,θ,N r u * ,θ,N -π r u r a r . (89i)
In Appendix C, it is shown that the following constraints for the acoustic impedance a l and a r :

if p r -p l + σκ[H(ϕ r ) -H(ϕ l )] ≥ 0,                a l = ρ l c l + Υ l ρ l { 2 p r -p l + σκ[H(ϕ r ) -H(ϕ l )] ρ r c r + ρ l c l + u l -u r } + , a r = ρ r c r + Υ r ρ r { 2 p l -p r -σκ[H(ϕ r ) -H(ϕ l )] a l + ρ r c r + u l -u r } + , (90)
if p r -p l + σκ[H(ϕ r ) -H(ϕ l )] < 0,                a r = ρ r c r + Υ r ρ r { 2 p l -p r -σκ[H(ϕ r ) -H(ϕ l )] c l ρ l + c r ρ r + u l -u r } + , a l = ρ l c l + Υ l ρ l { 2 p r -p l + σκ[H(ϕ r ) -H(ϕ l )] a r + ρ l c l + u l -u r } + , (91)
ensures the positivity of ρ * ,θ,N l and ρ * ,θ,N r . In the low Mach regime, the second terms on the RHS of equation ( 90) and (91) can reach a magnitude of O(M) with respect to the characteristic value ρĉ. With the definition in (90) and (91), the acoustic impedance constraints in [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF] are sufficient to preserve the density positiveness. In practice, we use the artificial acoustic impedance defined by [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF].

As the same acoustic impedance is employed for low Mach schemes with AW, PC or the new correction, the same constraints on the time step ∆t are adopted.

No significant difference in terms of CPU time can be observed between the three low Mach corrections. The accuracy of these three methods will be compared in the next section.

Numerical results of low Mach corrections

For the simulations that will be presented in this section, the time step ∆t in order to account for [START_REF] Rieper | A low-Mach number fix for Roe's approximate Riemann solver[END_REF], ( 20), ( 25) and ( 26) is set by following constraint ∆t = min(∆t acoustic , ∆t capillary , ∆t transport , ∆t viscosity ).

(92)

The analysis derived in the previous sections suggests that the acoustic impedance weighted or centered pressure type corrections may generate important errors in the vicinity of the interface while the new low Mach correction should preserve its accuracy. The performance of these schemes will be tested against two classic tests: a two-phase Gresho vortex and a static bubble simulation.

Two-phase Gresho vortex problem

The Gresho vortex [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 2:Implementation[END_REF][START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF] is a rotating flow with a time-independent solution of homogeneous Euler equations. Here we adapt this problem to two-phase flow and the surface tension effects are neglected. An incompressible fluid 1 is placed inside an incompressible fluid 2, with a vortex centered at (x, y) = (0.5, 0.5). By using polar coordinates (r, ϑ), the density distribution can be expressed as:

ρ(r, ϑ) =        ρ 0 1 r < 0.2, ρ 0 2 r > 0.2, (93)
We keep the same angular velocity distribution as proposed in [START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF]:

(u r , u ϑ ) =              (0, 5r) 0 ≤ r < 0.2, (0, 2 -5r) 0.2 ≤ r < 0.4, (0, 0) r ≥ 0.4. (94)
The velocity reaches its maximum value of u ϑ,max = 1 at the interface between two different fluids (r = 0.2). As the centrifugal force is exactly balanced by pressure gradient, the pressure distribution can be given by:

p(r, ϑ) =              p 0 + 12.5ρ 1 r 2 0 ≤ r < 0.2, p 0 + 0.5ρ 1 + 12.5ρ 2 (r 2 -0.04) + 4ρ 2 (1 -5r + ln(5r)) 0.2 ≤ r < 0.4, p 0 + 0.5ρ 1 -2.5ρ 2 + 4ρ 2 ln 2 r ≥ 0.4, (95)
where the reference pressure p 0 satisfies:

p 0 = ρ 0 1 γ 1 M 2 1,max -π ∞ 1 = ρ 0 2 γ 2 M 2 2,max -π ∞ 2 , (96)
where π ∞ 1 and π ∞ 2 are two constant pressures related to the EOS of fluids 1 and 2. As presented in Table . For the Problem #4, a low Mach scheme with AW or CP type correction provides the same results since the initial acoustic impedance ratio of two different fluids is set to 1. A global view of the vortex at t = 0.1 is presented in Fig. 7(b) or Fig. 7(c), we can find that the vortex is maintained. With the velocity magnitude profile on the horizontal line passing through the center of the vortex given in Fig. 8, we can see a more quantitative result: compared with the new correction, the velocity field of either AW or CP type correction is no more smooth, the vortex is distorted.

From the results of the Gresho vortex problem we conclude with a qualitative analysis that either the AW or the CP type low Mach scheme cannot preserve the two-phase vortex with large density ratio as a truncation error of magnitude of

O( 1 α m
) is involved. Let us mention that the results of this test case seem robust in the sense that even if the low Mach initial condition [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF] is not strictly respected, the new low Mach succeeds in preserving the vortex with good accuracy.

Static bubble

To indicate the influence of existing low Mach corrections on scheme temporal stability with capillary effects, the "Static bubble" test case is certainly the first case of rudimentary simulation to implement. The simplicity of this case of simulation makes it possible to well isolate the capillary phenomenon, the theoretical solution being simply dictated by Laplace's law. A 2D bubble with a radius of R = 0.4 containing fluid 1 is placed in fluid 2, which are both incompressible. The initial pressure into fluid 2 is p 0 2 = 1. To satisfy the Laplace's law, the initial pressure in fluid 1 is given by: p 0 1 = p 0 2 +σ/R 0 . The surface tension coefficient of σ is set to 0.6144, and the Laplace number is defined as La = 2ρRσ/µ 2 . For problems with different densities in each phase, we keep the same dynamic viscosity and different Laplace numbers in each fluid. Problems #5-#8 are designed with different fluid properties (see Table 3). With the same design as the Gresho vortex problem, Problem #5 considers two fluids with the same EOS and initial density. Problem #6 aims at checking the ability of low Mach correction that takes into account large density ratios. While Problem #7 takes into account the variation of Mach number while the initial density of each phase remains the same. Problem #8 is a more complicated case with a large density ratio and a large Mach number variation.

Fluid 1 Fluid 2 # ρ 1 π ∞ 1 γ 1 La 1 ρ 2 π ∞ 2 γ 2 La 2 α m β m Static bubble 5
The existence of discretization errors on the curvature act as an initial impulse that makes the bubble oscillate. With the effect of physical viscosity, the velocity amplitude of bubble oscillation decreases and the bubble should return to rest. In order to test the scheme numerical stability, configurations with high Laplace numbers should be considered, corresponding to low physical viscosity. In the present study, the Laplace is set to around 12000. As presented in [START_REF] Popinet | Numerical models of surface tension[END_REF], the redistancing step tends to perturb the curvature and prevents the system from reaching an exact balance. For this test case, this step is not activated.

Numerical results for Problems #5 -#8 are presented in Fig. 9. In this figure, we can find that the numerical results of some configurations are absent. From local zoom, we can find that these absences correspond to the divergence of the numerical result, the scheme stability is not ensured. With results presented in Fig. 9, we can get a global conclusion: numerical schemes without correction and with new low Mach correction reach equilibrium and the velocity fluctuations are reduced to machine precision. Low Mach scheme with AW or CP type correction can not always reach steady state and maintain the scheme stable. For configurations in which the velocity fluctuations are reduced to machine precision, the low Mach scheme can always reduce the velocity fluctuations at a lower rate than without correction, proving that the numerical dissipation is clearly weakened.

For two fluids with the same EOS and initial density (Problem #5), all of the present low Mach corrections show good behavior of recovering fluctuations with a machine precision magnitude. But the numerical results of these low Mach corrections are a little different as the pressure jump across the interface leads to different sound speeds.

For Problem #6 with a large density ratio, the acoustic impedance weighted type and the new low Mach corrections can reach low velocity fluctuations, as they provide the same truncation error for the momentum update which is independent of density ratio α m . While the CP type correction involves a truncation While for the Problem #7, we can observe similar results as Problem #6. At the beginning of the simulation, the numerical scheme with CP type and new low Mach correction can decrease velocity fluctuations. Compared to the CP type correction, the AW type correction provides an extra error of magnitude of O(1) related to the first term of RHS of π * ,θ,AW in (59) as presented in the second line of (B.2), the numerical equilibrium is not observed. For the CP type correction, the magnitude of dimensionless velocity fluctuations starts to decrease down to 10 -8 before suddenly increasing, as the necessary condition is not satisfied according to the asymptotic analysis in (75).

As for Problem #8, both AW and CP types corrections provides the truncation error of magnitude O( 1 α m ). As expected, these two low Mach corrections cannot give satisfactory numerical results.

Conclusion

In this work, we studied the simulation of a compressible two-phase flows model with sharp interface in the low Mach regime. We proposed a definition of the low Mach regime that accounts for possible great discrepancies between the orders of magnitude that are associated with the fluid properties of each phase. We showed that classic low Mach fixes of the literature for Godunovtype solvers provide accurate results on Cartesian grids far from the interface. However, focusing on the numerical method proposed in [START_REF] Zou | An accurate sharp interface method for two-phase compressible flows at low-mach regime, Flow, Turbulence and Combustion[END_REF], we analyzed the behavior of the scheme in the vicinity of the interface. We found out that the abrupt changes of the order of magnitude of the flow parameters across the interface dramatically impact the precision of the method. Indeed, one could see that for an initial condition that verifies the low Mach regime defined in this work, the computed solution will rapidly stop to satisfy the low Mach regime requirements in the vicinity of the interface. We then proposed a new numerical solver that can cope with these important variations and that preserve low Mach regime approximated solutions.

Numerical results were presented for the Gresho vortex test case that showed good performance of the proposed low Mach correction. The new methods were able to successfully capture static bubble equilibrium with a high Laplace number.

Future investigations will focus on improving the accuracy of the method by implementing higher-order discretization, deriving large time step stable methods by implementing an implicit update on the acoustic subsystem and enriching the physical settings of the model by accounting for additional effects like heat and mass transfer.

Appendix C. Subcharacteristic conditions for density positivity

We follow the work presented in [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws: and Well-Balanced schemes for Sources[END_REF]. With intermediate densities presented in (89), the acoustic impedances a l and a r should satisfy the following conditions:

Figure 2 :

 2 Figure 2: Ghost cell representation

Figure 3 :

 3 Figure 3: 1-D sketch of a cell j across the interface Γ, separating domains D i , with its full neighborhood V and partial neighborhoods N ± ( j). Level Set function ϕ, positive at cell j. Discontinuous pressure field p across the interface with non-dimensional values at cells j and k.

6. 1 .

 1 Low Mach flow across the interface: behavior of the pressure near the interface during the acoustic step Let us consider π * ,θ,CP jk the non-dimensional form of corrected flux defined by (60), for k ∈ N -( j) ∅ we have that

  with the same fluid property are considered. Problem #2 is dedicated to testing the ability of the low Mach correction in the case of large density ratio, while the initial sound speeds of different phases remain the same. Problem #3 is devoted to checking the ability of the low Mach correction that takes account of a large Mach variation, while the initial density remains the same at each phase. Problem #4 is a test case with both large density ratio and large Mach number variation.

Figure 4 :

 4 Figure 4: Magnitude of the velocity field of two-phase Gresho vortex problem #1 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

Fig. 4

 4 Fig. 4 shows snapshots of Problem #1 at t = 0.1 for the numerical scheme without and with low Mach correction respectively. Without low Mach correction, the effect of dissipation is obvious and the vortex has completely disappeared. All of the presented low Mach corrections give the same numerical results for problem #1 since the properties of the two fluids are the same. With

Figure 5 :

 5 Figure 5: Magnitude of the velocity field of two-phase Gresho vortex problem #2 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

Figure 6 :

 6 Figure 6: Magnitude of the velocity field of two-phase Gresho vortex problem #3 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

Fig. 6

 6 Fig. 6 shows the numerical results of Problem #3. All these low Mach corrections can well preserve the vortex. For two fluids of the same initial density and large sound speed ratio, all these low Mach corrections provide with errors of magnitude O(1), uniform with respect with M, α m and β m , but the asymptoticpreserving conditions are not strictly respected by AW or CP type correction.

Figure 7 :Figure 8 :

 78 Figure 7: Magnitude of the velocity field of two-phase Gresho vortex problem #4 with a resolution of 100 × 100. (a): without low Mach correction; (b) AW type correction; (c) CP type correction; (d) new low Mach correction.

1 0

 1 Problem #5: α m = 1 β m = 0.63 Problem #6: α m = 0.01 β m = Problem #7: α m = 1 β m = 0Problem #8: α m = 0.01 β m = 0.01

Figure 9 :

 9 Figure 9: Temporal evolution of dimensionless velocity fluctuations of a 2D bubble (La=12000) at resolution of R = 12.8∆x.

  error of magnitude O( 1 α m ), the numerical result shows an immediate growth in velocity fluctuations. When the velocity fluctuations amplitude is reduced to around machine precision, the low Mach scheme with AW type correction starts to diverge. Compared to the new low Mach correction, the AW type correction provides an extra truncation error of magnitude of O(1) on mass and energy update associated with u * , showing the impact of the correction for u * on preserving a constant thermodynamic pressure profile.

Usingα 2 β 2 )

 2 is issued from the first part of RHS of π * ,θ,AW jk in (59) and it may become significantly large. For example when α 2 β 2 ≤ O(1) then the error becomes of magnitude O(α 2 ).

0 . 2 .

 02 ∀ρ ∈ [ρ l , ρ * ,θ,N l ], ρ 2 p ′ (ρ) ≤ a 2 l , ∀ρ ∈ [ρ r , ρ * ,θ,N r ], ρ 2 p ′ (ρ) ≤ a 2 r , (C.1)then the approximate Riemann solver could preserve positiveness of densitiesρ l , ρ r , ρ * ,θ,N l , ρ * ,θ,N r .With the second law of thermodynamics, we haveT ds = c p dT -τd p, (C.2)where s, c p and τ are the entropy, specific heat capacity at constant pressure and specific volume respectively. According to the equation of state, c p can be expressed as: .2) and (C.3), we can obtain that:∆s = c p γ ln p + π ∞ p 0 + π ∞ + c p ln( fluid, we have:p + π ∞ = Eρ γ , (C.5)where E = p 0 +π ∞ ρ γ With (C.5), we can make the following assumptions: An inverse function ζ can be defined that:ρ √ p ′ (ρ) = a ⇔ ρ = ζ(a). (C.7)

Table 1 :

 1 1 and Table. 2, we define Problem #1 -#4 with different fluid properties. Problem #1 is the original test case from[START_REF] Miczek | New numerical solver for flows at various Mach numbers[END_REF] where two fluids Numerical parameters of Gresho vortex test problems.

				Fluid 1			Fluid 2
		#	ρ 1	M 1,max γ 1	ρ 2	M 2,max γ 2	α m	β m
		1	1	10 -3	1.666 1	10 -3	1.666 1	1
	Gresho	2	1	10 -3	1.666 100	10 -3	7.14	0.01	1
	vortex	3	1	10 -3	1.666 1	10 -5	7.14	1	0.01
		4	1	10 -5	1.666 100	10 -3	7.14	0.01	0.01
							# test case
				1		2		3	4
	Gresho vortex	π ∞ 1 π ∞ 2		0 0		0 13405362	0 1399959983 0 5988395358

Table 2 :

 2 Constant pressure π ∞ of Gresho vortex test problems.

Table 3 :

 3 Numerical parameters of static bubble test problems

Acknowledgments

Authors would greatly thank Pascal Tremblin and Thomas Padioleau from Maison de la Simulation for their fruitful discussions and helpful suggestions.

Appendix A. Acoustic relaxation approximation

We present an approximate Riemann solution to solve the acoustic subsystem within the Eulerian framework. The acoustic subsystem ( 6) is a quasilinear system. In order to derive the resolution of this system, we will perform several approximations.

We notice that for a smooth solution [START_REF] Liu | The ghost fluid method for compressible gaswater simulation[END_REF], the pressure verifies: ∂ t p + (ρc) 2 τ∇ • u = 0. We thus perform a Suliciu-type relaxation adapted to vacuum [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws: and Well-Balanced schemes for Sources[END_REF] of [START_REF] Liu | The ghost fluid method for compressible gaswater simulation[END_REF] by introducing a surrogate pressure π and considering following relaxed system by using non-conservative variables with τ = 1/ρ:

In regime χ → ∞, we formally recover [START_REF] Liu | The ghost fluid method for compressible gaswater simulation[END_REF]. In our numerical solver context, we classically mimic χ → ∞ regime by enforcing at each time step π n i = p EOS (τ n i , e n i ) (given by ( 1)) and then solving (A.1) with χ = 0.

Assuming a one-dimensional problem in x-direction, we have ρ(x, t)∂ t ≈ ρ(x, t n )∂ t , then if we define a mass variable m by dm(x) dx = ρ(x, t n ), we obtain up to a slight abuse of notation, a system that is written in vector form:

and its eigenvalues are: (-a, 0, 0, 0, 0, 0, 0, +a). Fields involved in this system are all linearly degenerated. Considering a discontinuity that propagates at celerity D, let ⟦b⟧ be the jump of the variable b across the interface. Following the Rankine-Hugoniot relationships we get the jump conditions:

with D taking the same value as eigenvalues of matrix A. The parameters q and q ′ are weight associated with the Dirac masses M and Mu, their definition will be given latter. D = 0 yields that u is an independent Riemann invariant for the wave λ = 0 and brings in the jump conditions ⟦π⟧ = q⟦ϕ⟧ and ⟦a⟧ ∈ R in (A.2c) and (A.2e). Thus a is not required to be continuous across the wave λ = 0, and the jump condition on π should conform the Laplace equation (3).

By ( 65), (66), (67), (69), ( 70) and (59), tedious calculations show that

Then we can see that when α 2 becomes large, it may produce important errors for the estimate of the pressure gradient in Ω j . Indeed, let us consider

With assumptions in (C.6), we can derive the following relationship:

According to the monotonicity of ζ and (89), the conditions (C.1) can be expressed as:

.

(C.10)

These conditions (C.10) include the positivity of ρ * ,θ,N l and ρ * ,θ,N r and can be satisfied with the following acoustic impedances:

(C.12)

Proof:

Here we give the proof for the constraint on acoustic impedance a l . The definition of a l always satisfies the first condition at the first line of (C.10). To derive required constraint on the second condition, we define a variable X

If X ≤ 0, the second condition at the first line of (C.10) is always satisfied. Assume that X ≥ 0, then a l can be rewritten as: a l = ρ l ( √ p ′ (ρ l ) + ΥX). Multiplying by ρ l , the second condition at first line of (C.10) reads:

By defining a variable ϖ that:

with (C.9), we can always have:

where 0 < ϖ ≤ 1 and Υ ≥ 1. Thus the definition of a l can always satisfy (C.14), the positivity of ρ * ,θ,N l is then guaranteed. With the same analysis, we can justify that constraints on a r presented in (C.11) and (C.12) can preserve the positivity of ρ * ,θ,N r .