TWO APPROACHES TO COMPUTE UNSTEADY COMPRESSIBLE TWO-PHASE FLOW MODELS WITH STIFF RELAXATION TERMS

. The paper deals with the numerical modeling of two-phase flows while using Baer–Nunziato type models. Focus is given here on the numerical treatment of source terms that involve three (or four) relaxation time scales. A new coupled approach relying on the continuous analysis of the system of ODEs is compared with a more widely used strategy grounded on the fractional step approach. Properties of schemes are given in both cases. Several numerical applications show that the coupled approach should be preferred for both stability and accuracy reasons.


Introduction
Many industrial studies urge the development of suitable models and reliable numerical tools, in order to predict two-phase or even multiphase flows.While restricting to two-phase flows, at least two different modelling strategies may be considered.
A first one basically assumes that inner processes involving various relaxation time scales are such that a full instantaneous equilibrium is reached everywhere in the flow.This has led to a rather broad class of so-called homogeneous two-phase flow models, among which we may at least cite [1,4,22,39,47].
When very fast transients are at stake or when droplet atomization occurs, a second strategy, which is grounded on full disequilibrium models, may be retained.The present contribution clearly lies within this framework.More precisely, focus will be given on a class of two-phase flow models that is now well-known, see [3,13,14,25,29,30,32,38,48,55,60] among others.It must be recalled that three (respectively four) relaxation time scales are embedded in these gas-liquid (respectively liquid-vapour, thus including mass transfer) flow models.The latter time scales obviously require suitable closure laws, that can be found in the literature, see [4,7,8,19,28,46,54] among others.These relaxation time scales may be quite distinct and may depend on the application.Assumptions regarding some relaxation time scales may lead to relaxation models such as [23,24,34,51] among others.In the present work, no hypothesis concerning the relaxation time scales is retained, in order to preserve the widest range of applications.
Moreover, the dynamics of the underlying relaxation process in Baer-Nunziato type models has not been thoroughly studied yet.This question is often (implicitly) addressed by seeing the total entropy as a Lyapunov function.In the homogeneous case, i.e. without convective terms, that ensures the stability close to the steady state, or equivalently here, close to the mechanical and thermodynamic equilibrium between the two-phases.However, in many cases, such as vapour explosion or loss-of-coolant accidents, the initial conditions may be set far from thermodynamic equilibrium.Thus, a better understanding of the effective relaxation process is crucial, in order to guarantee the return to equilibrium.In this paper, conditions of effective relaxation are exhibited and discussed for various equations of state (EOS).
Up to now, most of the numerical strategies that have been proposed in order to obtain approximate solutions of these non-equilibrium two-phase flow models, are grounded on the use of a two-step algorithm.The latter includes: (i) an explicit evolution step that treats all convective contributions together with help of Riemann solvers, relaxation solvers or Discontinuous Galerkin methods, see for example [2,14,16,17,25,35,55,58,59,61,63,64], and (ii) an implicit step that deals with source terms associated with the former relaxation time scales, as in [2,17,18,25,26,35,36,52,53,55,58] among others.
It seems worth emphasizing that, in the framework of immiscible three-phase flow models, such as those proposed in [40], the problem of the preservation of admissible states through the convective subsystem (step (i)) has been addressed in Section 2.2.4 of [6], while focusing on simple Stiffened Gas EOS.This result obviously applies in the framework of two-phase flow models considered here.The first explicit step (i) introduces some constraint on the time step for obvious stability issues.This, in turn, renders the implicit treatment of step (ii) mandatory.
Concerning the second implicit step (ii), a fractional step approach has been widely applied up to now, that treats separately the distinct relaxation time scales, see [2,17,18,26,35,55] among others.In addition, the numerical treatment of the source terms, using strong hypothesis on some relaxation time scales, has been investigated in the literature, see for example [36,52,53].However, the fractional step approach used for handling the source terms may suffer from deficiencies, and even lead to a blow-up of the code in some extreme situations, see for example [6] which tackles the problem of vapour explosion.Hence, this motivates to investigate further on the set of coupled ODEs accounting for source terms, and meanwhile derive relevant and more coupled schemes in order to tackle extreme situations.In the present work, a new robust numerical scheme without any assumption on the relaxation time scales is presented.
In the following, emphasis is given on the stability of the scheme.We will also examine whether the accuracy of the new scheme is improved, or not, on a given mesh size, by comparing it with the fractional step approach.
The paper is organised as follows.The governing set of equations is recalled first, and the main properties of the model are given.A few results are then provided and discussed, which concern the true inner relaxation process with respect to temperature, pressure, and velocity in gas-liquid flow models.Afterwards, two distinct algorithms will be considered.The first one is classical.It takes the three (respectively four) relaxation effects into account in a fractional step approach involving three steps when dealing with gas-liquid flows (respectively four in the case of liquid-vapour flows).The second one relies on the investigation of coupled relaxation effects, as discussed in Section 2. Some important properties of the latter schemes are detailed.Both schemes are tested against Chauvin experiment, and a more complex situation arising from the vapour explosion framework.As expected, the comparison of the two schemes is clearly in favour of the coupled algorithm.Appendices complete the paper.In particular, Appendix G discusses the influence of the ratio of pressure and thermal relaxation time scales, whereas Appendix H highlights the influence of the interfacial area.

Governing equations of the two-phase flow model and main properties
We first introduce the set of governing equations of the two-phase flow models examined in the sequel.These include the mass balance equations, the momentum and energy balance equations, together with the evolution of the statistical fractions.We also refer the reader to some companion references [3,13,26,27,29,30,32,34,37,38,48,55,60] that may help and provide additional details.
Both gas-liquid flows (without any mass transfer), and liquid-vapour flows, which involve a unique component (basically water in our framework), will be examined in the sequel.Thus, we introduce rather logical notations as follows.
The statistical fractions of the immiscible liquid phase and the gas/vapour phase are noted   p, q and  , p, q (for gas and vapour respectively).They are such that:   p, q `, p, q " 1. ( In order to ease notations, we will favour the liquid phase (indexed by ) in both cases.All variables in the gas (respectively the vapour) phase will be indexed by  (respectively by ).
Within the -phase, and for  P p, q or  P p, q,   ,   ,   ,   "     and   will respectively denote the phasic mean velocity, mean pressure, mean density, mass fraction and mean total energy, setting: where   p  ,   q denotes the internal energy within phase .The gas-liquid state variable   will be noted: while the liquid-vapour state variable will correspond to: In the sequel, in order to ease notations, the state variable is called  ;  can refer to   or   , depending on the current treated case.If no precision is given, then the results stand true for both cases.

Open set of equations
We may now write the governing set of PDE, which correspond to balance equations for mass, momentum and energy, within phase , and for the statistical fraction   .These are: B  p  q ` p q∇  "   p q; B  p  q `∇ ¨p    q " Γ  p q; B  p    q `∇ ¨p    b   `   ℐ q ´Π p q∇  "    p q; B  p    q `∇ ¨p    p  ` qq `Π p qB  p  q "    p q (5) for  P p, q, or  P p, q, and noting ℐ the identity matrix.
When focusing on gas-liquid flows, we obviously have: whereas for liquid-vapour flows, the following constraint holds: Γ  p q `Γ p q " 0.
For  P p, q, or  P p, q, interfacial transfer terms arising in momentum and energy balance equations, comply with: together with:

Entropy-consistent closure laws
We assume that the interfacial velocity   p q takes the form: with p q P r0, 1s.Hence, it satisfies Galilean invariance.We will specify later on some scalar functions p q that will guarantee unique field by field jump conditions (see also the reference paper [13]).The next definitions require introducing phasic entropies   p  ,   q and temperatures   p  ,   q, which are such that: where: Using the latter definitions, we get the admissible interfacial pressure Π  p q as: where the function p q is obtained straightforwardly: We emphasize that this enables to recover the standard Baer-Nunziato closure Π  p q "   ,   p q "  , .It also guarantees the Realisable Interfacial Pressure condition (see Appendix F and [43]).We recall now some classical results.
Property 1 (Structure of the one-dimensional convective subset).We restrict to the one-dimensional framework.We have the following results: -The homogeneous convective part (left-hand side) of system ( 5) is hyperbolic.The seven real eigenvalues read:  0 p q "   p q;  1 p q "   ´ ;  2 p q "   ;  3 p q "   ` ;  4 p q "  , ´, ;  5 p q "  , ;  6 p q "  , `, .

(16)
Right eigenvectors span the whole space away from the resonance state: -For  P p, q or  P p, q, waves associated with eigenvalues   ˘ are Genuinely Non-Linear, and waves associated with eigenvalues   are Linearly Degenerate.-The coupling wave associated with  0 p q is Linearly Degenerate if: p qp1 ´p qq " 0 (18) or if: -System (5) can be symmetrized away from resonant states (17).
The proof of the first three points can be found in [13], while the fourth one is available in [15].Moreover, details pertaining to Riemann invariants within fields associated with eigenvalues  1´6 can also be found in the latter references.Even more, explicit analytic forms of Riemann invariants in the coupling wave can be found in the particular case p q " 0 in [13,20].Owing to the LD structure of the coupling wave, when p q is chosen in a suitable way, unique jump conditions may be defined in GNL fields for system (5).This was pointed out in [13], and a direct consequence is that the computation of shocks is meaningful in that case, since nonconservative products are well-defined in shock waves (in the non-resonant case): convergent approximations of shocks are not scheme dependent (see [35]).Now, coming back to the three-dimensional setting, we introduce the entropy -entropy flux pair p,   q: " Considering smooth solutions, we can examine the time evolution of  in system (5).Straightforward calculations enable to get: with an explicit form of the right-hand side RHS  p q depending on the closure laws for   p q, Γ  p q,    p q and    p q.
In order to ease the presentation and calculations, we rewrite source terms in a slightly different form and define translated unknowns   p q and   p q such that: p q "    p q ´   p q  p q ´ , 2 Γ  p q.
Taking interfacial constraints on    p q and    p q into account, we must fulfil: p q `, p q " 0, (23) and also:   p q `, p q " 0.
We consider a consistent and Galilean invariant formulation for    p q, hence: where   p q must lie in r0, 1s.We classically note in the sequel: where the free enthalpy ℎ  writes: we obtain the following classical result ( [3,18,25,38,49,55]).
Property 2 (Entropy consistent source terms for a class of two-phase flow models).We consider the following source terms: p q " p q∆,   p q " ´p q∆,   p q " ´p q∆, Γ  p q " ´Λp q∆. ( Then smooth solutions of system (5) agree with: B  pq `∇ ¨pℱ  q " RHS  p q ě 0 (30) providing positive functions p q, p q, p q, Λp q.
Proof.The proof is simple.The right-hand side term RHS  is simply: p qp , ´ q `ˆ1 ´ p q   ` p q  , ˙ p qp , ´ q `ˆ1 ´p q   `p q  , ˙ p qp  ´, q `ˆ ,  , ´   ˙Γ p q (31) and thus positive, owing to (29).
Note that the second, third and fourth closure laws arising in (29) are also entropy-consistent when focusing on single-pressure two-fluid six-equation models [46].Closure laws for heat transfer and drag coefficients arising in p q and p q can be taken from the standard literature (see [46] among others).Besides, references [7,8,28,44] provide closure laws for the pressure relaxation time scales   p q involved in p q.

Relaxation process in a class of two-phase flow models 1
When focusing on gas-liquid flows, where no mass transfer occurs, a straightforward question arises, which concerns the (physically expected) decay of velocity, pressure and temperature gaps.Note that this problem has been investigated recently in [43], while focusing on the sole pressure gaps, though examining several twophase or multiphase flow models.While restricting to the present class of gas-liquid two-phase flow models, the problem at stake here is whether the whole relaxation process is active for the three quantities ∆ , ∆ and ∆ .For that purpose, we consider some homogeneous situation with no gradient of mean variables, thus considering an initial condition for system (5), such that: ∇p,  " 0q " 0 (33) whatever  stands for.Hence, system (5) reduces to: B  p  q "   p q; B  p  q " 0; B  p    q "    p q; B  p    q `Π p qB  p  q "    p q, (34) 1 Consider the vector Δ defined at each point p, q, the components of which are quantity discrepancies between phases, corresponding to states that make the source terms vanish.If the solution of the governing equations, without any convective term, complies with: then the relaxation process is said to be effective.
for  P p, q.This system may be rewritten in a slightly different form: B  p  q "   p q; B  p    q `Π p qB  p  q "    p q; B  p    q "    p q; B  p  q " B  p  q " 0; B  p    `   q " 0; B  p    `   q " 0. ( Equipped with this system (34), we may write the governing equations for the three quantities ∆ , ∆ and ∆ .
If we note: ∆  " p∆, ∆, ∆ q  (36) simple calculations enable to get: where the matrix ℛ  p q P ℳ 3 pRq is given by: Coefficients in matrix ℛ  p q read: together with p q " p q ´ ` `´p q´1 We have used the following notation here for  " , : Thus, we get: Property 3 (Relaxation effects in a class of gas-liquid flow models).We assume that equations of state within each phase are such that, for  " , : Considering positive functions p q, p q, p q, the relaxation process is guaranteed for solutions of (34), if eigenvalues of matrix ℛ  p q are real and positive, or if they are complex with a positive real part.This is guaranteed if the pressure gap ∆ is sufficiently small in the following sense: and if the following condition holds: Proof.The proof is simple and can be found in [42].We briefly recall it below.
-First, note that  "     p q is an obvious real eigenvalue of matrix ℛ  p q, and is positive.-The remaining two eigenvalues of ℛ  p q are the two solutions  ˘of the second-order polynomial: Owing to conditions ( 43) and ( 44), the sum     p q `   p q is positive, thus: ‚ In case of complex eigenvalues, these are complex conjugate, and their real part  is equal to p    p q à   p qq{2, thus positive; ‚ If the two eigenvalues  ˘are real, we have:  `` ´"     p q `   p q, thus the sum is positive.Moreover, the product  `´i s equal to: and we can conclude that both  `and  ´are positive, owing to (45).
For sake of clarity, conditions ( 43)-( 45) are specified and discussed below for various EOS.
To do so, we consider the specific case used in the sequel which is: p q " 0 and thus Π  "   .
-First, we note that condition (43) stands true in most EOS as it corresponds to a specific capacity at constant volume, which is expected to be positive.-For a mixture of perfect gases: Conditions ( 44) and ( 45) always hold true, whatever the state variable  is.-For a mixture of Stiffened-Gases we have: where   ą 1, Π ą 0 and  0 are constants.Admissible states are such that:   `Π  ą 0 and   ´0 ą 0.
Close to thermodynamic equilibrium, equation ( 44) is obviously satisfied.Condition (44) may be rewritten as: Equation (50) shows that, if Π ą Π , condition (44) holds true, even far from thermodynamic equilibrium, in the admissible range:   `Π  ą 0,  P t, u.Condition ( 45) is equivalent to: It must be emphasized that, close to the thermodynamic equilibrium, the right-hand side of (51) behaves as ´p Π ´Π  q 2 {  , and thus condition (45) is satisfied.This is an expected property as the entropy may be understood as a Lyapunov function of our system, owing to the closure of source terms (see Property 2).Condition (51) -and thus (45) -, no longer holds, for any admissible state  , far from equilibrium.-For a mixture of Noble-Abel Stiffened-Gases [50] we have: where   ą 1, Π ą 0,   ą 0 and ε0 are constants.Admissible states are such that:   `Π  ą 0, p1 ´   q ą 0 and   ´ε 0 ą 0. Condition (44) writes as: Condition ( 45) is identical to the one exhibited for a mixture of Stiffened Gases (51).
In practice, for complex EOS, conditions (44) and (45) have to be checked in computer codes, in particular far from the thermodynamic equilibrium.
Remark 1. -We may also note that the threshold effect on ∆ arising in condition (44) has already been pointed out in [6].It arises when taking energy balance into account, and it does not exist when restricting to the barotropic case (see [41]).Its counterpart in the framework of immiscible three-phase flows and miscible two-phase flows is discussed in [43].Straightforward numerical applications show that it can be hardly violated in practice.Thus, the latter constraint is indeed very weak.Actually, we may note that a simpler sufficient condition guarantees that it holds, whatever the values of statistical fractions are, which is: -Moreover, a glance at the precise form of condition (45) shows that the relative amplitudes of time scales   p q and   p q respectively involved in closure laws of p q and p q have no impact on the latter condition.-When accounting for mass transfer terms arising in liquid-vapour flow models, some additional constraints may arise, as emphasized in [42].

Relaxation time scales:
,   ,   and In order to close system (5), the functions p q, p q, Λp q and p q are to be given and thus, four positive relaxation time scales   ,   ,   and   must be introduced.
-Pressure relaxation coefficient: with  0 a positive reference pressure.

Finite volume techniques to compute system (5)
We restrict herein to the computation of unsteady two-phase flows, while applying the two-fluid approach, and considering (5).Before going further on, we rewrite in a quite formal way system (5) as follows: where  denotes the so-called conservative variable introduced in (3) or ( 4).Thus,  p q denotes the conservative flux accounting for convective effects, while p q collects non-conservative contributions arising from the left-hand side of (5).Eventually, p q accounts for the right-hand side source terms in (5).We emphasize that, in the sequel, we assume that shock relations are uniquely defined in model ( 5), and thus we will restrict for our applications to closure laws for the interfacial velocity   p q detailed in Property 1, or other suitable laws that comply with the LD structure of the coupling wave  0 p q "   p q.Of course, this provides a suitable framework in order to verify algorithms (while computing the error), when shocks occur in the flow.
Due to specific applications including water-hammer, loss of coolant accidents, steam explosion and other extreme situations including shock structures and high energy transfers between fluids or phases, robust algorithms are required, and for that purpose, focus is usually given first on low-order time-space Finite Volume schemes (see [21,33]).
In this framework, a simple way to compute approximate solutions of system (5) consists in using the following hybrid implicit/explicit time scheme: p∇ ¨p p qq `p q∇  q  d d (60) setting: ∆  "  `1 ´ , noting volpΩ  q the volume of cell Ω  , and using an explicit approximate Riemann solver associated with the hyperbolic system (evolution step): B  p q `∇ ¨p p qq `p q∇  " 0 (61) in order to get  #  solution of: in such a way that the entropy inequality holds true.Afterwards, the solution  `1  of the relaxation step: must be found within each cell.This strategy has been used by many authors, in order to compute two-phase or multiphase flows (see among others [6,18,25,55]).Actually, it is rather well-suited for transient flows including shock waves, owing to the fact that: -the implicit relaxation step enables to get rid of too heavy constraints linked with upper bounds on the time step; -the explicit scheme involved in the evolution step is in some sense optimal in terms of accuracy, when focusing on fast pressure waves impinging solid structures, and meanwhile it automatically provides the dynamical time stepping.
Obviously, as it occurs in the single phase framework for weakly compressible flows, other numerical strategies should be preferred when aiming at computing low velocity flows, as proposed for instance quite recently in [56].
Before going further on, we would like to recall some suitable algorithms that provide meaningful approximations in the evolution step (61).Among others, note first that the simple Rusanov solver [57] may be applied for that purpose.Moreover, four distinct schemes providing accurate approximations of the homogeneous part of the Baer-Nunziato model are: -an approximate Godunov scheme, which has first been proposed in [61]; -a relaxation scheme that was then introduced in [14,58] and extended in [59] for three-phase flow; -an HLLC approximate Riemann solver, which was proposed in [64]; -another relaxation scheme that was described in [2].
The reader is referred to the paper [16] which provides a detailed comparison of  1 norm of errors for the latter schemes, while focusing on test cases involving difficult Riemann problems.It also seems worth mentioning the recent Discontinuous Galerkin scheme described in [17].Eventually, we point out the recent work [62] dedicated to the analysis of the Riemann problem associated with the Baer-Nunziato model.
For all the simulations run in this paper, a robust Rusanov scheme [57] adapted for the handling of nonconservative products is used.
We will now focus on two different strategies in order to get approximations of the set of coupled ODE: In the sequel, we will consider: "   p q "   Π  p q "  , . (65)

Discrete source terms for gas-liquid flow models
The strategy consists in simulating the drag effects and thermodynamic effects separately.This approach can be justified by the block triangular structure of matrix ℛ  , see (38).Hence, the simulation of the source terms contains two steps: The same time step ∆  is used within each step.The velocity relaxation process is taken from [25] and recalled in Appendix B.
For sake of readability, in Sections 3.1 and 3.2, the state  `1´w ill be referred as   and ∆  will be mentioned as ∆.

Two approaches for the computation of the thermodynamic source terms
First, we recall that for a quantity Ψ, ∆Ψ is set as: Two algorithms are detailed in order to simulate the thermodynamic part of the model (5) in the case of a mix of liquid and gas (i.e.without mass transfer).

A -Fractional step algorithm
A first possible approach in order to account for source terms (34) is to use a fractional step scheme, which decouple (35).System (35), without the velocity relaxation (p q " 0), is simulated in two implicit steps: The time step ∆, given by the evolution step, is used for computing step I and step II.This idea of decoupling system (34) is quite standard and is already used in the literature (see [26,56] for example).

I -The pressure relaxation algorithm
This step is similar but different to the one used in [26].It simulates the solution of the following system: B  p  q " p q∆ ; B  p  q " 0; p P t, uq B  p    q " 0; p P t, uq B  p    `   q " 0   B  p  q ` B  p  q " 0. (67) The equation (67.4) allows writing: Solution of (67) also complies with: Thus: Using the definition of the sound speed, one can derive from (67): setting: It can be noted that the coefficient    is the one arising in ℛ  , see (38).
The pressure relaxation algorithm consists in three steps 2 .

Pressure relaxation algorithm
Step 1. Compute an approximate solution of (71) with an implicit Euler time discretization by considering the frozen coefficient    at time   : Step 2. Using the constraint (70) on the entropy   g "    and the conservation law of the sum of the internal energies (67.4), define: Then compute   l by inverting the immiscible constraint: Step 3. Update the statistical fractions   l and   g : and the total energies: Thus, we have: Property 4 (Pressure relaxation algorithm).
-The pressure relaxation process is effective during step 1 if the pressure gap satisfies condition (44).
-For any mixture of two generalized stiffened gas EOS, the solution   l of step 2, in the admissible range, exists and is unique.Moreover,   l P r0, 1s. Proof.
-If the pressure gap satisfies (44), then the coefficient    is positive and the discrete equation ( 73) is a contraction whatever the time step is.-Let us consider a mixture of two generalized stiffened gases: with Π a positive constant and  0 and  0 constants.Let us recall the constraint for the pressure in the case of a generalized stiffened gas: Thanks to (70), one can obtain: By enforcing (67.4) one can deduce: The immiscible constraint (1) can then be rewritten as follows: Which can be seen as: Standard calculations show that the function  is decreasing and, using (79), that the bounds are: which enables to conclude for   l , such that   l `Π ą 0 and   g `Π ą 0. Thus,   l P r0, 1s, owing to (80) and (82).

II -The temperature relaxation algorithm
It simulates the solutions of the following system: B  p  q " 0 B  p  q " 0 p P t, uq B  p    q " 0 p P t, uq B  p    `   q " 0   B  p  q ` B  p  q " ´p q∆. ( The evolution of the temperature gap ∆ can be deduced from (85), it reads: setting: With    arising in the matrix ℛ  , see (38).
The temperature relaxation algorithm consists in three steps.

Temperature relaxation algorithm
Step 1. Compute an approximate solution of the equation (86) with an implicit Euler time discretization by considering the coefficient    frozen at the time   ˚: Step 2. Since   is constant throughout this step, and considering the conservation law of the sum of the internal energies (85.4),  `1  is solution of: Step 3. Update the thermodynamic quantities  `1  "  `1  ´∆ `1 , and the total energies: We get: Property 5 (Temperature relaxation algorithm).
-Step 1 guaranties the relaxation process for the temperature throughout this algorithm, whatever the time step is.-Considering a mixture of two generalized stiffened gas EOS, the solution  `1  of step 2 in the admissible range exists and is unique.
Proof.-As the coefficient    is positive, the discretized equation of evolution of ∆ (88) is a contraction.
-For a mixture of two generalized stiffened gases, the temperature at time  `1 is:

B -Coupled algorithm
The basic idea of this new algorithm consists in simulating the thermodynamic relaxation effects in one step.The governing set of equations is as follows: B  p  q " 0 p P t, uq B  p    q " 0 p P t, uq B  p    `   q " 0 B  p  q " p q∆   B  p  q ` B  p  q " ´p q∆. (92) Using (92) one can obtain: with: Matrix ℛ   is a sub-matrix of matrix ℛ  (38): The coupled algorithm reads as follows.

Coupled (P-T) relaxation algorithm
Step 1. Compute an approximate solution of system (93), with the evolution matrix ℛ   frozen at time   , using an implicit Euler scheme: where ℐ 2 is the identity matrix in ℳ 2 pRq.
Step 2. Set: Compute  `1  and  `1  , solutions of the following system: Step 3. Update  `1  : and the total energies: Before going further on, we introduce the following lemma: Lemma 1.For a system of the form: B   p q " ´p q p q (101) with  P R  and  P ℳ  pRq,  P N, discretized by using an implicit Euler scheme as: with ℐ  the identity of ℳ  pRq.If the real part of every eigenvalue of   is positive, then system (102) ensures that: and that system (102) is invertible, whatever the time step is.
Proof.As   P ℳ  pRq, its eigenvalues  1 , . . .,   are in C: Then, the eigenvalues   1 , . . .,    of the matrix p  `∆  q ´1 are: Assuming that the real parts of the eigenvalues of   are positive: Whatever the time step is, this ensures that system (102) is invertible and that: Then, one can obtain: Property 6 (Coupled P-T relaxation algorithm).
-If the pressure gap satisfies condition (44), then equation (96) ensures the relaxation process of the pressure and the temperature over time whatever the time step is, if (45) holds.-Given a mixture of two perfect gas EOS, the solution of (98) in the admissible range, exists and is unique.
-If conditions ( 44) and ( 45) are verified, then the real parts of the eigenvalues of ℛ   are positive.Thus, Lemma 1 applies and equation (96) ensures the relaxation process over time, whatever the time step is.-Let us consider a mixture of two perfect gases: Then it can be deduced from the conservation law of the sum of the internal energies (92.3) that: Inserting ( 109) and ( 110) in (98.2), one can obtain: Standard calculations show that the function Θp  ) is decreasing and that its bounds are: Θp  q " ´1 lim   Ñmaxp0,Δ q Θp  q " `8 which enables to conclude for   l , such that   l `Π ą 0 and   g `Π ą 0.Moreover, as  P t, u,   " `Π     p  ´1q , temperatures   and   are in the admissible range.Thus,   l P r0, 1s, owing to (111).

Verification and comparison between the two approaches in a homogeneous case
The basic idea for the two approaches detailed before is to deduce from system (34) an equation of evolution of the gap of the thermodynamics quantities.The main advantage of this idea is that it can ensure the relaxation process between the phases over time, giving conditions that can be easily verified inside a code.The main difference is whether system ( 34) is simulated in one step or "decoupled".In this section, the two approaches are numerically compared.To do so, the flow is supposed to be homogeneous: @Ψp q, B  Ψp q " 0 (115) and velocities within each phase are assumed to be null.This simulation can be viewed as the return to thermodynamics equilibrium of a two-phase flow inside a box.The different simulations are performed with time steps ranging between 10 ´10 s and 10 ´2 s.All relaxation time scales are supposed to be constant.Their values, the initial conditions and the coefficients of the EOS for the different cases are given in Appendix A. Moreover, in order to ensure an effective relaxation time of each separated effect close to the relaxation time scale given by the user, the value of  0 has to be set at  " 0 to: The fractional step algorithm obviously does not capture the correct behaviour of the pressure   at the origin for coarse time steps in some cases (see Fig. 1).Actually, it can lead to a huge overestimation of the pressure, whereas the coupled algorithm better follows the monotony of the exact solution (see Fig. 2).For example, in Figure 1, the pressure   for large time steps can be ten times bigger than the one obtained with the coupled approach.The relaxation process overtime is effective in case 1 as it can be seen in Figures 3 and 4. Here, owing to (50), and since Π ą Π , condition ( 44) is automatically satisfied.Nonetheless, condition (45) must be checked inside the code -see Section 2.3-.Figures 5 and 6 show that the two numerical methods give similar results on the temperature profiles in this case.The main difference between the two numerical approaches here lays in the liquid pressure profiles.
In order to evaluate the performance of these two numerical schemes, a convergence study has been conducted.Let us introduce the error  Ψ p, ∆q of a thermodynamic quantity Ψ for a time step ∆ at time .Given the numerical approximate solution Ψ  of Ψ  , the error  Ψ p, ∆q is calculated as follows: with: In our cases we set   " 2.0 ˆ10 ´3 s.The convergence curves (see Figs. 7 and 8) show a speed of convergence close to 1 for the two approaches, as expected theoretically.It is worth noting that the value of the error can be quite large for the fractional step algorithm, as it can be two hundred times bigger than the error of the coupled algorithm for the pressure   , see Figure 8.The plateau of convergence for large time steps corresponds to simulations where ∆t ą t t .
Moreover, for a given value of ∆, the ratio of computational costs is between 2 and 8 in favour of the coupled algorithm, depending on the mesh size, see Table 1.       Figure 9. Sketch of the data settings in meters.

Application: shock wave through a two-phase gas-liquid medium [11]
A -Experimental set up and numerical settings This part aims at validating the coupled scheme by comparing its results with the experiment [11] (see [10] for further details).The experimental set-up is made up of a one-dimensional shock tube composed of air and a layer of liquid droplets.It contains, at  " 0 (see Fig. 9 for a sketch of this set up), a high pressure chamber (HP), followed by a zone at atmospheric pressure.Inside the low pressure chamber, a cloud of droplets (CD) with a unique diameter (  " 500 m) lies between  " 3.0 m and  " 3.40 m.
The initial conditions are set as follows [11]: In the experimental set-up, several pressure probes are set throughout the tube [11].They give the value of the total pressure.The total pressure  mix is defined as: Indeed, setting: and using ( 5), we have: Integrating (121) over a closed domain shows that the pressure applied on the walls is precisely  mix , i.e. the total pressure measured experimentally.The value of  mix is compared with the experimental data on three probes: - 1 located at  " 1.77 m.This station is located inside the initial gaseous area.
- 2 located at  " 3.08 m.This station is located at the beginning of the initial layer of droplets.
- 3 located at  " 3.19 m.This station is located in the middle of the initial cloud of droplets.
In this configuration, the fragmentation of the water droplets plays a substantial role in the dynamic, according to [31].That is why an interfacial area  is introduced for the liquid phase: The equation of evolution of the interfacial area (see Appendix C) and its numerical treatment are taken from [5].The velocity relaxation process has to be computed too.To do so, the algorithm detailed in [25]  This expression of   is derived from the Stokes formula [46].Pressure relaxation time scale: where   " 1.8 ˆ10 ´5 kg m ´1 s ´1 is the dynamic viscosity of the air at 1 bar and 293 K.It is the limit of the closure law proposed in [28] for small diameter droplets.Temperature relaxation time scale: where   " 10 is the Nusselt number and   " 0.6 (W m ´1 K ´1q is the thermal conductivity of the gas phase.This form is taken from [54].
The simulations have been performed using one-dimensional meshes including 1000, 10 000 and 20 000 cells.The CFL number is set to 0.45 and provides the time stepping.

B -Numerical results and comparison with experimental data
We recall that the Weber number,  , is a dimensionless number that quantifies the atomization of the droplets.Introducing   , a positive reference surface tension,   is defined as: More precisely, if the Weber number is smaller than a certain threshold called the critical Weber number    , then the atomization does not occur (see Appendix C).Two critical Weber numbers will be tested in the simulations:    " 3 and    " 12.These values, respectively suggested for liquid water and liquid aluminium droplets, are taken from [54].The total pressure  mix obtained through numerical simulations is compared with the experimental data.First, let us note that on Figures 10-12, the experimental data have been translated of a fixed value  translate " 0.0019 s, so that the experimental and simulated incident shock waves are synchronized on station 3.As the   experimental and numerical shock waves arrive almost at the same time at the station 1 and 2, it can be concluded that the simulation captures well the celerity of the shock wave.On fine meshes, Figures 10 and 11 show a pressure drop right after the incident shock wave of 2.4 bar.The amplitude of the incident shock wave is overestimated by the simulation by 10% at station 3 (Fig. 10), whereas there is almost no discrepancy between the simulation and the experimental data at stations 1 and 2. The pressure loss right after the shock wave is expected as it accounts for the atomization of the water droplets, see [31].A coarse mesh with 1000 cells is not sufficient to properly capture this phenomenon.In Figure 10, after the pressure loss, the total pressure increases and then plateaus around 3.5 bar, considering    " 3 and a 20 000 cells mesh.Focusing on this mesh, the numerical results overestimate the experimental data on  mix by about 7 % on station 3, 10% on station 2 and 2% on station 1.Eventually, the total pressure decreases due to the arrival of the reflected rarefaction wave coming from the left wall boundary.It is also worth noting that ∆ remains lower than 2 ˆ10 ´3% of the initial phasic pressure throughout the simulation.On the other hand, ∆ remains lower than 2% of the initial phasic temperature gap.
The numerical simulation of Chauvin experiment [11] using the fractional step approach can be found in Section 4.3 of [6].On a fine mesh, discrepancies between the two approaches on the total pressure profile are less than 1% on the incident shock wave amplitude, and about 5% on the pressure plateau.
In addition, these numerical simulations also show the influence of the critical Weber number -i.e.how the interfacial area is taken into account -on the overall behaviour of the solution after the shock wave, inside the cloud of collapsible droplets.Indeed, as shown in [12], an accurate modelling of the inter-facial area evolution is a key ingredient in order to capture the correct behaviour of the pressure profile after the shock.Herein, computational results are closer to the experimental data when using    " 3 on a refined mesh, as anticipated.
According to Figure 13, the velocity of the liquid phase (labelled one) and the gas phase (labelled two) are completely distinct, with a maximum gap of more than 200 m s ´1, right after the incident shock wave.The structural hypothesis of velocity disequilibrium between phases is thus retrieved in the numerical results.Figure 14 shows the impact of the value of the critical Weber number on the atomization process.Indeed, a  four times bigger critical Weber number results here in a twice larger droplet's diameter on the first plateau after the shock.
In Appendix H, the Chauvin experiment [11] is also computed, assuming rigid particles, as conducted by authors in [12].In that case, as the diameter of the droplets is constant, no pressure drop arises right after the incident shock wave, as expected (see Fig. 3 of [12], and also [6]), even for a refined mesh.

Discrete source terms for a mixture of liquid water and vapour
As in the liquid-gas section, the velocity relaxation is taken into account before the thermodynamic part of the source terms.

Two algorithms for the simulation of the thermodynamic source terms
We recall that for Ψ P t, ,  u, ∆Ψ is set as: And for clarity, we note: A -Fractional step approach For a mixture of water and vapour, the relaxation effects, which are presented for a mixture of liquid and gas, must be complemented with the mass transfer terms.This fractional step scheme is thus composed of three steps, one for each thermodynamic effect.The sequence is as follows, still using the same time step ∆ within each step: The first two steps are identical to the ones presented for a mixture of liquid and gas in Section 3.1.1.A. Besides, step III is taken from [18] and recalled in Appendix D.

B -Coupled approach
Once more, the basic idea of this approach is to account for all the thermodynamic effects in one step.First, the global system of source terms is recalled: B  p  q " p q∆ B  p  q " Γ  p q B  p    q "   ` 2 Γ  p q B  p    q ` B  p  q "    2 Γ  p q ´p q∆ B  p  ` q " 0 B  p    `   q " 0 B  p    `   q " 0. (129) From system (129), it can be obtained: B  p  q " p q∆ B  p  q " Γ  p q B  p    q "   ` 2 Γ  p q B  p    q ` B  p  q " ´p q∆ B  p  ` q " 0 B  p    `   q " 0 B  p    `   q " 0.
(130) It can also be deduced from system (130) that: since: System (130) can be recast in two subsystems, which are: and: " B  p    q "   `
From system (133), an equation of evolution for the thermodynamic quantity gaps can be derived, see [42]: The coefficients of ℛ    P ℳ 3 pRq are detailed in [42] and recalled in Appendix E. Using an implicit Euler method, system (135) is discretized as follows: The conservation laws included in (133) imply: and the one included in (134) gives: Coupled (P-T, ) relaxation algorithm Step 1. Compute ∆ `1 , ∆ `1 and ∆ Update  `1  and  `1  in agreement with (140).
Step 5. Then compute  `1  and  `1  solutions of: Step 6. Update the total energies as: Remark 4. -Equation ( 139) is deduced from the mass conservation (133.2), with ∆ frozen at time  `1 .
-System (140) is derived from the conservation law of the sum of the internal energies (131) and the immiscible constraint (1).-System (142) is derived from system (134) by using an implicit Euler scheme.Thus, we have: Proof.-Standard calculations give the determinant  of system (142): As  is strictly positive whatever the time step is, system (142) is always invertible.-System (140) is identical to (98), thus it benefits from the same properties.
-If the three fundamentals minors of ℛ    (135) are positive, then according to [42], Lemma 1 applies and ensures the relaxation process over time, whatever the time step is.

Verification and comparison between the two approaches in a homogeneous case
In this section, the two schemes presented in the previous part are compared on different test cases.To do so, the flow is supposed to be homogeneous: @Ψp q, B  Ψp q " 0 (145) and the velocities within each phase are assumed to be null.The different simulations are performed with time steps ranging between 10 ´10 s and 10 ´1 s.All the relaxation time scales are supposed to be constant.Their values, the initial conditions and the coefficients of the EOS for the different cases are given in Appendix A. Moreover, in order to ensure an effective relaxation time of each separated effect close to the relaxation time scale given by the user, the value of  0 and Γ 0 have to be set at  " 0 to: As it occurs in the liquid-gas framework, the fractional step algorithm can in some cases overestimate the pressure of the liquid phase for coarse time steps.For example, in case 3, (see Fig. 15), using a time step of 10 ´3 s, the latter method leads to pressure values almost one hundred times bigger than the converged solution.
For the same time stepping, the coupled algorithm overestimation is much lower (see Fig. 16).Even if it is not the case here, the overestimation of the pressure   by the fractional step algorithm could lead to the violation of the relaxation conditions detailed in [42].In addition, when case 3 is computed with a time step ∆ " 10 ´3 s, the fractional step algorithm is unable to return to pressure equilibrium, and even stops just before  " 0.02 s (see Figs. 17 and 18).The temperature profiles of the liquid water and vapour, computed with the fractional step algorithm and the coupled algorithm, are given in Figures 19 and 20 respectively.A convergence study is presented below.The error  Ψ p, ∆q is defined as in the liquid gas section (117).The speed of convergence is close to one (see Fig. 21) for the two approaches, as expected.Note that the value of the error for the fractional step algorithm is about one hundred times bigger for the pressure   than the one obtained with the coupled algorithm.
Besides, the ratio of computational costs is in the range r5, 12s in favour of the coupled algorithm, for a given ∆, see Table 2.Moreover, as mentioned before, the fractional step algorithm even blows up on coarse meshes.

Application: One dimensional shock tube in vapour with a cloud of water droplets
The aim of this section is to test the coupled algorithm for a mixture of liquid water and its vapour.The application test case is composed of a one dimensional shock tube in vapour.At  " 0, a cloud of water droplets lies inside a portion of the tube.The initial conditions along the tube are set as follows:     Three probes are set inside the tube: - 1 is located inside the initial vapour layer of the tube:  " 1.77 m.
- 2 is located at the beginning of the initial cloud of water droplets:  " 3.08 m.
- 3 is located in the middle of the initial zone of liquid -vapour mixture:  " 3.20 m.
Relaxation time scales and EOS coefficients are the same as those used for the convergence curves in Figure 21 and can be found in Appendix A, case 2. In this test case, the particles are supposed to be rigid.
As expected, the solution obtained in the homogeneous case is retrieved at station 3 in the first time steps, see Figures 22-24.Moreover, in this case, the pressure disequilibrium is substantial as ∆ can be four times larger than the initial pressure, see Figure 23.The temperature disequilibrium is also significant, as it is still about 70% of its initial value after the shock wave.The difference of velocities ∆ is also quite large, close to 600 m s ´1 right after the shock wave.Remark 5. Eventually, it is also worth noting that the fractional step algorithm fails to compute an approximate solution of this test case.

Conclusion
When focusing on the class of liquid-gas flow models considered here, conditions ( 43)-( 45) are the true conditions that guarantee the effective relaxation process.For most EOS, conditions (44) and (45) are not trivial and have to be numerically tested inside the code.In most cases, conditions (44) and ( 45) also involve initial conditions.For liquid-vapour flows, four conditions arise, which also involve EOS and initial conditions.

Liquid-vapour
See Tables A. 4 and A.5. B  p    q " ´p q   ` 2 ∆ B  p  ` q " 0 B  p    `   q " 0 B  p    `   q " 0.

Velocity relaxation algorithm
Step 1. Compute  `1ĺ  and  From (D.1), a conservation law of the internal energies is obtained during this step: " B  p    q " 0 B  p    q " 0 (D.3) which implies that: Γ  " Γ p  q. (D.4) The following algorithm is taken from [18].A detailed version in French is also available in [45].

Mass transfer algorithm
Step 1. Compute  `1 l , solution of: (D.7) Appendix E. Coefficients of the matrix ℛ    Matrix ℛ    p q P  3ˆ3 is as follow, see [42]: ℛ    p q " ¨   p q     p q     p q     p q     p q     p q    p q    p q    p q ‹ ‚. (E.1) The above mentioned coefficients are as follows: p q " p q ´ ` ´Δ (E.7) The relaxation process over time is guaranteed if the three fundamentals minors are positive [42].Chauvin experiment [11] with rigid particles of diameter   " 500 m.Evolution of  mix at station 3 using a 10 000 cells mesh and the coupled algorithm.
Appendix H. Chauvin experiment [11] simulated with rigid particles Here, Chauvin experiment [11] is simulated with rigid particles assuming constant interfacial area, associated with constant particle diameter:   " 500 m.The initial conditions and thermodynamic coefficients are recalled in Appendix A. The relaxation time scales are the same as those used in Section 3.1.3.
We recall that for collapsible droplets, as in the experiment, a pressure drop occurs right after the incident shock wave, as expected when droplet atomization arises, see [31].The simulation captures this phenomenon quite well, see Figure 10.
On the other hand, Figure H.1 shows that with rigid particles, the simulation exhibits a slow increase of the total pressure after the incident shock wave.Similar results have already been exhibited in [6,12].

Figure 1 .
Figure 1.Homogeneous case: evolution of the pressure   in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the fractional step algorithm.The initial conditions are given in Appendix A.

Figure 2 .
Figure 2. Homogeneous case: evolution of the pressure   in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the coupled algorithm.The initial conditions are given in Appendix A.

Figure 3 .
Figure 3. Homogeneous case: evolution of ∆ "   ´ in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the fractional step algorithm.The initial conditions are given in Appendix A.

Figure 4 .
Figure 4. Homogeneous case: evolution of ∆ "   ´ in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the coupled algorithm.The initial conditions are given in Appendix A.

Figure 5 .
Figure 5. Homogeneous case: evolution of   and   in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the fractional step algorithm.The initial conditions are given in Appendix A.

Figure 6 .
Figure 6.Homogeneous case: evolution of   and   in case 1 (  " 10 ´5 s and   " 10 ´3 s), computed with the coupled algorithm.The initial conditions are given in Appendix A.

Figure 7 .
Figure 7. Homogeneous case: convergence curve for case 1.The initial conditions are given in Appendix A.

Figure 8 .
Figure 8. Homogeneous case: convergence curve for case 5.The initial conditions are given in Appendix A.

Figure 11 .
Figure 11.Chauvin experiment: evolution of the pressure  mix (Pa) at station 2 for various meshes and    .

Figure 12 .
Figure 12.Chauvin experiment: evolution of the pressure  mix (Pa) at station 1 for various meshes.

Figure 13 .
Figure 13.Chauvin experiment: evolution of the velocities at station 3 for a mesh of 20 000 cells.Index number 1 corresponds to the liquid phase and index number 2 corresponds to the gaseous phase.

Figure 14 .
Figure 14.Chauvin experiment: evolution of the diameter   at station 3 for a mesh of 20 000 cells.

Property 7 (
Coupled P-T- algorithm).-System (142) is invertible whatever the time step is.-If the three fundamentals minors of matrix ℛ    (135) are positive, then step 1 ensures the thermodynamic relaxation process over time, whatever the time step is.-For a given couple of perfect gas EOS, solutions  `1  and  `1  of (140) in the admissible range, exist and are unique.Moreover,  `1  belongs to r0, 1s and the partial masses  `1  and  `1  remain positive.

Figure 15 .
Figure15.Evolution of   in case 3 computed with the fractional step algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 16 .
Figure16.Evolution of   in case 3 computed with the coupled algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 17 .
Figure17.Evolution of ∆ in case 3 computed with the fractional step algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 18 .
Figure 18.Evolution of ∆ in case 3 computed with the coupled algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 19 .
Figure19.Evolution of   and   in case 3 computed with the fractional step algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 20 .
Figure 20.Evolution of   and   in case 3 computed with the coupled algorithm.  " 10 ´5 s,   " 10 ´4 s,   " 10 ´6 s.The initial conditions are given in Appendix A.

Figure 22 .
Figure 22.Evolution of  mix for two meshes including 1000 cells and 10 000 cells.

Figure 23 .
Figure 23.Evolution of   and   at station 3 with the initial condition of case 2. A mesh of 10 000 cells has been used.

Figure H. 1 .
Figure H.1.Chauvin experiment[11] with rigid particles of diameter   " 500 m.Evolution of  mix at station 3 using a 10 000 cells mesh and the coupled algorithm.
p  p  ,   qq ˙{pB   p  p  ,   qqq   p  p  ,   qq B   p  p  ,   qq ¨(13)

Table 1 .
CPU normalised by the CPU of the fractional step algorithm using ∆ " 10 ´8 s for case 1.
´2 10 ´8 s 1.21 ˆ10 ´1 1 has been used and is recalled in Appendix B. In order to simulate this shock tube, stiffened gas EOS (see (78)) are used in each phase.The values of the thermodynamic coefficients are given in Appendix A.Moreover, the relaxation time scales   ,   and   are chosen as follows.
`1 , solution of system (136).The matrix ℛ    is considered frozen at time   .

Table 2 .
CPU normalised by the CPU of the fractional step algorithm using ∆ " 10 ´10 s for case 3.

Table A . 1 .
Simulation parameters for the homogeneous liquid gas cases and initial conditions.Table A.2. EOS for the homogeneous liquid gas cases.Table A.3.EOS for the Chauvin experiment [10].
Table A.4. Numerical parameters for the liquid gas simulations and initial conditions.Table A.5. EOS for all the liquid vapour simulations.   " 0 B    " 0   B    " ´p q∆ 1 l "   l ´∆Γ  and update  `1 v "   l ` v ´`1 l .Step 2. Compute  `1 l and  `1 v as: l ¯(D.5)˚" p    q l ¯`1  `1 l ` `1 v 2 p    q `1 ˚" p    q  ˚`p    q  ˚´p    q ˚" p    q l ¯`1  `1 l  `1 v 2 p    q `1 ˚" p    q  ˚`p    q  ˚´p    q `1 ˚.
B   p  q|  ¯,     p q " p q ´1   B   p  q|  `1 B  pq|  ´´  ` B  p  q|    B   p  q|  ´`B pq|  B  pq|   B  p  q|    B   p  q|  `´ 2  B  pq|  B  pq|  ¯,     p q " p q ´1   B   p  q|  `1 B  pq|  ´´  ´ B  p  q|    B   p  q|  `´´B  pq|  B  pq|    " ˆ1   B   p  q|   ´ℎ   B   p  q|   ˙(E.4)   " ℎ    B   p  q|   p  ` B   p  q|   q ´1   B   p  q|   p  ` B   p  q|   q (E.5)   " ℎ    B   p  q|   `´  `2  B   p  q|   ˘´1   B   p  q|   `´  `2  B   p  q|   ˘(E.6) ¯,   p q " Λp q  ¯.