Federica Olivero
email: fede.olivero@bristol.ac.uk

Domingo Paola
email: paola.domingo@mail.sirio.it

Ornella Robutti
email: robutti@dm.unito.it

Teaching proof in a dynamic geometry environment: what mediation?

This paper aims to discuss two issues. First we intend to analyse if and how dynamic geometry software prove to be a communication space for students and how they support students towards proving. Second, we reflect about the use of video-tapes as source of information both for Mathematics Education research and for the teaching practice. In particular we approach the question: how do video-tapes help in the construction of a shared experience in the classroom?

Resumé Ce travail se propose de réfléchir sur deux thèmes. En premier lieu nous voulons analyser si et comment l'emploi d'un logiciel en géométrie dynamique peut se proposer comme un espace de communication pour les étudiants et comment cela peut les aider dans l'activité de démonstration. En deuxième lieu nous réfléchissons sur l'emploi des vidéo -enregistrements en qualité de source d'information soit dans la recherche, soit dans la pratique didactique. En particulier on considère la question suivante : dans quelle façon les vidéo -enregistrements peuvent nous aider dans la construction d'une expérience partagée en classe.

Introduction

"To be mathematically literate means to understand mathematics as a part of the cultural heritage. Mathematical experiences have to be offered to all students, that include knowledge about the cultural and the historical background of mathematics." [START_REF] Manifesto | Where we are and where we go[END_REF].

Mathematical literacy involves both the aspects of understanding and communicating.

Moreover, it is important to consider the necessity of a literacy for all the students, so that everyone is given the possibility of participating into a mathematical experience. Nowadays, the presence of technology in our society influences the school and the educational change. The research that has been carried out over the years suggests that, within the problem of mathematical literacy for all, the main issues to be addressed, with respect to the use of new technologies, concern the potentialities they may offer, compared to more traditional tools, with respect to specific contexts of use. One of these potentialities may be to foster and widen the possibilities of understanding and communicating mathematical knowledge. The part of that knowledge we focus on in this paper concerns the theoretical organisation of mathematics, in its aspects of understanding and communicating. The chosen domain is Euclidean geometry and the technology we used is the microworld Cabri-Géomètre.

The "treasure hunting problem"

The research project we have been involved in for some years now within the Mathematics Education group at the University of Turin and the Graduate School of Education in Bristol, is aimed at investigating the potentialities of the dynamic geometry software Cabri-Géomètre [START_REF] Baulac | Cabri-Géomètre, un logiciel d'aide à l'apprentissage de la géomètrie[END_REF] in supporting students' proving process in geometry (e.g. [START_REF] Arzarello | Dragging in Cabri and modalities of transition from conjectures to proofs in geometry[END_REF]. By proving process we mean the process of exploring, conjecturing and proving in an open problem. The research methodology adopted involves the use of videotaping in the classroom, that provides a rich account of all actions and language of the students.

In this paper we present the analysis of the video of a classroom session which was part of a classroom experiment1 that involved five Further Mathematics A-level students in Bristol. It was aimed at introducing students to the proving process in geometry within a dynamic geometry environment. The teacher run most sessions. The students have a solid mathematical background.

They did geometry only at GCSE. They had never used Cabri before. In this experiment, we carried out a brief introduction of Cabri, but we did not cover all the possibilities it offers. The only thing we explained them was the dragging test. After that, other Cabri features were introduced during the sessions when needed, so to make them part of the classroom culture.

In the session we analyse in this paper, the students were asked to tackle the following problem (the problem was proposed and solved by George Gamow in: "Uno due tre…infinito" Arnoldo Mondadori Editore, 1952, Milano):

A man left instructions about how to discover the place where a treasure had been hidden in an island, knowing only the position of three trees: a pine, an oak and an apple tree. In order to discover the treasure a man has to position himself where the apple tree is and then: 1. walk till the oak counting the steps, then turn right2 (90 degrees), walk the same number of steps and put a mark where he arrives. 2. come back to the apple tree, walk to the pine counting the steps, then turn left (90 degrees), walk the same number of steps and put a mark where he arrives. The treasure will be in the midpoint of the line joining the two marks. The man who had these instructions goes to the island; when he gets there he finds out that the apple tree is missing. Will he be able to discover the treasure?

In the session, the students are divided in pairs. They are asked to work with paper and pencil first until they get a conjecture, and then they can work in Cabri. The session lasts 1 hour and 20 minutes. Almost all the students get to a conjecture. In the following session (1 hour and 30 minutes) they are asked to go on working on this problem and try to construct a proof. The process the observed 3 students followed in the solution of this problem can be divided into three phases. The first phase is carried out in paper and pencil. The students make the conjecture: "M (the treasure) is independent from A" (the apple tree) right at the beginning. In this first phase a logical conflict is present: does A (the apple tree) exist or not? The students make two different drawings, which represent two worlds: one with the apple tree and one without the apple tree. Then, they reconstruct the whole situation starting from two separate drawings, one containing the pine and one containing the oak. The second phase concerns the work in Cabri. The students start with a construction, which is not general, because the point corresponding to the apple tree is not constructed as a free point. After they realise this, they perform a new construction and carry out a lot of explorations, but they are not able to discover the procedure to construct M. This phase too is characterised by a logical conflict: all points can be moved, but what are the important relationships? The third phase concerns the attempts to proving that M is fixed. The students work with a static Cabri figure and they are not able to construct a proof.

In the following, we will analyse two short excerpts from the students' protocol, with respect to the following issues: does Cabri provide a communication space? Does it support the proving process?

Developing a communication space

In this first excerpt, the students (I and S) are at the end of the paper and pencil phase, they are exploring the situation and try to make sense of what happens when the apple tree is not there. 92. S: no, picking that (P on Fig1) as the point.

I: so he walks out (I points on fig1 from P to A)…

94. S: to there. 95. I: so this is point P (on Fig2), this is A, you go from A to P, from P to there, the same distance.

S: yeah, so if you make that (A in Fig1

) the centre, we have both (pointing O and P).

I: no, now you do the same for that point (O on Fig1

). 3 We video-recorded one pair of students for each session. In the second excerpt, which is in the second phase, the students are carrying out the first construction in Cabri.

S: you draw one circle as that…

Towards proving?

This extract is taken from the second phase, immediately after the students have finished the first Cabri construction.

211. S: So now we have to prove it stays exactly the same thing.??????????? Now what do you do?

212. As soon as he finished constructing I starts dragging the first circle to strange positions (Figure 5). Stuart wants to prove; he is at a theoretical level. He seems not to be surprised to see that the treasure stays fixed while moving the apple tree. On the contrary Ian is immersed in Cabri and keeps dragging and observing what happens to the treasure, being surprised (wow!). Stuart seems also disturbed by this continuous dragging, and keeps saying that they have to prove it. He knows that the treasure is fixed, but he wants to know why this is so.

Conclusions

The comparison of the first two excerpts suggests that the Cabri environment proves to be a communication space for students, that is a space in which the interaction between the students is more likely to converge towards a shared understanding than in the paper and pencil environment. Cabri provides "a setting in which the emerging knowledge [of the students] can be expressed, changed and explored … and in which the language that A and B can now use to communicate is the language of the medium" [START_REF] Noss | Windows on Mathematical Meanings[END_REF]. It seems that entering one another universe of thought in Cabri is easier than in paper and pencil. Seeing the dynamic variation of figures on the screen allows interactive participation of both students to the same experience. This dynamic variation is extremely stimulating from the perceptive point of view, but students can interpret this in different ways. In the third excerpt mentioned above, Ian remains at a perceptual level, as he keeps dragging the figure and looking at the variation on the screen, while Stuart immediately thinks of how to prove what he sees on the screen. We believe that Cabri is a good support, because it allows students to do explorations and experiences, which foster the production of conjectures and motivate to proving: proving why a certain proposition holds, within a theory, after they know that it is true, within the Cabri environment (see the second example). Moreover, the fact that when constructing a proof the students exploit their previous exploration in Cabri suggests that this tool may support the construction of a cognitive unity [START_REF] Mariotti | Approaching geometry theorems in contexts: from history and epistemology to cognition[END_REF] between conjectures and proofs.

Developing a Cabri shared experience may constitute a common reference point when constructing the proof. And, generally speaking, constructing a shared classroom experience may be very productive and meaningful. A means to develop this common experience, through the communication in the classroom, can be found in the use of videotapes. In fact, "by videotaping a session and watching the video with the students, we might reduce the risk that significant moments of the classroom life are neglected or forgotten by students over time. The teacher himself may hear voices would remain hidden in a normal classroom and may exploit these recovered voices to discuss problematic issues" [START_REF] Furinghetti | Students approaching proof through conjectures: snapshots in a classroom[END_REF]. From a research point of view, the use of video-tapes allows to carry out deep observation of dynamics and processes in the classroom. The implication of the analysis are both at cognitive level, through the development of frameworks about cognitive processes, and at practical level, through suggestions for classroom practice.

  87. I: a point. Then he walks there. So you have x and x. So that's.. I draws another figure (Fig.2) 88. S: A. 89. I: if you put a circle around there … 90. S: yeah, so what you are saying is that here is A (pointing on Fig1), picking A as a point. 91. I: no.

  98. S: yeah. 99. I: right. Now he walks from O (I draws Fig3 while saying). You're gonna get exactly the same thing, with the circle around again. So that's point O, so that's A.

  Figure 1 Figure 2 Figure 3 In this paper and pencil phase, it is the student who has got the pencil (Ian) to guide the actions, while the other one (Stuart) shows difficulty in understanding and participating in what Ian wants to do. Ian does most of the talking, while Stuart tries to explain to Ian what he sees in the drawings, but he is always wrong: Ian always answers No (#91, #97).

  No you don't do that, you draw a point on that circle and then you draw a circle which goes through that point. Now check whether that point is linked, check whether when you change one the other changes. 206. I moves the first circle 207. S: yeah, so here we go. So this is what you want. 208. I put letters, segments PA and AO, then moves circle then P 209. S: Now you need a perpendicular line and you need a midpoint. 210. I constructs Figure 4.

Figure 4

 4 Figure 4 In Cabri, Ian has got the mouse, but Stuart tells him what to do. Now Stuart does most of the talking. It seems that in Cabri it is easier to enter one another universe of thought. It is not necessary to have the mouse. The one who has got the mouse moves a point, but both of them see the variation of the figure on the screen; this is something dynamic and so it allows interactive participation of both students.

Figure 5

 5 Figure 5

The classroom experiment was carried out by G.Moenne, C.Mogetta, F.Olivero and R.Sutherland from the Graduate School of Education, University of Bristol, in collaboration with Churchill School, Bristol.

The position of the oak and the pine with respect to the apple tree was fixed, so that there was no ambiguity between turning right and left.


Michaël Baudin

Frédéric Coquel

Quang-Huy H Tran

© M Baudin

A RELAXATION SCHEME FOR TWO-PHASE MULTI-COMPONENT FLOWS

Keywords: 1991 Mathematics Subject Classification. 76T10, 76N15, 35L65, 65M06. August 7, 2008 Two-phase flows, multi-component fluids, finite volume schemes, relaxation methods, Whitham condition

Pursuing the program launched our previous works [Numer. Math. 99 (2005), 411-440] and [SIAM J. Sci. Comput. 27 (2005), 914-936], we propose a relaxation scheme for the numerical simulation of one-dimensional two-phase multi-component flows governed by a drift-flux model, the main features of which are a large number of components and a high degree of nonlinearity in the closure laws. In the explicit setting, the relaxation approach allows to ensure positivity for the densities and the mass fractions. The relaxation method is worked out further so as to fit into a hybrid explicitimplicit setting, where fast acoustic waves are treated implicitly to save computational time while slow kinematic waves are treated explicitly in order to maintain accuracy on the transportation of materials.

Introduction

As a sequel to [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], which were primarily concerned with the numerical simulation of a simplistic two-phase flow drift-flux model [START_REF] Masella | Transient simulation of two-phase flows in pipes[END_REF][START_REF] Pauchon | TACITE: A transient tool for multiphase pipeline and well simulation[END_REF] via a relaxation method, this contribution is devoted to a more realistic model, where several species -or components-may appear in each phase. Before elaborating on how the materials of [START_REF] Baudin | Méthodes de relaxation pour la simulation des écoulements polyphasiques dans les conduites pétrolières[END_REF][START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF] can be adapted and enriched by new ingredients, we would like to set out the scope which this project comes within, as well as the problems that await us.

In the transportation of oil and gas mixtures along a pipeline, it is capital to correctly model the mass exchange between the two phases. This cannot be achieved with the gas-liquid model considered in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], which is also referred to as the non-compositional model and where mass balances are written for each phase. Instead, we need a more physically relevant version, where mass balances are written for each specie and which is called multi-component or compositional model.

Due heed must be paid to the number of components involved, since it now equates the number of continuity equations in the model. Traditional techniques such as VFRoe schemes [START_REF] Gallouët | A rough Godunov scheme[END_REF][START_REF] Masella | On an approximate Godunov scheme[END_REF] turn out not to be efficient in terms computational cost, insofar as these require the numerical evaluation of the eigenvalues and the eigenvectors of a large-sized matrix at each edge and at each time-step. The problem is worsened by the high degree of nonlinearity exhibited by the algebraic closure laws. It is well-known [START_REF] Vidal | Thermodynamique : méthodes appliquées au raffinage et au génie chimique[END_REF] that the thermodynamics of a multi-component mixture is far more complex than that of a non-compositional one. At a practical level, the pressure, the partial fractions and the slip velocity are usually the outputs of a lengthy and expensive process consisting of nonlinear inversions and/or of look-up into tables.

In addition to the foregoing aspects, the novelty of which lies in the multi-component nature of the model, we must keep in mind two issues that we already encountered in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]. The first issue is positivity. In [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF], we put forward an explicit relaxation scheme that enables one to guarantee that the total density and the gas mass fraction remain positive from one time-step to another under an adequate tuning of the relaxation parameters. This physically sound requirement meets the need for robustness of the method. For the multicomponent model, positivity means that the total density and the various component mass fractions should remain positive.

The second issue we have to face is the design of a hybrid explicit-implicit time integration, the motivation of which comes from the following observations. On one hand, we want to use large time-steps in order to bring down the CPU cost to an acceptable level. This advocates the use of an implicit strategy. On the other hand, in the context of the flow regimes considered, there co-exist two kinds of waves that are clearly separated by the magnitude of their characteristic speeds: fast acoustic waves and slow kinematic waves. From the petroleum engineer's standpoint, however, only the kinematic ones are of interest since they represent mass transportation. In order to maintain accuracy on the slow waves, we want these slow waves to be treated as if the scheme were explicit. This seemingly paradoxical constraint gives rise to a hybrid time integration that is implicit with respect to fast waves while remaining explicit with respect to slow waves. In [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], we paved the way toward such a selectively implicit time integration, based on a reinterpretation of the explicit relaxation scheme as a Roe method [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference schemes[END_REF] and inspired from a paradigm by Faille and Heintzé [START_REF] Faille | A rough finite volume scheme for modeling two phase flow in a pipeline[END_REF]. Another attempt via the Arbitrary Lagrange-Euler formalism was recently made in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF].

The present work aims at taking up the challenge of the multi-component case by extending the strategies of [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]. The path we suggest differs from the one sketched out in [START_REF] Baudin | A relaxation method via the Born-Infeld system[END_REF]. Nevertheless, it is powerful enough to be valid for any slip law, unlike the latter. This new extension is made possible and tractable thanks to the versatility of the relaxation philosophy [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Chalons | Navier-Stokes equations with several independent pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF][START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF], from which we inherit valuable properties regarding stability, accuracy and simplicity on the basis of a virtual freedom from the closure laws. The feasibility of the extension also relies on a deepened and refreshed understanding of the methods formerly proposed for the gas-liquid model.

This paper is outlined as follows. We start, in §2, by introducing the two-phase multi-component model. The explicit framework of the relaxation scheme is presented in §3, where we focus on the ability to guarantee positivity with the help of only two parameters. The lower-bounds for these parameters are derived in an improved fashion, which is more incisive than in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF]. A thorough discussion on the relaxation parameters is provided in the Appendix. In §4, we turn to the hybrid explicit-implicit scheme and shed new light on the "correct" way to perform relaxation in an implicit context, as proposed by Chalons [START_REF] Chalons | Bilans d'entropie discrets dans l'approximation numérique des chocs non classiques. Application aux équations de Navier-Stockes multi-pression et à quelques systèmes visco-capillaires[END_REF] and ourselves [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]. We also revisit Roe's form of the relaxation scheme via the Osher-Solomon interpretation [START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF]. Numerical results are shown in §5 before we conclude.

Two-phase multi-component drift-flux model

Let K ≥ 2 be the number of components or species, indexed by k ∈ {1, 2, . . . , K}, each of which may appear in the two phases gas g and liquid ℓ. Let ξ k (resp. η k ) be mass fraction of component k in the gas (resp. liquid) phase. We also refer to (ξ k , η k ) ∈ [0, 1] 2 as partial fractions. Naturally, the partial fractions satisfy the consistency relations

ξ 1 + . . . + ξ K = η 1 + . . . + η K = 1. (1)
Furthermore, it is assumed that each component k must be present in at least one phase, that is,

(ξ k , η k ) ≠ (0, 0). (2)

Physical formulation

The isothermal model we consider consists of the K + 1 partial differential equations

∂ t (ρ g R g ξ k + ρ ℓ R ℓ η k ) + ∂ x (ρ g R g ξ k v g + ρ ℓ R ℓ η k v ℓ ) = 0 (3a) ∂ t (ρ g R g v g + ρ ℓ R ℓ v ℓ ) + ∂ x (ρ g R g v 2 g + ρ ℓ R ℓ v 2 ℓ + p) = 0, (3b)
where the first K equations (3a) express the mass balances of the K components, and the last equation (3b) express the total momentum balance of the mixture. The notation R g (resp. R ℓ ) stands for the volumetric fraction of the gas (resp. liquid) phase, and we have

R g + R ℓ = 1. ( 4
)
The symbol ρ g (resp. ρ ℓ ) denotes the gas (resp. liquid) density, assumed to be given functions of the pressure p. It is customary and convenient to define the apparent density of the kth component as

ρ k = ρ g R g ξ k + ρ ℓ R ℓ η k . (5)
The PDE system (3) must be accompanied by two sets of algebraic closure laws. The first set is concerned with thermodynamics and returns

p = p({ρ k } 1≤k≤K ), ξ h = ξ h ({ρ k } 1≤k≤K ), η h = η h ({ρ k } 1≤k≤K ), (6)
for all h ∈ {1, 2, . . . , K}. The second set comprises a single function returning the slip -or drift -velocity

v g -v ℓ = φ({ρ k } 1≤k≤K , ρ g R g v g + ρ ℓ R ℓ v ℓ ), (7)
that is, the difference between the gas velocity and the liquid velocity. Note that, for convenience, we have changed the sign convention from φ = v ℓv g in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF] to φ = v gv ℓ in this paper. Incidentally, the name drift-flux model originates from the choice of a single momentum balance (3b) and from the subsequent necessity to resort to the additional drift law [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] in order to close the system. The slip law encapsulates the knowledge of physicists about flow regimes based on experimental data. We shall not dwell on this first statement of the model, because it is in reality more helpful to work with an alternative formulation that bears more resemblance to gas dynamics.

Mathematical formulation

We define the total density ρ, the gas mass fraction Y , the mass fraction c k of the kth-component and the mass averaged velocity v by the equalities

(a) ρ = ρ g R g + ρ ℓ R ℓ , (b) ρY = ρ g R g , (c) ρc k = ρ k , (d) ρv = ρ g R g v g + ρ ℓ R ℓ v ℓ . (8)
With some abuse of notation, the closure laws ( 6)-( 7) read

p = p(ρ, {ρc k } 1≤k≤K-1 ), ξ h = ξ h (ρ, {ρc k } 1≤k≤K-1 ), η h = η h (ρ, {ρc k } 1≤k≤K-1 ), (9)
which accounts for the thermodynamics, with h ∈ {1, 2, . . . , K}, and

φ = φ(ρ, ρv, {ρc k } 1≤k≤K-1 ), (10)
which accounts for the hydrodynamics. Consider the natural phase space

Ω u = {u = (ρ, ρv, {ρc k } 1≤k≤K-1 ) ∈ R K+1 (11) ρ > 0, v ∈ R, c k ≥ 0, c 1 + . . . + c K-1 ≤ 1}.
We are now in a position to derive an equivalent system for (3).

Lemma 2.1. Weak solutions of (3) are equivalent to weak solution of the system

∂ t (ρ) + ∂ x (ρv) = 0 (12a) ∂ t (ρv) + ∂ x (ρv 2 + P (u)) = 0 (12b) ∂ t (ρc k ) + ∂ x (ρc k v+ σ k (u)) = 0, ( 12c
)
defined for u ∈ Ω u , where the subscript k runs from 1 to K -1 in (12c), and

(a) P (u) = p + ρY (1 -Y )φ 2 , (b) σ k (u) = (ξ k -η k )σ, (c) σ = ρY (1 -Y )φ. (13)
Proof. For smooth solutions, the first equation (12a) results from adding the K continuity equations (3a) together and from the equality

ρ = ρ 1 + . . . + ρ K , (14)
as well as ( 1) and (8d). By inverting the 2 × 2 linear system

Yv g + (1 -Y )v ℓ = v, v g -v ℓ = φ, (15)
we find

v g = v + (1 -Y )φ, v ℓ = v -Yφ, (16)
and substituting in (3a)-(3b) yields (12b)-(12b).

The converse follows the same steps. For weak solutions, equivalence still holds because the only combination we made is a mere sum or subtraction.

From now on, the new system [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] is credited the abstract but concise form

∂ t u + ∂ x f (u) = 0. (17)
At this juncture, two difficulties arise. First, as pointed out in the Introduction, the closure laws ( 9)-( 10) are highly nonlinear and very costly. Put another way, the flux f (u) in ( 17) often does not have any closed-form expression and is expensive. Second, because of nonlinearities, the hyperbolicity property cannot be ascertained for system [START_REF] Gallouët | A rough Godunov scheme[END_REF], i.e., we do not know beforehand whether the Jacobian matrix ∇ u f (u) has real eigenvalues and is R-diagonalizable for u ∈ Ω u , except for trivial instances such as when φ ≡ 0. Numerical computations reveal, however, that hyperbolicity usually holds within the scope of our application. Moreover, for the type of simulations under consideration, the increasingly ordered eigenvalues

µ 1 (u) < µ 2 (u) ≤ . . . ≤ µ K (u) < µ K+1 (u) (18)
always satisfy µ 1 (u)µ K+1 (u) < 0 and

µ 1 (u) ≈ µ K+1 (u)) ≫ µ 2 (u) ≈ . . . ≈ µ K (u) . (19)
We label this phenomenon as the separation of velocity scales: the 2 extreme characteristic speeds µ 1 (u) and µ K+1 (u) correspond to fast acoustic waves, while the K -1 remaining ones are associated with slow kinematic waves representing the actual transportation of the mixture.

Explicit relaxation scheme

Relaxation model

As explained in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF], instead of relaxing every nonlinearity like Jin and Xin [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF], we carry out a partial relaxation procedure in the same vein as Jin and Slemrod [START_REF] Jin | Regularization of the burnett equations via relaxation[END_REF] or Coquel and Perthame [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]. The advantage of this approach is to minimize dissipation by concentrating on the "true" nonlinearities of the system.

In the same spirit as in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF], we define the relaxation state space Ω U as

Ω U = {U = (ρ, ρv, ρΠ, {ρc k } 1≤k≤K-1 , {ρΣ k } 1≤k≤K-1 ) ∈ R 2K+1 (20)
ρ > 0, (v, Π) ∈ R 2 , Σ k ∈ R, c k ≥ 0, c 1 + . . . + c K-1 ≤ 1}.
For later use, let us introduce a few more notations. To every

G = (G 1 , G 2 , G 3 , {G} 4→K+2 , {G} K+3→2K+1 ) ∈ R 2K+1 , (21)
we associate its projected value

▾G = (G 1 , G 2 , {G} 4→K+2 ) ∈ R K+1 . ( 22
)
Since the projection operator ▾ merely deletes the components pertaining to the relaxation variables, it can act on any vector of R 2K+1 , including those which do not belong to Ω U . The subset

€ = {U ∈ Ω U Π = P (▾U) and Σ k = σ k (▾U), 1 ≤ k ≤ K -1} (23)
is called equilibrium manifold, and its elements are said to be at equilibrium.

In the opposite direction, we define

∧u = (ρ, ρv, ρP (u), {ρc k } 1≤k≤K-1 , {ρσ k (u)} 1≤k≤K-1 ) ∈ € (24)
to be the equilibrium extension of u ∈ Ω u . The extension operator ∧ will be of great help at the end of §3 and the beginning of §4.

Definition 3.1. For λ ≥ 0 and a > b > 0, the system

∂ t (ρ) λ + ∂ x (ρv) λ = 0 (25a) ∂ t (ρv) λ + ∂ x (ρv 2 + Π) λ = 0 (25b) ∂ t (ρΠ) λ + ∂ x (ρΠv + a 2 v) λ = λρ[P (▾U λ ) -Π λ ] (25c) ∂ t (ρc k ) λ + ∂ x (ρc k v + Σ k ) λ = 0 (25d) ∂ t (ρΣ k ) λ + ∂ x (ρΣ k v + b 2 c k ) λ = λρ[σ k (▾U λ ) -Σ λ k ], (25e)
defined for U λ ∈ Ω U , is said to be the relaxation model for the original model [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF].

A few remarks are in order. To begin with, this relaxation system is obtained through the same procedure as in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF]. We first switch to Lagrangian coordinates in order to highlight (P, {σ k }) as the nonlinear cores to "get rid of." A relaxation model is then proposed, replacing the nonlinear functions (P, {σ k }) by the new variables (Π, {Σ k }) to which evolution equations are imposed. Reverting back to Eulerian coordinates, we wind up with [START_REF] Masella | Transient simulation of two-phase flows in pipes[END_REF]. Noteworthy is the fact that the latter is split into an acoustic block (25a)-(25c) acting as a "nucleus" around which revolve K -1 kinematic blocks (25d)-(25e). For λ = 0, the blocks are totally unrelated, at least in Lagrangian coordinates, where they have been studied by several authors [START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic BGK models[END_REF][START_REF] Després | Lagrangian systems of conservation laws[END_REF]. For λ > 0, the blocks are coupled through the stiff relaxation terms appearing in right-hand side.

It is of paramount importance to mention that all the K -1 kinematic blocks are assigned the same relaxation parameter b > 0. This fundamental choice is justified not only by our desire for simplicity, but also by our will to be able to ensure positivity for the mass fractions in a predictive manner. Had we opted for a set of componentdependent parameters {b k }, we would have been overwhelmed with troubles. On the side of "simplicity," the solution of the Riemann problem (see §3.3) would have been more painful to write down, because the ordering of the eigenvalues v ± b k ρ -1 is not known in advance. The scheme would have been more expensive and would not have lent itself to a selectively implicit version. On the side of "positivity," no analysis leading to an adequate lower-bound for {b k } such as Theorem 3.3 would have been possible, because a preliminary consistency property stated in Lemma 3.1 would have been violated.

The relaxation system (25) can be given the convenient abstract form

∂ t U λ + ∂ x F(U λ ) = λS(U λ ), (26)
where F and S receive clear definitions. Let us recapitulate the classical properties that we will be using. From now on, we set

τ = ρ -1 (27)
to be the specific volume of the mixture.

Proposition 3.1. The first-order system (25) is hyperbolic over Ω U , i.e., for all U ∈ Ω U , the Jacobian matrix ∇F(U) has real eigenvalues and is R-diagonalizable. The increasingly arranged eigenvalues

v -aτ < v -bτ < v < v + bτ < v + aτ, (28)
where v ±bτ is of multiplicity K -1, all correspond to linearly degenerate fields and are associated with the strong Riemann invariants

Π -av, {Σ k -bc k } 1≤k≤K-1 , Π + a 2 τ, {Σ k + bc k } 1≤k≤K-1 , Π + av. ( 29
)
The right eigenvectors for the eigenvalues (28) is given by the matrix

P(U) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 . . . 0 1 0 . . . 0 1 v -aτ 0 . . . 0 v 0 . . . 0 v + aτ Π + a 2 τ 0 . . . 0 Π 0 . . . 0 Π + a 2 τ c 1 1 . . . 0 c 1 1 . . . 0 c 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . c K-1 0 . . . 1 c K-1 0 . . . 1 c K-1 Σ 1 -b . . . 0 Σ 1 b . . . 0 Σ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . Σ K-1 0 . . . -b Σ K-1 0 . . . b Σ K-1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 30
)
Proof. The calculations are easily adapted from [2, Lemma 3].

The two extreme fields v ± aτ represent fast acoustic waves, while the 2K -1 remaining ones represent slow kinematic waves. The separation of velocity scales is now reflected by a ≫ b, which implies

v -aτ ≈ v + aτ ≫ v -bτ ≈ v + bτ . (31)

Linear asymptotic stability via a Whitham condition

Since no entropy pair is known for the system under consideration, our approach cannot enter the general relaxation theory developed by Liu [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF] and Chen et al. [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. The only way for us to investigate the asymptotic stability of the relaxation system is the framework proposed by Whitham [START_REF] Whitham | Linear and Nonlinear Waves[END_REF] which relies on the Chapman-Enskog expansion in the limit λ → +∞.

Let v = (τ, v, {c k } 1≤k≤K-1 ) ∈ R * + × R × [0, 1] K-1 (32
) be the Lagrangian variables associated to u. Abusing notations once again, we write P = P (v) and σ k = σ k (v) whenever needed. Within this new set of unknowns, we define the partial derivatives P τ , P v , (σ k ) c h . . . the argument of which may be viewed as either v or u depending on the context. Theorem 3.1. At the first order approximation in λ -1 , the projected value u λ = ▾U λ of the solution U λ to the relaxation system (25) satisfies the equivalent equation

∂ t u λ + ∂ x f (u λ ) = λ -1 ∂ x {D a,b (u λ )∂ x u λ } ( 33
)
where the (K + 1) × (K + 1) matrix D(u) has the form

D a,b (u) = 1 ρ 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 . . . 0 × a 2 -A(v) × . . . × × × b 2 -B 1 (v) . . . × . . . . . . . . . . . . . . . × × × . . . b 2 -B K-1 (v) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (34)
in which

A(v) = -P τ + P 2 v + P c1 (σ 1 ) v + . . . + P c K-1 (σ K-1 ) v (35a) B k (v) = P c k (σ k ) v + (σ k ) c1 (σ 1 ) c k + . . . + [(σ k ) c k ] 2 + . . . + (σ k ) c K-1 (σ K-1 ) c k (35b)
and the possibly non-zero entries denoted by × do not depend on (a, b).

Proof. The idea is to plug the Chapman-Enskog expansions

Π λ = P (u λ ) + λ -1 P λ 1 + O(λ -2 ) (36a) Σ λ = σ(u λ ) + λ -1 Σ λ 1 + O(λ -2 ) (36b)
into the relaxation system and to keep the first-order terms. Barring from the opposite sign convention for σ and from the heaviness of the multi-component case, the calculations are similar to those of [2, Proposition 2].

Formally, as λ → +∞, the original model ( 12) is recovered. As for linear stability, the widely and commonly used practice consists in requiring that the matrix D a,b (u) be diffusive, i.e., its eigenvalues have positive real parts. This is what we did in [2, Proposition 1, Lemma 4] to determine the admissible region for (a, b). Here, such a task is out of reach, in view of the size of the matrix.

Instead of the eigenvalues, we are going to impose positivity to the diagonal entries of D a,b (u). This heuristic simplification is strongly supported by our experience with the gas-liquid model, as summarized in [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]Appendix]. It is neither more nor less "rigorous" than the requirement that the eigenvalues of D a,b (u) have positive real parts, which is in itself an approximation. We recall that the exact condition for L 2 stability, obtained by linearizing [START_REF] Zuber | Average volumetric concentration in two-phase flow systems[END_REF] in the neighborhood of any steady state solution u, is that all eigenvalues of the matrix

ıκ∇f (u) -λ -1 κ 2 D a,b (u), κ ∈ R, ı 2 = -1, (37)
should have negative real parts [START_REF] Majda | Stable viscosity matrices for systems of conservation laws[END_REF][START_REF] Serre | Systèmes de Lois de Conservation I & II, Collection Fondations[END_REF].

Definition 3.2. The pair (a, b) ∈ R 2 + , with a > b > 0, is said to satisfy the diagonal Whitham condition if a 2 > -P τ (u) + P 2 v (u) and b > max 1≤k≤K-1 (σ k ) c k (u) (38)
for all u ∈ Ω u under consideration.

Obviously, the double inequality (38) is somewhat different from the positivity of the diagonal entries of D a,b (u). It amounts to merely keeping the framed terms in (35). This choice is corroborated by three observations. First, a close inspection of the magnitudes of the summands in the right-hand side of (35) testifies to the fact that -at least for the numerical data we are working with-the leading terms coincide indeed with the framed ones.

Secondly, from a more theoretical point of view, the inequalities (38) could have been derived by the traditional procedure (requiring positivity for the eigenvalues of the allegedly diffusive matrix) with the additional assumptions that P does not depend on the c k 's and the σ k 's do not depend on (τ, v). The calculations are similar to those of [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]Theorem 2.3,Theorem 2.4]. This corresponds to the situation in which the acoustic block and the kinematic block are uncoupled from each other, whence the name of the condition. The last observation pleading in favor of (38) is that the numerical simulations using (38) work fine, behave stably while giving sharper results than those using the whole right-hand side of (35).

Riemann problem and positivity of intermediate states

As a prerequisite to the full numerical scheme in §3.4, the Riemann problem corresponding to the homogeneous (λ = 0) relaxation system needs to be solved. Besides the structure of the Riemann solution, we lay emphasis on the sufficient conditions for various quantities of the intermediate states to be positive. This will, in turn, be useful for the positivity of the updated variables.

Let U L ∈ Ω U and U R ∈ Ω U be the left and right states defining the initial data

U(t = 0, x) = U L 1 {x<0} + U R 1 {x>0} , (39)
where 1 {.} is the characteristic function. We introduce

(a) v * = v R + v L 2 - Π R -Π L 2a (b) Π * = Π R + Π L 2 -a v R -v L 2 (c) τ * L = v R -v L 2a - Π R -Π L 2a 2 + τ L (d) τ * R = v R -v L 2a + Π R -Π L 2a 2 + τ R (40)
and

(a) c * k = (c k ) R + (c k ) L 2 - (Σ k ) R -(Σ k ) L 2b (b) Σ * k = (Σ k ) R + (Σ k ) L 2 -b (c k ) R -(c k ) L 2 (41)
The above quantities are intended to appear in the formulae for the solution, as shown in the upcoming Proposition. Note that, in (41), the subscript k runs from 1 to K -1.

Proposition 3.2. If a > b > 0 and if a is large enough, e.g., in accordance with (47) of Theorem 3.2, then the solution to the homogeneous relaxation system [START_REF] Masella | Transient simulation of two-phase flows in pipes[END_REF] with the initial data (39) is the self-similar function

U(t, x) = U (x t; U L , U R ) made up of 6 constant states, given by U (x t; U L , U R ) = U L 1 { x t <µ L } + U ⊣ L 1 {µ L < x t <µ * L } + U L 1 {µ * L < x t <µ * } (42) + U ⊩ R 1 {µ * < x t <µ * R } + U ⊢ R 1 {µ * R < x t <µ R } + U R 1 { x t >µ R }
and defined for (t, x) ∈ R * + × R, where the characteristic speeds are

v L -aτ L < v * -bτ * L < v * < v * + bτ * R < v R + aτ R µ L µ * L µ * µ * R µ R ( 43
)
and the consecutively connected intermediate states are

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ * L ρ * L v * ρ * L Π * ρ * L {c k } L ρ * L {Σ k } L ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ↝ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ * L ρ * L v * ρ * L Π * ρ * L {c * k } ρ * L {Σ * k } ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ↝ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ * R ρ * R v * ρ * R Π * ρ * R {c * k } ρ * R {Σ * k } ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ↝ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ * R ρ * R v * ρ * R Π * ρ * R {c k } R ρ * R {Σ k } R ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ U ⊣ L U L U ⊩ R U ⊢ R (44)
Proof. The proof makes use of the Riemann invariants (29) of Proposition 3.1 and follows along the same lines as in [2, Proposition 3].

The Riemann problem can be solved in the most general case including configurations of eigenvalues other than (43). But in view of the separation of velocity scales [START_REF] Vidal | Thermodynamique : méthodes appliquées au raffinage et au génie chimique[END_REF], only the ordering (43) is of interest to us. However, in contrast with [START_REF] Pauchon | TACITE: A transient tool for multiphase pipeline and well simulation[END_REF], the assumption a > b > 0 alone is not sufficient to secure (43). Theorem 3.2. For a > b > 0, we have the ordering

v L -aτ L < v * -bτ * L < v * < v * + bτ * R < v R + aτ R (45)
if and only if the specific volumes of the intermediate states are positive, i.e.,

τ * L > 0 and τ * R > 0. ( 46
)
A sufficient condition for the equivalent conditions ( 45) and (46) to hold is

a > -(v R -v L ) + (v R -v L ) 2 + 8 min(τ L , τ R ) Π R -Π L 4 min(τ L , τ R ) . ( 47
)
Proof. Since the v ± aτ waves are linearly degenerate, we have

v L -aτ L = v * -aτ * L and v R + aτ R = v * + aτ * R . (48)
This can also be directly checked from (40). Therefore, the double inequality

v L -aτ L < v * < v R + aτ R (49)
is equivalent to the positivity property (46).

Inserting the value of v * from (40a) into (49) yields two second-degree polynomial inequalities in a, namely,

℘ L (a) = 2τ L a 2 + (v R -v L )a -(Π R -Π L ) > 0, (50a) ℘ R (a) = 2τ R a 2 + (v R -v L )a + (Π R -Π L ) > 0. ( 50b
)
For a > 0, we have ℘ L (a) ≥ ℘(a) and ℘ R (a) ≥ ℘(a) with the lower-bound

℘(a) = 2 min(τ L , τ R )a 2 + (v R -v L )a -Π R -Π L . (51)
Hence, it suffices to ask for ℘(a) > 0. As a second-degree polynomial, ℘ is convex and has two roots of opposite signs. The positive one is given by (47).

In practice, and as will be seen in §3.4, the left and right states usually belong to the equilibrium manifold €, so that Π L = P (▾U L ) and Π R = P (▾U R ). For short, we write P L and P R for the equilibrium values of Π L and Π R . Imposing the lower-bound

a > a ♯ (u L , u R ) = -(v R -v L ) + (v R -v L ) 2 + 8 min(τ L , τ R ) P R -P L 4 min(τ L , τ R ) , (52)
we may legitimately wonder about a possible logical connection between this lower-bound and the first part of the diagonal Whitham condition (38), i.e., a > -P τ + P 2 v .

(53) In other words, is (52) or (49) implied by some discrete version of (53)? As shown in the Appendix, the answer is in the affirmative for a pure acoustic system, in which P does not depend on (v, {c k }).

We now shift attention to the intermediate fractions c * k . As said earlier, the formulae (41a) are valid only for 1 ≤ k ≤ K -1. It is tempting to define the intermediate mass fraction for the Kth-component as

c * K = 1 -(c * 1 + . . . + c * K-1 ). ( 54
)
The good news is that this definition is consistent with formula (41a) for the Riemann solution.

Lemma 3.1. If U L ∈ € and U R ∈ €, then the intermediate fraction c * K defined by ( 54) is equal to the value obtained by formally putting k = K in (41), that is,

c ⋆ K = (c K ) R + (c K ) L 2 - (σ K ) R -(σ K ) L 2b , (55)
where

(σ K ) L = σ K (▾U L ) = (Σ K ) L and (σ K ) R = σ K (▾U R ) = (Σ K ) R .
Anticipating the proof, we underline that this miracle occurs solely when the same relaxation parameter b is used for every component. As will be clarified in Theorem 3.3, the very possibility of guaranteeing c * k ∈ [0, 1] for 1 ≤ k ≤ K relies on this essential property. This is a compelling argument in favor of a unique b for all components.

Proof. Since the left and right states are at equilibrium, we have Σ k = σ k for 1 ≤ k ≤ K -1. If we sum the equalities (41a) over 1 ≤ k ≤ K -1 and subtract the result to 1, we get

c ⋆ K = (c K ) R + (c K ) L 2 - ς R -ς L 2b , (56)
where

ς = -(σ 1 + . . . + σ K-1 ). (57)
Now, from definition [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF], it is readily verified that

σ 1 + . . . + σ K-1 + σ K = 0 (58)
because of the consistency relations [START_REF] Baudin | Méthodes de relaxation pour la simulation des écoulements polyphasiques dans les conduites pétrolières[END_REF]. As a consequence,

ς = σ K (59)
at the left and right states, which proves the claim.

We are now able to state a sufficient condition for the positivity of the intermediate mass fractions.

Theorem 3.3. If U L ∈ € and U R ∈ €, then a sufficient condition to secure c * k ∈ [0, 1], 1 ≤ k ≤ K, ( 60
)
is that b ≥ b ♯ (u L , u R ) = max { ρφ L , ρφ R } . ( 61
)
The reader's attention is drawn to the fact that the right-hand side of ( 61) is exactly equal to that of [2, Eq. ( 64)], notwithstanding a more complex model. As far as the proof is concerned, the version below is shorter than in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF].

Proof. Since c * 1 + . . . + c * K-1 + c * K = 1,
we are going to seek a sufficient condition to have

c * k ≥ 0, for 1 ≤ k ≤ K. (62)
From (41a) and (55), we infer that

c * k ≥ 0 is equivalent to b ≥ (σ k ) R -(σ k ) L (c k ) R + (c k ) L (63) if (c k ) L > 0 or (c k ) R > 0.
The case (c k ) L = (c k ) R = 0 will be settled later on. By the triangle inequality, we have

(σ k ) R -(σ k ) L (c k ) R + (c k ) L ≤ σ k R + σ k L (c k ) R + (c k ) L (64a) ≤ Y (1 -Y )θ k R + Y (1 -Y )θ k L (c k ) R + (c k ) L × max { ρφ L , ρφ R } (64b) with θ k = ξ k -η k .
(65) Combining (8c), ( 5) and (8b), we obtain

c k = Yξ k + (1 -Y )η k . ( 66
)
This equality can be cast under three different forms, namely,

c k = η k + Y θ k = ξ k -(1 -Y )θ k = γ k + Z k θ k ( 67
)
where the pair

(γ k , Z k ) = (η k , Y ) 1 {θ k ≥0} + (ξ k , 1 -Y ) 1 {θ k <0} (68)
still belongs to [0, 1] 2 . Therefore, dropping out the subscript k momentarily in c k , γ k , Z k , θ k so as to alleviate the notations, we have

Y (1 -Y )θ R + Y (1 -Y )θ L c R + c L = Z R (1 -Z R ) θ R + Z L (1 -Z L ) θ L γ R + Z R θ R + γ L + Z L θ L (69a) ≤ Z R θ R + Z L θ L γ R + Z R θ R + γ L + Z L θ L (69b) ≤ 1. ( 69c
)
Resuming the bounding process (64), we end up with

(σ k ) R -(σ k ) L (c k ) R + (c k ) L ≤ max { ρφ L , ρφ R } , (70)
from which b > b ♯ (u L , u R ) clearly appears to be a sufficient condition for c * k ≥ 0. To settle the case (c k ) L = (c k ) R = 0, we invoke the identity (66) to have

c k = 0 ⇒ Yξ k = (1 -Y )η k = 0. ( 71
)
But by virtue of assumption (2), ξ k and η k cannot vanish simultaneously. Thus, Y = 0 or 1 -Y = 0, and in any case

Y (1 -Y ) = 0. This implies (σ k ) L = (σ k ) R = 0, (72)
Again, as shown in the Appendix, the answer is in the affirmative for the pure scalar case, where σ k depends only on c k .

Numerical scheme and positivity of updated variables

The Godunov flux associated with the relaxation system is classically defined from the Riemann solution as

H(U L , U R ) = F(U (x t = 0 + ; U L , U R )). (74)
From now on, we make heavy use of the operators ▾ and ∧ previously defined in ( 22) and ( 24).

Definition 3.3. Let u L ∈ Ω u and u R ∈ Ω u be two states, and a > b > 0 be the relaxation parameters subject to

a > a ♯ (u L , u R ) and b ≥ b ♯ (u L , u R ) (75)
as in ( 52) and (61). The vector

h(u L , u R ) = ▾H(∧u L , ∧u R ) (76
) is said to be the relaxation numerical flux for the original problem [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF].

The domain is divided into cells of size ∆x. We look for an approximation u n i of u(x i , t n ) at the center x i of cell i and at time t n . Let ∆t = t n+1t n be the time step. The first-order explicit relaxation numerical scheme reads

u n+1 i = u n i - ∆t ∆x [h(u n i , u n i+1 ) -h(u n i-1 , u n i )], (77)
where h(., .) is given by Definition 3.3. Note that (a, b) can be chosen locally: at each edge i + 1 2, there are two relaxation parameters (a i+1 2 , b i+1 2 ) exerting an influence on the numerical flux h(u n i , u n i+1 ). The overall scheme (76)-( 77) can be thought of as the outcome of the following splitting procedure, by means of which we initially described the relaxation method in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF].

(1) Set to equilibrium. Put λ = +∞ and build the equilibrium manifold Π = P (u) and Σ = Σ(u) by applying the extension operator U n i = ∧u n i ;

(2) Relax evolution. Put λ = 0 and solve the relaxation system ∂ t U + ∂ x F(U) = O by the first-order explicit Godunov scheme

U n+1 i = U n i - ∆t ∆x [H(U n i , U n i+1 ) -H(U n i-1 , U n i )] (78)
from which deduce u n+1 i = ▾U n+1 i by projection. The benefit of the splitting interpretation is that it is valid at the continuous level. However, the conciseness of the new formulation (76)-(77) turns out to be a judicious starting point for the design of the selectively implicit scheme in §4. Before getting to this, let us examine the explicit scheme (77). 52) and (61), i.e.,

Theorem 3.4. Assume u n i ∈ Ω u for all i ∈ Z. If • the parameters a n i+1 2 > b n i+1 2 > 0 are in agreement with (
a n i+1 2 > a ♯ (u n i , u n i+1 ) and b n i+1 2 > b ♯ (u n i , u n i+1 ); (79)
• the time-step complies with the CFL condition

∆t ∆x max i max{ v * -aτ * L n i+1 2 , v * + aτ * R n i+1 2 } < 1 2 , (80)
then the first-order explicit relaxation scheme (77) satisfies

ρ n+1 i > 0 and (c k ) n+1 i ∈ [0, 1], 1 ≤ k ≤ K. ( 81
)
Proof. By Theorem 3.2 and Theorem 3.3, we are sure that the intermediate states appearing in the Riemann solution at each edge i + 1 2 satisfy

(τ * L ) i+1 2 > 0, (τ * R ) i+1 2 > 0, and (c * k ) i+1 2 ∈ [0, 1], 1 ≤ k ≤ K. (82)
The proof from this to (81) is similar to that of [2, Corollary 1]. In a nutshell, the positivity of ρ n+1 i stems from the fact that, as a result of the Godunov scheme (78) for the relaxation system, it is equal to the average of a piecewise-constant function at time n + 1, and each of the pieces returning a strictly positive value. As for

(c k ) n+1 i
, it is also equal to an average -albeit in the sense of a positive measure induced by the total densityof a piecewise-constant function bounded by [0, 1].

As explained in [3, §2.3.2], the relaxation parameters are in practice computed as

a n i+1 2 > max{ a ♭ (u n i , u n i+1 ), a ♯ (u n i , u n i+1 ) } (83a) b n i+1 2 ≥ max{ b ♭ (u n i , u n i+1 ), b ♯ (u n i , u n i+1 ) } (83b)
where the additional bounds

a ♭ (u L , u R ) = max -P τ (u L ) + P 2 v (u L ), -P τ (u R ) + P 2 v (u R ) (84a) b ♭ (u L , u R ) = max max 1≤k≤K-1 { (σ k ) c k (u L ) , (σ k ) c k (u R ) } (84b)
are intended to mimic the diagonal Whitham condition (38) at the discrete level. Numerical experiments corroborate to the fact that the reinforced choice ( 83) is a good one. A discussion about possible connections between the ♭ and ♯ bounds is supplied in the Appendix.

Hybrid explicit-implicit relaxation scheme

The design of a hybrid explicit-implicit time-integration in the context of relaxation methods is a delicate matter which requires several ingredients. We start by making right away some general but crucial observations about the fully and linearly implicit time integrations in such a context. Afterwards, we revisit the relaxation numerical flux in order to give it a "differentiable" form via the Roe interpretation. Finally, we propose a strategy for the selectively implicit time-integration.

Interplay between relaxation and implicit time integrations

Let us rewrite the first-order explicit scheme (76)-(77) as

u n+1 i = u n i - ∆t ∆x [▾H(∧u n i , ∧u n i+1 ) -▾H(∧u n i-1 , ∧u n i )]. (85)
Most naturally, the fully implicit counterpart of (85) reads

u n+1 i = u n i - ∆t ∆x [▾H(∧u n+1 i , ∧u n+1 i+1 ) -▾H(∧u n+1 i-1 , ∧u n+1 i )]. (86)
Introducing the enlarged variables

U n+1 i = ∧u n+1 i , U n i = ∧u n i , (87)
which are at equilibrium by construction of ∧, we can easily check that the implicit formulation (86) is equivalent to

▾U n+1 i = ▾U n i - ∆t ∆x [▾H(U n+1 i , U n+1 i+1 ) -▾H(U n+1 i-1 , U n+1 i )] (88a)
U n+1 i ∈ €. (88b)
The equilibrium condition (88b) compensates for the "loss" of information created by the projection operator ▾. In order to express it algebraically, let

⊥U = (ρΠ, {ρΣ k } 1≤k≤K-1 ) ∈ R × R K-1 . ( 89
)
The operator ⊥ is best seen as the map sending any vector of R 2K+1 to the components pertaining to the relaxation variables deleted by the operator ▾. The Cartesian equations for the equilibrium manifold € are abstractly condensed as ⊥U = £(▾U).

(90)

Proposition 4.1. The fully implicit scheme (86) in the original variables u is equivalent to the system

▾U n+1 i = ▾ U n i - ∆t ∆x [H(U n+1 i , U n+1 i+1 ) -H(U n+1 i-1 , U n+1 i )] (91a)
⊥U n+1 i = £(▾U n+1 i ) (91b)
in the enlarged variables U.

Proof. The first line (91a) results from (88a) and the linearity of the operator ▾. The second line (91b) is none other than (88b). This is a most fundamental result, insofar as that it tells us the correct form for the implicit equations in the relaxation variables U. The naive idea

U n+1 i = U n i - ∆t ∆x [H(U n+1 i , U n+1 i+1 ) -H(U n+1 i-1 , U n+1 i )], (92)
even followed by some projection step, leads to an erroneous discretization and disastrous results, as evidenced by Chalons [START_REF] Chalons | Bilans d'entropie discrets dans l'approximation numérique des chocs non classiques. Application aux équations de Navier-Stockes multi-pression et à quelques systèmes visco-capillaires[END_REF][START_REF] Chalons | Well-balanced time implicit formulation of relaxation schemes for the Euler equations[END_REF] for steady-state solutions and by Baudin et al. [START_REF] Baudin | Méthodes de relaxation pour la simulation des écoulements polyphasiques dans les conduites pétrolières[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF] for transient solutions. A partial version of the correct form (91) was suggested by the same authors under the name of implicit projection [3, §3.2.2] and justified by convincing but sophisticated arguments based on the splitting procedure relax evolution-set to equilibrium (see §3.4). To our knowledge, however, equation (86) has never been given a prominent role in the derivation the proposed modifications.

In industrial applications, a handy and cheaper substitute for the full implicit scheme is the so-called linearly implicit strategy. It consists in performing just one Newton iteration on the full implicit system. This amounts to using the Taylor expansion

H(U n+1 i , U n+1 i+1 ) ≈ H(U n i , U n i+1 ) + ∇ L H(U n i , U n i+1 ) ⋅ δU i + ∇ R H(U n i , U n i+1 ) ⋅ δU i+1 (93)
with the increments

δU = U n+1 -U n , (94)
provided of course that H is differentiable. In our case, the Godunov flux H(U L , U R ) is notoriously not differentiable, but we will find a way to circumvent this difficulty in §4.2. Therefore, assuming differentiability, we denote by ∇ L H and ∇ R H the (K + 1) × (K + 1) matrices representing the partial derivatives of H with respect to U L and U R . We introduce

B i = - ∆t ∆x ∇ L H(U n i , U n i+1 ), C i+1 = ∆t ∆x ∇ R H(U n i , U n i+1 ) (95)
and

A i = I 2K+1 -B i -C i . ( 96
)
We also define the Jacobian matrix

Λ i = ∇£(u n i ), (97)
the "explicit" update

U n i = U n i - ∆t ∆x [H(U n i , U n i+1 ) -H(U n i-1 , U n i )], (98)
which is out of equilibrium, and u n i = ▾U n i its projected value. The quotes are simply due to the fact that the time-step ∆t at hand generally does not comply with the usual CFL condition for a "true" explicit scheme.

For any matrix M having 2K +1 rows, we obtain ▾M from M by extracting the rows corresponding to the indices kept by the projector ▾. For any matrix M having 2K + 1 columns, we obtain M ▾ (resp. M ⊥ ) by selecting the columns corresponding to the indices preserved by the projector ▾ (resp. ⊥). The following Proposition gives the correct form for the linearly implicit scheme. Proposition 4.2. The linearly implicit scheme in the context of relaxation takes the form

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾B ▾ i-1 ▾B ⊥ i-1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾A ▾ i ▾A ⊥ i -Λ i I K ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾C ▾ i+1 ▾C ⊥ i+1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾δU i-1 ⊥δU i-1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾δU i ⊥δU i ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾δU i+1 ⊥δU i+1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ▾(U n i -U n i ) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ( 99
)
in the enlarged increments δU. It is algebraically equivalent to the form

▾B ▾ i-1 + ▾B ⊥ i-1 Λ i ▾A ▾ i + ▾A ⊥ i Λ i ▾C ▾ i+1 + ▾C ⊥ i+1 Λ i ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ δu i-1 δu i δu i+1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = u n i -u n i ( 100
)
in the original increments δu. The latter can also be obtained by linearizing (86) directly.

Proof. Plug the approximate value

H(U n+1 i , U n+1 i
) of ( 93), as well as the linear approximation

⊥δU i = Λ i ⋅ ▾δU i (101)
into the generic form (91). The calculations are similar to those in [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF]. Eliminating the ⊥δU's in (99) results in (100). Finally, it it not difficult to verify that a direct linearization of (86) yields (100).

In both formulations, we have to cope with a block-tridiagonal linear system. On the grounds of computational efficiency, it is recommended to implement (100), the size of which is almost twice as small -K + 1 versus 2K + 1-as (99). Nevertheless, the matrices in the enlarged formulation (99) are necessary for the design of a hybrid explicit-implicit time integration. Indeed, the matrices A i , B i and C i will eventually be altered in §4.3.

Roe's form of the relaxation flux

The Godunov flux H(U L , U R ) defined in (74) is only Lipschitz continuous and not differentiable. Moreover, the way it is expressed so far does not help us spot the non-smooth parts. Note that a Roe-or VFRoe-type flux, which is of the form

H(U L , U R ) = 1 2 [F(U L ) + F(U R )] -1 2 A (U L , U R ) (U R -U L ), (102)
would be easier to handle, because non-differentiability is concentrated in the absolute value of the matrix A (U L , U R ). In such a situation, the cure could be to resort to a partial linearization by freezing the nonsmooth parts, as advocated by Mulder and van Leer [START_REF] Mulder | Experiments with implicit upwind methods for the Euler equations[END_REF]. Concretely, we would decide that

∇ L H(U L , U R ) ∶= 1 2 [ ∇F(U L )+ A (U L , U R ) ], (103a) ∇ R H(U L , U R ) ∶= 1 2 [ ∇F(U R )-A (U L , U R ) ]. (103b)
This motivates our search for a Roe-type interpretation (102), the existence of which is ultimately a consequence of the linear degeneracy of all fields in the relaxation model. The success of such an undertaking relies on the following technical Lemma, which was taken for granted in [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], but which needs to be brought to light for the sake of rigor.

-(v R -v L ) + (v R -v L ) 2 + 8 min(τ L , τ R ) Π R -Π L 4 min(τ L , τ R ) , ( 104
)
then the (2K + 1) × (2K + 1) matrix R(U L , U R ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 . . . 0 1 0 . . . 0 1 v L -aτ L 0 . . . 0 v * 0 . . . 0 v R + aτ R Π L + a 2 τ L 0 . . . 0 Π * 0 . . . 0 Π R + a 2 τ R (c 1 ) L 1 . . . 0 c * 1 1 . . . 0 (c 1 ) R . . . . . . . . . . . . . . . . . . . . . . . . . . . (c K-1 ) L 0 . . . 1 c * K-1 0 . . . 1 (c K-1 ) R (Σ 1 ) L -b . . . 0 Σ * 1 b . . . 0 (Σ 1 ) R . . . . . . . . . . . . . . . . . . . . . . . . . . . (Σ K-1 ) L 0 . . . -b Σ * K-1 0 . . . b (Σ K-1 ) R ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 105
)
using the notations ( 40)-( 41), is invertible.

Proof. Let {R q }, 1 ≤ q ≤ 2K + 1, be the column vectors of R, and assume that there is a linear combination such that ∑ 1≤q≤2K+1 α q R q = 0, with α q ∈ R. We claim that

(α 1 , α K+1 , α 2K+1 ) = (0, 0, 0), (106)
since otherwise, the 3

× 3 matrix ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 v L -aτ L v * v R + aτ R Π L + a 2 τ L Π * Π R + a 2 τ R ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (107)
extracted from the first three rows of R, would have to be singular. By subtracting the middle column to the other columns and by expanding with respect to the new first row, we find for the -supposedly equal to zerodeterminant of this 3 × 3 matrix

-2a[v * -(v L -aτ L )][v * -(v R + aτ R )]. (108)
By virtue of Theorem 3.2, however, this quantity is strictly positive if the parameter a meets condition (47).

In turn, the equality (106) implies α 2 = α K+2 = 0 by examining rows 4 and K + 4, columns 2 and K + 2. We proceed likewise until α K = α 2K = 0.

The invertibility of R allows us to state the following Theorem, the most remarkable feature of which is, from a practical point of view, the fact that the matrix A can be computed by a closed formula. 47) is satisfied for all (U L , U R ) under consideration, then there exists a (2K + 1) × (2K + 1) matrix A (U L , U R ) such that the Godunov flux (74) can be expressed as

H(U L , U R ) = 1 2 [F(U L ) + F(U R )] -1 2 A (U L , U R ) (U R -U L ). ( 109
)
This matrix can be given, for instance, by the formula

A (U L , U R ) = R(U L , U R )D(U L , U R )R -1 (U L , U R ), (110)
in which the invertible matrix R(U L , U R ) is defined by (105) of Lemma 4.1, and

D(U L , U R ) = Diag(v L -aτ L , v * -bτ * L , . . . , v * -bτ * L , v * , v * + bτ * R , . . . , v * + bτ * R , v R + aτ R ) (111)
using the notations (40)-(41). Formula (110) defines a Roe-type linearization for the relaxation system [START_REF] Masella | Transient simulation of two-phase flows in pipes[END_REF], in the sense that

1. A (U, U) = ∇F(U); 2. A (U L , U R ) is R-diagonalizable; 3. A (U L , U R )(U R -U L ) = F(U R ) -F(U L ).
Proof. The proof is similar to [3, Proposition 3.2] regarding the claim that A is a Roe matrix. For the readers interested in the details, the starting point of this part is the set of equalities below, which are the counterparts of the jump relations in [3, Lemma 3.1]:

U ⊣ L -U L = [ρ * L -ρ L ]R 1 (112a) U L -U ⊣ L = [ρ * L c * 1 -ρ * L (c 1 ) L ]R 2 + . . . + [ρ * L c * K-1 -ρ * L (c K-1 ) L ]R K (112b) U ⊩ R -U L = [ρ * R -ρ * L ]R K+1 (112c)
U ⊢ R -U ⊩ R = [ρ * R (c 1 ) R -ρ * R c * 1 ]R K+2 + . . . + [ρ * R (c K-1 ) R -ρ * R c * K-1 ]R 2K (112d) U R -U ⊢ R = [ρ * R -ρ R ]R 2K+1 . (112e)
The vectors {R q }, 1 ≤ q ≤ 2K + 1, are the columns of R. The intermediate states

U ⊣ L , U L , U ⊩ L , U ⊢ L are those of the Riemann solution in Proposition 3.2.
As for checking that this matrix gives rise to the flux equality (109), the steps are identical to those in [3, Theorem 3.3]: the formalism by Harten et al. [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] is invoked in to express the flux at issue in two different ways, the half-sum of which finally yields the desired result. In comparison with [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], the only difference lies in the multiplicity of the eigenvalues v *bτ * L and v * + bτ * R , but no serious complication arises in any part. For informative purpose, we wish to briefly report a second way to establish Theorem 4.1 via the Osher-Solomon interpretation [START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF]. This theoretical tool not only makes the checking of (109) simpler, but also brings new insight into the relaxation model. 47) is satisfied for all (U L , U R ) under consideration, then there exists a path Γ = Γ(U L , U R ) ⊂ Ω U such that the Godunov flux (74) can be expressed as

H(U L , U R ) = 1 2 [F(U L ) + F(U R )] -1 2 Γ(U L ,U R ) ∇F(U) dU. ( 113
)
This path, which is made up of 2K + 1 straight line portions

Γ = Γ 1 ∪ (Γ 2 ∪ . . . ∪ Γ K ) ∪ Γ K+1 ∪ (Γ K+2 ∪ . . . ∪ Γ 2K ) ∪ Γ 2K+1 , (114)
each Γ q being directed by the eigenvector R q of the matrix R defined in (105), connects the intermediate states of the Riemann solution in accordance with the diagram

U L Γ1 → U ⊣ L Γ2∪...∪Γ K → U L Γ K+1 → U ⊩ R Γ K+2 ∪...∪Γ 2K → U ⊢ R Γ 2K+1 → U R.
Proof. Except for the claim that every Γ q is a straight line portion, the existence of such a path follows from the linear degeneracy of all fields. The Γ q 's are none other than the integral curves, along which the tangent vector is parallel to the right eigenvector, i.e.,

dU dζ (ζ) = R q (U(ζ)) (115)
for some parametrization ζ ↦ U(ζ) of Γ q . Again, the multiplicity of the eigenvalues v ± bτ is hardly a source of annoyance. This is exemplified by the the jump from U ⊣ L to U L , which is not a pure contact discontinuity because v *aτ * L is not a "simple" eigenvalue, but which can be decomposed into K -1 elementary "simple" jumps.

It comes as an amazing feature of the relaxation model ( 25) that along each portion Γ q , the right eigenvector R q (U(ζ)) is constant, which implies that Γ q is a straight line portion! Indeed, inspection of the matrix of eigenvectors [START_REF] Serre | Systèmes de Lois de Conservation I & II, Collection Fondations[END_REF] reveals that only columns 1, K + 1 and 2K + 1 may have non-constant entries. Using the strong Riemann invariants [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference schemes[END_REF], it can be verified that the components of R 1 are constant across the 1-contact discontinuity. We remind that a strong Riemann invariant for a given field is a weak Riemann invariants for all other fields. The same line of reasoning applies to the two other eigenvectors.

For short, let {ν q }, 1 ≤ q ≤ 2K + 1 be the diagonal entries of D(U L , U R ) defined in (111). As before, R q is the qth column of the matrix R(U L , U R ) given in (105). On one hand, by invariance of the qth-eigenvalue field and qth-eigenvector field along Γ q , we have

Γq ∇H(U) dU = ν q (ζ q+1 -ζ q )R q , (116)
assuming that the segment Γ q corresponds to the interval [ζ q , ζ q+1 ] within the parametrization ζ ↦ U(ζ) for the path Γ. Summing these equalities over q yields

Γ ∇H(U) dU = ν 1 (ζ 2 -ζ 1 )R 1 + . . . + ν 2K+1 (ζ 2K+2 -ζ 2K+1 )R 2K+1 = R D (ζ 2 -ζ 1 , . . . , ζ 2K+2 -ζ 2K+1 ) T (117)
On the other hand, from the plain equality

U R -U L = Γ1 dU + . . . + Γ 2K+1 dU = (ζ 2 -ζ 1 )R 1 + . . . + (ζ 2K+2 -ζ 2K+1 )R 2K+1 (118)
and from the invertibility of R proven in Lemma 4.1, we infer that

(ζ 2 -ζ 1 , . . . , ζ 2K+2 -ζ 2K+1 ) T = R -1 (U R -U L ). (119)
Therefore, combining (117) and (119), we end up with

Γ ∇H(U) dU = R D R -1 (U R -U L ) = RDR -1 (U R -U L ). (120)
This alternative proof of (109) relies on the property that the integral curves are straight lines. Had the eigenvector R q (U(ζ)) not been invariant along Γ q , we would have had to replace R by R, the matrix whose columns are

Rq = 1 ζ q+1 -ζ q ζq+1 ζq R q (U(ζ))dU, (121)
in equations ( 117)-(119). Thus the Roe matrix would have been RD R-1 . However, the invertibility of R is unclear.

From implicit to selectively implicit

After solving the non-differentiability problem for H(U L , U R ), we go back to the Taylor expansions (93). In order to become explicit with respect to kinematic waves, our modus operandi is to prevent the first-order terms

∇ L H(U n i , U n i+1 ) ⋅ δU i + ∇ R H(U n i , U n i+1 ) ⋅ δU i+1 (122)
from capturing any contribution due to slow waves. As initiated by Faille and Heintzé [START_REF] Faille | A rough finite volume scheme for modeling two phase flow in a pipeline[END_REF] and explained in [3, §3.2], we replace the approximations (103) by

∇ L H(U L , U R ) ≈ 1 2 [ ∇F(U L )+ Ã (U L , U R ) ] (123a) ∇ R H(U L , U R ) ≈ 1 2 [ ∇F(U R )-Ã (U L , U R ) ], (123b)
where ̃stands for the acoustic part operator, defined as follows. Let M be a (2K + 1) × (2K + 1) matrix representing any of the matrices ∇F or A , and of course supposed to be R-diagonalizable. From the diagonalized form

M = R Diag(ν 1 , ν 2 , . . . , ν K , ν K+1 , ν K+2 , . . . ν 2K , ν 2K+1 ) R -1 , (124)
we set

M = R Diag(ν 1 , 0, . . . , 0, 0, 0, . . . 0, ν 2K+1 ) R -1 . ( 125
)
We have to be aware of the dependence of the basis of eigenvectors R on the matrix M . This is why this selectively implicit procedure is best seen as a heuristic, although it is exact for linear hyperbolic systems and works fine in nonlinear cases. Anyhow, it allows us to work with a time-step subject to the CFL condition ∆t ∆x max

i∈Z max{ v * -bτ * L n i+1 2 , v * + bτ * R n i+1 2 } < 1 2 (126)
based on the slow characteristic speeds v ± bτ , thus saving CPU time and maintaining accuracy on kinematic waves. To get a good compromise between a large time-step and an acceptable amount of smearing for the acoustic part of the solution, it is advised to take into account the additional safety limitation ∆t ∆x max

i∈Z max{ v L -aτ L n i+1 2 , v R + aτ R n i+1 2 } < 20 (127)
based on fast characteristic speeds v ± aτ .

Numerical results

We carry out two numerical simulations corresponding to K = 2 components, endowed with an idealized but physically meaningful pressure law p = p(ρ, ρc 1 ). This thermodynamic law assumes a homogenized response of the two phase densities

ρ ℓ (p) = ρ 0 ℓ + p -p 0 a 2 ℓ and ρ g (p) = p a 2 g ( 128
)
with respect to the pressure p, but introduces the notion of dew and bubble curves in the (ρ, ρc 1 )-plane. We refer the reader to [START_REF] Bouvier | Algorithmes de résolution compostionnelle dans TACITE[END_REF][START_REF] Vidal | Thermodynamique : méthodes appliquées au raffinage et au génie chimique[END_REF] for details. In the test cases, we set The simulations are Riemann problems over a sufficiently long domain, discretized by the uniform space meshing ∆x = 0.5 m. The discontinuity in the initial data is located at x = 50 m. Actually, we use a secondorder versions of the first-order schemes presented in §3 and §4. The enhancement procedures are the same as in [3, §2.4 & §2.5].

Fast shock with zero-slip law

In the first test case, we consider the zero-slip law φ ≡ 0, which implies that the gas and liquid phases move at the same velocity. The reason why we start by taking such a simplified hydrodynamic law is that we first want to see whether or not the above two-phase two-component thermodynamic law, which is far more sophisticated than that used in [START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], is well-supported by the schemes. Another reason is that for φ ≡ 0, the original system ( 12) is known to be hyperbolic, and an exact solution to the Riemann problem associated with [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] is available for comparison.

The initial data

⎛ ⎜ ⎝ ρ c 1 v ⎞ ⎟ ⎠ L = ⎛ ⎜ ⎝
0.999 × 10 5 0.98 89.569

⎞ ⎟ ⎠ and ⎛ ⎜ ⎝ ρ c 1 v ⎞ ⎟ ⎠ R = ⎛ ⎜ ⎝
0.9 × 10 5 0.98 55.076

⎞ ⎟ ⎠ (130)
are tailored so that this exact solution is a pure acoustic shock propagating at 400 m/s. Here and contrarily to the common practice mentioned in the Introduction, we take a close look at the acoustic wave in order to assess the amount of dissipation. Figure 1 compares the pressure computed by the schemes of this paper to their VFRoe counterparts of the TACITE code [START_REF] Pauchon | TACITE: A transient tool for multiphase pipeline and well simulation[END_REF]. In the explicit setting, relaxation is less diffusive than VFRoe. In the hybrid explicit-implicit setting, relaxation is more diffusive than VFRoe. From this and other numerical runs, the general observation is that the hybrid explicit-implicit version of the relaxation scheme captures acoustic waves with a little more dissipation than VFRoe.

Slow transportation with Zuber-Finlday's law

We now switch to the more realistic slip law

φ(u) = - 0.2 v + 0.35 √ gD 1 + 0.2 κ(u) , and
κ(u) = 1 -Y (u) ρ ρ L (p(u)) -1 , (131)
where g = 9.81 m/s 2 is the gravity constant and D = 0.144 m is the diameter of the pipeline. The hydrodynamic relation ( 131) is known as the Zuber-Findlay law [START_REF] Zuber | Average volumetric concentration in two-phase flow systems[END_REF], and reflects the intermittent flow of medium-sized gas bubbles in a vertical duct. This is a difficult experiment, insofar as the initial data

⎛ ⎜ ⎝ ρ c 1 v ⎞ ⎟ ⎠ L = ⎛ ⎜ ⎝ 0.91827 0.8163265 -10 ⎞ ⎟ ⎠ and ⎛ ⎜ ⎝ ρ c 1 v ⎞ ⎟ ⎠ R = ⎛ ⎜ ⎝ 997.1 0.0816327 -10 ⎞ ⎟ ⎠ (132)
are adjusted so that we cover the whole two-phase domain from Y L = 0 to Y R = 1. Our goal is to put the schemes through the test of robustness at the boundaries of the admissible domain for thermodynamics. Both versions (explicit and selectively implicit) of the VFRoe method in the TACITE code failed this test and stopped after producing a negative mass fraction c 1 .

We do not know the exact solution to the Riemann problem equipped with this combination of closure laws. But on the basis of the analysis by Benzoni-Gavage [START_REF] Benzoni-Gavage | Analyse Numérique des Modèles Hydrodynamiques d' Écoulements Diphasiques Instationnaires dans les Réseaux de Production Pétrolière[END_REF] for the Zuber-Findlay law paired with a much simpler pressure law, and considering that v L = v R and p L ≈ p R ≈ 10 5 Pa, we can legitimately anticipate the "exact" solution as a slow shock propagating at a speed around -10 m/s. This guess is in good a agreement with the curves in Fig. 2, where we display the density and the gas mass fraction.

Conclusion

Once again, relaxation has proved to be a valuable tool in the design of an effective method for the simulation of two-phase flows governed by a drift-flux model. The explicit and the hybrid explicit-implicit versions bring a very satisfactory answer to the difficulties due to the multi-component nature of the mixture, namely: a higher degree of nonlinearity in the closure laws and an impressive increase in the number of unknowns to be updated. In trying to deploy the same strategies as in [START_REF] Baudin | A relaxation method for two-phase flow models with hydrodynamic closure law[END_REF][START_REF] Baudin | A semi-implicit relaxation scheme for modeling two-phase flow in a pipeline[END_REF], we have gained further insight into the theoretical foundations of the relaxation method. We hope these findings could be useful to other researches in the area. with σ ∈ C 1 ([0, 1]; R), by the Jin-Xin relaxation system [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimension[END_REF]

∂ t Y λ + ∂ m Σ λ =0 (134a) ∂ t Σ λ + b 2 ∂ m Y λ =λ[σ(Y λ ) -Σ λ ] (134b)
for λ > 0 and b > 0. By the Chapman-Enskog analysis, it is shown that the asymptotic equilibrium equation for (134) in the limit λ → +∞ is

∂ t Y λ + ∂ m σ(Y λ ) = λ -1 ∂ m {[b 2 -σ 2 Y (Y λ )]∂ m Y λ }, (135)
Y Σ (t = 0, m) = Y L σ(Y L ) 1 {m<0} + Y R σ(Y R ) 1 {m>0} (137)
is the self-similar function

Y Σ m t = Y L σ L 1 { m t <-b} + Y * Σ * 1 {-b< m t <b} + Y R σ R 1 { m t >b} , (138) where
(a) Y * = Y R + Y L 2 - R -σ L 2b (b) Σ * = σ R + σ L 2 -b Y R -Y L 2 (139)
with the shorthand notations σ L = σ(Y L ) and σ R = σ(Y R ). In imitation of (136), we consider

b ♮ (Y L , Y R ) = σ R -σ L Y R -Y L , (140)
which is well-defined for

Y L = Y R . Theorem A.1. Y * ∈ ⌊Y L , Y R ⌉ if and only if b ≥ b ♮ (Y L , Y R ).
Proof. From (139a), it is straightforward to check that

(Y * -Y L )(Y * -Y R ) = (σ R -σ L ) 2 b 2 -(Y R -Y L ) 2 . ( 141
Y n+1 j = Y n j - ∆t ∆m [Σ * (Y n j , Y n j+1 ) -Σ * (Y n j-1 , Y n j )], (142)
then by an averaging argument, we can show that

Y n+1 j ∈ ⌊Y n j-1 , Y n j , Y n j+1 ⌉ (143)
as soon as

b i-1 2 > b ♮ (Y n i-1 , Y n i ), b i+1 2 > b ♮ (Y n i , Y n i+1 ), ∆t ∆m max{b i-1 2 , b i+1 2 } < 1 2 . ( 144
)
In practice, the habit is to define

b ♭ (Y L , Y R ) = max{ σ Y (Y L ) , σ Y (Y R ) }. ( 145
(Y ) = Y (1 -Y )φ(Y ), a sufficient condition to secure Y * ∈ [0, 1] is that b ≥ b ♯ (Y L , Y R ) = max{ φ(Y L ) , φ(Y R ) }. (146)
Proof. By formally setting ρ = 1, the proof follows the same lines as -and is even easier than-Theorem 3.3. Things change, however, when instead of σ(Y ) we have σ(τ, Y, v), as is formally the case of (12c) in Lagrangian coordinates. In such a situation, b ♮ may no longer be finite for Y L = Y R , and even if it is, there is no longer any clear inequality between b ♮ and b ♭ . The advantage of b ♭ over b ♮ is that it remains well-defined, using the partial derivative with respect to Y . This is why we heuristically have to look for b ≥ max{b ♯ , b ♭ }, as indicated in (83b).

Plainly, Y * ∈ ⌊Y L , Y R ⌉ ⇒ Y * ∈ [0, 1

A.2. Euler's isentropic model

Let us approximate the p-system [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]

∂ t τ -∂ m v =0 (147a) ∂ t v+∂ m p(τ )=0, (147b)
where the pressure law is assumed to satisfy p τ < 0 and p τ τ > 0, (148) by the relaxation system

∂ t τ λ -∂ m v λ =0 (149a) ∂ t v λ + ∂ m Π λ =0 (149b) ∂ t Π λ +a 2 ∂ m Π λ =λ[p(τ λ ) -Π λ ] (149c)
for λ > 0 and a > 0. By the Chapman-Enskog analysis, it is shown that the asymptotic equilibrium equation for (149) in the limit λ → +∞ is

∂ t τ λ -∂ m v λ =0 (150a) ∂ t v λ +∂ m p(τ λ )=λ -1 ∂ m {[a 2 + p τ (τ λ )]∂ m v}, (150b)
so that a sufficient condition for linear stability is the Whitham condition

a > -p τ (τ ) (151)
at the continuous level. Introduce w = (τ, v) and let w L ∈ R * + × R and w R ∈ R * + × R. The solution to the Riemann problem of the homogeneous system (149) corresponding to the initial data

w Π (t = 0, m) = w L p(τ L ) 1 {m<0} + w R p(τ R ) 1 {m>0} (152)
is the self-similar function

w Π m t = w L p L 1 { m t <-a} + (τ * L , v * ) Π * 1 {-a< m t <0} + (τ * R , v * ) Π * 1 {0< m t <a} + w R p R 1 { m t >a} , (153)
where

(a) v * = v R + v L 2 - p R -p L 2a (b) Π * = p R + p L 2 -a v R -v L 2 (c) τ * L = v R -v L 2a - p R -p L 2a 2 + τ L (d) τ * R = v R -v L 2a + p R -p L 2a 2 + τ R (154)
with the shorthand notations p L = p(τ L ) and p R = p(τ R ). In imitation of (151), we consider

a ♮ (τ L , τ R ) = - p R -p L τ R -τ L 1 2 (155)
which is well-defined for τ L = τ R .

Proposition A.1. (τ * Rτ * L )(τ Rτ L ) ≥ 0 if and only if a ≥ a ♮ (τ L , τ R ). Proof. This equivalence is a direct consequence of (154c)-(154c).

Hence, for the p-system, the discrete version a > a ♮ (τ L , τ R ) of the continuous Whitham condition a > √ -p τ allows the initial ordering of the specific volumes to be carried over to the intermediate states.

What we want to investigate above all is the positivity of the intermediate volumes τ * L , τ * R and the extent to which this positivity could be implied by a Whitham-like condition. Let us recall a preliminary result which we already encountered in Theorem 3.2.

Proposition A.2. The intermediate specific volumes τ * L , τ * R are positive if and only if we have the ordering

v L -aτ L < v * < v R + aτ R (156)
for the characteristic speeds of the Eulerian version of (149).

Proof. See Theorem 3.2 or directly check using (154).

The main result we wish to put forward is the following.

Theorem A.3.

(1) If v R ≥ v L , then for all a ≥ a ♮ (τ L , τ R ) we have τ * L > 0 and τ * R > 0.

(

) If v R < v L but v R -v L < a ♮ (τ L , τ R )(τ R + τ L ), 2
) then the choice a = a ♮ (τ L , τ R ) ensures τ * L > 0 and τ * R > 0. Condition (157) holds in particular for a subsonic regime where

v L < a ♮ (τ L , τ R )τ L and v R < a ♮ (τ L , τ R )τ R . (158)
If, in addition to (157), the ratio in the densities does not exceed 3, i.e.,

max{τ L , τ R } ≤ 3 min{τ L , τ R }, (159)
then for all a ≥ a ♮ (τ L , τ R ) we have τ * L > 0 and τ * R > 0. Proof. It is convenient to use the notations

Ψ = Ψ R -Ψ L and Ψ = Ψ R + Ψ L 2 (160)
for any quantity Ψ. Plugging the expression (154a) for v * into (156) yields two quadratic inequations

℘ L (a) = 2τ L a 2 + v a -p > 0, (161a) ℘ R (a) = 2τ R a 2 + v a + p > 0. (161b)
On putting a = a ♮ (τ L , τ R ) = (p τ ) 1 2 , we obtain

℘ L (a ♮ ) = ℘ R (a ♮ ) = a ♮ ( v + 2a ♮ τ ). ( 162
)
Case 1. If v ≥ 0, then a = a ♮ is a suitable choice because ℘ L (a ♮ ) = ℘ R (a ♮ ) > 0. Furthermore, both functions ℘ L and ℘ R are increasing with respect to a ∈ R * + . Consequently, every a > a ♮ also satisfies (161). Case 2. If v < 0 but within the smallness constraint (157), we still have ℘ L (a ♮ ) = ℘ R (a ♮ ) > 0, and a = a ♮ remains a suitable choice. To deduce (157) from the subsonic hypothesis (158), we add up the inequalities

v R + a ♮ τ R > 0 and -v L + a ♮ τ L > 0. ( 163
)
The delicate point with v < 0 is that since ℘ L and ℘ R are not increasing with respect to a ∈ R * + , it may happen that as a → +∞, the quantities ℘ L (a) and/or ℘ R (a) becomes negative for a while before turning positive again.

To fix ideas, assume τ ≤ 0. Then, p ≥ 0. The polynomial ℘ L has two real roots of opposite signs. Between the roots, it takes negative values. Therefore, a ♮ must be larger than the positive root. As a consequence, ℘ L (a) > 0 for all a > a ♮ . As far as ℘ R is concerned, it may have two positive roots or may no real roots, according to the sign of the discriminant v 2 -8τ R p . If there is no real roots, the conclusion is ℘ R (a) > 0 for all a > a ♮ . If there are two real roots, their half-sumv 4τ L is also the critical point of ℘ R , where a minimum occurs. Assume that

a ♮ < - v 4τ L . ( 164
)
This implies 4τ L a ♮ + v < 0 < v + 2τ a ♮ , (165) from which it follows that 4τ L < 2τ , whence 3τ R < τ L . This contradicts (159), and establishes that a lies on the right branch

a ♮ > - v 4τ L (166)
of the parabola. The other case τ ≥ 0 is similar.

If the ratio condition (159) is violated, counter-examples can be found so as to invalidate the statement about positivity. We see that the Whitham condition a > a ♮ still has strong connections with the positivity of the intermediate densities, even though it is not geared as neatly as in the scalar model. Like the scalar equation, an averaging argument shows that the positivity of the intermediate specific volumes at each edge implies that of the specific volume at the center of the cells, if we use the first-order explicit scheme

τ n+1 j = τ n j - ∆t ∆m [v * (w n j , w n j+1 ) -v * (w n j-1 , w n j )] (167)
under the CFL condition ∆t ∆m max{a i-1 2 , a i+1 2 } < 1 2 . A direct control of positivity is possible via another lower-bound for a.

Theorem A.4. A sufficient condition to secure τ * L > 0 and τ * R > 0 is that

a > a ♯ (w L , w R ) = -(v R -v L ) + (v R -v L ) 2 + 8 min(τ L , τ R ) p R -p L 4 min(τ L , τ R ) . ( 168
)
Proof. See Theorem 3.2.

For v ≥ 0 or v < 0 subject to (157)-(159), we actually do not to worry about a > a ♯ . The one and only condition to be taken care of is a > a ♮ . Since p is convex, we can even work with the stronger requirement a > a ♭ , where a ♭ (τ L , τ R ) = max{ -p τ (τ L ), -p τ (τ R )}.

(169) Things change, however, when instead of p(τ ) we have P (τ, Y, v), as is formally the case of (12a)-(12b) in Lagrangian coordinates. In such a situation, a ♮ may no longer be finite for τ L = τ R , and even if it is, there is no longer any clear inequality between a ♮ and a ♭ . The advantage of a ♭ over a ♮ is that it remains well-defined, using the partial derivative with respect to τ . This is why we heuristically have to look for a ≥ max{a ♯ , a ♭ }, as indicated in (83a).

  and by applying (41a) or (55), we arrive at the equality c * k = 0. The proof is now completed. Again, the question naturally arises as to whether or not there is a link between condition (61) or the positivity of c * k and some discrete version of the second part of the diagonal Whitham condition (38), i.e., b ≥ (σ k ) c k .

Lemma 4 . 1 .

 41 If the parameter a fulfills condition (47), i.e., a >

Theorem 4 . 1 .

 41 If a > b > 0 and if condition (

Theorem 4 . 2 .

 42 If a > b > 0 and if condition (

p 0 = 10 5

 5 Pa, ρ 0 ℓ = 997.10 kg/m 3 , a g = 330 m/s, a ℓ = 500 m/s. (129)

Figure 1 .

 1 Figure 1. Explicit (left) and hybrid explicit-implicit (right) schemes for experiment 1

Figure 2 .

 2 Figure 2. Density (left) and the gas mass fraction (right) for experiment 2

  so that a sufficient condition for linear stability is the Whitham conditionb > σ Y (Y ) (136) at the continuous level. Let Y L ∈ [0, 1] and Y R ∈ [0, 1]. The solution to the Riemann problem of the homogeneous system (134) corresponding to the initial data

  ) Therefore, the left-hand side is negative if and only if b > b ♮ (Y L , Y R ). Hence, for a scalar conservation law, the discrete version b > b ♮ (Y L , Y R ) of the continuous Whitham condition b > σ Y yields a min-max principle on the intermediate state Y * of the Riemann problem when the left and right data are at equilibrium. This min-max principle on the intermediate state in turn gives rise to the min-max principle on the updated value at the cell-centers. Indeed, if we use the first-order explicit scheme

)

  If σ is a convex or concave function over the interval ⌊Y L , Y R ⌉, then b ♭ ≥ b ♮ and we are justified in choosing b larger than b ♭ (Y L , Y R ). Theorem A.2. For flux functions of the form σ

  ]. Put another way, for a pure scalar conservation law, we do not have to worry about b ≥ b ♯ . The Whitham-like condition b ≥ b ♮ is the one and only condition to be taken care of. By the way, it is likely that b ♮ ≥ b ♯ .

Appendix A. Whitham-like conditions and positivity principles at the discrete level

This Appendix aims at carefully examining the connection between various requirements on the relaxation parameters and their effects at the discrete level for two basic models. It seems to us that the results gathered here in a unifying view are not as widely known as they deserve.

A.1. Scalar conservation law

Let us approximate the scalar hyperbolic equation