PRESSURE RELAXATION IN SOME MULTIPHASE FLOW MODELS

. We consider in this paper multiphase flow models involving two, three or four fields, in total mechanical and thermodynamical disequilibrium. Thus several pressure fields arise, and we precisely focus here on the pressure relaxation process, while restricting to four distinct multiphase flow models. The first two models only involve immiscible compressible components, while the last two hybrid models involve both miscible and immiscible components. It is shown that some -weak- restrictions may occur on pressure gaps, which are unlikely to appear in practice. Evenmore, three-phase flow models may involve a non-monotone behaviour in the return to pressure equilibrium.


INTRODUCTION
We examine herein the pressure relaxation process in some two-phase and three-phase flow models involving either immiscible compressible components, or a mixture of miscible and immiscible components.
Actually, the main problem at stake here is whether the pressure relaxation process is guaranteed by solutions of the latter models, or in other words, whether the pressure gaps between fields/phases decay in time, or not. According to authors, this problem has seldomly been adressed in the multiphase flow literature, and the present contribution aims at giving some better understanding of this expected behaviour. The four models discussed in the sequel were introduced in [10, 26, 28-30, 33, 34], and more details concerning the derivation and the main properties of the associated PDEs can be found in the latter references. All these models comply with the following basic specifications: • A physically relevant entropy inequality holds for smooth solutions of PDEs; • Jump conditions are uniquely defined, field by field; • The governing set of PDEs can be symmetrized, even in the multi-dimensional framework, which of course implies that the convective subset is hyperbolic. Before going further on, we also refer the reader to [15-18, 20, 22], among other references, which give another insight of the whole modeling concept, while retaining the powerful tools relying on Hamilton's principle. In particular, while focusing on two-phase flows involving a cloud of bubbles in a liquid, the early paper [20] provides a very promising approach. Bridges between the two approaches still remain to be built.
The paper is organized as follows. In section 2, a class of gas-liquid and liquid-vapor two-phase flow models is examined. These are of the Baer-Nunziato (BN) type (see [3,10]), and Property 1 will give some first conditions pertaining to initial conditions on the pressure gap in order to fulfill the pressure relaxation process. In the following section, some three-phase flow models proposed in [26,29] are considered, which aim at describing immiscible three-phase flow mixtures. Property 2 characterizes interfacial pressures that are consistent with the decay of the mixture entropy, while Property 3 gives some conditions on initial pressure gaps in order to guarantee the decay of the latter. Two distinct hybrid multiphase flow models are then investigated. The first one, which was first introduced in [34], enables to describe a mixture of liquid water, water vapor, and a non condensable gas which is assumed to be perfectly miscible with the water vapor. Property 4 provides the counterpart of Property 1 for this hybrid three-field two-phase flow model. Eventually, a four-field three phase flow model [33] is discussed, which was derived in order to cope with mixtures of liquid metal, liquid water and its vapor, together with some non condensable gas. It is briefly shown in Property 6 that this model also requires some conditions on initial pressure gaps in order to guarantee their decay in time. Once more, it should be noted that these conditions are again very unlikely to happen in practice.
Three appendices complete the whole paper. The first one concerns the definition of the Realizable Interfacial Pressure (RIP) condition, while the second one focuses on some algorithms in order to account for pressure relaxation effects. The last one discusses the velocity relaxation effects in a three-phase flow model taken from [28].

GAS-WATER AND LIQUID-VAPOR TWO-PHASE FLOW MODELS
We start our discussion while focusing on two-phase flow models of BN-type, and thus refer the reader to [3] and also to some additional references with slightly different visions, including [9,10,21,37]. These models give the time-space variations of the state variable, and the governing equations include the mass balance, the mean momentum balance, and the energy budgets for each phase. The structure of these models will be briefly recalled below.
For that purpose, we first emphasize that the two immiscible phases (for instance liquid water and water vapor), which are indexed by l, v, are in full desequilibrium, and characterized by their -positive-statistical fractions α l , α v which are such that: We need first to define the state variable W BN which will be the following: where m k = α k ρ k , ρ k , U k and E k respectively denote the -positive-mass fraction, the -positive-density, the mean velocity and the mean total energy, within phase k, setting: The function ε k (P k , ρ k ), stands for the internal energy of phase k, and must be given by user, in terms of the pressure P k and density ρ k , through an equation of state (EOS). Hence the whole BN-type model, including energy budgets, writes: where: and : The left hand side in (4) contains all convective terms, whereas source -relaxation terms lie on the right hand side. For a given closure law of the interfacial velocity: with β BN (W BN ) ∈ [0, 1] to be given in a suitable way, we recall that the interfacial pressure P I (W BN ) is chosen so that the following entropy inequality holds for smooth solutions of (4): where the entropy η BN simply denotes the mixture entropy: thus relying on specific entropies S k (P k , ρ k ), while F BN η denotes the mixture entropy flux: A crucial point to note is that the interfacial pressure P I (W BN ) is defined in a unique way, for a given interface velocity (7), or in other words for a given function β BN (W BN ) : where T l , T v respectively stand for the temperature in phase l, v (see [10]).
It also seems worthwhile recalling at this stage that the basic strategy retained to fix the function β BN (W BN ) simply consists in enforcing the Linearly Degenerate structure for the field associated with the eigenvalue λ = V I (W BN ). This is indeed a major ingredient, since it results in the fact that non-conservative products are well-defined, which means that uniqueness of field-by-field jump conditions is reached. Again, we refer the reader to [10] and [27] for this important feature, however it will not interfer with the sequel.
Eventually, the right hand side term S BN (W BN ), which complies with the entropy inequality (8), is described in [13], and we refer to the latter references for more details. In the sequel we will only need the exact definition of the first and last components of the right-hand side term.
From now on, we consider some homogeneous situation, where: ∂ x (ψ) = 0 whatever ψ stands for, and we first restrict to gas-liquid flows without any mass transfer. We also focus on very large temperature and velocity relaxation time scales, which amounts to get rid of velocity and temperature relaxation effects. Thus the governing set of equations reduces to the following simple system: where the so-called (positive) pressure relaxation time scale τ P is a given function of the state variable W BN . Two slightly distinct forms of τ P are proposed in [6,7,14], which involve molecular viscosities of both components l, v. We recall here that the entropy increases throughout this step (11), since: In order to simplify the presentation, we restrict here to the classical BN-like choice of the interfacial velocity: which in turn implies that : We introduce the variable : ∆P vl = P v − P l , and also, for k = l, v, celerities c k such that: Very simple calculations enable to rewrite (11) so that we get: Hence we obtain: Property 1: (Pressure relaxation effects in a two-phase liquid-gas model) We consider system (11) together with the closure law (14). Then the pressure relaxation process is guaranteed for smooth solutions of (11) if the pressure gap ∆P vl is sufficiently small in the following sense: □ This threshold effect, which was pointed out in [5], arises when taking energy balance into account. Otherwise, in the barotropic case, such a constraint on initial conditions does not arise. In order to fix ideas, an estimation of the upper bound can be found in realistic situations in [33], considering a mixture of liquid metal, liquid water and vapor, and one can easily check that the latter cannot be reached in practice. Note also that Property 1 still holds when considering liquid-vapor water flows.
We will see later on that we retrieve similar conditions when focusing on three-phase flows with immiscible components.

FLOW MODELS INVOLVING THREE IMMISCIBLE COMPONENTS
Very few models have been proposed in the literature till now, in order to tackle three-phase flows involving immiscible components. These are indeed mandatory in order to predict complex flows such as those occuring in vapor explosion for instance. In the latter case, a mixture of hot liquid metal droplets interacts with liquid water and its vapor, and huge pressure waves occur which could create damage when they hit surrounding solid structures. As noted before, the positive statistical fractions for liquid water and its vapor will be noted α l , α v , while α m will stand for the liquid metal statistical fraction. The three of them must satisfy: Obviously, this unsteady application framework urges to consider full disequibrium models. Hence we will focus here on the early model introduced in [26], and also refer the reader to some companion work [28,39]. The state variable is now: with k ∈ (l, v, m). The basic strategy applied in [26] is still the same as before. The governing closed set of equations is expected to comply with the following main requirements: • Smooth solutions of the closed model should verify a physically relevant entropy inequality ; • Unique field-by-field jump conditions must clearly arise; • A symmetric form of the governing set of equations should exist, considering a three-dimensional framework.
The first and third specifications are actually structurant tools when deriving the closed set of equations.
We also take the opportunity here to underline that the two-phase flow framework discussed in section 2 will be retrieved, at least formally, by enforcing α m = 0 in all equations of the three-field model.
We will now briefly describe the whole set of equations and meanwhile comment the key ingredients.
The three components indexed by l, v, m possess their own description in terms of mean density ρ k , mean velocity U k , mean pressure P k , for k ∈ (l, v, m). Internal energy functions ε k (P k , ρ k ) must be prescribed, and the mass fractions and total energies E k within each phase have been defined in the previous section (see (3)). The interfacial velocity is again assumed to be a convex combination of phasic velocities, thus: with the Galilean invariance (GI) constraint arising for The governing set of equations simply reads, for k ∈ (l, v, m): The source terms φ t p f m k (W t p f m ) on the right hand side, which depend on the sole local variable W t p f m , should fufill the following constraint: since the three phases are immiscible. Of course, the latter is mandatory in order to fulfill the maximum principle for the three statistical fractions. Since we focus here on the modeling of interfacial transfer terms, we will also need to enforce the three constraints: in the mass balance equation, but also: for the momentum balance, and eventually: Source terms must also be such that they comply with the entropy inequality (24). The main point at this stage first consists in finding relevant closure laws for the six interfacial pressures Π t p f m i j . This is achieved assuming that an entropy inequality holds for smooth solutions of (19): where the entropy η t p f m simply denotes the mixture entropy: while the entropy flux F t p f m η is: recalling that: (19), the derivation of the time derivative of the entropy (25) is straightforward, and it yields: Thus we get (see [26,29]): Property 2 (Interfacial pressure closure laws in system (19)) • There exists a unique set of six interfacial pressures Π t p f m i j such that smooth solutions of (19) comply with the entropy inequality (24). These read: where k ̸ = i, k ̸ = j and k ∈ (l, v, m), and: The latter interfacial pressures listed in (28) comply with the Realizable Interfacial Pressure (RIP) condition recalled in appendix A.
□ Proof: The first part of the proof is detailed in appendix G of [26], which provides existence and uniqueness in the general case. It simply requires to account for the identity: whatever the state variable W t p f m is. Thus it only remains to check that the formula (28) satify the latter constraint, which is easy though cumbersome. Eventually, since (28) simply stands for some average of phasic pressures, the RIP condition is obviously guaranteed.

□
We wonder now whether the steady states involving equilibrated pressures remain stable, more precisely whether the solutions of the following system : relax towards pressure equilibrium.
As already mentioned in section 2, we emphasize that solutions of (29) are such that: when considering above-mentioned admissible closure laws for Π t p f m k j (W t p f m ), and source terms φ t p f m k (W t p f m ) in agreement with (35), see [29].
Hence, we define the vector of unknowns involving the pressure gaps: for (i, j) ∈ (l, v, m) 2 , as follows: Thanks to (29), it comes: where the matrix A P (W t p f m ) ∈ R 2×2 reads: Note that in the above mentionned matrix, the following notations have been used: but also : It remains to insert entropy-consistent closure laws for the source terms φ t p f m k (W t p f m ), which should be such that: or equivalently: Introducing three independent positive time scales τ P i j , we apply for the following closure laws ( [5]): which are admissible, since (35) is equivalent to: and the latter quadratic form is positive since the matrix D P (W T PFM ): , is symmetric positive. Eventually, we may write:  (19)) We consider system (29) equipped with closure laws (28) and (36), and we use notations introduced in (34) and (31). If we assume that all pressure gaps ∆P i j , for (i, j) ∈ (l, v, m) 2 , are such that: then the pressure relaxation process is guaranteed for solutions of (29).

Sketch of proof:
The two eigenvalues of the real matrix R P = A P (W t p f m )D P (W t p f m ) arising from (38) are either real or complex conjugate. We also know that : det(D P (W t p f m )) > 0.
In the case of real eigenvalues, we must check that both the trace of R P and the determinant of A P (W t p f m ) are positive. Now we have: Moreover, a straightforward calculation also yields: Hence we may conclude that both trace(R P ) and det(A P (W t p f m )) are positive if the following sufficient condition: holds. The latter condition coincides with (39).
In the complex case, we only need to check that the trace of R P is positive. □ When eigenvalues are real, the decay is uniform with respect to time, whereas some oscillations may occur in the relaxation process if complex eigenvalues arise. This was already mentioned in [4] when considering a barotropic three-phase flow model with immiscible components. To the knowlegde of authors, this non-monotone behaviour has not been reported in practical studies up to now.
Eventually, and for computational purposes, the reader is refered to appendix B which provides a straightforward application of the latter results , while focusing on the hybrid model discussed in the following section.

An hybrid two-phase water-vapor model including incondensable gas
The model considered in this section is somewhat different. We focus here on a slightly distinct framework, where the mixture of liquid water and its vapor -respectively indexed by l, valso includes some incondensable gas -indexed by g-(air for instance). We assume that the gas and the water vapor are perfectly miscible, which yields: while we get of course, as it occured in section 2: Owing to the latter two constraints, we use the sole statistical fraction α g in the remaining. We still note the partial mass m k = α k ρ k for k ∈ (l, v, g), and U k , P k , ρ k and E k again stand for the phasic mean velocity, the mean pressure, the mean density and the mean total energy respectively ; the latter phasic total energy E k stands for the sum of the internal energy and the kinetic energy, as recalled in (3). The state variable now stands for: The governing set of equations for mass balance, momentum and energy balance, within phase k, and for the statistical fraction α g are (see [34]): for k ∈ (l, v, g). In order to ease the presentation, we restrict here to the specific closure law: Hence, we get the unique set of interfacial pressures: which is such that smooth solutions of (40) comply with the entropy inequality: where η lvg , F lvg η denotes the mixture entropy -entropy flux pair: Classical constraints still hold for the source terms detailed in [34]: Meanwhile, the source term arising on the right hand side of the statistical fraction equation reads: with a positive function K g (W lvg ).
At this stage, it is worthwhile noting that again, the RIP condition is satisfied, considering a flow initially at rest, with uniform initial pressure and temperature fields such that: . Now we wish to examine whether the pressure relaxation process holds in a natural way for solutions of (40). Thus we consider an homogeneous flow, and evenmore get rid of velocity and temperature effects embedded in source terms, retaining very large temperature and velocity relaxation time scales. This means that we focus on solutions of the set of coupled ODE: Once more, we note that solutions of (45) comply with the entropy inequality, since: Again, introducing : and still setting A k = ρ k c 2 k α k for k ∈ (l, v, g), simple calculations enable to derive:
If the pressure gap ∆P defined in (46) complies with the condition: then the pressure relaxation is guaranteed for solutions of (45).

□
Obviously, we retrieve the same threshold effect as pointed in (16) for the liquid-vapor two-phase flow model, which was expected and hoped. Once more, easy estimations of the upper bound in practical situations can be obtained, which confirm that the latter effect is unlikely to happen. We note that the pressure relaxation acts uniformly with respect to time, unlike in the previous model including three immiscible components.
As mentioned in the previous section, appendix B provides some way to use the latter results in the building of stable schemes to cope with stiff pressure relaxation time scales.

An hybrid four-field three-phase model with incondensable gas
We conclude this section with another hybrid model recently proposed in [33], which is expected to be relevant in order to represent vapor explosion. Thus we focus on a mixture of liquid metal droplets (with index m), liquid water and its vapor (indexed by l, v), together with an incondensable gas (with index g). The derivation of the model was achieved in [33], and the reader is refered to the latter reference, which provides details on the whole procedure, together with the main properties of the closed set of PDEs. Actually, we will only provide here some hints and show that this model belongs to the same class of compressible multiphase flow models.
Since the water vapor and the incondensable gas are assumed to be fully miscible, the statistical fractions will agree with: α v (x,t) = α g (x,t); and we also need to guarantee: α v (x,t) + α l (x,t) + α m (x,t) = 1. since components indexed by v, l, m are immiscible. This means that one should only consider two independent statistical fractions, which are -arbitrary-chosen to be α m , α l in the sequel. With some abuse of notations, the state variable W lvgm now stands for: Using similar notations as in the previous sections, we define the interfacial velocity V lvgm i (W lvgm ) as: with β k (W lvgm ) ≥ 0, and still keeping the Galilean Invariance constraint: Now, the whole model will read: for: k ∈ (l, m, v, g).
Let us now consider the following entropy-entropy flux pair η lvgm , F lvgm η defined by: and: Standard computations enable to derive: We have (see [30,33]): Property 5 (Interfacial pressure closure laws in system (48)) • There exists a unique set of eight interfacial pressures Π lvgm i j (W lvgm ), with i ∈ (l, v, g, m), j ∈ (l, m), such that: Π lvgm i j (W lvgm )∇α j = 0 • The latter interfacial pressures Π lvgm i j (W lvgm ) comply with the Realizable Interfacial Pressure (RIP) condition. □ We refer to [30,33] which give a proof for the first item. Owing to the form of admissible closure laws for the source terms φ lvgm k (W lvgm ) , S lvgm Q k (W lvgm ) and S lvgm E k (W lvgm ) (see [33]), the preservation of the RIP condition can be easily checked, considering initial null velocities, initial pressure fields such that: P l (x,t = 0) = P m (x,t = 0) = P 0 ; P v (x,t = 0) = P 1 ; P g (x,t = 0) = P 0 − P 1 and uniform initial temperature fields : We focus now on the following closure law for the interfacial velocity: and for the statistical fraction evolution: where both functions K l (W lvgm ) and K m (W lvgm ) take positive values, and we note: We consider the subsystem: and we still note A k = ρ k c 2 k /α k . Then we get: Property 6 (Pressure relaxation in the hybrid model (53)) Consider system (53) with closure laws for interfacial pressures and source terms respectively introduced in Property  5 and (51). Assume that the pressure gap ∆P defined in (52) complies with the conditions: (54) with some abuse of notation K k = K k (W lvgm ), and also: Then the pressure relaxation is guaranteed for solutions of (53).

□
We refer to section 4.3 in [33] for more details.

CONCLUDING REMARKS
• We first recall that the four models recalled in (4), (19), (40), or introduced in [33], which have been discussed in the present paper, comply with the general specifications first used in [10], which means that: -Smooth solutions of the global sets of PDEs, including source terms, comply with the entropy inequality for the mixture, and this in turn not only provides a unique definition of interfacial pressures, but it also gives a relevant framework for closure laws of source terms. Evenmore, the former interfacial pressures were shown to agree with the RIP condition; -Owing to the LD structure of the coupling wave associated with the eigenvalue λ = V I , unique field-by-field jump conditions can be derived for these models. This is a crucial point, and of course it is mandatory if one aims at computing flows where shock patterns arise. Defficiencies of other strategies were pointed out some time ago in [23], while restricting to compressible two-phase flow models; -Though all models contain non-conservative first-order contributions in the convective subset, associated systems of PDE can be symmetrized, far from resonance states. This can be simply shown, using two different tools (see [11,29,33,34] and also [19] for a different proof). Of course a straighforward basic consequence is that these models are hyperbolic in the same framework. Another one concerns the existence of smooth solutions [38,43].
• We have focused herein on the pressure relaxation process, thus omitting the velocity and temperature relaxation process, which are of course an intrinsic part of all models considered in the three previous sections. We refer to [29][30][31] that provide more information on the latter two process, and also to appendix C which briefly summarizes the relaxation behaviour of velocity variables in the three-phase flow model (19).
• It has been shown that some threshold effects may arise in the pressure relaxation procedure, unless one restricts to barotropic models (see [4,28]). More over, it has been emphasized that the return towards pressure equilibrium is monotone for the two-phase flow models (4) and (40) ; meanwhile, some -stable-oscillatory behaviour may be observed in three-phase flow models such as (19) or the one introduced in [33].
• We pointed out that the RIP condition, which is recalled in appendix A, is indeed another mandatory building block in order to derive relevant models and closure laws for so-called interfacial pressures. This is not only true for two-phase or three-phase flow models with immiscible components, where it may be understood as a "consistency" condition (see section 2 and 3), but it is also valid when tackling models that aim at describing fast transient multiphase flows including non condensables gases (see section 4).
• The latter properties concerning pressure, velocity and temperature relaxation process at the continuous level may be used when deriving numerical methods. Restricting first to the two-phase flow framework, and applying for the fractional step technique, this has been first used in [13] and [1]. It has also been improved in [12,32] ; one difference between the latter two and the former is that the exact conservation of the total energy of the mixture is enforced, together with the preservation of the min-max principle for the statistical fractions, solving a non-linear 3x3 system with discrete unknowns α n+1 l , P n+1 l , P n+1 v , resulting from a simple implicit Euler discretization of the total energy balance within each phase, and of the governing equation for the statistical fraction.
• This simple strategy was also extended to the three-phase framework in [4] in the barotropic case first, and then in [5] when accounting for energy balances, though choosing slightly distinct algorithms. The leading idea here consists in computing the pressure relaxation process with an implicit scheme, considering an approximate solution of (38), and then enforcing the saturation condition: while ensuring the preservation of the total energy for the mixture, and keeping two frozen entropies, in order to comply with the continuous requirements. This ends up with the resolution of a scalar non-linear equation with one unknown P k 0 . Nonetheless, in some extreme cases corresponding to the vapor explosion situation, a lack of robustness was still observed (see [5]), which suggests to investigate even more stable techniques. In order to cure this weakness, an extension of the latter has been recently proposed (see appendix B), which is currently investigated, while focusing on the three-phase flow model (19) of [26], and also on the hybrid two-phase flow model (40) introduced in [34]. Handling this pressure relaxation step in the most severe cases is actually mandatory in order to cope with the break-up of liquid metal droplets, which has major consequences when simulating vapor explosions, since the heat transfer drastically increases with the interfacial area, which basically modifies predictions of the resulting mixture pressure ∑ k α k P k .

APPENDIX A : THE RIP CONDITION
We define in this appendix the Realizable Interfacial Pressure (RIP) condition, For that purpose, we focus here on the case of immiscible components, but this principle may be easily extended to any multiphase framework, such as those discussed in the present paper. Evenmore, for the sake of clarity, we restrict to the three-phase framework (19), without any mass transfer.
We define some initial condition: W 0 (x) = W (x,t = 0), such that the flow is still: with a uniform pressure field: P k (x,t = 0) = P 0 and temperature profiles within phase k in agreement with: One expects in that case that the fluid will remain steady, that is: Obviously, the source term S(W 0 (x)) vanishes (see [5]), and plugging W 0 (x) in system (19), while using the closure law (18), we get at once: The main objective here is to define two pressure relaxation algorithms in order to obtain approximate solutions of ODE arising in the three-field two-phase flow model (40) associated with [34], while restricting to the sole pressure effects.
In a third part we give some numerical results that are obtained while focusing on the three-phase flow model [26] involving immiscible components, the governing equations of which are recalled in (19).
(1) We focus first on the hybrid two-phase flow model (40), still using the closure law : for the interfacial velocity. Hence, starting from the pressure relaxation model: with : Π 0 τ P , standard calculations yield: for k = g, v, and also: In order to account for pressure relaxation effects associated with (57), a first algorithm consists in taking all steady states into account. This can be achieved by eliminating all variables but α g . Eventually, this ends up in looking for the solution α g of the following scalar equation: where : setting : with ρ k (α g ) = m 0 k α g for k = v, g. The latter equation (58) simply stands for the application of an implicit Euler scheme to the first equation in (57), denoting ψ 0 the initial condition, and ∆t the time step.
Turning now to the existence of solutions α g ∈ [0, 1] for equation (58), we assume that both the EOS for the gas phase and the vapor are perfect gas EOS, and meanwhile, that a stiffened gas EOS holds for the liquid phase, so that: P k = (γ k − 1)ρ k ε k for k = v, g ; P l + γ l Π l = (γ l − 1)ρ l ε l with γ k > 1, for k ∈ (l, v, g). It may be checked that: whatever the time step ∆t is. This implies that (58) admits at least one solution in the admissible range. The proof of uniqueness can be achieved since the function H (X) is monotone.
Once α g has been computed, it only remains to update the remaining variables.
(2) A second strategy may be considered, which is more focused on the pressure relaxation process. We consider now the equivalent form of (57), which writes: where ∆P = P v + P g − P l , remembering that: The algorithm is now twofold: • Compute an approximate value ∆P n+1 of the pressure gaps ∆P at time t n+1 , by solving the last equation of (59).
• Define Y = p n+1 g and note that: Owing to the constraint: α v − α g = 0, we may deduce p n+1 v as a function of Y , solving: which is noted p n+1 v (Y ) in the sequel. By the way, we note that p n+1 v (Y ) is increasing wrt Y . Due to the invariance of the sum of internal energies, we get at once: Eventually, it remains to solve the following scalar equation of unknown Y : where: Considering similar assumptions as before for the three EOS for liquid, vapor and gas phases, it is easy to verify that: lim Y →0 + G (Y ) = +∞ ; lim Y →+∞ G (Y ) = −1 > 0 and also that G (Y ) is decreasing.
Thus we may conclude that the solution Y = P n+1 g ∈ [0, +∞[ of (60) exists and is unique.
• Eventually update all remaining components.
(3) To conclude, we show now an application of the latter techniques in the framework of the three-phase flow model (19).
The exact solution is not recalled herein, and is taken from appendix 4 in [5], which also provides the EOS within each phase, and the initial conditions of the run.
We use the previous pressure relaxation algorithm, which is adapted to the three-phase flow model (19).
We first plot in figure (1) the L 1 norm of the error on the statistical fraction α 1 , at time t = τ P /2, considering successively two different values of τ P = 10 −2 and τ P = 10 −5 . We obviously retrieve the expected first order rate of convergence. The two curves are almost exact translations of one another.
In figure (2), the time evolution of the three pressures is given, considering two different time steps.

COMPONENTS
We turn now to the velocity relaxation process. For sake of simplicity, we concentrate on the barotropic three-phase flow model taken from [28], and focus again on the sole velocity relaxation effects.
with V jk = (U j +U k )/2, and with symmetric positive functions e jk (W t p f m ) involving velocity relaxation time scales τ U jk (see [29] ). We have the following result: Property 7 (Velocity relaxation process) : We consider solutions ∆u(t) of (62). These guarantee the return towards velocity equilibrium, whatever the initial condition ∆u(0) is. □ Starting from (62), we get at once: ∂ t (∆u) = −R u ∆u (63) where the matrix R u ∈ R 2×2 is: are both strictly positive. Hence the two eigenvalues µ u ± of R u are: µ u ± = trace(R u ) ± (∆ u ) 1/2 /2 setting : ∆ u = (trace(R u )) 2 −4det(R u ). These guarantee the -not necessarily uniform-decay of relative velocities u k −u l . □