Laura Gastaldo
email: [laura.gastaldo@irsn.fr

Rapha Èle Herbin

Jean-Claude Latch
email: jean-claude.latche]@irsn.fr

A discretization of phases mass balance in fractional step algorithms for the drift-flux model

Keywords: Two-phase flows, drift-flux model, finite volume methods, monotone schemes

 that the discretization of the convection operator must be such that it vanishes for constant y, we prove the existence and uniqueness of the solution, together with the fact that it remains within its physical bounds, i.e. within the interval [0, 1]. Then this scheme is combined with a pressure correction method to obtain a semi-implicit fractional-step scheme for the so-called drift-flux model. To satisfy the above-mentioned assumption, a specific time-stepping algorithm with particular approximations for the density terms is developed. Numerical tests are performed to assess the convergence and stability properties of this scheme.

Introduction

This paper adresses a class of physical problems which can be set under the form of the Navier-Stokes equations, supplemented by the balance equation of an independent unknown field y:

∂ t ρ + ∇ • (ρu) = 0, (1.1a) ∂ t (ρ u) + ∇ • (ρu ⊗ u) + ∇p -∇ • τ = f , (1.1b) ∂ t (∂ ρ y) + ∇ • (ρ y u) + ∇ • (ρ φ (y) u r ) = ∇ • (D ∇y),
(1.1c) ρ = η(p, y), (1.1d) where t stands for the time, u for the fluid velocity, p for the pressure, ρ for the fluid density. The tensor τ is the viscous part of the stress tensor, given by the following expression:

τ = µ (∇u + ∇ t u) - 2 3 µ (∇ • u) I, ( 1.2)
IMA Journal of Numerical Analysis c Institute of Mathematics and its Applications 2005; all rights reserved.

where µ is the fluid viscosity and I stands for the identity tensor. For a constant viscosity, this relation yields:

∇ • τ = µ ∆ u + 1 3 ∇∇ • u . (1.3)
and, in this case, this term is coercive. The function η, which gives the density as an explicit function of y and the pressure, is obtained from the equation of state of the considered fluid. The nonlinear function φ is such that φ ∈ C 1 ([0, 1], R) and φ (0) = φ (1) = 0; for physical reasons, y is supposed to satisfy 0 y 1, so that φ (•) can be extended by continuity to R \ [0, 1] by 0 without altering the model. The volumic diffusion coefficient D and the velocity field u r are known quantities. The problem is posed over an open polygonal bounded connected subset Ω of R d , d 3, and over a finite time interval (0, T ).

It must be supplemented by suitable boundary conditions, and initial conditions for ρ, u and y. In the sequel, we assume that f ∈ L 2 (Ω ).

Several physical problems enter this framework. For instance, taking for y the gas mass fraction, and for ρ the mixture density of dispersed two-phase flows yields the so-called drift-flux model, in the isothermal case. In this case, u r is the relative velocity between the two phases and φ is given by φ (y) = y(1y). Dispersed two-phase flows and, in particular, bubbly flows are widely encountered in industrial applications; one may think, in particular, about bubble columns and airlift reactors, where the agitation due the gaseous phase is used to promote the contact and consequently the chemical reactions between chemical species in the flow. They are also of wide concern in the framework of nuclear safety studies, either for the modelling of boiling of water in the primary coolant circuit in case of an accidental depressurization or for the simulation of the late phases of a core-melt accident, when the flow of molten core and vessel structures comes to chemically interact with the concrete of the containment floor. This is the context of the present work.

When designing a numerical scheme for the solution of system (1.1), one faces at least two difficulties. First, the unknown y is expected, from both physical and mathematical reasons, to remain in the interval [0, 1]. Similarly, the unknown ρ should remain positive at all times. This suggests to build a numerical scheme which reproduces these features at the discrete level. This is performed by discretizing the mass balance equation (1.1a) with an upwind finite volume scheme and combining a monotone flux approach [13, section 21] for the term ∇ • (ρ φ (y) u r ) in the advection diffusion equation (1.1c) with the argument introduced by Larrouturou [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF]: the discrete counterpart of the advection operator ∂ ρ y/∂t + ∇ • (ρ y u) satisfies a maximum principle provided that this operator applied to a constant value of y vanishes, i.e. that a discrete version of the mass balance is satisfied. We first prove that, with the proposed scheme, the variable y is kept within its expected range [0, 1]. By a topological degree argument, this yields the existence of a discrete solution, which is then shown to be unique by a duality argument. Now even if the model at hand represents a compressible flow, in fact the liquid is almost incompressible, so that zones may appear in the flow where the velocity of acoustic waves is very large, and the Mach number accordingly very small. We thus need to design a numerical method which is stable in the low Mach number limit, and therefore able to deal with incompressible flows. To this purpose, we use a fractional step algorithm of the class of finite element projection methods, which are widely used for incompressible flows, see e.g. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Marion | Navier-Stokes equations: Theory and approximation[END_REF] and references herein. An extension to the barotropic Navier-Stokes equations close to the scheme developed here can be found in [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF], together with references to related works. For stability reasons, the spatial discretization must preferably be based on pairs of velocity and pressure approximation spaces satisfying the so-called inf-sup or Babuska-Brezzi condition (e.g. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]). Among these elements, nonconforming approximations with degrees of freedom for the ve-locity located at the center of the faces seem to be well suited to a coupling with a finite volume method for the advection-diffusion of y; this is the choice made here. The fractional step approach is extended to the entire scheme, and the whole set of equations is thus solved in sequence. As a consequence, to ensure both conservativity and the above-mentioned monotonicity condition for the computation of y, a particular time stepping must be developed. This paper is organized as follows. The finite volume scheme for transport equations (1.1a) and (1.1c) are described and analysed in section 2. The fractional step algorithm for the solution of the whole problem is presented in section 3, along with the non conforming finite element discretization of Equation (1.1b) . We show how the compatibility between the discretization of the different steps of the algorithm allows to prove the well-posedness and the stability properties of the fully discrete algorithm. Numerical tests are reported in section 4; first a problem exhibiting an analytical solution allows to assess convergence properties of the discretization, then a physical situation consisting of a phase separation problem under gravity is addressed.

A finite volume scheme for the nonlinear advection-diffusion equation

In this section, we present a finite volume discretization of the advection diffusion (1.1c). Precisely speaking, the problem that we address here is the following: supposing that the velocity field u and the density ρ are known and satisfy a discrete mass balance equation of a particular form, we build a scheme for the computation of y which enjoys the property that the unknown y stays in the interval [0, 1]. We thus obtain a "brick" of a fractional step algorithm, which may be implemented for the solution of a wide range of problems, and presents the high practical interest to keep y within its physical bounds. This section is built as follows. We present in Paragraph 2.1 the considered finite volume mesh and discretization space. In Paragraph 2.2, we give the discretization of the mass balance equation (1.1a) which is used in Section 3 for the solution of the complete problem (1.1). It also serves as an example of a general form of the mass balance equation for ρ and u which is required for the stability analysis of the nonlinear advection diffusion equation discretization (1.1c); in particular, we show that it preserves the positivity of the density. We then address in Paragraph 2.3 the stability analysis of the finite volume discretization of the advection diffusion (1.1c) provided a specified general discrete mass balance is satisfied.

Discretization mesh and spaces

In this section, we denote by T a set of non intersecting convex subdomains of Ω , such that:

(i) Ω = K∈T K.
(i.e. E \ E ext ) is denoted by E int . If K, L ∈ T , we suppose that either K ∩ L = / 0, K ∩ L is a vertex or K ∩ L ∈ E int ,
and, in the latter case, this common edge or face of K and L is denoted by K|L.

For the discretization of the advection diffusion (1.1c), the following additional orthogonality condition is assumed to hold:

(iii) There exists a family P = (x K ) K∈T of points of Ω such that x K ∈ K for all K ∈ T and, if σ = K|L, x K = x L and the straight line going through x K and x L is orthogonal to σ .

Even though this condition is not needed in the present analysis because we are only concerned here with stability issues, it enables to perform a consistent approximation of the normal diffusion fluxes to the boundaries of the cells, which in turn is known, in simpler cases, to yield the convergence of the scheme (see e.g. [START_REF] Eymard | Finite volume methods[END_REF]Chapter 2] for the convergence proof for a linear steady diffusion problem).

The set of edges (in 2D) or faces (in 3D) of a cell K of T is denoted by E (K). For each internal edge or face of the mesh σ = K|L, n KL stands for the normal vector to σ , oriented from K to L (so n KL = -n LK ). By |K| and |σ |, we denote the measure of the control volume K and of the edge or face σ , respectively. For any K ∈ T and σ ∈ E (K), we denote by d K,σ the Euclidean distance between x K and σ . For any σ ∈ E , we define

d σ = d K,σ + d L,σ , if σ ∈ E int (in which case d σ is the Euclidean distance between x K and x L ) and d σ = d K,σ if σ ∈ E ext .
We denote by X T the space of piecewise constant functions on each control volume K ∈ T :

X T = q ∈ L 2 (Ω ) : q| K = constant, ∀K ∈ T . (2.1)
Let N = card(T ); any function q ∈ X T can be defined by the data of the N values of q over the elements of T , and hereafter we sowewhat improperly identify the function itself and this family of real numbers, therefore allowing such expressions as q ∈ X T , q = (q K ) K∈T . The gas mass fraction y is approximated by functions of X T , so y = (y K ) K∈T .

Finally, we define M = card(E int ) and, throughout this paper, for any real number a, we define a + = max(a, 0) and a -=min(a, 0), so that a = a +a -with a + 0 and a -0.

Space discretization of the mass balance equation (1.1a)

We consider here a first order backward Euler time discretization of the mass balance equation (1.1a), which reads in the semi-discrete form:

ρ -ρ * δt + ∇ • (ρu) = 0. (2.2)
Assuming a known velocity field u and initial density ρ * , we discretize the density field ρ by a finite volume method on the above described mesh T . We assume that from the discrete velocity field u, we are able to derive a family of values representing the velocity on each internal edge or face (u σ ) σ ∈E int ∈ R M×d (for the sake of simplicity, we suppose that the velocity obeys homogeneous Dirichlet boundary conditions, and thus that its value over the external edges vanishes); note that these quantities are readily obtained for instance when using a Crouzeix Raviart or Rannacher-Turek discretization as in Section 3.2. Hence, for the known families ρ * ∈ X T and (u σ ) σ ∈E int ∈ R M×d , we look for ρ ∈ X T , solution of the following upwind finite volume scheme:

∀K ∈ T , |K| δt [ρ K -ρ * K ] + ∑ σ =K|L F up σ ,K (ρ, u) = 0, (2.3)
where F up σ ,K (ρ, u) is defined as the upwind mass flux with respect to u through the interface σ and is defined by:

∀u = (u σ ) σ ∈E int ∈ (R M ) d , ∀ρ = (ρ K ) K∈T ∈ R N , F up σ ,K (ρ, u) = (|σ | u σ • n KL ) + ρ K -(|σ | u σ • n KL ) -ρ L , ∀K ∈ T , ∀σ = K|L. (2.4)
Note that (2.4) may also be written

F up σ ,K (ρ, u) = |σ | u σ • n KL ρ up σ , with ρ up σ =    ρ K if u σ • n KL 0, ρ L otherwise.
(2.5)

The flux thus defined satisfies the so-called local conservativity, which is one of the finite volume method's most classical features, and which we can express here as:

F up σ ,K (ρ, u) + F up σ ,L (ρ, u) = 0. (2.6)
Moreover, this scheme ensures that condition (2.8a) holds, i.e. that the density stays positive at all times, as we now show:

LEMMA 2.1 Let (ρ * K ) K∈T ∈ R N and (u σ ) σ ∈E int ∈ R M×d be two given families such that ρ * K > 0 for any K ∈ T . Then the linear system (2.3) has a unique solution which satisfies ρ K > 0 for any K ∈ T .

Proof. With a natural equation ordering, the matrix of the linear system (2.3) is of the form I + A where I denotes the (N, N) identity matrix and A = (a i, j ) 1 i, j N . It is easy to check that thanks to the upwind choice (2.4), one as: a i,i 0 for 1 i N, a i, j 0 for 1 i, j N, i = j, and ∑ N i=1 a i, j = 0 for 1 j N (note that the sum is over the lines and not the columns). Hence I + A is an M-matrix. This yields that the matrix I + A is invertible and that if ρ * K > 0 for all K ∈ T then ρ K > 0.

Discretization of the nonlinear advection-diffusion equation (1.1c)

We now turn to the discretization of the balance equation (1.1c); in a semi-discrete form obtained by a first order backward Euler time discretization, it reads:

ρy -ρ * y * δt + ∇ • (ρyu) + ∇ • (ρ φ (y) u * r ) = ∇ • (D∇y), in Ω × (0, T ), (2.7)
where the density fields ρ and ρ * , the beginning-of-step mass fraction y * and the velocity fields u and u * r are here supposed to be known quantities. For the sake of simplicity, we assume in this section that both u and u * r vanish on the boundary ∂ Ω of the computational domain, and that y obeys a homogeneous Neumann condition on the whole boundary; however, it is clear from the subsequent developments that similar results would hold for a Dirichlet boundary condition, provided it lies in the interval [0, 1]. Now at the fully discrete level, we assume that (ρ * K ) K∈T ∈ R N and (ρ K ) K∈T ∈ R N and a family of fluxes (F σ ,K ) K∈T , σ =K|L ∈ R 2M are given, and satisfy the following relations:

∀K ∈ T , ρ * K > 0, ρ K > 0, (2.8a) ∀K ∈ T , |K| δt [ρ K -ρ * K ] + ∑ σ =K|L F σ ,K = 0, (2.8b) ∀σ = K|L ∈ E int , F σ ,K = -F σ ,L . (2.8c)
This set of relations may for instance be derived by the upwind finite volume discretization of the mass balance equation (2.3)-(2.5) with F σ ,K = F up σ ,K (ρ, u), as in Section 3 below, thanks to (2.6) and Lemma 2.1. However, this is by far not the only way it may be obtained. For instance, the density ρ may be not a natural unknown for the problem: indeed, if instead of considering a compressible (gaseous) dispersed phase, we consider an incompressible (solid supension) phase, the density is then a function of y alone, and its positivity is ensured from the expression of the equation of state, provided that 0 y 1; the complete problem essentially becomes an incompressible one, and the mass balance may be seen as a constraint which requires the presence of a (dynamic) pressure in the momentum balance to be satisfied, and for which a centered discretization is natural (see [START_REF] Babik | On some fractional step schemes for combustion problems[END_REF] for the treatment of such a system). In some cases (and especially for incompressible flow problems where the mass balance, i.e. the divergence free constraint, must have a specific discretization to ensure the stability of the scheme), the finite volume mesh for the computation of y may even be different from the mesh used for the discretization of the mass balance; for instance, one will find in [START_REF] Calgaro | An hybrid finite-volume-finite element method for variable density incompressible flows[END_REF] a way to derive from a (possibly high order) mixed finite element solution of the momentum and mass balance equations an approximation for the velocity satisfying the mass conservation over a dual vertex-centered mesh.

With the previous notations, the discrete problem with unknown y = (y K ) K∈T ∈ R N considered in this section reads:

∀K ∈ T , |K| δt (ρ K y K -ρ * K y * K ) + ∑ σ =K|L F σ ,K y σ + ∑ σ =K|L G σ ,K Φ σ (y K , y L ) + D ∑ σ =K|L |σ | d σ (y K -y L ) = 0, (2.9)
where

• y σ = y K if F σ ,K 0 and y σ = y L otherwise (upwind choice),

• the quantity G σ ,K stands for the mass flux (i.e. the analogue to F σ ,K ) associated to the relative velocity u * r :

G σ ,K = ρ σ σ u * r • n KL
with, for ρ σ , any reasonable approximation for the density on σ , for instance ρ σ = 1 2 (ρ K + ρ L ) (centered choice), or ρ σ = ρ K if F σ ,K 0 and ρ σ = ρ L otherwise (upwind choice with respect to F σ ,K ),

• Φ σ (y K , y L ) stands for g(y K , y L ) if G σ ,K 0 and for g(y L , y K ) otherwise, g being a numerical monotone flux function with respect to φ in the sense of the following definition (see [START_REF] Eymard | Finite volume methods[END_REF] for the theory and some examples). 1. g(a, b) is non-decreasing with respect to a and non-increasing with respect to b, for any real numbers a and b, 2. g is Lipschitz-continuous with respect to both variables over R,

3. g(a, a) = φ (a), for any a ∈ R.
Then g is said to be a numerical monotone flux function with respect to the function φ .

Note that thanks to the definition of G σ ,K Φ σ (y K , y L ), the conservativity property holds, that is:

∀σ = K|L ∈ E int , G σ ,K Φ σ (y K , y L ) = -G σ ,L Φ σ (y K , y L ). (2.10)
Moreover, it is easily seen that G σ ,K Φ σ (y K , y L ) is non-decreasing with respect to y K and non-increasing with respect to y L , Lipschitz-continuous with respect to both variables over R, and that, for any real number a:

G σ ,K Φ σ (a, a) = ρ σ σ u * r • n KL φ (a). (2.11)
The result proven in this section is the following.

THEOREM 2.2 (EXISTENCE, UNIQUENESS AND L ∞ BOUNDS FOR A DISCRETE SOLUTION) Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (Φ σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying (2.8
). Let g be a numerical monotone flux function such that φ (x) = g(x, x) vanishes for x 0 and x 1. Then, if y * K ∈ [0, 1], ∀K ∈ T , there exists a unique solution to the discrete problem (2.9), which verifies

y K ∈ [0, 1], ∀K ∈ T .
This theorem summarizes a series of lemmata, which are detailed in the following subsections: first (section 2.3.1), we prove an a priori L ∞ estimate for y, precisely speaking the inequalities 0 y(x) 1, ∀x ∈ Ω ; then, on the basis of this bound, we apply a topological degree technique to obtain the existence of a solution (section 2.3.2); finally, this latter is shown to be unique (section 2.3.3).

2.3.1

An L ∞ stability property From a physical point of view, for instance thinking of the field y as a mass fraction, it seems natural for y to satisfy an "L ∞ stability property", more specifically to remain at any time in the [0, 1] interval. The aim of this section is to prove that this property holds for the solution of the scheme (2.9), provided that (2.8) holds at each time step and that the initial condition for y takes its values in [0, 1].

Let us first review the proof for the continuous problem, assuming all functions to be regular enough for the following calculations to make sense. Starting from the equation satisfied by y:

∂ t (ρy) + ∇ • (ρyu) + ∇ • (ρ φ (y) u r ) = ∇ • (D ∇y),
we first prove that y 0. Multiplying the previous equation by -y -and integrating over Ω yields:

-

Ω ∂ t (ρy) y -- Ω ∇ • (ρyu) y -- Ω ∇ • (ρ φ (y) u r ) y -- Ω D ∇y • ∇y -= 0.
(2.12)

Consider the first two terms of the previous relation, i.e. the terms associated to the advection operator:

T adv = - Ω ∂ t (ρy) y -- Ω ∇ • (ρyu) y -.
Expanding the derivatives, we obtain the so-called non-conservative form of the equation, and using the fact that when y -is non-zero, y = -y -, we obtain:

T adv = Ω [∂ t ρ + ∇ • (ρu)] [y -] 2 - Ω ρ y -∂ t y + (ρuy -) • ∇y .
The first term vanishes because of the mass balance equation, and the second one reads:

-

Ω ρ y -∂ t y + (ρuy -) • ∇y = 1 2 Ω ρ ∂ t ((y -) 2 ) + (ρu) • ∇(y -) 2 .
Hence, integrating by parts and using once again the mass balance equation, we get:

T adv = 1 2 Ω ρ ∂ t ((y -) 2 ) -(y -) 2 ∇ • (ρu) = 1 2 Ω ρ ∂ t ((y -) 2 ) + (y -) 2 ∂ ρ ∂t = 1 2 ∂ t Ω ρ (y -) 2 .
(2.13)

Substituting the term T adv in the relation (2.12) yields:

1 2 ∂ t Ω ρ (y -) 2 - Ω ∇ • (ρ φ (y) u r ) y -+ Ω D |∇y -| 2 = 0.
Since φ (x) vanishes for x 0, the second integral vanishes and we have:

1 2 ∂ t Ω ρ (y -) 2 = -D Ω |∇y -| 2 0.
Thus y is non-negative, provided that the initial condition for y is non-negative. Considering the equation satisfied by ỹ = 1y, one may prove similarly that y 1.

The proof we give in the discrete setting closely follows this calculation. The first step is thus to obtain an estimate for the terms related to the advection operator, which can be seen as a discrete counterpart to relation (2.13); this is achieved by the following lemma. As the mass balance plays a central role in the continuous setting, it is natural that the discrete mass balance (2.8b) also does in the discrete one, which is indeed the case.

LEMMA 2.2 Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (F σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying (2.8).
Then, for any family (y K ) K∈T ∈ R N , the following property holds:

-∑ K∈T y - K |K| δt (ρ K y K -ρ * K y * K ) + ∑ σ =K|L F σ ,K y σ 1 2 ∑ K∈T |K| δt ρ K (y - K ) 2 -ρ * K ((y * K ) -) 2 , (2.14)
where y σ = y K if F σ ,K 0 and y σ = y L otherwise.

Proof. The left-hand side of (2.14) may be written as:

T = T 1 + T 2 with T 1 = -∑ K∈T |K| δt y - K (ρ K y K -ρ * K y * K ) and T 2 = -∑ K∈T y - K ∑ σ =K|L F σ ,K y σ .
We first remark that when y - K is non-zero, y K = -y - K ; hence T 1 may be written as:

T 1 = ∑ K∈T |K| δt (y - K ) 2 (ρ K -ρ * K ) + ∑ K∈T |K| δt ρ * K y - K (y - K + y * K ).
Using the fact that

ρ * K y - K y * K -ρ * K y - K (y * K ) -,
we then obtain:

T 1 ∑ K∈T |K| δt (y - K ) 2 (ρ K -ρ * K ) + ∑ K∈T |K| δt ρ * K y - K [y - K -(y * K ) -];
thanks to the fact that ab 1 2 (a 2 + b 2 ), this yields:

T 1 ∑ K∈T |K| δt (y - K ) 2 (ρ K -ρ * K ) + 1 2 ∑ K∈T |K| δt ρ * K [(y - K ) 2 -((y * K ) -) 2 ],
which in turn gives:

T 1 1 2 ∑ K∈T |K| δt (y - K ) 2 (ρ K -ρ * K ) T 1,1 + 1 2 ∑ K∈T |K| δt ρ K (y - K ) 2 -ρ * K ((y * K ) -) 2 T 1,2
.

(2.15)

We now turn to T 2 . By conservativity (Relation (2.8c)), we have:

T 2 = - ∑ σ ∈E int (σ =K|L) (y - K -y - L ) F σ ,K y σ .
Let us assume, without loss of generality that the orientation of the edges σ = K|L is chosen such that F σ ,K 0. Then, for any σ = K|L, we have y σ = y K , and we get that:

T 2 = - ∑ σ ∈E int (σ =K|L) (y - K -y - L ) y K F σ ,K = ∑ σ ∈E int (σ =K|L) (y - K ) 2 F σ ,K + ∑ σ =K|L y - L y K F σ ,K .
Let us then write y -

L y K = 1 2 (y - L + y K ) 2 -1 2 (y - L ) 2 -1 2 (y K ) 2
, which yields:

T 2 = 1 2 ∑ σ ∈E int (σ =K|L) ((y - K ) 2 -(y - L ) 2 ) F σ ,K + 1 2 ∑ σ ∈E int (σ =K|L) ((y - K ) 2 -y 2 K + (y - L + y K ) 2 ) F σ ,K . (2.16)
The term (y - K ) 2y 2 K + (y - L + y K ) 2 is always non negative; indeed:

(y - K ) 2 -y 2 K + (y - L + y K ) 2 =    (y - L + y K ) 2 if y K 0, (y - L ) 2 + 2 y K y - L otherwise.
Hence, reordering the first summation in (2.16) and using (2.8b), we get:

T 2 1 2 ∑ K∈T ∑ σ =K|L (y - K ) 2 F σ ,K = - 1 2 ∑ K∈T |K| δt (y - K ) 2 (ρ K -ρ * K ) = -T 1,1 .
(2.17) Summing (2.15) and (2.17) then yields (2.14). We prove in the next tow lemmas that y remains in the interval (2.8). Let g be a numerical monotone flux function such that φ (x) = g(x, x) vanishes for x 0. Then, if y * K 0, ∀K ∈ T , the discrete solution of (2.9) also verifies y K 0, ∀K ∈ T .

[0, 1]. LEMMA 2.3 (LOWER BOUND ON y) Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (F σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying
Proof. As in the continuous case, the starting point is to multiply the equation by -y -, which, in the discrete case, consists in multiplying relation (2.9) by -y - K and summing over the control volumes. We get T adv + T nl + T dif = 0 with:

T adv = ∑ K∈T -y - K |K| δt (ρ K y K -ρ * K y * K ) + ∑ σ =K|L F σ ,K y σ , T nl = ∑ K∈T -y - K ∑ σ =K|L G σ ,K Φ σ (y K , y L ) , T dif = D ∑ K∈T -y - K ∑ σ =K|L |σ | d σ (y K -y L ) .
By Lemma 2.2:

T adv 1 2 ∑ K∈T |K| δt ρ K (y - K ) 2 -ρ * K ((y * K ) -) 2 .
Reordering the sum in the term T nl by conservativity (2.10), we have:

T nl = ∑ σ ∈E int (σ =K|L) T nl,K|L with T nl,K|L = -G σ ,K Φ σ (y K , y L ) (y - K -y - L ).
Let K and L be two neighbouring control volumes. If both y K and y L are non-negative, the term T nl,K|L vanishes. If y L 0, we get from (2.11) that Φ σ (y L , y L ) = φ (y L ) = 0, and thus:

T nl,K|L = -G σ ,K [Φ σ (y K , y L ) -Φ σ (y L , y L )] (y - K -y - L ).
Since the function G σ ,K Φ σ (•, •) is non-decreasing with respect to the first argument and the function x → x -is non-increasing, we obtain that T nl,K|L 0. Otherwise, y K is necessarily negative and we have:

T nl,K|L = -G σ ,K [Φ σ (y K , y L ) -Φ σ (y K , y K )] (y - K -y - L ),
which is also non-negative, since G σ ,K Φ σ (•, •) is non-increasing with respect to the second argument.

Let us now turn to the third term. Reordering the sum, we have:

T dif = - ∑ σ ∈E int (σ =K|L) D |σ | d σ (y - K -y - L ) (y K -y L ) 0.
Finally, we have:

1 2 ∑ K∈T |K| δt ρ K (y - K ) 2 -ρ * K [(y * K ) -] 2 0,
and thus, if (2.8). Let g be a numerical monotone flux function such that φ (x) = g(x, x) vanishes for x 1. Then, if y * K 1, ∀K ∈ T , the discrete solution of (2.9) also verifies y K 1, ∀K ∈ T .

(y * K ) K 0 ∀K ∈ T , then 1 2 ∑ K∈T |K| δt ρ K (y - K ) 2 0. LEMMA 2.4 (UPPER BOUND ON y) Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (F σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying
Proof. Thanks to (2.8b) and (2.9), we get the following discrete equation on the variable 1y:

|K| δt [ρ K (1 -y K ) -ρ * K (1 -y * K )] + ∑ σ =K|L F σ ,K (1 -y σ ) -G σ ,K Φ σ 1 -(1 -y K ), 1 -(1 -y L ) + D |σ | d σ [(1 -y K ) -(1 -y L )] = 0. Let g(•, •) be the function defined by g(a, b) = -g(1 -a, 1 -b).
This function is non-decreasing with respect to the first variable and non-increasing with respect to the second one. Moreover, g(x, x) = φ (x) = φ (1x) vanishes for x 0, as φ (x) vanishes for x 1. Thus the assumptions of Lemma 2.3 hold and, ∀K ∈ T , (1y) K is non-negative, which concludes the proof.

REMARK 2.1 (STRICT BOUNDS ON y) Strict bounds bounds on y may also be obtained. Indeed, if ∀K ∈ T , y * K > 0, then y satisfies the strict inequality ∀K ∈ T , y K > 0. To prove this result, let us assume that there exists K ∈ T such that y K =0. Replacing y K by zero in the equation (2.9) of the scheme, we get:

|K| δt (-ρ * K y * K ) + ∑ σ =K|L -F - σ ,K y L + G σ ,K Φ σ (0, y L ) + D |σ | d σ (-y L ) = 0.
The first term is by assumption negative, while, since y L 0, the second and last ones are non-positive. Since the function s → G σ ,K Φ σ (0, s) is non-increasing and Φ σ (0, 0) = 0, G σ ,K Φ σ (0, y L ) 0, and the third term also is non-positive, which contradicts the fact that the whole sum vanishes. By applying this result to 1y, we similarly prove that, if y * K < 1 ∀K ∈ T , then y K < 1 ∀K ∈ T . Returning to the initial physical problem, this result shows that, when using this scheme for the computation of the gas mass fraction y, monophasic zones cannot appear in the flow if they are not present at the initial time.

Existence for the approximate solution

The existence of a solution to the scheme (2.9) is obtained through a topological degree argument. We recall this result in the following theorem and refer to [10, chapter 5] for the general theory and [START_REF] Eymard | Error estimates for the approximate solutions of a nonlinear hyperbolic equation given by finite volume schemes[END_REF][START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] for its use in the case of other non-linear numerical scheme).

THEOREM 2.3 (APPLICATION OF THE TOPOLOGICAL DEGREE, FINITE DIMENSIONAL CASE)

Let (V, • ) be a normed finite dimensional vector space on R, let f be a continuous function from V to V and let b ∈ V . Let us assume that there exists a continuous function F from V × [0, 1] to V , and R > 0 satisfying:

(i) F (•, 1) = f ; (ii) For all α ∈ [0, 1], if v is such that F (v, α) = b then v ∈ B(0, R) (that is v < R);
(iii) the topological degree of F (•, 0) with respect to b and B(0, R) is equal to d 0 = 0.

Then the topological degree of F (•, 1) with respect to b and to B(0, R) is also equal to d 0 = 0; consequently, there exists at least a solution v ∈ B(0, R) such that f (v) = 0. (2.8). Let g be a numerical monotone flux function such that φ (x) = g(x, x) vanishes for x 0 and x 1. Then, if y * K ∈ [0, 1], ∀K ∈ T , there exists a solution to the considered discrete problem (2.9).

LEMMA 2.5 (EXISTENCE OF A DISCRETE SOLUTION

) Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (F σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying
Proof. In order to apply Theorem 2.3, we consider the space V = X T (recall that X T is the space of functions which are piecewise constant on each cell K) equipped with the max norm, which we denote by | • | ∞ ; we then introduce a function F : X T × [0, 1] → X T , such that F (., 1) = f where f is a function from X T to X T such that any solution of the nonlinear system (2.9) is a zero of f . For y = (y K ) K∈T ∈ X T and α ∈ [0, 1], this function F is defined by F (y, α) = (q K ) K∈T ∈ X T , with:

q K = |K| δt (ρ K y K -ρ * K y * K ) + ∑ σ =K|L F σ ,K y σ + α G σ ,K Φ σ (y K , y L ) + D |σ | d σ (y K -y L ) .
The function F is continuous from X T × [0, 1] to X T . We then remark that the lemmata 2.3 and 2.4 apply to the solution to the equation F (y, α) = 0, for 0 α 1. Hence, any solution to this equation belongs to:

B(0, 2) = {y ∈ X T , such that |y| ∞ < 2}.
Moreover, thanks to the estimate y K 1, ∀K ∈ T , the linear system F (y, 0) = 0 has a unique solution, which belongs to B(0, 2) (indeed, if it had two solutions a classical argument allows to conclude that there would be an unbounded set of solutions; existence follows from uniqueness since the dimension of the space is finite). From the existence of a solution to the linear system F (y, 0) = 0, we get that the topological degree of F (•, 0) with respect to B(0, 2) and 0 is non zero. Applying Theorem 2.3, we then get that the topological degree of F (•, 1) with respect to B(0, 2) and 0 is non zero, which in turn yields that there exists at least one solution to the nonlinear system (2.9).

Uniqueness of the approximate solution

The uniqueness of the solution to the scheme (2.9) is obtained through a duality argument. First, we introduce this technique in the semi-discrete time setting. Let y and ỹ be two smooth functions satisfying (2.7). Then the difference δ y = yỹ satisfies:

ρδ y δt + ∇ • (ρ δ y u) + ∇ • ρ φ (y) -φ ( ỹ) δ y δ y u r = ∇ • (D∇δ y).
Multiplying by a (smooth) test function ψ and integrating on Ω , we get:

1 δt Ω ρ δ y ψ + Ω ∇ • (ρ δ y u) ψ + Ω ∇ • ρ φ (y) -φ ( ỹ) δ y δ y u r ψ + D Ω ∇δ y • ∇ψ = 0. (2.18)
We then consider the following dual problem (with unknown ȳ and data y, ỹ and δ y, and which is consequently linear):

∀ψ, 1 δt Ω ρ ȳ ψ + Ω ∇ • (ρ ψ u) ȳ + Ω ∇ • ρ φ (y) -φ ( ỹ) δ y ψ u r ȳ + D Ω ∇ ȳ • ∇ψ = Ω δ yψ. (2.19)
Under some regularity assumptions, the dual problem (2. [START_REF] Marion | Navier-Stokes equations: Theory and approximation[END_REF]) is known to satisfy the maximum principle (e.g. [16, chapter 8]), and then, by an application of the Fredholm alternative, to admit a unique solution.

Taking as test function ψ = δ y in the dual problem, we get:

1 δt Ω ρ ȳ δ y + Ω ∇ • (ρ δ y u) ȳ + Ω ∇ • ρ φ (y) -φ ( ỹ) δ y δ y u r ȳ + D Ω ∇ ȳ • ∇δ y = Ω (δ y) 2 .
However, thanks to (2.18), the left-hand side of the previous relation is equal to zero, and therefore Ω (δ y) 2 = 0 so that δ y = 0, which proves the uniqueness of the solution. The proof that we give for the discrete problem is adapted from this technique. (2.8). Let g be a numerical monotone flux function. Then there exists at most one solution y ∈ X T to the discrete equation (2.9).

LEMMA 2.6 (UNIQUENESS OF THE DISCRETE SOLUTION

) Let (ρ K ) K∈T ∈ R N , (ρ * K ) K∈T ∈ R N and (F σ ,K ) K∈T , σ =K|L ∈ R 2M satisfying
Proof. Let y = (y K ) K∈T ∈ R N and ỹ = ( ỹK ) K∈T ∈ R N be two solutions of (2.9); then δ y = (y K -ỹK ) K∈T ∈ R N satisfies, for all K ∈ T :

|K| δt ρ K δ y K + ∑ σ =K|L F σ ,K δ y σ + G σ ,K Φ σ (y K , y L ) -Φ σ ( ỹK , ỹL ) + D ∑ σ =K|L |σ | d σ (δ y K -δ y L ) = 0. The quotients (Φ σ (a, •) -Φ σ (b, •)) / (a -b) and (Φ σ (•, a) -Φ σ (•, b)) / (a -b
) may be extended to the case a = b thanks to the fact that Φ σ is Lipschitz-continuous; we may thus write the above system as:

|K| δt ρ K δ y K + ∑ σ =K|L F σ ,K δ y σ + G σ ,K Φ σ (y K , y L ) -Φ σ ( ỹK , y L ) δ y K δ y K +G σ ,K Φ σ ( ỹK , y L ) -Φ σ ( ỹK , ỹL ) δ y L δ y L + D |σ | d σ (δ y K -δ y L ) = 0.
(2.20)

We now introduce the following discrete dual problem (with unknown ȳ ∈ X T and data y, ỹ and δ y ∈ X T ):

Find ȳ ∈ R N such that, ∀ψ ∈ R N , ∑ K∈T |K| δt ρ K ȳK ψ K + ȳK F σ ,K ψ σ + G σ ,K Φ σ (y K , y L ) -Φ σ ( ỹK , y L ) δ y K ψ K + Φ σ ( ỹK , y L ) -Φ σ ( ỹK , ỹL ) δ y L ψ L + D ȳK ∑ σ =K|L |σ | d σ (ψ K -ψ L ) = ∑ K∈T |K| (δ y K ) ψ K .
(2.21)

If there exists ȳ ∈ X T satisfying (2.21), then taking as a test function ψ = δ y we get from (2.20):

∑ K∈T |K| (δ y K ) 2 = 0,
which yields y = ỹ. In order to conclude the proof of the lemma, there only remains to show that there exists a (unique) solution to the dual problem (2.21). Let A be the matrix of the linear system (2.21), obtained with the natural ordering: the line of A associated to the control volume K is obtained by taking ψ = 1 K , where 1 K is the characteristic function of the element K. In this line, the diagonal entry is given by the term of the sum associated to K, and extra-diagonal entries are given by the terms of the sum corresponding to control volumes sharing an edge with K. Denoting by A KK this diagonal entry, we have:

A KK = |K| δt ρ K + ∑ σ =K|L F + σ ,K + G + σ ,K g(y K , y L ) -g( ỹK , y L ) y K -ỹK -G - σ ,K g( ỹL , y K ) -g( ỹL , ỹK ) y K -ỹK + D |σ | d σ .
Since the density is positive and the function g is non-decreasing with respect to its first argument and non-increasing with respect to its second argument, A KK is positive. Non-diagonal entries on the same line are given by:

A KL = -F + σ ,K -G + σ ,K g(y K , y L ) -g( ỹK , y L ) y K -ỹK + G - σ ,K g( ỹL , y K ) -g( ỹL , ỹK ) y K -ỹK -D |σ | d σ .
where, in this relation, L is a neighbouring control volume of K and σ = K|L. By the same arguments, these terms are non-positive. Moreover, the sum of all coefficient on a line reads:

∑ L∈T A KL = |K| δt ρ K ,
which is positive, since ρ K > 0. Thus A is a strictly diagonal dominant matrix; therefore it is invertible and problem (2.21) admits a unique solution, which completes the proof.

A fractional step algorithm for dispersed two-phase flows

In this section, we address the solution of the full system (1.1). To this purpose, we build a numerical scheme by complementing the discrete nonlinear advection-diffusion equation (2.9) with an incremental-projection-like algorithm.

Let us consider a partition 0 = t 0 < t 1 < . . . < t N = T of the time interval (0, T ), which, for the sake of simplicity, we suppose uniform. Let δt be the constant time step δt = t n+1 -t n for n = 0, 1, . . . , N -1. In a a semi-discrete time setting, the proposed algorithm consists in the following three steps scheme:

1 -Solve for y n+1 :

ρ n y n+1 -ρ n-1 y n δt +∇ • (ρ n y n+1 u n ) +∇ • ρ n y n+1 (1 -y n+1 ) u n r -∇ • (D ∇y n+1 ) = 0. (3.1)
2 -Solve for ũn+1 :

ρ n ũn+1 -ρ n-1 u n δt + ∇ • ( ũn+1 ⊗ ρ n u n ) + ∇p n -∇ • τ( ũn+1 ) = f n+1 . (3.2)
3 -Solve for p n+1 , u n+1 and ρ n+1 :

ρ n u n+1 -ũn+1 δt + ∇(p n+1 -p n ) = 0, (3.3a) ρ n+1 -ρ n δt + ∇ • (ρ n+1 u n+1 ) = 0, (3.3b
)

ρ n+1 = η(p n+1 , y n+1 ). (3.3c)
After a computation of the unknown y (step 1), step 2 consists in a semi-implicit solution of the momentum balance equation to obtain a predicted velocity.

Step 3 is a nonlinear pressure correction step, which degenerates in the usual projection step used in incompressible flow solvers when the density is constant (e.g. [START_REF] Marion | Navier-Stokes equations: Theory and approximation[END_REF]). Taking the divergence of (3.3a) and using (3.3b) to eliminate the unknown velocity u n+1 yields a non-linear elliptic problem for the pressure. This computation is formal in the semi-discrete formulation, but, of course, is necessarily made clear at the algebraic level, as described in section 3.3.

Once the pressure is computed, the first relation yields the updated velocity and the third one gives the end-of-step density.

The main difficulty to design such an algorithm lies in the approximation of the density. Indeed, we have to meet two requirements: first, to satisfy the compatibility condition (2.8b) when computing y (i.e. at Step 1 of the algorithm), second to ensure the conservativity of the scheme. The first point has been shown in Section 2.3 to be necessary to a reliable computation of the unknown y, and, from our experience, a violation of this condition may be at the origin of strong instabilities, for the estimation of y itself, but also for the whole algorithm. Still from our experience, using a non-conservative scheme for the approximation of y leads to large errors. This is especially important when the equation of state is strongly non-linear, as for flows involving phases of very different densities. To meet both requirements, we use here a time-shift of the density: in the advection terms of both the computation of y (Equation (3.1)) and the prediction of the velocity (Equation (3.2)), the density is taken one time step before the unknown y. This technique shows remarkable stability properties, but is of course limited to first order in time; this convergence property is assessed by numerical experiments.

REMARK 3.1 (ON ANOTHER TIME DISCRETIZATION OF THE DENSITY)

To satisfy condition (2.8b), another way to proceed has already been proposed [START_REF] Babik | On some fractional step schemes for combustion problems[END_REF][START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF]. The idea is to use, as a preliminary stage of the time step, the mass balance equation with a known value for the velocity (for instance, the velocity at the previous time step, or any extrapolation of it) to obtain a prediction of the density. If, for stability reasons, the discretization of this equation is chosen to be implicit, this preliminary step reads:

ρ -ρ n δt + ∇ • ( ρ ū) = 0,
where ū and ρ are the velocity used and the density obtained in this step, respectively. The first two terms of the balance equation for y (step 3.1 of the present algrotithm) are now:

ρ y n+1 -ρ n y n δt + ∇ • ( ρ y n+1 ũ) • • •
This approach can be easily modified to obtain a (formally) second order scheme [START_REF] Babik | On some fractional step schemes for combustion problems[END_REF]. Unfortunately, it seems difficult to consider ρ as the end-of-step value for the density, as its computation does not make use of the equation of state; this scheme thus cannot be conservative. In the present context, it may however be used at the first time step (and only at the first time step), to initialize the density by the following prediction step:

ρ 0 -ρ -1 δt + ∇ • (ρ 0 u -1 ) = 0,
where ρ -1 and u -1 are suitable approximations for the initial density and the velocity, respectively.

In order to obtain the full space-time discrete algorithm, we discretize (3.1) and (3.3b) by the finite volume method as described in Section 2. There remains to discretize the steps (3.2) (Paragraph 3.2 below) and (3.3a) (Paragraph 3.3) This is performed with the mixed Crouzeix-Raviart or Rannacher-Turek finite elements, which we now describe. This choice is quite convenient here since the discrete velocities are located on the edges of the mesh while the discrete pressures are piecewise constant on the cells; it is therefore compatible with the finite volume discretization described in Section 2. Furthermore, these elements are known to satisfy the inf-sup inequality, which contributes to the stability of the scheme, especially in the incompressible limit.

The Crouzeix Raviart and Rannacher-Turek finite elements

Let us briefly describe the Crouzeix-Raviart element for simplicial meshes (see [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF] for the seminal paper and, for instance, [11, p. 83-85] for a synthetic presentation), and the so-called "rotated bilinear element" introduced by Rannacher and Turek for quadrilateral or hexahedric meshes [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF].

In the following, we consider that the mesh of Ω which was introduced in Section 2 is either made of simplices (triangles in 2D or tetrahedra in 3D) or, in the case where the shape of Ω allows it, of rectangles or rectangular parallelepipeds. Note that conditions (i), (ii) and (iii) of Section 2.1 hold for both types of mesh, taking for x K the intersection of the orthogonal bisectors of the edges of K. The reference element for the Crouzeix-Raviart is the unit d-simplex and the discrete functional space is the space P 1 of affine polynomials. The reference element K for the rotated bilinear element is the unit d-cube (with edges parallel to the coordinate axes); the discrete functional space on K is Q1 ( K) d , where Q1 ( K) is defined as follows:

Q1 ( K) = span 1, (x i ) i=1,...,d , (x 2 i -x 2 i+1 ) i=1,...,d-1 .
For both velocity elements used here, the degrees of freedom are determined by the following set of nodal functionals:

{m σ ,i , σ ∈ E (K), i = 1, . . . , d} , m σ ,i (v) = |σ | -1 σ v i . (3.4)
The mapping from the reference element to the actual discretization cell is the standard affine mapping for the Crouzeix-Raviart element, and the standard Q 1 mapping for the Rannacher-Turek element. Finally, in both cases, the continuity of the average value of discrete velocities (i.e., for a discrete velocity field v, m σ ,i (v), 1 i d) across each face of the mesh is required, thus the discrete space V T is defined as follows:

V T = { v ∈ L 2 (Ω ) d : v| K ∈ Q1 (K) d , ∀K ∈ T ; m σ ,i (v) continuous across each edge σ ∈ E int , for 1 i d ; m σ ,i (v) = 0, ∀σ ∈ E ext , 1 i d }.
(3.5)

For both the Crouzeix-Raviart and the Rannacher-Turek discretizations, the approximation space for the pressure is the space X T of piecewise constant functions defined by (2.1) while the approximation space for the velocity is V T . Since only the continuity of the integral over each edge of the mesh is imposed, the velocity is discontinuous through each edge; the discretization is thus nonconforming in H 1 (Ω ) d . These pairs of approximation spaces for the velocity and the pressure are inf-sup stable, in the usual sense for "piecewise H 1 " discrete velocities, i.e. there exists c i > 0 which does not depend on the mesh such that:

∀p ∈ X T , sup v∈V T Ω ,T p ∇ • v ||v|| 1,b c i ||p -m(p)|| L 2 (Ω ) ,
where m(p) is the mean value of p over Ω , the symbol

Ω ,T stands for ∑ K∈T K
and || • || 1,b stands for the broken Sobolev H 1 semi-norm:

||v|| 2 1,b = Ω ,T |∇v| 2 = ∑ K∈T K |∇v| 2 .
From the definition (3.4), each velocity degree of freedom can be uniquely associated to an element edge. Hence, the velocity degrees of freedom may be indexed by the number of the component and the associated edge, and the set of velocity degrees of freedom reads:

{v σ ,i , σ ∈ E int , 1 i d}. We define v σ = ∑ d i=1 v σ ,i e (i)
where e (i) is the i th vector of the canonical basis of R d . We denote by ϕ (i) σ the vector shape function associated to v σ ,i , which, by the definition of the Crouzeix-Raviart and Rannacher-Turek finite elements, reads:

ϕ (i) σ = ϕ σ e (i) ,
where ϕ σ is the scalar basis function.

Space discretization of the momentum equation (1.1b)

As previously mentioned, we seek discrete approximations of u and p in the Crouzeix-Raviart or Rannacher-Turek finite element spaces, which we again denote u ∈ V T and p ∈ X T . With the notations of the previous section, we write the discrete functions as linear combinations of the basis functions:

u = ∑ d i=1 ∑ σ ∈E int u σ ,i ϕ (i)
σ and p = ∑ p K 1 K , where 1 K is the characteristic function K.

A blunt discretization of (3.2) with this finite element pair is in fact not quite satisfying for stability reasons. Indeed, we slightly modify it in order to obtain a finite volume type scheme on a dual mesh, which allows to respect a discrete counterpart of the L 2 -stability of the advection operator for the velocity, i.e. the discrete counterpart of the following relation:

Ω ∂ ρu ∂t + ∇ • (ρu) • u = d dt Ω 1 2 ρ|u| 2 , (3.6)
which is satisfied for any smooth functions ρ, u satisfying (1.1a)-(1.1b). This discrete counterpart is stated in theorem 3.1 below (see [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] for the proof). It is central in the proof of a priori estimates (e.g. the kinetic energy conservation theorem for incompressible flows) for the solution of the overall system, and has been found to be essential for convection dominant flows [START_REF] Babik | An L 2 -stable approximation of the Navier-Stokes advective operator for non conforming finite elements[END_REF]. It is also one of the ingredients used in [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] to derive a pressure correction scheme for compressible barotropic flows which conserves the entropy of the system. We state it below for the dual diamond mesh M which is constructed as follows: for each internal edge σ = K|L, let D K,σ be the cone with basis σ and the mass center of the cell K as opposite vertex. The volume D σ = D K,σ ∪ D L,σ is referred to as the "diamond cell" associated to σ and D K,σ is the half-diamond cell associated to σ and K (see Figure 1). The diamond cells D σ are naturally refered to by the primal edge σ ∈ E int , and their edges by the letter ε.

THEOREM 3.1 (STABILITY OF FINITE-VOLUME ADVECTION OPERATORS) Let M = (D σ ) σ ∈E int . Let (ρ σ ) σ ∈ E int ∈ R M , (ρ * σ ) σ ∈ E int ∈ R M and (F ε,σ ) σ ∈ E int , ε∈ E (D σ ) ∈ R kM (
where k is the number of sides of the diamond cells) satisfying:

∀σ ∈ E int , ρ * σ > 0, ρ σ > 0, (3.7a
)

∀σ ∈ E int , |D σ | δt (ρ σ -ρ * σ ) + ∑ ε∈E (D σ ) F ε,σ = 0 = 0, (3.7b
)

∀ε = σ |σ ′ , F ε,σ = -F ε,σ ′ . (3.7c)
(note that this is the dual mesh equivalent of the primal mesh assumptions (2.8)). Let (v * σ ) σ ∈E int and (v σ ) σ ∈E int be two families of real numbers. For any internal edge ε = D σ |D σ ′ , we define v ε either by

v σ = 1 2 (v σ + v σ ′ ) (centred choice), or by v ε = v σ if F ε,σ 0 and v ε = v σ ′ otherwise (upwind choice).
In both cases, the following inequality holds:

∑ D σ ∈M v σ   |D σ | δt (ρ σ v σ -ρ * σ v * σ ) + ∑ ε=D σ |D σ ′ F ε,σ v ε   1 2 ∑ D σ ∈M |D σ | δt ρ σ v 2 σ -ρ * σ v * σ 2 . (3.8)
Let us then derive a discretization scheme for (3.2) for which we are able to apply Theorem 3.1, taking for v the velocity components. Let us first note that, for the Crouzeix-Raviart element and, for the rectangular (in two dimensions) or cubic (in three dimensions) Rannacher-Turek element, one has K ϕ σ = |D K,σ |, where |D K,σ | is the measure of the half diamond cell D K,σ . Thus, a mass lumping of the finite element discretization of the term ρ n u n+1 in Equation (3.2) associated with σ leads to an expression of the form ρ n σ u n+1 σ , where ρ n σ is defined by:

∀σ ∈ E int , |D σ | ρ n σ = |D K,σ | ρ n K + |D L,σ | ρ n L . (3.9)
We also note that when the diffusion term ∇ • τ(u) reduces to the Laplace operator, its 2D Crouzeix-Raviart finite element discretization is identical to the finite volume discretization on the dual mesh consisting of the diamond cells D σ [START_REF] Chou | Analysis and convergence of a covolume method for the generalized Stokes problem[END_REF][START_REF] Chatzipantelidis | A finite volume method based on the Crouzeix-Raviart element for elliptic PDE's in two dimensions[END_REF]. This property readily extends to the Rannacher-Turek element and suggests that the advection term ∇

• ( ũn+1 ⊗ ρ n u n ) in (3.2) be discretized on each edge σ ∈ E int by the term ∑ ε∈E (D σ ) F n ε,σ ũn+1 ε , where E (D σ ) is the set of the edges of D σ , ũn+1
ε is a centred approximation of ũn+1 on ε and F n ε,σ = |ε| q n ε • n ε , where q n ε denotes an approximation of the momentum ρ n u n on the edge ε, |ε| is the measure of ε and n ε is the normal to ε outward D σ .

FIG. 2. sub-volume of K.

We then need to express q n ε • n ε in such a way that the discrete mass balance (3.7) holds. To this goal, we use the following result [START_REF] Babik | An L 2 -stable approximation of the Navier-Stokes advective operator for non conforming finite elements[END_REF] and give its (elementary) proof for the sake oompleteness.

LEMMA 3.1 (MASS BALANCE IN A SUB-VOLUME OF

A MESH) Let K ∈ T , let (ρ * , ρ) ∈ (R + ) 2 , and consider a family (F σ ,K ) σ ∈E (K) ⊂ R k such that: |K| δt (ρ -ρ * ) + ∑ σ ∈E (K) F σ ,K = 0. ( 3.10)
Let w be a momentum field on K, such that ∇ • w is constant over K and satisfying:

K ∇ • w = ∂ K w • n ∂ K = ∑ σ ∈E (K) F σ ,K , (3.11)
where ∂ K and n ∂ K stand for the boundary of K and the normal vector to ∂ K outward to K, respectively. Let D be a subset of K with boundary ∂ D (see Figure 2). Then the following property holds :

|D| δt (ρ -ρ * ) + ∂ D w • n ∂ D = 0,
where n ∂ D stands for the normal vector to ∂ D outward to D.

Proof. Multiplying (3.10) by |D|/|K| and using (3.11) yields that:

|D| |δt| (ρ -ρ * ) + |D| |K| K ∇ • w = 0,
which concludes the proof thanks to the fact that ∇ • w is constant over K. Now let us apply this lemma to obtain the mass conservation property on the diamond mesh M . At step n, the discrete mass balance for ρ n is obtained from the solution of (3.3b) at step n -1; its (finite volume) discretization is (2.3) with ρ = ρ n , ρ * = ρ n-1 and u = u n :

K δt (ρ n K -ρ n-1 K ) + ∑ σ ∈E (K) F up σ ,K (ρ n , u n ) = 0, ∀K ∈ T . (3.12)
We construct the momenta on the edges ε ∈ E (D σ ) by building on each cell K a field w n with constant divergence and such that:

∀σ ∈ E (K), σ w n • n σ = F up σ ,K (ρ n , u n ). (3.13)
Assuming the field w n to be constructed, the value F n ε,σ of the mass flux on the edge or face ε of the diamond cell is then computed by integrating w n • n ε over ε, i.e.:

F n ε,σ = ε w n • n ε , (3.14)
which yields a discrete mass balance over both half-diamond cells; thanks to the fact that ρ n σ is defined by (3.9), summing the two discrete balance equations on the half diamonds D K,σ and D L,σ gives the discrete mass balance over D σ = D K,σ ∪ D L,σ , that is:

|D σ | δt (ρ n σ -ρ n-1 σ ) + ∑ ε∈E (D σ )
F n ε,σ = 0; Condition (3.7b) of Theorem 3.1 is therefore satisfied.

Let us then turn to he construction of the field w n ; such a field is derived for the Crouzeix-Raviart element by direct (i.e. using the standard expansion of the Crouzeix-Raviart elements) interpolation of quantities ((ρ u) n σ ) σ ∈E (K) :

w n (x) = ∑ σ ∈E (K) ϕ σ (x)(ρ u) n σ where (ρ u) n σ is such that σ (ρ u) n σ • n σ = F up σ ,K (ρ n , u n ).
With the chosen discretization for the mass balance (see (2.4)), the natural choice for (ρ u) n σ reads (ρ u) n σ = ρ n up,σ u n σ where ρ n up,σ is the upwind density on σ with respect to u n . Indeed, thanks to the fact that the Crouzeix Raviart basis functions satisfy σ ϕ σ = |σ |, one has σ w n • n σ = F up σ ,K (ρ n , u n ). For the Rannacher-Turek element, the divergence of discrete functions is no longer constant, but when the mesh is rectangular or cubic, we may use the following interpolation formula:

w n (x) = ∑ σ ∈E (K) α σ (x • n σ ) F up σ ,K (ρ n , u n ) n σ
where the α σ are affine interpolation functions which are determined in such a way that the relations (3.13) hold. The extension to more general grids is underway.

Finally, the standard (Crouzeix-Raviart) finite element expansion is used to discretize the terms ∇p n -∇ • τ(u n+1 ) of (3.2), and we thus obtain the following discrete momentum balance equation:

∀σ ∈ E int , for 1 i d, |D σ | δt (ρ n σ ũn+1 σ ,i -ρ n-1 σ u n σ ,i ) + ∑ ε∈E (D σ ), ε=D σ |D σ ′ 1 2 F n ε,σ ( ũn+1 σ ,i + ũn+1 σ ′ ,i ) - Ω ,T p n ∇ • ϕ (i) σ + a d ( ũn+1 , ϕ (i) σ ) = Ω f n+1 • ϕ (i) σ (3.15)
where

a d (v, w) = µ Ω ,T ∇v : ∇w + 1 3 ∇ • v ∇ • w if (1.
3) holds (case of constant viscosity), Ω ,T τ(v) : ∇w with τ given by (1.2) otherwise.

Spatial discretization of the projection step

The discretization of the first projection step (3.3a) is consistent with that of the momentum equation(3.2); a mass lumping is performed for the unsteady term and a standard finite element formulation is used for the gradient of the pressure increment, yielding:

∀σ ∈ E int , for 1 i d, |D σ | δt ρ n σ (u n+1 σ ,i -ũn+1 σ ,i ) - Ω ,T (p n+1 -p n ) ∇ • ϕ (i) σ dx = 0,
where ρ n σ is defined by (3.9). Since the pressure is piecewise constant, the discrete gradient operator takes the form of the transposed of the finite volume standard discretization of the divergence (based on the primal mesh T , and not on the diamond mesh M ) and can be rewritten as follows:

∀σ ∈ E int , σ = K|L, |D σ | δt ρ n σ (u n+1 σ -ũn+1 σ ) + |σ | (p n+1 L -p n L ) -(p n+1 K -p n K ) n KL = 0. (3.16)
The mass balance equation (3.3a) is discretized as described in Section 2, which ensures that the density stays positive at all times; taking ρ = ρ n+1 , ρ * = ρ n and u = u n+1 in the finite volume scheme (2.3) yields:

∀K ∈ T , |K| δt ρ n+1 K -ρ n K + ∑ σ =K|L F up σ ,K (ρ n+1 , u n+1 ) = 0, (3.17a)
ρ n+1 K = η(p n+1 K , y n+1 K ). (3.17b)
The projection step therefore consists in finding (u n+1 σ ) σ ∈E int and (p n+1 K ) K∈T solution of the non linear system (3.16)- (3.17). Under some assumptions for the function η, we now prove that this projection step admits one solution. LEMMA 3.2 Let us suppose that the equation of state η is such that for any y ∈ [0, 1], the function p → η(p, y) is defined and increasing on [0, +∞), η(0, y) = 0 and lim p→+∞ η(p, y) = +∞. Then the non-linear algebraic system (3.16)-(3.17) admits at least a solution.

Proof. The proof of this lemma consists in an application of the Brouwer fixed point theorem (see e.g. [10, chapter 5]

). Let ũn+1 = ( ũn+1 σ ) σ ∈E int ∈ R M×d , y n+1 = (y n+1 K ) K∈T ∈ R N , ρ n = (ρ n K ) K∈T ∈ R N and p n = (p n K ) K∈T ∈ R N be
known families of velocities, mass fractions, densities and pressures (calculated at the previous iteration or at the prediction step), such that ρ n K > 0 for all K ∈ T . Let H be the function from

V T × X T to itself defined by (u, p) = ((u σ ) σ ∈E int , (p K ) K∈T ) → H(u, p) = (v, q) = ((v σ ) σ ∈E int , (q K ) K∈T ), with (v, q) satisfying: ∀K ∈ T , |K| δt [ρ K -ρ n K ] + ∑ σ =K|L F up σ ,K (ρ, u) = 0, (3.18a) ∀K ∈ T , ρ K = η(q K , y n+1 K ), (3.18b) ∀σ = K|L ∈ E int , for 1 i d, |D σ | δt ρ n σ (v σ ,i -ũn+1 σ ,i ) + Ω ,T (q -p n ) ∇ • ϕ (i) σ dx = 0, (3.18c)
where F up σ ,K is defined by (2.4) and ρ n σ by (3.9); note that ρ n σ > 0 for any σ ∈ E int . Equation (3.18a) is linear in ρ, and by Lemma 2.1, it has a unique solution (ρ K ) K∈T satisfying ρ K > 0 for all K ∈ T . Thanks to the assumption on η, for each K ∈ T , there exists a unique q K such that (3.18b) is satisfied. Finally, when q is computed, Equation (3.18c) yields a diagonal system for v. Hence the system (3.18) has a unique solution, so that the function H is well defined. We then note that any fixed point of the function H is a solution to the system (3.16)- (3.17). By conservativity of the finite volume scheme (3.17a), we easily see that ∑ K∈T |K| ρ K = ∑ K∈T |K| ρ n K , and therefore, there exists c ρ ∈ R such that 0 < max K∈T ρ K < c ρ and again from the assumption on η, also for the pressure: max K∈T q K c p . Multiplying (3.18c) by v σ ,i and summing over σ ∈ E int and i = 1, . . . , d, we get that:

∑ σ ∈E int |D σ | ρ n σ |v σ | 2 ∑ σ ∈E int |D σ | ρ n σ v σ ũn+1 σ + δt Ω ,T (q -p n ) ∇ • v.
Using the positivity of ρ n σ and the fact that all norms are equivalent on the finite dimensional space, we get that there exists c u 0 such that v

2 2 = ∑ σ ∈E int |D σ | |v σ | 2 c 2 u .
By the Brouwer theorem, the mapping H therefore admits a fixed point in the convex set C = {(v, q) ∈ X T × V T such that v 2 c u and max K∈T q K c p }. THEOREM 3.2 (PROPERTIES OF THE NUMERICAL SCHEME) Under the assumptions for the equation of state of Lemma 3.2, there exists a set of families (u n ) n 0 , (p n ) n 0 , (ρ n ) n -1 and (y n ) n 1 given by the proposed algorithm. Moreover, the following properties hold, for all n N: (i) positivity of the density:

ρ n K > 0, ∀K ∈ T . (ii) L ∞ stability property: y n K ∈ [0, 1], ∀K ∈ T .
(iv) conservativity property:

∑ K∈T |K| ρ n-1 K y n K = M 0 y and ∑ K∈T |K| ρ n = M 0 , ∀n 1,
where M 0 y (resp. M 0 ) denotes the initial gas mass (resp. total mass) in Ω .

Proof. The existence, uniqueness and positivity of ρ 0 satisfying (3.22) follows from Lemma 2.1, which also gives that ρ 0 K > 0 for any K ∈ T . Thanks to the assumptions on η and since ρ 0 > 0 and y 0 ∈ [0, 1], p 0 is uniquely determined by the equation of state ρ 0 K = η(p 0 K , (y 0 K ). Then, we get from Theorem 2.2 that there exists a unique family (y 1 K ) K∈T satisfiying (3.23), and that y 1 K ∈ [0, 1]. Summing Equation (3.23) over K ∈ T and using the conservativity of the fluxes (2.6) and (2.10) yields that ∑ K∈T |K| ρ 0 K y 1 K = M 0 y . The predicted velocity ũ1 satisfies the set of linear equation (3.15). Multiplying (3.15) by ũ1 σ , summing over the edges σ ∈ E int and using Theorem 3.1 (which we may use thanks to the careful discretization of the convection term) yields that there exists C ∈ R + such that ũ1 C 1 f L 2 (Ω ) , which in turn yields the existence and uniqueness of ũ1 . We then obtain from Lemma 3.2 that there exists u 1 , p 1 and ρ 1 satisfying (3.16)-(3.17), and from Lemma 2.1, that ρ 1 > 0. Summing Equation (3.17a) over K ∈ T and using the conservativity of the fluxes (2.6) yields that ∑ K∈T |K| ρ 1 K = M 0 . The proof of theorem is then completed by an easy induction.

Numerical results

In this section, we present two numerical tests performed to assess the behaviour of the above described fractional step scheme. In these tests, we compute the flow of an isothermal two-phase mixture of nonmiscible liquid and gas. In this case, the model (1.1) is the so-called drift-flux model, that is a mixture model that takes into account the relative velocity u r between the liquid and the gas phase (the so-called drift velocity), for which a phenomenologic relation must be supplied. Then, ρ stands for the mixture density and takes the general form ρ = (1 -α g )ρ ℓ + α g ρ g (p) where α g stands for the void fraction and ρ g (p) yields the gas density as a function of the pressure. Introducing the mass gas fraction y and using the relation α g ρ g = ρ y leads to the equation of state: the L 2 norm for the velocity and in the discrete L 2 norms for the pressure and the gas mass fraction. For large time steps, these curves show a decrease corresponding to approximately a first order convergence in time, until a plateau is reached, due to the fact that errors are bounded by below by the residual spatial discretization error. The value of the errors on this plateau then shows a spatial convergence order close to one, which is consistent with the choice of an upwind discretization for the advection terms in the gas mass fraction and mass balance equations. Finally, results seem to be significantly more accurate with the Rannacher-Turek element.

η(p, y) = 1 y/ρ g (p) + (1 -y)/ρ ℓ . ( 4

A phase separation problem

We now present numerical results obtained for a phase separation problem, with data inspired from a classical benchmark test for the simulation of two-phase flows [START_REF] Coquel | A numerical method using upwind schemes for the resolution of two-phase flows[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Staedtke | Advanced threedimensional two-phase flow simulation tools for application to reactor safety (ASTAR)[END_REF] with two-fields models (i.e. models considering separate balance equations for each phase). The considered physical domain is a vertical tube of length L = 7.5 m, filled at initial time with a two-phase mixture of air and water with α = 0.5, u = 0 and p = p 0 where p 0 = 10 5 Pa is the ambiant pressure. Under the action of gravity (with g = 9.81 m.s -2 ), phases separate and the solution at t = +∞ is the superposition of a zone of pure water and a zone of pure air, both at rest. In the original problem, the interactions between both phases are neglected; instead, we assume here that the relative velocity is constant and given by u r = 1 m.s -1 , which is clearly non-physical (at small times, water droplets just fall with a constant acceleration g). However, even under this assumption, the solution of the problem qualitatively reproduces the original phase separation phenomenon.

The equation of state for the mixture is the same as in the previous test case and the densities, for water and air respectively, are ρ ℓ = 1000 kg.m -3 and ρ g = p/RT where RT is such that ρ g = 1.2 kg.m -3 at p = 10 5 Pa. The diffusion coefficient D and the viscosity µ are set to zero. At the top and the bottom boundaries, both the velocity and the relative velocity are prescribed to zero.

For this test case, we use a regular mesh composed of rectangular cells (with the Rannacher-Turek element); since this problem is one-dimensional, only one cell is used in the horizontal direction, and 200 in the vertical one. Calculations with time steps up to δt = 10 -1 s have been performed without observing any instability. With respect to the time discretization, the convergence for the void fraction and the density is readily achieved, and profiles obtained with δt 10 -2 s are all similar; the results with δt = 10 -2 s are reported on Figure 10. However, probably because of the very steep variations of the density near y = 1, themselves due to the large difference of the densities of the two phases, convergence for y is more difficult to reach, and variations of the obtained profiles are observed when decreasing the time step down to δt = 5.10 -4 s.

Conclusion

In this paper, we first address a parabolic equation which models the phases mass balance in twophase flows. This equation differs from the mass balance for chemical species in compressible multicomponent flows studied in [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF] by the addition of a non-linear term of the form ∇ • ρϕ(y) u r , where y is the unknown, ϕ(•) is a regular function such that ϕ(0) = ϕ(1) = 0 and u r is a general (in particular non-necessarily divergence free) velocity field. We prove the existence and uniqueness of the finite volume approximation together with the fact that it remains within physical bounds, i.e. within the interval [0, 1]. As in [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF], the necessary condition for this L ∞ stability is that the discretization of the convection operator must be such that it vanishes for constant y, which amounts to demanding that a particular discrete mass balance equation be satisfied. The second ingredient of this scheme is a discretization of the non-linear term based on the notion of monotone flux functions [START_REF] Eymard | Finite volume methods[END_REF]. This work extends the theory developed in [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF] in two directions: it copes with a new non-linear term, and introduces different techniques well-suited for non-linear problem, first to prove L ∞ a priori estimates for the solution (recasting the equation under variational form and choosing the possible negative part of the solution as test function), to prove its existence (by a topological degree argument) and its uniqueness (introducing an auxiliary "linear dual problem").

In a second part, we propose a discretization by a fractional step method for the set of equations composed of Navier-Stokes equations (i.e. mass and momentum balance) and the phases mass balance. This formulation decouples the resolution of the phases mass balance from the resolution of the Navier-Stokes equations, which is performed by a pressure correction method based on low degree nonconforming finite elements. This algorithm meets two essential requirements: it is conservative, and the discrete mass balance needed for a stable computation of y is satisfied. To achieve this goal, the key ingredient is a particular time-discretization of the density terms, which unfortunately limits the time accuracy of the scheme to first order. This technique is now routinely used in the ISIS computer code [START_REF] Babik | On some fractional step schemes for combustion problems[END_REF] developed at IRSN and devoted to the modelling of reacting flows; in this context, it demonstrates very satisfactory stability properties, even for cases where very steep variations of the density appear, and for which numerical difficulties are often reported in the literature. Finally, let us also mention that the proposed numerical scheme degenerates to classical projection methods in the incompressible limit.

As far as extensions of this work are concerned, first, (formally) second order in space discretizations (typically using MUSCL-like techniques) should be developed. Second, the proposed algorithm is based on the simplest fractional step approach (all the equations are decoupled) and, consequently, the less time-consuming one, and this choice should probably be retained as far as it works. Unfortunately, in the specific case of two-phase compressible flows involving phases of very different densities, instabilities are observed, the cure of which seems to need a drastic reduction of the time step. These instabilities appear to be linked to the fact that the present algorithm does not preserve a constant pressure through moving interfaces between phases (i.e. contact discontinuities of the underlying hyperbolic system); a solution to this problem, still based on the same essential ingredients for the evaluation of the density terms and the discretization of the phases mass balance but coupling this latter equation to the projection step, is now under development and shows promising results.

(

  ii) There exists a set E of bounded subsets of hyperplanes of R d included in Ω , which are the edges (in 2D) or faces (in 3D) of the cells K ∈ T . The set of boundary edges or faces (i.e. the edges or faces included in the boundary ∂ Ω of the domain Ω ) is denoted by E ext and the set of internal ones

Definition 2 . 1 [

 21 Numerical monotone flux function] Let the function g ∈ C(R 2 , R) satisfy the following assumptions:

FIG. 1 .

 1 FIG. 1. Diamond-cells for the Crouzeix-Raviart and Rannacher-Turek element.

FIG. 3 .

 3 FIG. 3. Pattern used to split square cells into acute-angle triangles, for the Crouzeix-Raviart element.

FIG. 5 . 32 FIG. 6 .

 5326 FIG. 5. Crouzeix-Raviart element -Velocity error as a function of the time step.

FIG. 8 .

 8 FIG.8. Rannacher-Turek element -Gas mass fraction error as a function of the time step.

Aknowledgment

The authors would like to thank the referees for thier useful suggestions which lead to a significant improvement of the paper.

 REMARK 3.2We know from Remark 2.1 that, if y n > 0 (y n < 1), then y n+1 > 0 (resp. y n+1 < 1), and thus, by induction, that pure monophasic zones does not appear if not already present at initial time. Hence, if no monophasic zone exists at initial time, it is sufficient for the existence of the solution that assumptions on the equation of state of lemma 3.2 hold for y ∈ (0, 1); this allows to deal with cases where one phase is supposed to be incompressible.

Let us now combine the two algebraic relations (3.16) and (3.17) to build a discrete elliptic problem for the pressure. To this purpose, let us introduce the algebraic formulation of this system: .19b) where:

• M ρ n stands for the diagonal mass matrix weighted by the density ρ n σ at t n (at edges or faces center) (i.e. the d diagonal entries of M(ρ n ) associated to the edge σ are given by |D σ | ρ n σ );

and B is the opposite of the divergence operator matrix.

• Q up ρ n+1 ,u n+1 is a diagonal matrix; its entry corresponding to an edge σ = K|L ∈ E int is ρ n+1 up,σ , i.e. the upwind density with respect to u n+1 (see (2.5)). The matrix R is diagonal and, for any

The elliptic problem for the pressure is obtained by multiplying (3.19a) by B Q up ρ n+1 ,u n+1 (M ρ n ) -1 and using (3.19b). The resulting equation reads:

where L up

be evaluated through a finite volume type expression of the Laplace operator [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF]:

Note that, even if ∀σ ∈ E int , ρ n+1 up,σ = ρ n σ , the operator L up ρ n+1 ,u n+1 differs from the usual finite volume Laplace operator (in the case of rectangles or cubes, by a factor d = 2 or 3 respectively [START_REF] Gallouët | An unconditionnally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF]); this is linked to the fact that the non conforming (V T , X T ) finite element approximation is not consistent for the mixed form approximation of the Laplace equation, neither in the Crouzeix-Raviart case nor if in the Rannacher-Turek case. Provided that p n+1 is known, (3.19a) gives the updated value of the velocity:

(3.21)

Note that contrary to usual projection methods, equations (3.20) and (3.21) are not decoupled, because of the upwinding of the density with respect to u n+1 . We thus implement an iterative algorithm, which reads when the equation of state is linear with respect to the pressure:

Step 4.1 -Solve for p n+1 k+1 :

Step 4.2 -Compute u n+1 k+1 as :

When the equation of state is nonlinear with respect to p, which is in general the case, step 4.1 is replaced by one iteration of a quasi Newton algorithm where only the diagonal term ρ(p n+1 k+1 , y n+1 ) is differentiated with respect to p n+1 k+1 .

The fully discrete algorithm

To sum up, we consider the following algorithm:

with F up σ ,K defined by (2.4), and let (p 0 K ) K∈T be given by ρ 0 K = η(p 0 K , y 0 K ). Then, for n = 0, 1, 2, . . .:

1. Computation of y -Compute (y n+1 K ) K∈T by an upwind finite volume discretization, as explained in Section 2: The following theorem gathers the properties of the proposed numerical scheme.

In the perfect gas approximation and for a constant temperature, ρ g (p) is simply proportional to the pressure:

where R is the gas constant and T is the absolute temperature. We further assume that the liquid phase can be considered as incompressible so that the liquid density ρ ℓ is constant. In this section, the nonlinear function φ is given by φ (y) = y(1y) and the flux function by g(a, b) = ab 2 with, ∀x ∈ R, x = max(0, min(x, 1)).

Assessing the convergence against an analytic solution

We first assess the convergence rate of the proposed scheme with respect to space and time discretizations, by dealing with a case where an analytic solution can be exhibited.

We choose for the computational domain Ω = (0, 1) × (-0.5, 0.5), and for the momentum and density the following expressions:

The pressure and the gas mass fraction are linked to the density by the equation of state (4.1), where the liquid density ρ ℓ is set at ρ ℓ = 5 and the product RT in the equation of state of the gas (4.2) is given by RT = 1 (so ρ g = p). We choose the following expression for the unknown y: y(x,t) = 2.5 -0.5 ρ(x,t) 4.5 ρ(x,t)

The relative velocity is constant and given by u r = (0, 1) t and the diffusive coefficient D is equal to 0.1. The analytical expression for the pressure is obtained from the equation of state. These functions satisfy the mass balance equation; for the gas mass fraction and momentum balance, we add the corresponding right-hand side. In this latter equation, we suppose that the divergence of the stress tensor is given by: ∇

and we use the corresponding expression for the bilinear form a d (•, •).

For the Rannacher-Turek element, computations are made with 20×20, 40×40 and 80×80 uniform meshes. For the Crouzeix-Raviart one, the meshes are built as follows: the computational domain is first split in square subdomains, then each subdomain is split in 26 simplices, all having angles of at most 80 • , according to the pattern given in [3, Figure 5 -bbbb], see Figure 3 Velocity, pressure and gas mass fraction errors obtained at t = 0.5 as a function of the time step are drawn, for the Rannacher-Turek element, on Figure 4, Figure 6 and Figure 8, respectively and, for the Crouzeix-Raviart element, on Figure 5, Figure 7 and Figure 9, respectively. These errors are evaluated in