S Del Pino

E Labourasse
email: emmanuel.labourasse@cea.fr

G Morel

An asymptotic preserving multi-dimensional ALE method for a system of two compressible flows coupled with friction

Keywords: Compressible gas dynamics, multi-fluid, finite-volumes, unstructured meshes, asymptotic preserving, arbitrary-Lagrangian-Eulerian (ALE)

We present a multi-dimensional asymptotic preserving scheme for the approximation of a mixture of compressible flows. Fluids are modeled by two Euler systems of equations coupled with a friction term.

The asymptotic preserving property is mandatory for this kind of model, to derive a scheme that behaves well in all regimes (i.e. whatever the friction parameter value is). The method we propose is defined in ALE coordinates, using a Lagrange plus remap approach. This imposes a multi-dimensional definition and analysis the scheme.

Introduction

A multifluid model is a model for a fluid mixture for which each fluid is described by is own full set of variables (for instance density, velocity and energy). The model is generally closed in a way that defines interactions between the constituents, depending on the envolved physic. These models are widely used in different communities. One very popular model of this kind is the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] for deflagration-to-detonation transition of reactive flows.

Many numerical methods to approximate this model have been designed, let us just cite a few of them [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF][START_REF] Abgrall | An Asymptotic Preserving Scheme for the Barotropic Baer-Nunziato Model[END_REF][START_REF] Ambroso | Interface model coupling via prescribed local flux balance[END_REF]. Such kind of models is also used in plasma physics to account for plasmas collision or Non-Local-Thermodynamic-Equilibrium (NLTE) Ion-Electron interactions [START_REF] Sentis | Mathematical Models and Methods for Plasma Physics[END_REF]. Scannapieco and Cheng [8] also derive similar kind of model for turbulent flows and apply it to describe a mixing zone driven by Rayleigh-Taylor or Richtmyer-Meshkov instabilities [START_REF] Cheng | Buoyancy-drag mix model obtained by multifluid interpenetration equations[END_REF].

In this paper, we present a multi-dimensional scheme to approximate solutions of two compressible inviscid fluids coupled with friction, refer to equation [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF]. This model is a slightly simplified version of the Scannapieco-Cheng [START_REF] Scannapieco | A multifluid interpenetration mix model[END_REF] model where friction is considered uniform in space. It can also be viewed as a simplification of the model proposed in [START_REF] Sentis | Mathematical Models and Methods for Plasma Physics[END_REF] for which the electron effect is neglected or of the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], neglecting the interfacial terms and in the case where there is no phase transition. scribes the Lagrangian step of the proposed scheme. It is demonstrated that the scheme preserves the properties of conservation, stability and consistency with respect to the continuous model for all regimes (independantly of the value of the friction parameter). Then in section 5, our ALE strategy is described. Finally, section 6 is devoted to numerical experiments on several problems (Sod shock tube, triple point and Rayleigh-Taylor). Some comparisons with a non-AP scheme are provided.

A two fluids model with friction

Let us consider a mixture of two fluids f 1 and f 2 . In the following, we will denote by "multifluid model", a model for which each fluid α ∈ { f 1 , f 2 } is represented by its own set of variables: (ρ α , u α , E α ). Conversely, we will refer as "mono-fluid model", a model describing a mixture where mean quantities are considered (ρ, u, E), each fluid position being precised by an additional equation on the concentration (e.g. c α = ρ α ρ α +ρ β ).

In this part, we present a simplified version of Scannapieco-Cheng's model where the interaction between the two constituents reduces to a friction term. In semi-Lagrangian coordinates, for each fluid α ∈ { f 1 , f 2 } (β denoting the other fluid), the model writes

ρ α D α t τ α = ∇ • u α , ρ α D α t u α = -∇p α -νρδu α , ρ α D α t E α = -∇ • (p α u α ) -νρδu α • u, (1)
where ρ α , u α and E α respectively denote the mass density, the velocity and the total energy density of fluid α. Also, τ α = 1 ρ α denotes the specific volume. The pressure p α satisfies the equation of state p α := p α (ρ α , e α ), where e α , the internal energy density, is defined by e α := E α -1 2 u α 2 . The total density ρ and the mean velocity u are defined as ρ := ρ α + ρ β and ρu := ρ α u α + ρ β u β . The term δu α is the velocity difference, the δ(•) α operator being defined by δφ α = -δφ β = φ αφ β . Finally, ν is the friction parameter. Also, remark that the Lagrangian derivative D α t := ∂ t + u α • ∇, is obviously not the same for each fluid.

The entropy η α defined by Gibbs formula T α dη α = de α +p α dτ α satisfies the following entropy inequality

T α D α t η α ≥ ν τ α τ β δu α • δu α ≥ 0. ( 2
)
Prior to establishing a numerical scheme that discretizes this set of six equations, we recall some properties of the model itself.

Property 1 (Conservation). The model ( 1) is conservative in volume and mass for each fluid.

Also, it is conservative in the sum of momenta and in the sum of the total energies of the two fluids.

Proof. Conservation of mass and volume is obvious since the first equation of ( 1) is the continuity equation written for each fluid.

Conservation of momenta sum and total energies sum requires more cautiousness, since Lagrangian derivative are not the same for each fluid. To establish them one rewrites [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] in an Eulerian framework.

Developing Lagrangian derivatives and using the identity ∂ t (ρ α τ α ) = 0 elementary calculations allow to rewrite [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] as

∂ t ρ α + ∇ • (ρ α u α ) = 0, ∂ t (ρ α u α ) + ∇ • (ρ α u α ⊗ u α ) + ∇p α + νρδu α = 0, ∂ t (ρ α E α ) + ∇ • (ρ α E α u α ) + ∇ • (p α u α ) + νρδu α • u = 0. (3)
Summing the two later equations over α gives a system of the conservative form ∂ t U + ∇ • F(U) = 0, where

U = ρ α u α + ρ β u β ρ α E α + ρ β E β , and
F(U) = ρ α u α ⊗ u α + ρ β u β ⊗ u β + p α + p β I ρ α E α u α + ρ β E β u β + p α u α + p β u β ,
where I is the identity matrix of R 2×2 .

Property 2 (Hyperbolicity). The model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] is hyperbolic.

Proof. The proof is straightforward but calculatory, see [START_REF] Enaux | Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles[END_REF] for details.

Asymptotic model. When ν → +∞, (1) behaves has the following five equations model

ρD t u = -∇ p α + p β , (4)
while, for each fluid α ∈ { f 1 , f 2 }, β denoting the other one, one has

ρ α D t τ α = ∇ • u, ρ α D t E α = - ρ α ρ u • ∇ p α + p β -p α ∇ • u, ( 5
)
where u is the same velocity for both fluids, and thus the Lagrangian derivative is also the same.

Formal derivation (established in [START_REF] Enaux | Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles[END_REF]). Let = ν -1 so that (1) rewrites

ρ α D α t τ α = ∇ • u α , ρ α D α t u α = -∇p α - 1 ρδu α , ρ α D α t E α = -∇ • (p α u α ) - 1 ρδu α • u. (6)
We will now study its limit while → 0 + focusing first on the momentum equations since the friction term's goal is to impose that δu 0 →0 -→ 0.

Developing the Lagrangian derivatives and dividing each momentum equation by ρ α > 0, one has

∂ t u α + (∇u α ) u α = - ∇p α ρ α - 1 ρ ρ α δu α .
Since fluid β satisfies the same equation and recalling that δφ α = -δφ β = φ αφ β , one gets

∂ t (δu α ) + δ ((∇u) u) α = -δ ∇p ρ α - 1 λδu α , where λ = ρ 2 ρ α ρ β .
We now perform an Hilbert expansion for all variables in the equation, that is

φ = φ 0 + φ 1 + O( 2 ). One has ∂ t (δu α,0 ) + δ ((∇u)u) α,0 = -δ ∇p ρ α,0 -λ 0 1 δu α,0 + δu α,1 -λ 1 δu α,0 + O( ). (7)
Multiplying this equation by one has λ 0 δu α,0 = O( ), which gives δu α,0 = 0 when → 0 since λ > 0.

So, when → 0, formula [START_REF] Sentis | Mathematical Models and Methods for Plasma Physics[END_REF] recasts

δu α,1 = - 1 λ 0 δ ∇p ρ α,0 . (8)
Now, we perform an Hilbert expansion for the whole system (6), neglecting the non negative powers of . Choosing α ∈ { f 1 , f 2 }, β being the other one, it reads

ρ α,0 D α t τ α,0 =∇ • u α,0 , ρ α,0 D α t u α,0 = -∇p α,0 -ρ 0 1 δu α,0 + δu α,1 -ρ 1 δu α,0 , ρ α,0 D α t E α,0 = -∇ • p α,0 u α,0 -ρ 0 1 δu α,0 • u α,0 + δu α,1 • u α,0 + δu α,0 • u α,1 -ρ 1 δu α,0 • u α,0 ,
Since we just established δu α,0 = 0, one has u 0 = u 0 = u α,0 = u β,0 . Also, since D α t φ = ∂ t φ + u α,0 • ∇φ + O( ), Lagrangian derivatives are the same when → 0, so that using (8) the system simplifies to

ρ α,0 D t τ α,0 =∇ • u 0 , ρ α,0 D t u 0 = -∇p α,0 + ρ 0 1 λ 0 δ ∇p ρ α,0 , ρ α,0 D t E α,0 = -∇ • p α,0 u 0 + ρ 0 1 λ 0 δ ∇p ρ α,0 • u 0 , Recalling λ = ρ 2 ρ α ρ β and developing δ ∇p ρ α,0
, momentum equation rewrites

ρ α,0 D t u 0 = -∇p α,0 + ρ α,0 ρ β,0 ρ 0 ∇p α,0 ρ α,0 - ∇p β,0 ρ β,0 , = - ρ α,0 ρ 0 ∇ p α,0 + p β,0 .
Proceeding the same way with total energy equation, one gets

ρ α,0 D t E α,0 = -∇ • p α,0 u 0 + ρ α,0 ρ β,0 ρ 0 ∇p α,0 ρ α,0 - ∇p β,0 ρ β,0 • u 0 , = - ρ α,0 ρ 0 ∇p α,0 + ∇p β,0 • u 0 -p α,0 ∇ • u 0 , Remark 1. Defining E := ρ α E α +ρ β E β ρ
and τ := ρ -1 , it is easy to check that if (ρ α , ρ β , u, E α , E β ) is a solution of the asymptotic model ( 4)-( 5), one has

ρD t τ = ∇ • u, ρD t u = -∇ p α + p β , ρD t E = -∇ • p α + p β u .
One recognizes Euler equations for the mixture. The mixing pressure follows Dalton's law as one could have expected since we consider here non-reactive gases.

However, notice that unless each fluid follows a barotropic equation of state (p α = p α (ρ α )), equation ( 5) must be solved to determine e α .

Cell-centered schemes

We recall briefly the muti-dimensional finite volume schemes [START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique lagrangienne bidimensionnelle[END_REF][START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF], since it is the basis of this work. For convinience, we use the notations defined in [START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]. In the following, for all cell j, and for any quantity φ, one defines its mean value φ j := 1 V j j φ, where V j := j 1 is the cell volume. Also, let us denote the cell's mass as m j := j ρ = ρ j V j , which is constant in time in semi-Lagrangian coordinates (d t m j = 0). We consider first-order schemes, so that one has the following relations

d dt j 1 = m j d t τ j , d dt j ρ j = 0, d dt j ρ j u j = m j d t u j , d dt j ρ j E j = m j d t E j .
Let J r denote the set of cells connected to node r and let R j the set of nodes of cell j. Also, let us introduce C jr := ∇ x r V j , the gradient of the volume of the polygonal cell j, according to the position of one of its vertices r. Then, the cell-centered schemes we consider in this paper have the following structure: for any cell j of the mesh one has

m j d t τ j = r∈R j C jr • u r , d t m j =0, m j d t u j = - r∈R j F jr , m j d t E j = - r∈R j F jr • u r , (9)
where the fluxes u r and F jr are defined for any node r

∀ j ∈ J r , F jr = C jr p j -A jr (u r -u j ), (10)
and

j∈J r F jr = 0. (11)
In one hand, relation [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF] is the matrix form of the acoustic Riemann solver (see for instance [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF][START_REF] Maire | Contribution à la modélisation numérique de la fusion par confinement inertiel[END_REF]), while in the other hand [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] imposes conservation.

In the following to simplify notations, we omit sets R j and J r when there is no confusion.

• If A jr := ρ j c j C jr ⊗C jr C jr , then ( 9)-( 11) defines the Glace scheme [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF].

• Let N + jr = -1 2 (x r+1 -x r ) ⊥ and N - jr = -1 2 (x r -x r-1 ) ⊥ . If A jr := ρ j c j N + jr ⊗N + jr N + jr + N - jr ⊗N - jr N -

jr

, the scheme ( 9)-( 11) is Eucclhyd [START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF][START_REF] Maire | Contribution à la modélisation numérique de la fusion par confinement inertiel[END_REF]. One has N + jr + N - jr = C jr , see figure 1.

These schemes are conservative in volume, mass, momentum and total energy. One easily shows that they are entropy stable. These results can be found in [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF][START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | Contribution à la modélisation numérique de la fusion par confinement inertiel[END_REF], for instance.

Also, a consistency result has been established in [START_REF] Després | Weak consistency of the cell-centered lagrangian glace scheme on general meshes in any dimension[END_REF].

Asymptotic Preserving scheme in semi-Lagrangian coordinates

We shall now present a multi-dimensional finite-volume scheme written in semi-Lagrangian coordinates that preserves the asymptotic.

In this section, we present the Lagrangian step of our ALE method. In this step, each fluid is associated to its own mesh. If the meshes may evolve differently, we assume that they coincide at the begining of the Lagrangian step. The rezoning/remapping procedure that is detailed in section 5 is used to ensure that the meshes will coincide for the next Lagrangian step.

We first focus on the semi-discrete continuous in time scheme. Most of the properties of the scheme are proved using this simpler formulation without any lost of generality. In paragraph 4.2, we describe the fully discrete scheme. It is analysed in the remaining of this section.

Continuous in time semi-discrete scheme

Let ω ∈ [0, 2], for each fluid α ∈ { f 1 , f 2 }, β denoting the other fluid, we define the scheme

m α j d t τ α j = r C jr • u α r , d t m α j =0, m α j d t u α j = - r F α jr -ω r νρ r B jr δu α j -(1 -ω) r νρ r B jr δu α r , m α j d t E α j = - r F α jr • u α r - r νρ r u T r B jr δu α r + ω r νρ r u T jr B jr (δu α r -δu α j ), (12)
where the fluxes are given by

F α jr = C jr p α j -A α jr (u α r -u α j ) -νρ r B jr δu α r , and (13)
j F α jr = 0. (14)
In order to write [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF], we introduced ρ α r := 1 #J r j∈J r ρ α j and ρ r := ρ α r + ρ β r . Also, we set

u r := ρ α r u α r +ρ β r u β r ρ α r +ρ β r
and u jr :=

ρ α r u α j +ρ β r u β j ρ α r +ρ β r
. B jr are symmetric and positive definite matrices such that r∈R j B jr = V j I. Matrices A α jr are the standard "hydro-matrices" as defined in section 3.

Remark 2. One can choose B jr := V jr I, where V jr is the volume of the subcell associated to vertex r of cell j. Another obvious choice could be for instance B jr := 1 #R j V j I.

Observe that simple calculations allow to write

ρ r u r = ρ r u α r -ρ β r δu α r and ρ r u jr = ρ r u α j -ρ β r δu α j . (15)
Injecting ( 13) in [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF], and using [START_REF] Enaux | Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles[END_REF], one gets the alternative form

m α j d t τ α j = r C jr • u α r , d t m α j =0, m α j d t u α j = r A α jr (u α r -u α j ) + ων r ρ r B jr δu α r -δu α j , m α j d t E α j = - r C jr p α j • u α r + r u α r T A α jr (u α r -u α j ) + ν r ρ β r t δu α r B jr δu α r -ων r ρ β r t δu α j B jr (δu α r -δu α j ) + ων r ρ r t u α j B jr (δu α r -δu α j ). ( 16
)
This form enlightens the fact that knowing the fluxes (u α r , u β r ) at any vertex r is enough to define the scheme. We shall now show that these nodal velocities are well defined.

Injecting [START_REF] Colella | A conservative finite difference method for the numerical solution of plasma fluid equations[END_REF] in ( 14) allows to calculate (u α r , u β r ). Obviously, as soon as ν 0, both nodal velocities are coupled at vertex r. Omitting boundary conditions in the sake of simplicity, for each vertex of the mesh (u α r , u β r ) is the unique solution of the following linear system:

j       A α jr + νρ r B jr -νρ r B jr -νρ r B jr A β jr + νρ r B jr       A ν r := u α r u β r = j       A α jr u α j + C jr p α j A β jr u β j + C jr p β j       b r := .
Proof. Since matrices A α jr and B jr are symmetric, A ν r is symmetric. To prove that (u α r , u β r ) is unique, it remains to show that it is positive definite. Elementary calculation gives,

∀(v α , v β ) ∈ R 2 × R 2 , (v α , v β ) T A ν r (v α , v β ) =v αT         j A α jr         v α + v β T         j A β jr         v β + (v α -v β ) T         j νρ r B jr         (v α -v β ),
which is strictly positive if (v α , v β ) (0, 0) since matrices j A α jr and j νρ r B jr are positive definite.

The scheme being well-defined, we now establish its properties.

Nodal velocities a priori estimates

Here, we establish estimates for the nodal velocities with regard to the frictionless case. These are actually some instantaneous stability results with regard to the mono-fluid schemes [START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF],

i.e. velocity fluxes are controled by the frictionless ones.

Property 3 (A priori estimates). For each fluid α ∈ { f 1 , f 2 }, let u α,ν r denote the nodal velocities at vertex r. Let A α r := j A α jr and B r := j B jr . Let β denote the other fluid, then one has the following relations, ∀ν ≥ 0

u α,ν r T A α r u α,ν r + u β,ν r T A β r u β,ν r ≤ u α,0 r T A α r u α,0 r + u β,0 r T A β r u β,0 r , (17)
u α,ν r -u β,ν r T B r u α,ν r -u β,ν r ≤ 1 2νρ r u α,0 r T A α r u α,0 r + u β,0 r T A β r u β,0 r , (18)
and

u α,ν r -u β,ν r T B r u α,ν r -u β,ν r ≤ u α,0 r -u β,0 r T B r u α,0 r -u β,0 r . ( 19
)
Let us first comment these estimates. The estimate ( 17) is a stability results. It shows that the nodal velocity (u α,ν r , u β,ν r ) A 0 r is bounded by (u α,0 r , u β,0 r ) A 0 r independently of ν. It shows that friction nodal velocities are stable with regard to the classic frictionless case for a given state.

The second estimate [START_REF] Franck | Construction et analyse numérique de schéma asymptotic preserving sur maillages non structurés. application au transport linéaire et aux systèmes de friedrichs[END_REF] shows that the nodal velocity difference δu α,ν

r B r is at most O(ν -1/2 )
according to (u α,0 r , u β,0 r ) A 0 r .

The last inequality [START_REF] Franck | Modified Finite Volume Nodal Scheme for Euler Equations with Gravity and Friction[END_REF] states that the nodal velocity difference is bounded by the frictionless case independently of ν in the • B r norm, which is purely geometric.

Proof of Property 3. ∀ν ≥ 0, (u α,ν r , u β,ν r ) is the unique solution of

A α r + νρ r B r -νρ r B r -νρ r B r A β r + νρ r B r u α,ν r u β,ν r = b r , with b r :=       j C jr p α j j C jr p β j       . So, since b r is independent of ν, one has ∀ν ≥ 0, A 0 r + νρ r ∆ r u ν r = A 0 r u 0 r , (20)
where

A 0 r := A α r 0 0 A β r , ∆ r := B r -B r -B r B r and u ν r := u α,ν r u β,ν r .
Multiplying on the left by

u ν r yields u ν r T A 0 r u ν r + νρ r u ν r T ∆ r u ν r = u ν r T A 0 r u 0 r .
Since B r is a positive matrix ∆ r is also positive, and since νρ r ≥ 0, one gets

∀ν ≥ 0, u ν r T A 0 r u ν r ≤ u ν r T A 0 r u 0 r .
Finally, A 0 r being symmetric and positive definite, the simple following Youngs inequality,

u ν r T A 0 r u 0 r ≤ 1 2 u ν r T A 0 r u ν r + 1 2 u 0 r T A 0 r u 0 r ,
allows to prove [START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF].

The proof of ( 18) follows the same way. Multiplying [START_REF] Franck | Finite volume scheme with local high order discretization of the hydrostatic equilibrium for the euler equations with external forces[END_REF] on the left by u ν r , one has

∀ν ≥ 0, νρ r u ν r T ∆ r u ν r + u ν r T A 0 r u ν r = u ν r T A 0 r u 0 r .
Then, using the same Youngs inequality, one gets after a few arrangements

∀ν ≥ 0, νρ r u ν r T ∆ r u ν r + 1 2 u ν r T A 0 r u ν r ≤ 1 2 u 0 r T A 0 r u 0 r ,
which yeilds to (18) since A 0 r is positive.

The third inequality is a bit more difficult to establish. Let us introduce the quadratic form

J ν v := 1 2 v T A 0 r + νρ r ∆ r v -b r • v.
So, since u ν r is the unique solution of the linear system, one has

∀ν ≥ 0, ∀v, J ν u ν r ≤ J ν v .
In the particular case v = u 0 r , one gets J ν

u ν r ≤ J ν u 0 r
. It is then easy to check that

J ν u 0 r = 1 2 u 0 r T A 0 r + νρ r ∆ r u 0 r -b r • u 0 r = J 0 u 0 r + νρ r 2 u 0 r T ∆ r u 0 r .
So, one has established a first inequality

J ν u ν r ≤ J 0 u 0 r + νρ r 2 u 0 r T ∆ r u 0 r . (21)
Similarly, since u 0 r is the unique solution of the linear system in the case ν = 0, one has J 0

u 0 r ≤ J 0 u ν
r , which can be written as

J 0 u 0 r ≤ J ν u ν r - νρ r 2 u ν r T ∆ r u ν r .
This actually gives a lower bound to J ν u ν r which combined with its upper bound (21) yields

J 0 u 0 r + νρ r 2 u ν r T ∆ r u ν r ≤ J 0 u 0 r + νρ r 2 u 0 r T ∆ r u 0 r .
Since νρ r is positive, elementary calculations allow to write [START_REF] Franck | Modified Finite Volume Nodal Scheme for Euler Equations with Gravity and Friction[END_REF].

Conservativity

Property 4 (Conservation). The scheme defined by ( 12)-( 14) ensures conservation of mass and volume for each fluid α or β. It also ensures that the sum of the fluids' momenta and total energies are conserved.

Proof. Conservations of mass and volume for each fluid are obvious since the associated balance equations are unchanged with regard to the mono-fluid schemes (see for instance [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Carré | A cell-centered lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF][START_REF] Maire | Contribution à la modélisation numérique de la fusion par confinement inertiel[END_REF]).

Summing momenta equations in [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF] for both fluids gives

m α j d t u α j + m β j d t u β j = - r F α jr - r F β jr -ω r νρ r B jr (δu α j + δu β j ) -(1 -ω) r νρ r B jr (δu α r + δu β r ).
Recalling that by definition, δu α j + δu β j = 0, one has

m α j d t u α j + m β j d t u β j = - r F α jr - r F β jr .
The conservativity proof is ended in a standard way. One now sums these equations over the cells which gives

j m α j d t u α j + j m β j d t u β j = - j r∈R j F α jr - j r∈R j F β jr , which rewrites, j m α j d t u α j + j m β j d t u β j = - r j∈J r F α jr - r j∈J r F β jr .
This proves that momenta sum is conserved using ( 14) and recalling that cell masses are Lagrangian.

Conservation of total energies sum is obtained in the exact same way.

Stability

Before proving this result, let us recall that the fully discrete scheme's stability is presented bellow (see paragraph 4.2).

Property 5 (Entropy). The first-order continuous in time scheme defined by ( 12)-( 14) satisfies, ∀ω ∈ [0, 2], the following entropy inequality

∀α ∈ { f 1 , f 2 } m α j T α j d t η α j ≥ 1 - ω 2 r νρ β r t δu α r B jr δu β r + ω 2 r νρ β r t δu α j B jr δu α j ≥ 0.
This inequality is consistent with (2).

Let us establish a simple technical Lemma that will be useful in the following and to demonstrate Property 5.

Lemma 1. Let M denote a symmetric matrix of R d×d . Let ω ∈ R, then ∀v, w ∈ R d , v T Mv -ωw T M(v -w) = 1 - ω 2 v T Mv + ω 2 w T Mw + ω 2 (w -v) T M(w -v). Proof. Let ξ := v T Mv -ωw T M(v -w).
Obviously, one has

ξ = v T Mv + ωw T Mw -ωw T Mv. Since M is symmetric, one has -2w T Mv = (v -w) T M(v -w) -v T Mv -w T Mw. Injecting this
equality in the expression of ξ ends the demonstration.

Corollary 1. Let M denote a symmetric and positive matrix of R d×d . Let ω ≥ 0, then

∀v, w ∈ R d , v T Mv -ωw T M(v -w) ≥ 1 - ω 2 v T Mv + ω 2 w T Mw.
Proof. This is a direct consequence of Lemma 1, since ωM is a positive matrix.

We can now give the proof of Property 5.

Proof of Property 5. Gibbs formula reads T dη = de + pdτ, so that one has

T α j d t η α j = d t e α j + p α j d t τ α j ,
which rewrites also

m α j T α j d t η α j = m α j d t E α j -u α j • m α j d t u α j + p α j m α j d t τ α j .
Using [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF], one gets

m α j T α j d t η α j = - r C jr p α j • u α r + r u α r T A α jr (u α r -u α j ) + ν r ρ β r t δu α r B jr δu α r -ων r ρ β r t δu α j B jr (δu α r -δu α j ) + ων r ρ r t u α j B jr (δu α r -δu α j ) + u α j •        r A α jr (u α r -u α j ) + ων r ρ r B jr δu α r -δu α j        + r C jr • u α r p α j ,
which simplifies as

m α j T α j d t η α j = r (u α r -u α j ) T A α jr (u α r -u α j )+ν r ρ β r δu α r T B jr δu α r -ων r ρ β r δu α j T B jr (δu α r -δu α j ).
Since B jr matrices are symmetric and positive and since ω ≥ 0, one can apply Corollary 1 to obtain

m α j T α j d t η α j ≥ r (u α r -u α j ) T A α jr (u α r -u α j ) + 1 - 1 2 ω ν r ρ β r δu α r T B jr δu α r + 1 2 ων r ρ β r δu α j T B jr δu α j ,
Matrix A α jr being positive, one finally has

m α j T α j d t η α j ≥ 1 - 1 2 ω ν r ρ β r δu α r T B jr δu α r + 1 2 ων r ρ β r δu α j T B jr δu α j ,
which is positive as soon as ω ∈ [0, 2].

Asymptotic preserving

We now establish the main result of this paper. It consists in stating that when the friction parameter ν tends to infinity, the scheme ( 12)-( 14) behaves asymptotically as a scheme that is consistent with the asymptotic model ( 4)- [START_REF] Abgrall | An Asymptotic Preserving Scheme for the Barotropic Baer-Nunziato Model[END_REF].

To this end, we first compute the asymptotic scheme by means of Hilbert expansions, then we show its consistency with the asymptotic model. This later result relies strongly on B. Després's work [START_REF] Després | Weak consistency of the cell-centered lagrangian glace scheme on general meshes in any dimension[END_REF].

Asymptotic scheme. Let ω 0. If ∀α ∈ { f 1 , f 2 }, ∀ j, (ρ α j , u α j , E α j )
are constant cell data, then the scheme ( 12)-( 14), behaves asymptotically as

(m α j + m β j )d t u j = - r F α jr - r F β jr , (22)
d t V j = m α j d t τ α j = r C jr • u r , ( 23
)
d t m α j = 0, m α j d t E α j = - r C jr p α j • u r + r u T r A α jr (u r -u j ) - ρ α j ρ β j ρ j r u T j δ A jr ρ j α (u r -u j ), ( 24
)
where u j = u α j = u β j , and where nodal velocities u r = u α r = u β r satisfy

F α jr + F β jr = C jr p α j + p β j -A α jr + A β jr (u r -u j ), and j F α jr = 0. ( 25
)
Formal derivation. Let α ∈ { f 1 , f 2 }, β denoting the other fluid. Let us introduce := ν -1 . One rewrites [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] as

m α j d t τ α j = r C jr • u α r , (26)
d t m α j =0, m α j d t u α j = r A α jr (u α r -u α j ) - 1 ω r ρ r B jr (δu α j -δu α r ), (27)
m α j d t E α j = - r C jr p α j • u α r + r t u α r A α jr (u α r -u α j ) + 1 r ρ β r (δu α r ) T B jr δu α r + 1 ω r (ρ r u α j T -ρ β r δu α j T )B jr (δu α r -δu α j ), (28) and j
A α jr u α r + j 1 ρ r B jr δu α r = j A α jr u α j + j C jr p α j . (29)
Following the analysis of the asymptotic model, we perform an Hilbert expansion.

The first information one gets is from equation ( 29) that writes

j A α,0 jr u α,0 r + j 1 ρ 0 r B jr δu α,0 r + j ρ 0 r B jr δu α,1 r + j ρ 1 r B jr δu α,0 r = j A α,0 jr u α,0 j + j C jr p α,0 j + O( ),
so that multiplying this equation by leads to ρ 0 r ( j B jr )δu α,0 r = 0 that is

δu α,0 r = 0, (30)
since j B jr is symmetric positive definite and ρ r = ρ 0 r + O( ) > 0 so that ρ 0 r > 0 when → 0.

One gets volume conservation equation [START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique lagrangienne bidimensionnelle[END_REF]. Now, the momentum equation ( 27) is considered, using [START_REF] Després | Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques[END_REF], one has

m α j d t u α,0 j = r A α,0 jr (u α,0 r -u α,0 j ) - 1 ω r ρ 0 r B jr δu α,0 j -ω r ρ 1 r B jr δu α,0 j -ω r ρ 0 r B jr (δu α,1 j -δu α,1 r ) + O( ),
which gives

δu α,0 j = 0. (31)
Using, ( 31) and ( 30), one defines u 0 j := u α,0 j = u β,0 j and u 0 r := u α,0 r = u β,0 r .

So, Hilbert expansions of equations ( 26), ( 27) and ( 28) simplify as

m α j d t τ α,0 j = r C jr • u 0 r , m α j d t u 0 j = r A α,0 jr (u 0 r -u 0 j ) -ω r ρ 0 r B jr (δu α,1 j -δu α,1 r ), (32)
m α j d t E α,0 j = - r (C jr p α,0 j -A α,0 jr (u 0 r -u 0 j )) • u 0 r + ω r ρ 0 r u α,0 j T B jr (δu α,1 r -δu α,1 j ), ( 33
)
Our aim is now to evaluate the term ω r ρ 0 r u α,0 j T B jr (δu α,1 r -δu α,1 j ). To do so, we divide momentum equation ( 32) by ρ α j (> 0), which gives

V j d t u 0 j = 1 ρ α j r A α,0 jr (u 0 r -u 0 j ) -ω r ρ 0 r ρ α j B jr (δu α,1 j -δu α,1 r ).
The same relation can be written for fluid β. The difference of these two equations writes, recalling that δφ α = -δφ β ,

0 = r δ        A 0 jr ρ j        α (u 0 r -u 0 j ) - ρ j ρ α j ρ β j ω r ρ 0 r B jr (δu α,1 j -δu α,1 r ).
Injecting this relation in [START_REF] Brull | Asymptotic-preserving scheme for a bi-fluid euler-lorentz model[END_REF] gives the limit scheme total energy balance equation [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF]. The momentum equation ( 22) is obtained the same way or by simply summing equations [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] for both fluids α and β.

In order to establish that the scheme is asymptotic preserving, it remains to show that the limit scheme ( 22)-( 25) is consistent with the asymptotic model ( 4)- [START_REF] Abgrall | An Asymptotic Preserving Scheme for the Barotropic Baer-Nunziato Model[END_REF].

Before establishing this result, we recall the fundamental result by B. Després [START_REF] Després | Weak consistency of the cell-centered lagrangian glace scheme on general meshes in any dimension[END_REF], that we adapt to the present context.

Property 6 (B. Després). Let m j := m α j + m β j , ρ j := ρ α j + ρ β j , τ j = ρ -1 j and E j :=

ρ α j E α j +ρ β j E β j ρ j
. Then, the scheme

d t m j = 0, m j d t τ j = r C jr • u r , m j d t u j = - r F jr , m j d t E j = - r F jr • u r , where F jr = C jr (p α j + p β j ) -(A α jr + A β jr )(u r -u j ), and j (A α jr + A β jr )u r = j (A α jr + A β jr )u j + j C jr (p α j + p β j ),
is weakly consistent with the following system of equations

ρD t τ = ∇ • u, ρD t u = -∇(p α + p β ), ρD t E = -∇ • (p α + p β )u.
Proof. The proof can be found in [START_REF] Després | Weak consistency of the cell-centered lagrangian glace scheme on general meshes in any dimension[END_REF].

Property 7. The limit scheme ( 22)-( 25) is weakly consistent with the asymptotic model ( 4)- [START_REF] Abgrall | An Asymptotic Preserving Scheme for the Barotropic Baer-Nunziato Model[END_REF].

Proof. Consistency for volume, mass and momentum is a direct consequence of Property 6, it remains to show the consistency for total energy.

We rewrite equation ( 5) using a more convenient form

ρ α D t E α = -∇ • (p α + p β )u + p β ∇ • u + ρ β ρ ∇(p α + p β ) • u.
As a starting point we recall [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF] for fluid α

m α j d t E α j = - r C j p α j • u r + r u T r A α jr (u r -u j ) - ρ α j ρ β j ρ j r u T j δ A jr ρ j α (u r -u j ),
that we rewrite

m α j d t E α j = - r C j (p α j + p β j ) • u r + r u T r (A α jr + A β jr )(u r -u j ) + r C j p β j • u r - r u T r A β jr (u r -u j ) - ρ α j ρ β j ρ j r u T j         A α jr ρ α j - A β jr ρ β j         (u r -u j ).
Simple algebraic manipulations on the later term allow to write

m α j d t E α j = - r C j (p α j + p β j ) • u r + r u T r (A α jr + A β jr )(u r -u j ) + r C j p β j • u r - r (u r -u j ) T A β jr (u r -u j ) - ρ β j ρ j u T j r A α jr + A β jr (u r -u j ).
• According to Property 6 the term

1 V j        - r C j (p α j + p β j ) • u r + r u T r (A α jr + A β jr )(u r -u j )        , is weakly consistent with -∇ • (p α + p β )u x j . • Also since 1 V j ( r C j • u r ) is weakly consistent with ∇ • u, 1 V j        p β j r C j • u r        ≈ p β ∇ • u x j .
• Now, since r C jr = 0, one has

- r F jr = r (A α jr + A β jr )(u r -u j ),
so that Property 6 implies that

1 V j         - ρ β j ρ j u T j r A α jr + A β jr (u r -u j )         ≈ ρ β ρ ∇(p α + p β ) • u x j .
To conclude, it remains to prove for the remaining term

1 V j        - r (u r -u j ) T A β jr (u r -u j )        ≈ 0.
Let ζ α denote its limit:

1 V j        - r (u r -u j ) T A β jr (u r -u j )        -→ V j →0 ζ α .
We have shown

ρ α j d t E α j ≈ -∇ • (p α + p β )u + p β ∇ • u + ρ β ρ ∇(p α + p β ) • u x j + ζ α .
Since the same result holds for fluid β, simple calculations lead to

ρ α j d t E α j + ρ β j d t E β j = ρ j d t E j ≈ -∇ • (p α + p β )u x j + ζ α + ζ β .
According to Property 6

ρ j d t E j ≈ -∇ • (p α + p β )u x j , so that ζ α + ζ β ≈ 0. 199 Actually, one has 1 V j        r (u r -u j ) T A β jr (u r -u j )        + 1 V j        r (u r -u j ) T A α jr (u r -u j )        → 0,
since A α jr and A β jr are positive matrices, one has finally

1 V j        r (u r -u j ) T A α jr (u r -u j )        → 0 and 1 V j        r (u r -u j ) T A β jr (u r -u j )        → 0,
which ends the proof.

Discrete scheme

We now describe the fully discrete scheme. According to the previously established results, let ω ∈]0, 2]. One defines the following scheme for each fluid α ∈ { f 1 , f 2 }, β denoting the other one.

m α j τ α j n+1 -τ α j n ∆t = r C n jr • u α r n , (34)
m α j u α j n+1 -u α j n ∆t = - r F α,n jr -ω r νρ n r B n jr δu α j n+1 -(1 -ω) r νρ n r B n jr δu α r n , (35)
m α j E α j n+1 -E α j n ∆t = - r F α,n jr • u α r n - r νρ n r t u n r B n jr δu α r n + ω r νρ n r t u n+1 jr B n jr δu α r n -δu α j n+1 , (36)
where the fluxes are computed explicitly as

F α,n jr = C n jr p α j n -A α,n jr (u α r n -u α j n ) -νρ n r B n jr δu α r n , (37) and j
A α,n jr u α r n + j νρ n r B n jr δu α r n = j A α,n jr u α j n + j C n jr p α j n . (38)
To complete the scheme definition, observe that we introduced the following mean velocities

u n+1 jr = ρ α r n u α j n+1 +ρ β r n u β j n+1 ρ α r n +ρ β
Similarly to the semi-discrete case, for convinience, we inject the fluxes expression into momentum and total energy balance equation and use (39)

m α j u α j n+1 -u α j n ∆t = r A α,n jr (u α r n -u α j n ) + ων r ρ n r B n jr (δu α r n -δu α j n+1 ), (40)
m α j E α j n+1 -E α j n ∆t = - r C n jr p α j n • u α r n + r t u α r n A α,n jr (u α r n -u α j n ) + ν

Stability of the discrete scheme

In this section we establish that the scheme is stable for arbitrary equation of state: there

exists ∆t > 0 such that for each fluid α ∈ { f 1 , f 2 }, τ α j n+1 > 0, e α j n+1 > e(T = 0) and η α j n+1 ≥ η α j n .
For the sake of simplicity, and without loss of generality, we will consider in the following the case e α j n+1 > 0.

Actually, we will provide explicit timesteps for the positivity of density and internal energy, but we will only show that the increasing physical entropy timestep will be greater that the one of the mono-fluid case for given velocity fluxes, for which we established Property 3. The main reason is that there only exists existence results for entropy stability for cell-centered semi-Lagrangian schemes (even in 1D), see [START_REF] Després | Lagrangian systems of conservation laws[END_REF][START_REF] Gallice | Positive and Entropy Stable Godunov-type Schemes for Gas Dynamics and MHD Equations in Lagrangian or Eulerian Coordinates[END_REF].

Positivity of density

Since p = p(ρ, e) one has to ensure that density cannot be made negative.

Property 8 (Positivity of density).

Assuming that ∀α ∈ { f 1 , f 2 }, ∀ j ∈ M, ρ α j n > 0. Denoting
C αn the set of compressive cells for each fluid α,

C αn := j ∈ M/ r C n jr • u α r n < 0 , there exists ∆t ρ > 0 such that, ∀α ∈ { f 1 , f 2 }, ∀ j ∈ C αn , ∆t ρ < V α j n -r C n jr • u α r n .
Then, the scheme (34)-(38) defined by ∆t ∈ ]0, ∆t ρ ] ensures that

∀ω ∈ [0, 2], ∀α ∈ { f 1 , f 2 }, ∀ j ∈ M, ρ α j n+1 > 0.
Observe that, as expected, only compressive cells ( j ∈ C αn ) can lead to negative densities, so in the case of non-compressive flows, ∆t ρ may be arbitrarly large. Also, in the case of trianglular meshes, this constrain implies that no cell will tangle during the timestep.

Proof. Obviously, this is equivalent to show that τ α

j n+1 = 1 ρ α j n+1 > 0.
According to (34), one has

τ α j n+1 = τ α j n + ∆t m α j r C n jr • u α r n .
So, one has the following alternative:

• if j C αn that is r C n jr • u α r n ≥ 0, then ∀∆t > 0 one has τ α j n+1 > 0, • else if j ∈ C αn , one has r C n jr • u α r n < 0, then ∀∆t < τ α j n m α j -r C n jr •u α r n , one has τ α j n+1 > 0. Since m α j -r C n jr •u α
r n > 0, the existence of such a ∆t > 0 is obvious.

Positivity of internal energy

First, as a primary result, we give internal energy variation for fluid α ∈ { f 1 , f 2 }, β denoting the other one. Internal energy is updated as

e α j n+1 = e α j n + ∆t m α j        r t (u α j n -u α r n )A α jr n (u α j n -u α r n ) - r p α j n C n jr • u α r n        + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n+1 B n jr (δu α j n+1 -δu α r n )        - ∆t 2 2m α j 2        r A α jr n (u α j n -u α r n )        2 + ∆t 2 2m α j 2        ων r ρ n r B n jr (δu α j n+1 -δu α r n )        2 . (42)
Proof. Rewriting e α j n+1 = -1 2 u α j n+1 2 + E α j n+1 and using (41), one gets after a few arrangements

e α j n+1 = 1 2 u α j n 2 - 1 2 u α j n+1 2 + e α j n - ∆t m α j        r p α j n C n jr • u α r n + r t u α r n A α jr n (u α j n -u α r n )        + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n -ω r ρ β r n t δu α r n+1 B n jr δu α r n -δu α j n+1 +ω r ρ n r t u n+1 j B n jr δu α r n -δu α j n+1        . ( 43
)
As a first step one estimates kinetic energy variation

-∆K α j = 1 2 u α j n 2 - 1 2 u α j n+1 2 = u α j n + u α j n+1 2 • u α j n -u α j n+1 ,
which rewrites using (40)

-∆K α j =        u α j n - ∆t 2m α j        r A α jr n (u α j n -u α r n ) + ων r ρ n r B n jr (δu α j n+1 -δu α r n )               • ∆t m α j        r A α jr n (u α j n -u α r n ) + ων r ρ n r B n jr (δu α j n+1 -δu α r n )        , that is -∆K α j = ∆t m α j        r t u α j n A α jr n (u α j n -u α r n ) + ων r t u α j n ρ n r B n jr (δu α j n+1 -δu α r n )        - ∆t 2 2m α j 2        r A α jr n (u α j n -u α r n ) + ων r ρ n r B n jr (δu α j n+1 -δu α r n )        2 .
So, one has

e α j n+1 = e α j n + ∆t m α j        r t (u α j n -u α r n )A α jr n (u α j n -u α r n ) - r p α j n C n jr • u α r n - ∆t 2m α j        r A α jr n (u α j n -u α r n ) + ων r ρ n r B n jr (δu α j n+1 -δu α r n )        2         + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n+1 B n jr (δu α j n+1 -δu α r n )        + ων ∆t m α j r ρ n r t (u α j n -u α j n+1 )B n jr (δu α j n+1 -δu α r n ),
which using (38) is nothing but (42).

Actually, (42) can be rewritten as

e α j n+1 = e h α j n+1 + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n+1 B n jr (δu α j n+1 -δu α r n )        + ∆t 2 2m α j 2        ων r ρ n r B n jr (δu α j n+1 -δu α r n )        2 , (44)
where e h α j n+1 denotes the obtained internal energy without friction: i.e. injecting nodal velocities u α r n into the classic mono-fluid scheme. The remaining terms can be viewed as the heating do to friction.

Since ω ∈]0, 2], using Corollary 1 allows to minorate e α j n+1 e α j n+1 ≥ e h

α j n+1 + ν ∆t m α j        1 - ω 2 r ρ β r n t δu α r n B n jr δu α r n + ω 2 r ρ β r n t δu α j n+1 B n jr δu α j n+1        + ∆t 2 2m α j 2        ων r ρ n r B n jr (δu α j n+1 -δu α r n )        2 , (45)
which implies e α

Entropy stability for general equations of state

In the previous paragraph, we provided explicitly a choice of ∆t > 0 that ensures positivity of internal energy and density for the proposed scheme, but this is not sufficient for stability. In this section, we give an existence result of a strictly positive timestep ∆t that ensures production of physical entropy for arbitrary physical equation of state.

Let U = τ, u T , E T and let η be the entropy of the fluid. Gibbs formula reads T dη = de+ pdτ. Following [START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique lagrangienne bidimensionnelle[END_REF][START_REF] Després | Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques[END_REF], we estimate the entropy change, by means of a third-order Taylor expansion, due to the proposed scheme:

η(U α j n+1 ) -η(U α j n ) = (U α j n+1 -U α j n ) T ∂η ∂U U α j n + 1 2 (U α j n+1 -U α j n ) T ∂ 2 η ∂U 2 U α j n (U α j n+1 -U α j n ) + O (U α j n+1 -U α j n ) 3 . 20 One has ∂η ∂U U α j n = 1 T α j n (p α j n , -u α j n , 1) T and the variable change V = (p, -u, η) T reads (U α j n+1 -U α j n ) T ∂ 2 η ∂U 2 U α j n (U α j n+1 -U α j n ) = (V α j n+1 -V α j n ) T ∂ 2 η ∂V 2 V α j n (V α j n+1 -V α j n ) + O (U α j n+1 -U α j n ) 3 ,
where, see [START_REF] Després | Lagrangian systems of conservation laws[END_REF][START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique lagrangienne bidimensionnelle[END_REF] for instance,

∂ 2 η ∂V 2 V α j n = - 1 T α j n            (ρc) α j n -2 0 0 0 1 0 0 0 0            . Let O 1 := (U α j n+1 -U α j n ) T ∂η ∂U U α j n
, using (34), ( 40) and ( 41), one gets

O 1 = 1 T α j n ∆t m α j        p α j n r C n jr • u α r n -t u α j n        r A α,n jr (u α r n -u α j n ) + ων r ρ n r B n jr (δu α r n -δu α j n+1 )        - r C n jr p α j n • u α r n + r t u α r n A α,n jr (u α r n -u α j n ) + ν
O 1 = 1 T α j n ∆t m α j        r t (u α r n -u α j n )A α,n jr (u α r n -u α j n ) + ν r ρ β r n t δu α r n B n jr δu α r n -ων r ρ β r n t δu α r n+1 B n jr δu α r n -δu α j n+1        + 1 T α j n ∆t m α j        ων r ρ n r t (u α j n+1 -u α j n )B n jr δu α r n -δu α j n+1       
.

Now using Lemma 1, one gets

O 1 = 1 T α j n ∆t m α j        1 - 1 2 ω ν r ρ β r t δu α r n B n jr δu α r n + 1 2 ων r ρ β r t δu α j n+1 B n jr δu α j n+1        + 1 T α j n ∆t m α j        r t (u α r n -u α j n )A α,n jr (u α r n -u α j n ) + ων r ρ β r n t δu α r n -δu α j n+1 B n jr δu α r n -δu α j n+1        + 1 T α j n ∆t m α j        ων r ρ n r t (u α j n+1 -u α j n )B n jr δu α r n -δu α j n+1       
. production of the continuous in time scheme established in Property 5 page 11.

241

One now focuses on the second-order term of the entropy variation

O 2 := 1 2 (V α j n+1 -V α j n ) T ∂ 2 η ∂V 2 V α j n (V α j n+1 -V α j n ),
which rewrites

O 2 = 1 2 (∆Ψ) T        (ρc) α j n -2 0 0 1        ∆Ψ, with ∆Ψ = p α j n+1 -p α j n -u α j n+1 + u α j n .
One has to estimate p α j n+1p α j n . Assuming that the equation of state p : (τ, e) → p(τ, e) is regular enough, one has

p α j n+1 -p α j n = (τ α j n+1 -τ α j n ) ∂p ∂τ jn + (e α j n+1 -e α j n ) ∂p ∂e jn + O(∆t 2 ).
Using ( 34) and ( 42) and keeping only first-order terms, one has

p α j n+1 -p α j n = ∆t m α j        r C n jr • u α r n        ∂p ∂τ jn + ∆t m α j               r t (u α j n -u α r n )A α jr n (u α j n -u α r n ) - r p α j n C n jr • u α r n        + ν        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n+1 B n jr (δu α j n+1 -δu α r n )        +ων        r ρ n r t (u α j n -u α j n+1 )B n jr (δu α j n+1 -δu α r n )               ∂p ∂e jn + O(∆t 2 ).
Then, using (40), one gets

O 2 = - 1 T α j n ∆t 2 2m α j 2        r C n jr • u α r n        ∂p ∂τ jn +        r t (u α j n -u α r n )A α jr n (u α j n -u α r n ) - r p α j n C n jr • u α r n + ν r ρ β r n t δu α r n B n jr δu α r n + ων r ρ β r n t δu α j n+1 B n jr (δu α j n+1 -δu α r n ) +ων r ρ n r t (u α j n -u α j n+1 )B n jr (δu α j n+1 -δu α r n )        ∂p ∂e jn 2 (ρc) α j n -2 +        r A α,n jr (u α r n -u α j n ) + ων r ρ n r B n jr (δu α r n -δu α j n+1 )        2 + O(∆t 3 ).
Finally, putting all the pieces together, one has

η(U α j n+1 ) -η(U α j n ) = 1 T α j n ∆t m α j        1 - 1 2 ω ν r ρ β r t δu α r n B n jr δu α r n + 1 2 ων r ρ β r t δu α j n+1 B n jr δu α j n+1        + 1 T α j n ∆t m α j        a - ∆t m α j (b + c) + O(∆t 2 )        ,
with a ≥ 0 and b ≥ 0.

Thus it remains to study the positiveness of a-∆t m α j (b+c)+O(∆t 2 ). There are two possibilities.

Case a > 0. In that case, there obviously exists ∆t > 0 such that

T α j n m α j η(U α j n+1 ) -η(U α j n ) ∆t ≥ 1 - 1 2 ω ν r ρ β r t δu α r n B n jr δu α r n + 1 2 ων r ρ β r t δu α j n+1 B n jr δu α j n+1 . Case a = 0. If a = 0, one has r t (u α r n -u α j n )A α,n jr (u α r n -u α j n ) + ων r ρ β r n t δu α r n -δu α j n+1 B n jr δu α r n -δu α j n+1 = 0.
Since ω ≥ 0 and since A α,n jr and B n jr are positive matrices, all the terms in the sum are zeros. Let us first focus on t (u α r nu α j n )A α,n jr (u α r nu α j n ) = 0 terms. Two cases occur. In case of Eucclhyd scheme, A α,n jr is positive definite so that one has u α r n = u α j n . For Glace scheme

t (u α r n -u α j n )A α,n jr (u α r n -u α j n ) = (ρc) α j n C n jr C n jr • (u α r n -u α j n ) 2 = 0.
So, for both scheme, one has Finally, if a = 0, one has

C n jr • u α r n = C n jr • u α j n and A α,n jr (u α r n -u α j n ) = 0. Recalling that r C n jr = 0,
T α j n m α j η(U α j n+1 ) -η(U α j n ) ∆t = ν r ρ β r t δu α r n B n jr δu α r n - ∆t 2m α j        ν r ρ β r n t δu α r n B n jr δu α r n        2 ∂p ∂e 2 jn (ρc) α j n -2 + O(∆t 2 ).
Before enunciating the result, one should remark that in the general case, one has

η(U α j n+1 ) -η(U α j n ) = 1 T α j n ∆t m α j        (a + a ν ) - ∆t m α j (b + c) + O(∆t 2 )        ,
with a ≥ 0, a ν ≥ 0 and b ≥ 0. Again, one has two alternatives a + a ν > 0 or a + a ν = 0. In the first case, there exists ∆t such that η(U α j n+1 ) -η(U α j n ) > 0. In the second case, one has a = a ν = 0 so as previously,

a = 0 =⇒ C n jr • u α r n = C n jr • u α j n , A α,n jr (u α r n -u α j n ) = 0 and δu α r n -δu α j n+1 = 0.
Also, since a ν = 0 and since B n jr is positive definite one has δu α r n = 0, so the scheme (34)-(

gives U α j n+1 = U α j n and finally one has ∀∆t > 0, η(U α

j n+1 ) = η(U α j n ).
The obtained results are summarized in the following Property.

Property 10 (Entropy). Let U := τ, u T , E T and let η the entropy. There exists ∆t η > 0 , such that ∀α, β ∈ { f 1 , f 2 }, such that α β, if the pressure law p α : (ρ, e) → p α (ρ, e) is a differentiable function, then the scheme (34)-(38) defined by ∆t = ∆t η and ∀ω ∈ [0, 2], ensures that, 1. the scheme is entropy stable:

∀ j ∈ M, η U α j n+1 ≥ η U α j n ,
2. and ∀ j ∈ M, one has the following alternative. In this case, the scheme acts simply as a first-order ODE solver. Since then dη = de and since η is strictly convex, a time integration error is to be expected.

If ∀r ∈ R j , C n jr • u α r n = C n jr • u α j n and δu α r n -δu α j n+1 = 0, then T α j n m α j η(U α j n+1 ) -η(U α j n ) ∆t ≥ ν r ρ β r t δu α r n B n jr δu α r n + O(∆t), else T α j n m α j η(U α j n+1 ) -η(U α j n ) ∆t ≥ ν r ρ β r t δu α
To sum up, we proved that the proposed scheme is stable, meaning that there exists 0 < ∆t ≤ min(∆t ρ , ∆t e , ∆t η ) such that the scheme is entropy stable and preserves positivity of density and internal energy. Moreover, it is consistent with (2).

Minoration of ∆t η

As stated before, to prove that the scheme is asymptotic preserving, it remains to show that lim ν→+∞ ∆t η 0. Even if we will not provide here an explicit value, we will give a lower bound independent of ν. Choosing ν 1 and solving (38) implies that δu α r is arbitrary small, and we write u r = u α r = u β r , and the sum equations (38) can be rewritten has j

(A α jr + A β jr )u r = j C jr p β j , (46)
since u α = u β = 0 and p α is constant (recalling that j C jr = 0). Since ρ α c α = , A α jr is neglectable with regard to A β jr . So, one has j A β jr u r = j C jr p β j , that is the timestep for fluid α is the same as one imposed by the fluid β which is much smaller than the arbitrary one ontained in the mono-fluid case ν = 0.

Example 2 (∆t η0 < ∆t ην ). Let us now consider a similar example. The two fluids in ]0, 1[ are now described as follows. Let fluid α be in the initial state of a Sod shock tube. Fluid β is this time a very heavy fluid at rest: ρ β = 1 with 0 < 1, u β = 0 and e β being set such that the sound speed c β = 1.

If ν = 0, stability for the fluid α is the one of the mono-fluid Sod shock tube.

Choosing ν 1, (46) holds again. Since, ρ β = 1 is very large, one gets in the limit u r = 0, which provides unconditionally stability for both fluids. Using this scheme, one easily checks that internal energy variation reads

e α j n+1 = e α j n + ∆t m α j        r t (u α j n -u α r n )A α jr n (u α j n -u α r n ) - r p α j n C n jr • u α r n        + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n B n jr (δu α j n -δu α r n )        - ∆t 2 2m α j 2        r A α jr n (u α j n -u α r n ) + ων r ρ n r B n jr (δu α j n -δu α r n )        2 .
that is

e α j n+1 = e h α j n+1 + ν ∆t m α j        r ρ β r n t δu α r n B n jr δu α r n + ω r ρ β r n t δu α j n B n jr (δu α j n -δu α r n )        - ∆t 2 m α j 2        ων r ρ n r B n jr (δu α j n -δu α r n )        •        r A α jr n (u α j n -u α r n )        - ∆t 2 2m α j 2        ων r ρ n r B n jr (δu α j n -δu α r n )        2 ,
where e h

ALE scheme

The semi-Lagrangian scheme presented in this paper is defined assuming that both fluid meshes are identical at the begining of the timestep. One understands easily that this is of huge help in the construction of an asymptotic preserving scheme. One could imagine a purely Lagrangian approach, but even dealing with a non-AP approach seems very difficult since one would have to consider meshes intersections and complex geometrical calculations.

Thus, the algorithm, we propose in this paper, consists in ensuring that for each timestep both fluids meshes coincide. To do so an ALE formulation is mandatory.

Figure 2 depicts the general ALE case. Our ALE method is a Lagrange-rezoning-advection procedure which ensures that the solution is defined at time t n+1 on a unique mesh. • Finally, for both fluids, the numerical solution is computed on the common mesh by remapping the conservative variables (ρ α , ρ α u α , ρ α E α ) T at velocity -v α,n+1 r , with a second-order accurate scheme. One can then compute another timestep.

Numerical tests

Reference scheme

We introduce the following scheme that will be used as a "non-AP" scheme to illustrate the advantages of scheme ( 22)- [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF]. For each fluid α ∈ { f 1 , f 2 }, one writes

d t m α j = 0, m α j d t τ α j = r C jr • u α r , m α j d t u α j = - r F α jr - r νρ r B jr δu α j , m α j d t e α j = - r F α jr • u α r - r νρ r u T jr B jr δu α j , (47)
where u jr and ρ r are defined as in section 4.1 and the fluxes are given by

F α jr = C jr p α j -A α jr (u α r -u α j ) j A α jr u α r = j A α jr u α j + j C jr p α j . (48)
It can be showed that this scheme is entropic, conservative in volume and mass for each fluid and in the sum of momentums and total energies. Also, the scheme is weakly consistent with (1).

One can moreover show that its associated discrete in time scheme, where only u j terms are implicit, is stable in the same way as scheme (34)-(38).

However, this scheme does not a priori preserve the asymptotic. For these reasons this scheme is a very good candidate for the comparisons that we perform in this section.

Tests conditions

In all the following tests, we choose A α jr = ρ α j c α j i N i jr ⊗N i jr N i jr (Eucclhyd scheme) and B jr = V jr I 2 , with V jr = 1 #R j V j . Also, for each test one chooses γ α = γ β = 1.4.

Results are compared with the non-AP scheme (47)-(48). Also for the 2D tests, we compare our results (ν >> 1) to the mono-fluid case, where mass fraction ρ α ρ α +ρ β is treated as a passive scalar.

As it is often the case for multi-velocity models [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF], the scheme is only defined in regions where both fluids are present. Thus in regions where a fluid should be absent, one keeps a neglectable amount of it.

Sod shock tube

The computational domain we consider is Ω :=]0, 1[×]0, 0.1[. Initial data is given as U := (ρ, u, p) T , so that one defines U L := (1, 0, 1) T , U R := (0.125, 0, 0.1) T and U L := (ε, 0, ε) T . For both fluids initial states are then The same test is performed for a friction parameter ν = 10 6 . The density sum is presented in figure 4 at time t = 0.14.

U α = 1 ]0,0.5[ (U L -U ) + 1 ]0.5,1[ U and U β = 1 ]0,0.5[ U + 1 ]0.5,1[ (U R -U ),
One retrieves the results presented in [START_REF] Enaux | Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles[END_REF], even if the scheme does not degenerate in 1D to the scheme proposed in [START_REF] Enaux | Analyse mathématique et numérique d'un modèle multifluide multivitesse pour l'interpénétration de fluides miscibles[END_REF].

Finally we show the importance of choosing ω ∈]0, 2] by comparing the velocity differences in the cases ω = 0 and ω = 1 for ν = 10 6 using the scheme (34)-(38). One will note the oscillations obtained for the choice ω = 0, for which we could not show the asymptotic preserving property.

Triple-point problem

The triple-point problem is a standard benchmark [START_REF] Loubère | Validation test case suite for compressible hydrodynamics computation[END_REF]. It is a multi-dimensional Riemann problem whose data are close to the Sod shock tube. The self-similarity of the problem yields an notices the nice agreement of the solution for the proposed scheme with regard to the mono-fluid case, even for this small amount of cells, whereas the non-AP scheme is not even able to compute the large structures of the flow at this grid resolution.

Then we studdy the effect of the friction parameter. Figure 8 presents the obtained solutions, on a finer 210 × 90 grid, for ν ∈ {10, 100, 10 6 }.

A Rayleigh Taylor instability

For this test, we modify the scheme in order to incorporate the gravity treatment. Obviously, we use a well-balanced approach [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF] to take this term into account. For this modified scheme, the properties we established for (34)-(39) remain true. We did not take the gravity term into account in section 4 to avoid a more complex presentation, since there is no additional difficulties to overcome.

The interface perturbation is defined by the function f (y) = 0.05 cos(8πy) and centered at Initially, velocities are set to 0 in Ω, and densities are defined as ρ α = 1 Ω α (0.8 -ε) + 1 Ω β ε, and ρ β = 1 Ω α ε + 1 Ω β (0.25 -ε).

Choosing the gravity acceleration as g = 9.8 e x , we define the pressure in the whole domain at a quasi-equilibrium state (omitting the y dependancy), that is

p(x) = x 0 ρ α + ρ β g • e x .
Again, symmetry boundary conditions are imposed all over ∂Ω. We represent the mass fraction of fluid α that is ρ α ρ α +ρ β . One takes ω = 1. We use the same ALE strategy as in the previous test:

a barycentric remapping is performed on the mesh of fluid α and we set M n+1 At first, we validate the approach by comparing the obtained result to the mono-fluid scheme.

The results are presented on figure 6.5, one observes again a very good agreement even on a 112 × 40 coarse grid. As expected, the non-AP scheme clearly shows lack of convergence. Finally, we studdy the influence of the friction parameter ν for successive values of 100, 1000 and 10 6 . A slightly finer grid (224 × 80) is used for it.

Conclusion

In this paper, we presented a multi-dimensionnal asymptotic preserving scheme to solve a bi-fluid model defined as a set of two Euler systems coupled with a friction term. The originality of the approach is that the scheme is ALE: the only constrain being that meshes must coincide at the begining of each timestep.

The scheme is conservative and weakly-consistent by construction. Moreover, we showed that it is at least as stable as the underlying hydro-scheme in the sense that the timestep required to increase entropy does not tend to zero when friction increases. We showed consistency of the limit scheme (ν → +∞) to the limit model. So, we proved that the scheme is asymptotic preserving. On the way we proved some stability results with regard to the fluxes u r , which give some bounds independently of ν (Property 3), which complete the numerical analysis of the scheme.

The numerical results show that the scheme behaves as expected and appears to be a good candidate to studdy interpenetration mixing [START_REF] Scannapieco | A multifluid interpenetration mix model[END_REF], which is the goal of this studdy. Actually, all the results1 can be established with a varying positive friction ν. In the paper we kept ν constant for the sake of simplicity. The numerical analysis and tests are performed in 2D, however the analysis in 3D is completely unchanged.

On the numerical point of view, a second-order accurate version of the scheme would be of interest. However, this is not an easy task for two main reasons. First, on the theoretical point of view, establishing properly the asymptotic preserving property would be challenging. Second, using a Runge-Kutta-like approach to get second-order accuracy in time would probably impose to incorporate the remeshing into the time integration or to consider a one-step approach.

Another extension is to introduce more Physics in the model. The friction coupling is a very simple approach, one could use more appropriate closures based on the presented work, for instance [START_REF] Scannapieco | A multifluid interpenetration mix model[END_REF] or [START_REF] Brull | Asymptotic-preserving scheme for a bi-fluid euler-lorentz model[END_REF] where Lorentz forces are taken into account in a ion-electron mixture.

Figure 1 :

 1 Figure 1: Illustration of C jr and N i jr vectors at vertex r for a polygonal cell j.

Property 9 (

 9 Positivity of internal energy). Assuming that ∀α∈ { f 1 , f 2 }, ∀ j ∈ M, e αj n > 0, there exists ∆t e > 0 such that the scheme (34)-(38) ensures that ∀ω ∈ [0, 2], ∀∆t ∈]0, ∆t e [, ∀α ∈ { f 1 , f 2 }, ∀ j ∈ M, e α j n+1 > 0. Proof. The proof is obvious since e α j n+1 ≥ e h α j n+1 and since e h α j n+1 (∆t) is a polynomial of degree 2 satisfying e h α j n+1 (0) = e α j n > 0. ∆t e is nothing but the smallest root of these polynomial for each cells of each fluid.

one also has r p α

  0. Here again two cases occur ω = 0 or ω > 0. In that second case, since B n jr are positive definite, this implies δu α r n -δu α j n+1 = 0.

Remark 3 .

 3 Let us comment point 2 of property 10. Actually, this is a consistency result with regard to (2). In the first case (if ∀r ∈ R j , C n jr •

Property 11 ., 1 T p 1 ,Remark 4 .

 11114 Let ω ∈]0, 2]. ∀ j ∈ M, let τ n j , state of fluid α ∈ { f 1 , f 2 }. Let {u r } r∈R j , be an arbitrary set of nodal velocities (or velocity fluxes). Then, if ∀ν ≥ 0, τ ν,n+1 j , e ν,n+1 j denotes the thermodynamic state obtained by scheme (34)-(36), one has where η := η(τ, e) is the physical entropy expressed according to the independent variables τ and e.Proof. Gibbs formula reads ∇ τ,e η = where T := T (τ, e) is a positive function. So, for any τ, η(τ, •) is an increasing function.Since (34) is independent of ν and according to (45), one has ∀{u r } r∈R j , ∀ν ≥ 0, ∀∆t τ Property 11 establishes that, for a given set of velocities {u r } r∈R j , the maximum timestep required for the scheme to be entropy stable is greater than the mono-fluid timestep independently of ν.However, one emphisises that velocities {u r } r∈R j are actually functions of ν as expressed in (38). It turns out that entropy stability timestep depends on ν though the velocity fluxes and can be either bigger or smaller than the mono-fluid timestep. Example 1 (∆t η0 > ∆t ην ). Let us consider two fluides in the monodimensional ]0, 1[ domain. Let the first fluid α be a very light fluid at rest ρ α = with 0 < 1, u α = 0 and e α being set such that the sound speed c α = 1. Let the second fluid β be the initial state of a Sod shock tube. If ν = 0, one has obviously ∀r, u α r = 0. So, (34)-(36) implies U α

4. 3 . 5 .

 35 On the importance of the implicit velocities in (34)-(36)Using the notations defined in section 4.2, let us consider the fully explicit scheme that consists in replacing momentum and total energy updates in (34)-(38) by their explicit counterparts

  the obtained internal energy without friction. The later term being a negative factor of ν 2 , in the explicit case, ∀∆t > 0 for large values of ν, one can have e α even if a similar result to Property 10 can be established (existence of an entropy stable timestep), one cannot prove an equivalent of Property 11. If cell velocities are explicit, one eventually gets lim ν→+∞ ∆t e = lim ν→+∞ ∆t η = 0 for a given set of nodal velocities {u r } r∈R j .

•••Figure 2 :

 2 Figure2: Left: at time t = t n , both fluid share the same mesh. Middle: at the end of the Lagrangian phase, one gets two different meshes, one for each fluid. Right: meshes are displaced so that they coincide. Solution is remapped and a new timestep can be performed.

where 1 O

 1 denotes the characteritic function of the set O and where we take = 10 -3 . Choosing ω = 1, one figure 3, we compare the solution at time t = 0.14 obtained by the proposed scheme (34)-(38) to the reference scheme (47)-(48). One plots the density sum: ρ α +ρ β .The grid is 200x3 cells and the solution is compared to a reference solution obtained using a 10 5 × 3 grid. The simulation is Eulerian: the smoothing strategy consists in a return to the initial mesh.

Figure 3 :

 3 Figure 3: ν = 100. ρ α + ρ β profile. AP-scheme (left) gives a much better solution than the non-AP scheme (right).

Figure 4 :

 4 Figure 4: ν = 10 6 . ρ α + ρ β profile. AP-scheme (left) gives a much better solution than the non-AP scheme (right). The expected solution is close to the classical mono-fluid case.

Figure 5 :ρ

 5 Figure 5: ν = 10 6 . Comparison of the δu j obtained for ω = 0 (red) and ω = 1 (blue) at time t = 1.4 for a 200 × 3 grid.

Figure 6 :

 6 Figure 6: Geometry, pressures and densities for the triple-point problem at time t = 0.

Figure 7 :

 7 Figure 7: 91 × 40 mesh. Mass fraction of fluid α at time t = 5. Left: mono-fluid solution. Right: bi-fluid solution with ν = 10 6 . Bottom: bi-fluid solution with non-AP scheme with ν = 10 6 .

x

  = 0.35 in the computational domain Ω =]0, 0.7[×]0, 0.25[. Thus, two regions are defined:Ω α = {(x, y) ∈ Ω/x < 0.35 + f (y)} and Ω β = Ω \ Ω α .

Figure 8 :

 8 Figure 8: 210 × 90 mesh. Time t = 5. Mass fraction of fluid α. Effect of the friction parameter ν. Left: ν = 10. Right: ν = 100. Bottom: ν = 10 6 .

Figure 9 :

 9 Figure 9: Rayleigh-Taylor test initial geometry. Fluid α being heavier than fluid β, instability will grow.

β = M n+1 α to allow the calculation of timestep n + 1 .

 1

Figure 10 :

 10 Figure 10: 112 × 40 mesh. Mass fraction of fluid α. Time t = 0.7s. Top: mono-fluid solution. Middle: bi-fluid solution with ν = 10 6 . Bottom: bi-fluid solution with non-AP scheme with ν = 10 6 .

Figure 11 :

 11 Figure 11: 224 × 80 mesh. Time t = 0.7. Mass fraction of fluid α. Influence of the friction parameter. Top: ν = 100. Middle: ν = 1000. Bottom: ν = 10 6 .

If friction parameter depends on the cell data (ν = ν j ), Property 3 takes a slightly different form.

Acknowledgement

Authors are grateful to C. Buet, H. Jourdren and F. Lagoutière for their valuable comments and remarks about this work.

Symmetry boundary conditions are set at each straight boundary of the computational domain.

The ALE strategy we use for this test consists in a barycentric smoother for the grid of the fluid α and then to impose x β r = x α r .

We run the test on a 91 × 40 grid and setting ω = 1. Choosing the friction parameter ν = 10 6 , we compare the obtained result to the solution of the mono-fluid solver and to the non-AP scheme, see figure 7. For the comparison, we plot the mass fraction in each case: ρ α ρ α +ρ β . One