<figure><img src="image_1.png" /></figure>

HAL
open science

# Simulation of an homogeneous relaxation model for a three-phase mixture with miscible phases

Jean Bussac, Hélène Mathis

► To cite this version:

Jean Bussac, Hélène Mathis. Simulation of an homogeneous relaxation model for a three-phase mixture with miscible phases. 2023. hal-03566970v3

HAL Id: hal-03566970

https://hal.science/hal-03566970v3

Preprint submitted on 20 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# SIMULATION OF A HOMOGENEOUS RELAXATION MODEL FOR A THREE-PHASE MIXTURE WITH MISCIBLE PHASES

J. BUSSAC* AND H. MATHIS

**ABSTRACT.** This paper addresses the numerical approximation of a compressible three-phase flow. The mixture is composed of a liquid, its vapor and an inert gas, with realistic equations of state. Since the liquid is immiscible and the gaseous phases are miscible, it yields constrictive volume and mass constraints. The fluid dynamics is depicted by a one-velocity hyperbolic system with relaxation terms toward the thermodynamical equilibrium of the mixture. The core of the paper is the comparison and the analysis of two possible choices of relaxation terms, accounting for a detailed analysis of mixture entropy and the associated thermodynamical equilibrium. As numerics are concerned, a fractional step approach is proposed, based on a finite volume approximation of the convective part. The complexity of the source term requires an implicit Broyden resolution. Numerical results allow to compare both the performance of the chosen numerical schemes and the impact of the source term choices.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The model</td><td>3</td></tr><tr><td>2.1. Thermodynamics of the mixture</td><td>3</td></tr><tr><td>2.2. Mixture dynamics</td><td>8</td></tr><tr><td>2.3. Source terms</td><td>10</td></tr><tr><td>3. Numerical approximation by a fractional step method</td><td>15</td></tr><tr><td>3.1. Convective part</td><td>16</td></tr><tr><td>3.2. Computation of the internal energy</td><td>18</td></tr><tr><td>3.3. Source term approximation</td><td>21</td></tr><tr><td>4. Numerical results</td><td>21</td></tr><tr><td>4.1. Infinitely slow relaxation</td><td>22</td></tr><tr><td>4.2. Global out-of-equilibrium simulation</td><td>28</td></tr><tr><td>5. Conclusion</td><td>30</td></tr><tr><td>Appendix A. Equilibrium computation for 3 Stiffened Gas</td><td>31</td></tr><tr><td>References</td><td>34</td></tr></table>

## 1. INTRODUCTION

The modelling of compressible multiphase flows is of major importance for many industrial applications, notably in the context of safety of nuclear power plants. Loss of coolant accident in pressurized water reactors and vapor explosions are two typical examples of compressible multiphase flows in this framework [12, 13]. The

*Date:* March 20, 2023.
literature about immiscible multiphase configurations is huge; the reader may re-
fer to recent publications [17, 14] for the modelling and analysis and to [29, 32]
for numerical approximations of such immiscible models. Recently models of com-
pressible flows with miscible phases have been proposed [28, 23, 24, 19, 18], see
also [2]. They depict the dynamical and thermodynamical behaviour of a mixture
composed of a liquid, its vapor and an inert gas ([18] considers in addition an-
other liquid phase). As the fluid dynamics is concerned, two distinct approaches
are available. In [19, 18] the phases evolve with their own velocities, leading to a
system of nonconservative evolution equations with so-called interfacial quantities
(velocity and pressure terms) with relaxation source terms towards equilibrium. In
[28, 23, 24] all the phases evolve with a common velocity, leading to a so-called
Homogeneous Relaxation Model (HRM) similar to an Euler system. It is endowed
with additional evolution equations on fractions of volume, mass and energy of
the phases which contain relaxation source terms towards the thermodynamical
equilibrium. The derivation of this model is proposed in [28], based on the pre-
liminar study [2]. The HRM model corresponds to a relaxation of a Homogeneous
Equilibrium Model (HEM), namely an Euler system with a complex pressure law
which depicts the thermodynamical behaviour of the three-phase mixture. When
considering realistic equations of state for each phase, the exact computation of
this mixture pressure is out of reach and precludes the use of the HEM model. The
relaxation method, used for instance in [10, 25, 16], provides an HRM model with
a simpler Euler-type structure with an explicit pressure law, which inherits from
the hyperbolic and entropic properties of the HEM model. However, it remains
to provide relaxation source terms to close the system. These source terms have
not only to comply with the second law but also to capture the correct thermody-
namical equilibrium of the mixture. Among all possible choices, two sets of source
terms have been proposed in [28], namely BGK-like source terms, corresponding
to a linearization around the equilibrium state, and source term corresponding to
the gradient of the entropy. When considering realistic equations of state for each
phase, Quibel and Hurisse propose in [23, 24] numerical simulations of the HRM
model with BGK-like source term. The aim of the present paper is, among other
things, to compare with the gradient of entropy source term.

The paper is organised as follows.

In a first part, the thermodynamics of the mixture is investigated, when con-
sidering that the three phases are depicted by complete equations of state (EoS),
focusing on Stiffened Gas (SG) and Noble-Able-Stiffened-Gas-Chemkin (NASG-
CK) laws. The mixture entropy is presented in its intensive form (see [28] for an
extensive formulation) and its concavity is investigated. The main difficulties of the
model are, first, that the gaseous phases are miscible while the liquid is immiscible
with them, and second, that the gas is inert, that is it cannot exchange mass with
the two remaining phases. This results in constrictive constraints on the volume
and mass fractions. We define the thermodynamical equilibrium as a maximiza-
tion of the mixture entropy under these constraints and addresses the question of
disappearance of phases.

In a second part, the derivation of the evolution equations is proposed, leading
to the expected HRM model. Note that this derivation differs from the relaxation
method proposed in [28]. A brief study of the hyperbolic structure is given.
Section 2 ends with a comparison of the BGK-like source term and the gradient of entropy one. Properties of the associated dynamical systems are given. If the BGK-like source term can be solved explicitly, it requires the knowledge of the thermodynamical equilibrium for any given state. For realistic EoS, namely NASG-CK or tabulated laws, this is merely out of reach. On the other hand, the gradient of entropy source term cannot be solved explicitly but it enjoys nice properties as well: maximum principle, existence of solutions, asymptotic state in agreement with the expected thermodynamical equilibrium.

Section 3 addresses the numerical approximation of the HRM model by means of a fractional step method. The approximation of the convective part is performed by a finite volume scheme with different numerical schemes. Due to the considered EoS, the computations of internal energy or temperature can be really tedious. For the sake of transparency, we provide in this section detailed algorithms of the computation of internal energy when dealing with either mixture of three SG or with a NASG-CK liquid phase. As the source terms are concerned, their nonlinearity (due to the use of complex EoS) require to use an implicit Broyden algorithm, presented in Section 3.3.

Finally, Section 4 contains the numerical results. First verification test cases of the convective part are presented, corresponding to infinitely slow relaxation. Exact Riemann problems are computed for a mixture of three SG and a mixture of a NASG-CK liquid with two SG. In both configurations a low and high pressure test cases are proposed. Convergence results are given which assess the good behaviour of the relaxation and the VFRoe-ncv schemes. However, the VFRoe-ncv scheme is shown to be much more CPU time-consuming than the relaxation scheme, and has a lack of robustness due to numerous internal energy computations, see Section 4.1.3. To finish a global simulation is proposed, accounting for the relaxation source terms, corresponding to a non-equilibrium test case. Numerical illustrations highlight the impact of the choice of the source terms.

## 2. THE MODEL

This Section presents the modelling of the three-phase mixture. First its thermodynamics is investigated. Compared to [28], complementary results about the mixture entropy concavity and the thermodynamical equilibrium are given. Notably the question of disappearance of phases is addressed. The dynamics of the fluid is depicted by a hyperbolic system with relaxation terms, whose derivation is detailed. The core of the section is the comparison and the analysis of two possible choices of source terms proposed in [28].

### 2.1. Thermodynamics of the mixture.
The mixture is composed of three phases: a liquid (labelled *l*), its vapor (labelled *v*) and an inert gas (labelled *g*). The gas is assumed to be inert in the sense that it cannot exchange mass with the two other phases (hence it is not *stricto sensu* a phase but the terminology is conserved in the sequel). We assume that no vacuum can occur. The main feature of this model relies on the fact that the two gaseous phases are miscible, while the liquid phase is immiscible with them.

In order to describe the thermodynamical behavior of the mixture, we consider a given state (τ, *e*), with *τ* the specific volume of the mixture and *e* its specific internal energy. The question addressed in the section 2.1.2 is the proper definition of the mixture entropy function which entirely depicts the thermodynamical behaviour of
the mixture. To do so, Section 2.1.1 presents first the complete equation of state for each phase $k \in \{l, g, v\}$. Section 2.1.2 starts by the constraints of mass and energy conservation as well as a volumic constraint which traduces the miscibility of the gaseous phases. It allows to express the mixture entropy as a function of fraction of mass, volume and energy of the phases. According to the second principle of thermodynamics, the thermodynamical equilibrium is reached at the maximum of the mixture entropy. This maximization process is realized under the abovementioned constraints. We present in Section 2.1.3 the associated optimality conditions, that give several constraints on phase disappearance.

2.1.1. *Phase setting.* Each phase $k = l, g, v$ is depicted by its specific volume $\tau_k > 0$ and its specific internal energy $e_k > 0$. The thermodynamical behaviour of each phase $k$ is fully described by its intensive entropy function $(\tau_k, e_k) \mapsto s_k(\tau_k, e_k)$ defined on $\Omega_k \subset (\mathbb{R}_+^*)^2$, that is a complete equation of state. By adopting the Gibbs formalism, each entropy function $s_k$ complies with the following differential form:

$$ (1) \qquad T_k ds_k = de_k + p_k d\tau_k, $$

where the temperature $T_k$ and pressure $p_k$ are defined by:

$$ (2) \qquad \frac{1}{T_k} = \left. \frac{\partial s_k}{\partial e_k} \right|_{\tau_k}, \qquad p_k = T_k \left. \frac{\partial s_k}{\partial \tau_k} \right|_{e_k}, $$

and the chemical potential by the relation:

$$ (3) \qquad \mu_k = -T_k s_k + p_k \tau_k + e_k. $$

Following hypotheses that are usually made on the extensive entropies [26], the phase entropy functions are supposed to be:

* strictly concave functions on $\Omega_k$,

* of class $C^2$ on $\Omega_k$, such that $\forall (\tau_k, e_k) \in \Omega_k$, $T_k(\tau_k, e_k) > 0$.

In practical applications, we focus on two equations of state, namely the Stiffened Gas (SG) and the Noble-Able Chemkin Stiffened Gas (NASG-CK) equations of state (EoS), whose accuracy on a large domain have been studied in [30, chap. 2]. We remind here their definitions.

For a given state $(\tau_k, e_k) \in \Omega_k$, the complete SG reads

$$ (4) \qquad s_k(\tau_k, e_k) = C_{v,k} \ln \left( (e_k - Q_k - \Pi_k \tau_k) \tau_k^{\gamma_k - 1} \right) + s_k^0, $$

where $C_{v,k} > 0$ is the calorific capacity at constant volume, $-\Pi_k$ is the minimal pressure, $Q_k$ is a reference enthalpy, $\gamma_k > 1$ is the adiabatic coefficient and $s_k^0 > 0$ is a reference specific entropy. The case $\Pi_k = Q_k = 0$ corresponds to a perfect gas law [27]. The associated pressure and temperature of phase $k$ are

$$ (5) \qquad p_k(\tau_k, e_k) = (\gamma_k - 1) \frac{e_k - Q_k}{\tau_k} - \gamma_k \Pi_k, \qquad T_k = \frac{e_k - Q_k - \Pi_k \tau_k}{C_{v,k}}. $$

In [30, chap. 2] are given sets of parameters for the vapor and liquid phases, which have been fitted with respect to experimental data around a reference point ($p_{ref} = 80\text{bar}, T_{ref} = 425\text{K}$). We use them for further simulations, and a perfect gas law for the gaseous phase, see Tables 1 and 2.

The Noble-Able Stiffened Gas Chemkin (NASG-CK) EoS is an extension of the SG where the heat capacity now depends on the temperature. According to [30,
<table><caption>TABLE 1. Stiffened Gas parameters</caption><thead><tr><th>Parameters</th><th>Liquid phase</th><th>Vapor phase</th></tr></thead><tbody><tr><td>γ<sub>k</sub> (dimensionless)</td><td>1.39864082368510</td><td>1.15442237458290</td></tr><tr><td>C<sub>k,l</sub> (JK<sup>-1</sup>kg<sup>-1</sup>)</td><td>3.19641035947920×10<sup>3</sup></td><td>2.91668522329726×10<sup>3</sup></td></tr><tr><td>Q<sub>k</sub> (Jkg<sup>-1</sup>)</td><td>-1.24606074764184 × 10<sup>6</sup></td><td>1.25942536895827 × 10<sup>6</sup></td></tr><tr><td>Π<sub>k</sub> Pa</td><td>4.79690712132593 × 10<sup>8</sup></td><td>-3.24993579473092 × 10<sup>2</sup></td></tr><tr><td>s<sub>k</sub><sup>0</sup> (JK<sup>-1</sup>kg<sup>-1</sup>)</td><td>-34597.52986978335</td><td>-33792.18353359583</td></tr></tbody></table>

<table><caption>TABLE 2. Perfect Gas parameters for the gaseous phase</caption><thead><tr><th>Parameters</th><th>Gaseous phase</th></tr></thead><tbody><tr><td>γ<sub>k</sub> (dimensionless)</td><td>1.4</td></tr><tr><td>C<sub>k,l</sub> (JK<sup>-1</sup>kg<sup>-1</sup>)</td><td>719</td></tr><tr><td>Q<sub>k</sub> (Jkg<sup>-1</sup>)</td><td>0</td></tr><tr><td>Π<sub>k</sub> Pa</td><td>0</td></tr></tbody></table>

chap. 2], the NASG-CK EoS provides a more accurate description of the liquid thermodynamical behaviour on a larger domain than the SG EoS. This accuracy has a cost since it requires nonlinear computations, the NASG-CK law being semi-implicit. Indeed, the chemical potential (which is a complete EoS) is classically given in the pressure-temperature plane as follows

$$ (6) \qquad \mu_l(p_l, T_l) = \mu_l^0(T_l) + b_l p_l + C_l(T_l) \ln(p_l + \Pi_l), $$

where

$$ (7) \qquad \begin{aligned} \mu_l^0(T_l) &= RT_l \left( A_l(1 - \ln(T_l)) - \frac{B_l}{2}T_l - \frac{C_l}{6}T_l^2 - \frac{D_l}{12}T_l^3 - \frac{E_l}{20}T_l^4 + \frac{R_l}{T_l} - G_l \right), \\ C_l(T_l) &= C_{v,l}(\gamma_l - 1)T_l, \end{aligned} $$

with $R = 0.4615228083134561$ the perfect gas constant for water and $b_l$ a translation of the specific volume coming from the NASG law, see [5]. For further numerical experiments, we focus on parameters given in [24] for the liquid water

<table><caption>TABLE 3. NASG-CK parameters for the liquid phase</caption><thead><tr><th>Parameters</th><th>Value</th></tr></thead><tbody><tr><td>γ<sub>l</sub> (dimensionless)</td><td>3.27113568773712</td></tr><tr><td>C<sub>v,l</sub> (JK<sup>-1</sup>kg<sup>-1</sup>)</td><td>7.24509640448929 × 10<sup>2</sup></td></tr><tr><td>Π<sub>k</sub> (Pa)</td><td>1.24425779880749 × 10<sup>9</sup></td></tr><tr><td>A<sub>l</sub> (dimensionless)</td><td>4.69738865636393 × 10<sup>1</sup></td></tr><tr><td>B<sub>l</sub> (K<sup>-1</sup>)</td><td>-4.19269571479452 × 10<sup>-1</sup></td></tr><tr><td>C<sub>l</sub> (K<sup>-2</sup>)</td><td>1.70702143968620 × 10<sup>-3</sup></td></tr><tr><td>D<sub>l</sub> (K<sup>-3</sup>)</td><td>-3.04805662517983 × 10<sup>-6</sup></td></tr><tr><td>E<sub>l</sub> (K<sup>-4</sup>)</td><td>2.02814588067819 × 10<sup>-9</sup></td></tr><tr><td>F<sub>l</sub> (K<sup>-5</sup>)</td><td>-2.519604765 × 10<sup>6</sup>/R</td></tr><tr><td>G<sub>l</sub> (dimensionless)</td><td>-5.573536947 × 10<sup>4</sup>/R</td></tr><tr><td>b<sub>l</sub> (m<sup>3</sup>kg<sup>-1</sup>)</td><td>5.66559849022606 × 10<sup>-4</sup></td></tr></tbody></table>
The computation difficulties come from the implicit calculation of the pressure $p_l$ and the temperature $T_l$ in the $(\tau_l, e_l)$ plane. According to the following derivative in the $(p_l, T_l)$ plane

$$ (8) \qquad \tau(p, T) = \left. \frac{\partial \mu}{\partial P} \right|_T, \quad s(p, T) = \left. \frac{\partial \mu}{\partial T} \right|_p. $$

that can be obtained by using the Gibbs relation (1) and the Legendre transform, the pressure $p_l$ and the internal energy $e_l$ are given by

$$ (9) \qquad p_l(\tau_l, T_l) = \frac{C_{v,l}(\gamma_l - 1)T_l}{\tau_l - b_l} - \Pi_l, $$

$$ (10) \qquad e_l(\tau_l, T_l) = RT_l \left( A_l + \frac{B_l}{2}T_l + \frac{C_l}{3}T_l^2 + \frac{D_l}{4}T_l^3 + \frac{E_l}{5}T_l^4 + \frac{F_l}{T_l} \right) - C_{v,l}(\gamma_l - 1)T_l + \Pi_l(\tau_l - b_l). $$

Hence for a given state $(\tau_l, e_l)$, first the temperature is computed by solving (10) in an implicit way, and then, the pressure is deduced using (9).

### 2.1.2. Intensive constraints and mixture entropy.
We now turn to the description of the thermodynamical behaviour of the mixture. Let $w = (\tau, e)$ be a state vector of the fluid mixture, which is composed of a liquid $l$, its vapor $v$ and an inert gas $g$. We denote $y_k \in [0, 1]$ the mass fraction of the phase $k$, $\alpha_k \in [0, 1]$ the volume fraction and $z_k \in [0, 1]$ the energy fraction. The conservation of total mass and energy provides that

$$ (11) \qquad \left\{ \begin{array}{l} y_l + y_v = 1 - y_g, \\ z_l + z_v + z_g = 1, \end{array} \right. $$

while the miscibility of the gaseous phases and the immiscible behaviour of the liquid impose that

$$ (12) \qquad \left\{ \begin{array}{l} \alpha_g = \alpha_v, \\ \alpha_l + \alpha_g = 1. \end{array} \right. $$

Note that the mass fraction $y_g$ of the gas is fixed, since the gas cannot exchange mass with the two remaining phases. The knowledge of the fractions and the mixture state $(\tau, e)$ allows to recover the state $(\tau_k, e_k)$ since

$$ (13) \qquad \tau_k = \frac{\alpha_k}{y_k} \tau, \quad e_k = \frac{z_k}{y_k} e. $$

For a given mixture state $(\tau, e)$, the mixture entropy is defined as a convex combination of the phase entropies with weights $y_k \in [0, 1]$. Denoting $Y = (y_l, \alpha_l, z_l, z_g)$, it reads

$$ (14) \qquad \sigma(\tau, e, Y, y_g) = \sum_{k=l,g,v} y_k s_k \left( \frac{\alpha_k}{y_k} \tau, \frac{z_k}{y_k} e \right). $$

The definition domain $\Omega$ of the mixture entropy is a convex subset of $(\mathbb{R}_+)^2 \times ]0; 1]^5$ [26].
**Remark 1.** (Strict concavity of the mixture entropy) We can show that $(\tau, e) \mapsto \sigma(\tau, e, Y, y_g)$ and $Y \mapsto \sigma(\tau, e, Y, y_g)$ are concave functions [22, 28]. Some comments concerning this strict concavity of the mixture entropy have been provided for instance in [22], in the case of an immiscible mixture. They rely on the study of the concavity of the extensive mixture.

The idea is to show, for a given extensive state $X = (M, V, E)$ of mass $M$, volume $V$ and energy $E$, that the extensive mixture entropy $S$ is a strictly concave function when restricted on the set $\mathcal{H}(M)$ of the states whose total mass is equal to $M$. If this property is true, then the intensive mixture entropy is strictly concave.

In order to prove the strict concavity of $S$ on $\mathcal{H}(M)$, the idea is to show that the degeneracy manifold of $S$, that is $\ker \nabla_X^2 S(X)$, coincides with $\mathrm{Vect}(X) := \{aX, a \in \mathbb{R}\}$ [26]:

$$ \forall X, \ker \nabla_X^2 S(X) = \mathrm{Vect}(X). $$

If so, then the restriction of $S$ on $\mathcal{H}(M)$ is strictly concave since

$$ \mathcal{H}(M) \cap \ker \nabla_X^2 S(X) = \{X\}. $$

However, only the inclusion $\mathrm{Vect}(X) \subset \ker \nabla_X^2 S(X)$ is shown [26]. One remarks that $\ker \nabla_X^2 S(X)$ is a vectorial subset in the 9-dimensional space in the case of a three-phase flow.

Following [8], few numerical verifications on the sign of $\det \nabla_Y^2 \sigma$ have been done for some given $(\tau, e, y_g)$ and no default of strict concavity were found. We emphasize that these tests help us to know if this latter hypothesis is reasonable, but in complete numerical applications, the diversity of the states could lead to difficulties. A more complete study on the whole domain of the entropy would be welcomed.

In this paper, we will consider that $Y \mapsto \sigma(\tau, e, Y, y_g)$ is strictly concave.

Accounting for the Gibbs relations (1), identities (13) and the definition (3) of the chemical potential $\mu_k$, the differential of the mixture entropy is

$$ d\sigma = \sum_{k=l,g,v} \mu_k dy_k + \tau \sum_{k=l,g,v} \frac{p_k}{T_k} d\alpha_k + e \sum_{k=l,g,v} \frac{1}{T_k} dz_k + \sum_{k=l,g,v} z_k \frac{1}{T_k} de + \sum_{k=l,g,v} \alpha_k \frac{p_k}{T_k} d\tau. $$

Now using the intensive constraints (11) and (12), it holds

$$ \begin{aligned} d\sigma = & (\mu_l - \mu_v)dy_l + \tau \left( \frac{p_l}{T_l} - \left( \frac{p_v}{T_v} + \frac{p_g}{T_g} \right) \right) d\alpha_l \\ & + e \left( \frac{1}{T_l} - \frac{1}{T_g} \right) dz_l + e \left( \frac{1}{T_v} - \frac{1}{T_g} \right) dz_v \\ & + \left( \frac{z_l}{T_l} + \frac{z_v}{T_v} + \frac{z_g}{T_g} \right) de + \left( \alpha_l \frac{p_l}{T_l} + (1 - \alpha_l) \left( \frac{p_v}{T_v} + \frac{p_g}{T_g} \right) \right) d\tau. \end{aligned} \tag{15} $$

Hence it follows that the mixture pressure and temperature are defined by:

$$ \begin{aligned} \frac{1}{T}(\tau, e, Y, y_g) &= \frac{z_l}{T_l} + \frac{z_v}{T_v} + \frac{z_g}{T_g}, \\ p(\tau, e, Y, y_g) &= \frac{\alpha_l \frac{p_l}{T_l} + (1 - \alpha_l) \frac{p_v}{T_v} + (1 - \alpha_l) \frac{p_g}{T_g}}{\frac{z_l}{T_l} + \frac{z_v}{T_v} + \frac{z_g}{T_g}}, \end{aligned} \tag{16} $$

where the notations $p_k$ (resp. $T_k$) stands for $p_k(\tau_k, e_k)$ (resp. $T_k(\tau_k, e_k)$) with relations (13).
2.1.3. **Thermodynamical equilibrium and optimality conditions.** As thermodynamical equilibrium is reached, the mixture entropy achieves its maximum. This maximization process is performed under the constraints (11)-(12), and we remind here that we suppose that $Y \mapsto \sigma(\tau, e, Y, y_g)$ is strictly concave. Thus we can define the equilibrium fractions by

**Definition 1. (Equilibrium fractions)** Let us fix $(\tau, e, y_g)$. Assuming that $Y \mapsto \sigma(\tau, e, Y, y_g)$ is strictly concave, we define the equilibrium fractions $Y_{eq}$ by

$$ (17) \qquad Y_{eq}(\tau, e, y_g) = \underset{Y \in C_w}{\operatorname{argmax}} \sigma(\tau, e, Y, y_g), $$

where $C_w$ is the definition domain of $Y \mapsto \sigma(\tau, e, Y, y_g)$ for a given $(\tau, e, y_g)$.

Let us mention phase disappearances. The model does not seem to be able to handle strict phase disappearance, i.e. a phase whose fractions $y_k$ and $\alpha_k$ are equal to zero. This question is very complex and would require to study how the objects (entropy function, thermodynamical potentials...) evolve towards these one or two phase situations.

However, it is necessary to consider the disappearance cases for the equilibrium fractions. Indeed, a three-phase mixture that admits a two-phase equilibrium is a quite standard situation (at least in our numerical tests). Since the gas is inert, there are only two particular equilibrium cases: the disappearance of liquid or the disappearance of vapor.

* If the equilibrium state corresponds to a liquid-inert gas mixture, the equilibrium vector should be $Y = (1 - y_g, \alpha_l, 1 - z_g, z_g)$. This case will be illustrated on Figure 2-B.

* If the equilibrium state corresponds to a vapor-gas mixture, the equilibrium vector should be $Y = (0, 0, 0, z_g)$. This particular value does not belong to the definition domain of $\sigma$, but it is mandatory. This case will be illustrated on Figure 3.

The case of a three-phase equilibrium corresponds to a saturation state characterized by the following relations, see [28]

$$ (18) \qquad \left\{ \begin{array}{l} p_l = p_v + p_g, \\ T_l = T_v = T_g, \\ \mu_l = \mu_v. \end{array} \right. $$

A practical approach that can handle phase disappearance can be found in [24].

2.2. **Mixture dynamics.** In this section are introduced the evolution equations of the mixture, assuming that each phase evolves with a common velocity field $u$. Their derivation consists in studying the evolution of an isolated fluid element along a streamline. Note that the resulting set of equations has already been obtained in [28] as a relaxed version of a three-phase equilibrium model.

2.2.1. *Derivation of the evolution equations of the mixture.* The purpose of this section is to derive a simple system of evolution equations of the observable quantities of the mixture in time $t$ and space $x$, considering that each phase evolves with a common velocity field $u$. The derivation is performed for $x \in \mathbb{R}^d$, even if we focus
in the sequel on the one-dimensional setting. To do so, one considers an isolated fluid element, composed of the three phases, along a streamline. Some standard assumptions have to be stated to perform the derivation, following the formalism depicted in [22]. First, the total mass of the element is conserved along time, as well as the mass fraction of gas $y_g$. Moreover, we assume that the volume variation of the element is solely due to the velocity field divergence. Hence one obtains that

$$d\tau = \tau \mathrm{div}(u)dt,$$

where the derivative notation $d\varphi$ corresponds to the derivative along a streamline of the flow

$$d\varphi = \partial_t\varphi + u\nabla_x\varphi.$$

Considering now that the only force acting on the fluid element is due to the pressure gradient, the fundamental principle of mechanics states that

$$du = -\tau \nabla_x p dt,$$

where the mixture pressure $p$ corresponds to the mixture pressure (16). The first law of thermodynamics then relates the internal energy to the pressure assuming that the fluid element is not submitted to other heat sources

$$de = -pd\tau.$$

Following [10, 4, 22], we assume that the time-evolution of the phase quantities are of the form

$$dY = \lambda \Gamma(\tau, e, Y, y_g)dt,$$

where the vector $\lambda$ corresponds to the time relaxation parameters.

Combining the above equations with the definition of the fluid density $\rho = 1/\tau$, one ends up with a compressible Euler-type system with five additional equations on the fractions with possible relaxation terms. The conservative one-dimensional set of equations reads

$$ (19) \qquad \left\{ \begin{array}{l} \partial_t(\rho Y) + \partial_x(\rho Y u) = \rho \lambda \Gamma(\tau, e, Y, y_g), \\ \partial_t(\rho y_g) + \partial_x(\rho y_g u) = 0, \\ \partial_t \rho + \partial_x(\rho u) = 0, \\ \partial_t(\rho u) + \partial_x(\rho u^2 + p) = 0, \\ \partial_t(\rho E) + \partial_x((\rho E + p)u) = 0, \end{array} \right. $$

where $E = u^2/2 + e$ denotes the total energy of the system and $p$ is the mixture pressure (16).

### 2.2.2. Properties of the homogeneous model.
Focusing on the homogeneous part of the model, its hyperbolicity and entropic structure have been studied in [28]. The concavity of $\sigma$ with respect to $(\tau, e, Y, y_g)$ is equivalent to the convexity of $H = -\rho\sigma$ with respect to the conservative variables $W = (\rho Y, \rho u, \rho E)$. Thus $H = -\rho\sigma$ is a Lax entropy for the system (19) and weak solutions satisfy

$$\partial_t H + \partial_x(Hu) \le 0.$$

Again hyperbolicity is ensured as soon as the three phases are present and depicted by concave entropies leading to a strictly positive sound speed $c$ defined by

$$ (20) \qquad \frac{c^2}{\tau^2} = p \partial_e p - \partial_\tau p = -T(p^2 \partial_{ee} \sigma - 2p \partial_{\tau e} \sigma + \partial_{\tau \tau} \sigma). $$
The eigenstructure of the system (19) is composed of one linearly degenerate wave, associated to the eigenvalue $u$, and two genuinely non-linear waves associated to the eigenvalues $u \pm c$, where $c$ is the mixture sound speed defined by (20)

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Structure: ghost-wave, contact, shock.</figcaption></figure>

Figure 1 represents the wave structure in the $(x, t)$ plane. As the contact wave is conserved (see Figure 1- dashed line), mixture pressure $p$ and velocity $u$ are preserved across it contrary to the fractions $Y$. For two given states $W_L$ and $W_R$ separated by a genuinely nonlinear wave (plain lines of Figure 1), the jump conditions read

$$ (21) \qquad \begin{cases} [Y] = 0, \\ J = -\frac{[p]}{[u]}, \\ J^2 = -\frac{[p]}{[\tau]}, \\ [e] + [\tau] \frac{p_L + p_R}{2} = 0, \end{cases}. $$

where $J = \rho u - \sigma$. It means in particular that the fractions are constant through nonlinear waves. In numerical verifications, we provide analytical solutions of Riemann problems when the system (19) is endowed with an initial condition $W(0, x) = W_L \mathbf{1}_{x<0} + W_R \mathbf{1}_{x \ge 0}$, see Section 4.1. To do so we make use of the jump relations (21) to compute solutions composed of four constant states separated by a ghost wave, a contact and a shock wave (see also [24] for detailed exact computations).

2.3. **Source terms.** The thermodynamical behaviour of the mixture is not only driven by the mixture pressure $p$ involved in the momentum and total energy conservation laws, but mostly depicted by the relaxation terms acting on the advection equations on the fraction vector $Y$. Indeed, all the thermodynamical disequilibrium is depicted by the dynamical system

$$ (22) \qquad \frac{dY}{dt} = \lambda\Gamma(\tau, e, Y, y_g). $$

Even if there is no consensus on the choice of relaxation source terms, they have to fulfill some physical requirements. First a thermodynamic evolution of the system must lead to an entropy growth and the equilibrium state corresponds to its
maximum. It turns out that the mixture entropy $\sigma$ verifies the following transport equation:

$$ (23) \qquad \partial_t \sigma + u \partial_x \sigma = \lambda \Gamma(Y, \tau, e) \nabla_Y \sigma. $$

Hence eligible source terms $\Gamma$ are such that

$$ (24) \qquad \lambda \Gamma(Y, \tau, e) \nabla_Y \sigma \geq 0. $$

Moreover, the relaxation terms have to preserve the asymptotic states. For a given state $(\tau, e)$ of the mixture, the source terms have to capture the correct thermodynamical equilibrium, that is

$$ \lim_{t \to +\infty} Y(t) = Y_{eq}. $$

The core of this section is to analyze two possible choices of source terms. The first one, referred in [21] as BGK source term, corresponds to a linearization around the equilibrium vector $Y_{eq}$. Introduced in [10, 16, 4], it is used in multiphase configurations with complex equations of state [24, 23]. The alternative source terms, proposed in [28], corresponds to the gradient of the mixture entropy which directly satisfies (23).

A similar analysis of these two choices in the immiscible (resp. miscible) two-phase framework is performed in [8]. The advantage of the two-phase case is that the dynamical system admits three (resp. two) equations leading to a convenient vizualisation of the trajectories in the phase plane. For the present problem, the dimension 4 of the system prevents such nice vizualisation.

**Remark 2.** Source terms $\Gamma$ contain time scales $\lambda = (\lambda_1, \lambda_2, \lambda_3, \lambda_4) \in (\mathbb{R}_+)^4$ that model the relaxation speed of the system. We will not consider any time dependence for $\lambda$ in this paper.

More precisely, each $\lambda_i$ is a factor contained by $\Gamma_i$, and concerns a specific fraction in system (19). On a physical point of view, $\lambda_1$ correspond to the mass transfer, $\lambda_2$ to the mechanical, $\lambda_3$ and $\lambda_4$ to thermal relaxation times.

2.3.1. BGK-like source terms. For a given state $(\tau, e, Y, y_g) \in \Omega$, assuming that $Y \mapsto \sigma(\tau, e, Y)$ is strictly concave, there exists a unique equilibrium fraction state $Y_{eq}$ which maximises the entropy $\sigma$, where $Y_{eq}$ only depends on $(\tau, e, Y)$, see Proposition 2.1.3.

**Definition 2.** For any state $(\tau, e, Y, y_g) \in \Omega$, we define:

$$ (25) \qquad \Gamma_1(\tau, e, Y, y_g) = \lambda(Y_{eq}(\tau, e, y_g) - Y), $$

where $\lambda \in (\mathbb{R}_+^4)$ is the time scale vector.

**Proposition 1.** (Eligibility) Assume that the time scale vector $\lambda > 0$ is such that $\lambda_1 = \lambda_2 = \lambda_3 = \lambda_4 = \lambda$, possibly time dependent. For any given state $(\tau, e) \in (\mathbb{R}_+^*)^2$, assume that the entropy function $(\tau, e, Y, y_g) \mapsto \sigma(\tau, e, Y)$ is strictly concave with respect to $Y$. Then, the source terms $\Gamma_1$ complies with the entropy growth criterion (24). Moreover, the stationary states of the dynamical system (22)-(25) coincide with the possible thermodynamical equilibria (18).

*Proof.* It is obvious that the stationary states of the source term $\Gamma_1$ are vector $Y_{eq}$ of Proposition 18. It remains to check the entropy growth. Let $(\tau, e, Y, y_g) \in \Omega$ be
a given state. It holds:

$$
\begin{align*}
\Gamma_1(\tau, e, Y, y_g) \nabla_Y \sigma(\tau, e, Y, y_g) &= (\lambda(Y_{eq} - Y)) \nabla_Y \sigma(\tau, e, Y, y_g) \\
&= \lambda((Y_{eq} - Y) \nabla_Y \sigma(\tau, e, Y, y_g)) \\
&\geq \lambda(\sigma(\tau, e, Y_{eq}, y_g) - \sigma(\tau, e, Y, y_g)) \\
&\geq 0,
\end{align*}
$$

thanks to the concavity of $Y \mapsto \sigma(\tau, e, Y, y_g)$.

<figure><img src="image_3.png" /></figure>

Note that invoking the concavity argument requires to consider a *unique* relaxation parameter $\lambda$. It corresponds to impose equal relaxation times for mass, mechanical and thermal transfers. This is a very strong constraint that seems not to be relevant on a physical point of view. Recently a modification of BGK-like source terms which allows distinct time scales has been introduced in [21].

**Proposition 2.** Let $(\tau, e, Y_0, y_g) \in \Omega$ be an initial state. The Cauchy problem

$$
\begin{cases}
\frac{dY}{dt}(t) = \Gamma_1(\tau, e, Y(t), y_g), & t > 0, \\
Y(0) = Y_0,
\end{cases}
$$

admits a global solution $Y$ given by

$$
(26) \quad \forall t \ge 0, \quad Y(t) = e^{-\Lambda(t)} (e^{\Lambda(0)} Y_0 + Y_{eq} \int_0^t \lambda(s) e^{\Lambda(s)} ds),
$$

where $\Lambda$ is a primitive of $\lambda$. In the case of a constant in time relaxation parameter $\lambda$, it holds

$$
(27) \quad \forall t \ge 0, \quad Y(t) = e^{-\lambda t} Y(0) + (1 - e^{-\lambda t}) Y_{eq}.
$$

The proof is direct and this last expression can be interpreted as a barycenter of the initial state $Y(0)$ and the equilibrium state $Y_{eq}$. Hence it guarantees that $Y(t) \in [0, 1]^4$ for all $t \ge 0$.

One major advantage of the BGK-like source term is that it can be explicitly integrated, since trajectories are exponential towards equilibrium fractions $Y_{eq}$. However, it requires to compute the equilibrium state $Y_{eq}$ for every considered state $(\tau, e, Y, y_g) \in \Omega$, that is to say maximise the mixture entropy (14) under the constraints (11)-(12).

For two-phase mixtures, the maximization can be done explicitly when classical EoS are considered [20]. In our three-phase context, the computations are much more complicated, especially when complex EoS are considered. For sake of illustration, the detailed computation of the equilibrium state $Y_{eq}$ for a mixture of three Stiffened Gases is given in Appendix A.

Figure 2 presents two evolutions in time of the fractions $Y(t)$ for a mixture of three SG with parameters 1. The time scale is constant in time and set to $\lambda = 1$.

On Figure 2-A, one observes the evolution in time of the fractions $Y(t)$ for an initial state $Y(0) = (0.1, 0.1, 0.1, 0.05)$, $y_g = 0.05$, with specific volume $\tau = 5.10^{-3}$ and internal energy $e = 10^8$. The associated equilibrium state is a three-phase equilibrium.

On Figure 2-B, one observes the evolution in time of the fractions $Y(t)$ for an initial state $Y(0) = (0.1, 0.1, 0.1, 0.05)$, $y_g = 0.3$, with specific volume $\tau = 4.54 \cdot 10^{-2}$ and internal energy $e = 10^8$. The associated equilibrium is near to a two-phase mixture, where $y_{v,eq} \sim 10^{-2}$ and $z_{v,eq} \sim 10^{-2}$.
<figure><img src="image_4.png" /><figcaption>(A) Trajectories toward a three-phase equilibrium. (B) Trajectories toward a liquid-gas equilibrium.</figcaption></figure>

<caption>FIGURE 2. BGK-like source term $\Gamma_1$ for a mixture of three SG.</caption>

Another important point is that the time when the asymptotic state is reached does not depend on the initial state of the fluid. Hence, the equilibrium fractions (that can be computed following Appendix A) are always reached at time $t = 5s$ (unless $\lambda \neq 1$).

When considering a NASG-CK EoS for the liquid, the computation of the equilibrium state vector $Y_{eq}$ requires an implicit solver, for instance using a Broyden algorithm, see [24] for a hint of the method.

### 2.3.2. Entropy gradient source terms.
This second type of source terms corresponds to the gradient of the mixture entropy.

**Definition 3.** For all $(\tau, e, Y, y_g) \in \Omega$, let $\Gamma_2$ be the source terms function defined by

$$ \Gamma_2(Y) = Y(1-Y)\nabla_Y\sigma(\tau, e, Y, y_g) = \begin{pmatrix} y_l(1-y_l)\partial_{y_l}\sigma \\ \alpha(1-\alpha)\partial_{\alpha}\sigma \\ z_l(1-z_l)\partial_{z_l}\sigma \\ z_g(1-z_g)\partial_{z_g}\sigma \end{pmatrix}. $$

Using intensive potentials (2)-(3), $\Gamma_2$ rewrites:

$$ (28) \qquad \Gamma_2 = \begin{pmatrix} y_l(1-y_l)\left(\frac{\mu_v}{T_v} - \frac{\mu_l}{T_l}\right) \\ \alpha_l(1-\alpha_l)\tau\left(\frac{P_l}{T_l} - \frac{P_v}{T_v} - \frac{P_g}{T_g}\right) \\ z_l(1-z_l)\left(\frac{1}{T_l} - \frac{1}{T_v}\right) \\ z_g(1-z_g)\left(\frac{1}{T_g} - \frac{1}{T_v}\right) \end{pmatrix}. $$

This function corresponds to the entropy gradient times $Y(1-Y)$ understood as a term by term product. Indeed, the set of definition of the mixture entropy $Y \mapsto \sigma(\tau, e, Y, y_g)$ is $\Omega$ which is included in $]0; 1]^4$. Thus a solution $Y$ of the dynamical system $Y'(t) = \nabla_Y\sigma(\tau, e, Y, y_g)$ is such that $Y(t) \in ]0; 1]^4$, for all $t > 0$. However, in numerical applications, this maximum principle is not guaranteed, as illustrated in Figure 3. Trajectories of the fraction vector in time are represented for a mixture of three SG with parameters 1, with initial state $Y(0) = (0.1, 0.1, 0.1, 0.05)$, $y_g = 0.05$, with specific volume $\tau = 5.5 \cdot 10^{-3}$ and internal energy $e = 10^8$ and $\lambda = (1, 1, 1, 1)$.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Gradient of entropy source term for a three SG mixture. Trajectories of the fractions in time without the correction (dashed lines) and with the correction (plain lines).</figcaption></figure>

Dashed lines correspond to trajectories when no correction is considered. One observes that the volume fraction becomes nonpositive at time $t < 0.0005s$. The correction term $Y(1-Y)$ ensures the maximum principle as illustrated by the plain lines. Of course, this correction term modifies the trajectories, in particular the return to equilibrium seems to be reached at a longer time. We emphasize that it directly depends on the initial state.

**Remark 3.** The source term $\Gamma_2$ involves the difference of phase pressures which acts on the volume fraction equation. Such an expression is rather classical. It can be found in the two-fluid literature, for instance in Baer-Nunziato like model [3] and the derivation of such mechanical transfer term is now well understood [11].

**Proposition 3.** (Eligibility) Let $\lambda_0 = \min \lambda_i$ such as $\lambda_0 > 0$. Then, for any state $(\tau, e, Y, y_g) \in \Omega$, the source term $\Gamma_2$ satisfies the eligibility criterion (24) for all $t > 0$.

*Proof.* Let $(\tau, e, Y, y_g) \in \Omega$ a given state. It holds

$$ 
\begin{aligned}
\Gamma_2(\tau, e, Y, y_g)\nabla_Y\sigma(\tau, e, Y, y_g) &= \lambda Y(1-Y)\nabla_Y\sigma(\tau, e, Y, y_g)\nabla_Y\sigma(\tau, e, Y, y_g) \\
&= \sum_{k=l,g,v} \lambda_k y_k(1-y_k)(\partial_{y_k}\sigma)^2(\tau, e, Y, y_g) \\
&\geq \lambda_0 \sum_{k=l,g,v} y_k(1-y_k)(\partial_{y_k}\sigma)^2(\tau, e, Y, y_g) \\
&\geq \lambda_0 \inf_{k=l,g,v} y_k(1-y_k) \sum_{k=l,g,v} (\partial_{y_k}\sigma)^2(\tau, e, Y, y_g) \\
&\geq \lambda_0 \inf_{k=l,g,v} y_k(1-y_k) ||\nabla_Y\sigma||^2 \\
&\geq 0.
\end{aligned}
$$

□

Contrary to the BGK-like source term, the source term $\Gamma_2$ does not require conditions on different time scales, thus we can choose distinct ones. We observed
that for equal time scales, the fractions evolve towards the equilibrium at the same speed. It would be interesting to be able to estimate the time taken to reach equilibrium for a given initial state, which is always the same for the BGK-like source term.

Unlike the $\Gamma_1$ source term, solving exactly the ODE system is out of reach. However, existence en uniqueness of local solution to Cauchy problem is guaranteed by the Cauchy-Lipschitz theorem.

**Proposition 4.** Consider $(\tau, e) \in (\mathbb{R}_*)^2$. For all $Y_0 \in (\mathbb{R}_*)^4$ such that $(\tau, e, Y_0, y_g) \in \Omega$, there exists a unique maximal solution $(I, Y) \in \mathbb{R}_+ \times ]0; 1]^4$, with $Y: I \to ]0; 1]^4$, of the Cauchy problem

$$ (29) \qquad \begin{cases} \frac{dY}{dt}(t) = \Gamma_2(\tau, e, Y(t), y_g), & t > 0, \\ Y(0) = Y_0. \end{cases} $$

The following purpose is to ensure that the stationary states of the Cauchy problem (29), namely vector $\bar{Y}$ such that $\Gamma_2(\bar{Y}) = 0$, correspond to physically relevant thermodynamical equilibria of the three-phase mixture.

For a given state $(\tau, e) \in (\mathbb{R}_*)^2$, if $\bar{Y} \in ]0; 1]^4$, then $\Gamma_2(\bar{Y}) = 0$ implies that $\nabla_Y \sigma(\tau, e, \bar{Y}, y_g) = 0$. Hence $\bar{Y}$ is such that (18) holds and the thermodynamical equilibrium is a stationary state.

On the other hand, since $\Gamma_2$ is not defined on the border of the cube $[0, 1]^4$, one cannot determine if single phase states $\bar{Y} = 0_{\mathbb{R}^4}$ or $\bar{Y} = 1_{\mathbb{R}^4}$ are stationary states of the dynamical system and if these states are attractive (using standard arguments like the Jacobian linearization). Similarly, a state $Y = (1 - y_g, 0, 0, 0)$ such that $y_v = 0$ could be a stationary state of the dynamical system. The following result rejects this possibility in the case of a mixture of three Stiffened Gases with parameters 1.

**Proposition 5 (Exclusion of the border).** Let $(\tau, e) \in (\mathbb{R}_*)^2$ be a given state. For any initial data $Y(0) \in ]0, 1]^3$, a solution $Y(t)$ of the system (29) cannot converge towards the asymptotic state $Y_0 = (1 - y_g, 0, 0, 0)$.

*Proof.* The study of the vector field around $Y_0 = (1 - y_g, 0, 0, 0)$ allows to exclude any possibility of convergence toward this point. Indeed, one can show that $\lim_{Y \to Y_0} \Gamma_{2,1}(Y) = -\infty$, where $\Gamma_{2,1}$ is the first component of $\Gamma_2$. We refer to [8] for a detailed proof in the case of a two-phase mixture. $\square$

The high nonlinearity of the source term $\Gamma_2$ requires a robust numerical approximation, especially when considering complex EoS such as NASG-CK EoS. The numerical method we consider is depicted in Section 3.3.

### 3. NUMERICAL APPROXIMATION BY A FRACTIONAL STEP METHOD

The fractional step method [34, 4] consists in approximating the convective part and the source terms in two distinct steps. The approximation of the convective part of the model (19) is performed by a Finite Volume method set on regular meshes of size $\Delta x = x_{i+1/2} - x_{i-1/2}$, $i \in \mathbb{Z}$. Time steps are calculated under CFL condition, with the notation $\Delta t^n = t^{n+1} - t^n$. Denoting $W = (\rho Y, \rho u, \rho E)$ the vector of conserved quantities, the first step of the method corresponds to the
approximation of

$$ \begin{cases} \partial_t W + \partial_x \mathcal{F}(W) = 0, \\ W(0, x) = W_0(x), \end{cases} $$

where $\mathcal{F}(W) = uW + pD$, with $D = (0, 0, 1, u)^T$ and $W_0(x)$ is a given initial state. Let $W(t^n, w)$ be approximated by

$$ (30) \qquad W_i^n = \int_{x_{i-1/2}}^{x_{i+1/2}} W(t^n, x)dx. $$

Integrating the convective part over the space-time domain $[x_{i-1/2}, x_{i+1/2}] \times [t^n, t^{n+1,*})$ gives

$$ \int_{x_{i-1/2}}^{x_{i+1/2}} W(t^{n+1,*}, x) - W(t^n, x)dx + \int_{t^n}^{t^{n+1,*}} \mathcal{F}(W(t, x_{i+1/2})) - \mathcal{F}(W(t, x_{i-1/2}))dt = 0. $$

Using the notation (30), the Finite Volume formulation reads

$$ (31) \qquad W_i^{n+1,*} = W_i^n - \frac{\Delta t}{\Delta x} (F_{i+1/2}^n - F_{i-1/2}^n), $$

where the notation $F_{i+1/2}^n$ stands for the two-point approximation of the flux through the interface $x_{i+1/2}$, namely $F_{i+1/2}^n = F(W_i^n, W_{i+1}^n)$. As the second step of the splitting method is concerned, it remains to provide a discretization of the source term with

$$ (32) \qquad \begin{cases} \frac{d}{dt}Y(t) = \Gamma(\tau, e, Y, y_g), \\ \frac{d}{dt}\tilde{W}(t) = 0, \end{cases} $$

where $\tilde{W} = (\rho y_g, \rho, \rho u, \rho E)$. Finally, the numerical approximation $W_i^{n+1}$ is an approximated solution of (32) at time $t = \Delta t$ with the initial condition $W^{n+1,*}$, deduced from the convection step.

### 3.1. Convective part.
Several two-point numerical fluxes have been considered and will be compared in the following section. We focus on standard numerical schemes which are recalled hereafter for the sake of completeness.

#### 3.1.1. Rusanov scheme.
Considering two neighbouring states $W_L$ and $W_R$, the Rusanov [31] flux reads

$$ F(W_L, W_R) = \frac{1}{2} (\mathcal{F}(W_L) + \mathcal{F}(W_R)) - \frac{\max(\Lambda_R, \Lambda_L)}{2} (W_R - W_L), $$

where $\Lambda_k = \Lambda_k(W_k)$ is the spectral radius of the convection matrix $\nabla_W \mathcal{F}(W_k)$, that is

$$ \max(\Lambda_R, \Lambda_L) = \max_{k=L,R} (|u_k \pm c_k|), $$

where the sound speed $c_k$ is given by (20). Every time step $t^n$ is computed according to the CFL constraint: $\Delta t^n = C \frac{\Delta x}{\max(\Lambda_R, \Lambda_L)}$, with $C$ such that $0 < C < \frac{1}{2}$.
### 3.1.2. VFRoe-ncv scheme.
The VFRoe-ncv scheme [7] is an approximation of the Godunov scheme. It consists in considering a linearization of the model. Considering a change of variable for regular solutions

$$W \mapsto V = \Psi^{-1}(W),$$

with $\Psi$ a $C^1$-diffeomorphism, the system is rewritten under a nonconservative form

$$(33) \qquad \partial_t V + A(V) \partial_x V = 0,$$

with $A(V) = (\nabla_V \Psi)^{-1} \nabla_W \mathcal{F}(\nabla_V \Psi)$ the new convection matrix. For two given states $V_L$ and $V_R$, the linearization is performed around the arithmetic mean $\bar{V} = \frac{V_L+V_R}{2}$.

Thus, the numerical flux is defined with the exact solution of the Riemann problem associated to the following linear system

$$\partial_t V + A(\bar{V}) \partial_x V = 0.$$

Let $\mathcal{V}(\frac{x}{t} = 0; V_L, V_R)$ be this solution. We get the numerical flux:

$$F(W_L, W_R) = \mathcal{F}(\Psi(\mathcal{V}(x/t = 0, V_L; V_R))),$$

and the linear problem solution is given by

$$\mathcal{V}(x/t = 0, V_L; V_R) = V_L + \sum_{\bar{\beta}_i \le 0} \bar{\gamma}_i \bar{r}_i,$$

where $(\bar{\beta}_i)_i$ are the eigenvalues of $A(\bar{V})$, $(\bar{r}_i)_i$ its eigenvectors and $(\bar{\gamma}_i)_i$ the coefficients which are solutions of the system $\Delta V = \sum_i \bar{\gamma}_i \bar{r}_i$, with $\Delta V = V_R - V_L$. To avoid the appearance of any non-physical shock, the numerical flux is corrected by an additional diffusive term [15], leading to

$$(34) \qquad F(W_L, W_R) = \mathcal{F}(\Psi(\mathcal{V}(x/t = 0, V_L; V_R))) - \frac{\delta}{2}(W_R - W_L),$$

where

$$(35) \qquad \delta = \begin{cases} \min(-(u+c)_L, (u+c)_R), & \text{if } (u+c)_L < 0 \text{ and } (u+c)_R > 0, \\ \min(-(u-c)_L, (u-c)_R), & \text{if } (u-c)_L < 0 \text{ and } (u-c)_R > 0, \\ 0, & \text{else.} \end{cases}$$

### 3.1.3. Relaxation scheme.
The principle of relaxation schemes is to enlarge the system and to relax the nonlinearities (coming from the complex pressure and energy terms), in a way such as all the characteristic fields are linearly degenerate. Doing so, the computations of the solution are explicit, and a proper relaxation term allows to recover the initial system solutions. Following for instance [9, Theorem 6], the enlarged system reads

$$(36) \qquad \left\{ \begin{aligned} & \partial_t(\rho Y) + \partial_x(\rho Y u) = 0, \\ & \partial_t(\rho y_g) + \partial_x(\rho y_g u) = 0, \\ & \partial_t \rho + \partial_x(\rho u) = 0, \\ & \partial_t(\rho u) + \partial_x(\rho u^2 + \Pi) = 0, \\ & \partial_t(\rho \Sigma) + \partial_x(\rho u \Sigma + u \Pi) = 0, \\ & \partial_t(\rho T) + \partial_x(\rho T u) = \frac{1}{\epsilon} \rho (\tau - T). \end{aligned} \right.$$

Here the quantity $T$ is such that for all $x$, $T(0, x) = \tau(0, x)$ and $\epsilon$ is a relaxation parameter. The pressure term is relaxed in the following sense

$$(37) \qquad \Pi = p(T, e, Y, y_g) + a^2(T - \tau),$$
while the relaxed energy reads

$$ (38) \qquad \Sigma = e + \frac{u^2}{2} + \frac{\Pi^2 - (p(\mathcal{T}, e, Y, y_g))^2}{2a^2}, $$

with a a parameter chosen such that

$$ (39) \qquad a > \max \left( \frac{c_l(\tau_l, e_l, Y_l, y_g)}{\tau_l}, \frac{c_r(\tau_r, e_r, Y_r, y_g)}{\tau_r} \right). $$

This above condition, referred as Whitham condition, ensures the hyperbolicity of the relaxed system and well-ordered waves. Finally, considering two states $W_L$ and $W_R$ at a given interface, the two-point numerical flux is

$$ (40) \qquad F(W_L, W_R) = \left( \frac{Yu}{\tau}, \frac{u}{\tau}, \frac{u^2}{\tau} + \Pi, \frac{u\Sigma}{\tau} + u\Pi \right), $$

where $Y, \tau, u, \Sigma$ are the components of the solution of the enlarged Riemann problem at the considered interface.

## 3.2. Computation of the internal energy.
The implementation of the above numerical fluxes requires to handle transits from conservative variables $W$ to observable quantities e.g. pressure and temperature fields and vice versa. Depending on the considered equations of state, these changes of variables can be non-explicit and require most of the time an implicit solver. The aim of this section is to provide detailed algorithms for the computation of the internal energy $e$ of the mixture, starting from a mixture state with specific volume $\bar{\tau}$, fraction vector $\bar{Y}$ and pressure $\bar{p}$. Similar algorithms can be adapted when considering an initial temperature $\bar{T}$.

Depending on the equations of state, the algorithms differ. We focus on two configurations: first, all the phases are depicted by Stiffened Gas laws, second, the gaseous phases follow SG laws while the liquid is modeled by a NASG-CK law.

### 3.2.1. Computation for a 3 SG mixture.
Let $(\bar{\tau}, \bar{p}, \bar{Y})$ be a given state. The question is to compute the associated internal energy $e$, which is the root of the function $\mathcal{P}$, defined by

$$ (41) \qquad \mathcal{P}(e) = p(\bar{\tau}, e, \bar{Y}, y_g) - \bar{p}. $$

Using the mixture pressure definition (16) and the expressions of the Stiffened Gas pressure and temperature, this function can be written as a third-degree polynomial

$$ (42) \qquad \mathcal{P}(e) = \bar{A}e^3 + \bar{B}e^2 + \bar{C}e + \bar{D}. $$
We precise here the exact coefficients, for $k, k', k'' \in \{l, g, v\}$ pairwise distinct (note that they slightly differ from those of [20]):

$$
\begin{align*}
\bar{A} &= \sum_{k,k',k''} z_k y_k (\gamma_k - 1) \frac{z_{k'}}{C_{k'}} \frac{z_{k''}}{C_{k''}} \frac{1}{\tau}, \\
\bar{B} &= \sum_{k,k',k''} -\alpha_k y_k \frac{1}{C_{k'} C_{k''}} \left( z_{k'} z_{k''} \left( \Pi_k \gamma_k + \frac{y_k Q_k (\gamma_k - 1)}{\alpha_k \tau} \right) \right. \\
&\qquad \left. + \frac{z_k (\gamma_k - 1)}{\alpha_k} \frac{1}{\tau} (z_{k'} \omega_{k''} + z_{k''} \omega_{k'}) \right) - z_k y_k \frac{z_{k'} z_{k''}}{C_{k'} C_{k''}}, \\
\bar{C} &= \sum_{k,k',k''} \alpha_k y_k \frac{1}{C_{k'} C_{k''}} \left( \frac{z_k (\gamma_k - 1)}{y_k} \frac{1}{\tau} \omega_{k'} \omega_{k''} \right. \\
&\qquad \left. + \left( \frac{y_k Q_k (\gamma_k - 1)}{\alpha_k} \frac{1}{\tau} + \Pi_k \gamma_k \right) (z_{k'} \omega_{k''} + z_{k''} \omega_{k'}) \right) \\
&\qquad + \left( z_k y_k \frac{1}{C_{k'} C_{k''}} (z_{k'} \omega_{k''} + z_{k''} \omega_{k'}) \right), \\
\bar{D} &= \sum_{k,k',k''} -\alpha_k y_k \frac{1}{C_{k'} C_{k''}} \left( \frac{y_k Q_k (\gamma_k - 1)}{\alpha_k} \frac{1}{\tau} + \Pi_k \gamma_k \right) \omega_{k'} \omega_{k''} \\
&\qquad - y_k z_k \frac{1}{C_{k'} C_{k''}} \omega_{k'} \omega_{k''},
\end{align*}
$$

where $\omega_k = y_k Q_k + \alpha_k \Pi_k \tau$. Remark that the coefficients $\bar{A}$ to $\bar{D}$ are explicitly defined by the data $(\bar{\tau}, \bar{p}, \bar{Y})$, the notation $\bar{x}$ being omitted for sake of readability. The algorithm reads as follows.

* For a given state $(\bar{\tau}, \bar{p}, \bar{Y})$, compute the coefficients (43) of $\mathcal{P}$.

* Compute the discriminant $\Delta$ of the derivative of $\mathcal{P}$.

I) If $\Delta \le 0$, the equation (42) admits a unique solution, which is computed with a Newton's method.

II) If $\Delta > 0$, let $X_1, X_2$ be the ordered roots of $\mathcal{P}'$, such as $X_1 < X_2$. They correspond to the local extrema of $\mathcal{P}$, see Figure 4.

Let $\delta = 10^9$ and $\gamma = 100$ two parameters.

i) If $\mathcal{P}(X_2)\mathcal{P}(X_2 + \gamma) \le 0$, there exists a unique root $e_0$ such as $e_0 > X_2$, see Figure 4-right. We make use of a Newton method with initial condition $X_2 + \delta$ and verify that the solution is greater than $X_2$. In case of non convergence or if the solution is lower than $X_2$, a bisection method is launched in the range $[X_2, X_2 + \gamma]$ to ensure to find the correct root.

ii) If $\mathcal{P}(X_2)\mathcal{P}(X_2 + \gamma) > 0$, $\mathcal{P}$ admits only one root $e_0$ such as $e_0 < X_1$, see Figure 4-left. Indeed, we can't have $X_1 < e_0 < X_2$ because it wouldn't be a unique root. Then we proceed like previously with a Newton method starting from $X_2 - \delta$.

The parameters $\delta$ and $\gamma$ are set arbitrarily and may be modified to optimize the algorithm. Note that the case II)-ii) excludes the possibility of $\mathcal{P}(X_i) = 0, i = 1$ or 2. Indeed, this situation would correspond to the existence of two roots which have not been encountered in practical simulations.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Computation of the internal energy for a mixture of three SG. Typical graph of $\mathcal{P}$ in the two cases where $\Delta > 0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Computation of the internal energy for a mixture of two SG and a NASGK-CK liquid. Typical graph of $\mathcal{P}$.</figcaption></figure>

### 3.2.2. Computation for a mixture of 2 SG and a NASG-CK liquid phase.
In this configuration, the expression $\mathcal{P}(e) = P(\bar{\tau}, e, \bar{Y}) - \bar{p}$ is no more explicit due to the NASG-CK EoS. Computing the root of $\mathcal{P}(e)$ is tedious because the graph of $\mathcal{P}$ is basically a hyperbola, see Figure 5. To handle this difficulty, a bisection method is preferred (compared to a secant or a Newton method to ensure to find the root in the desired range). The algorithm reads as follows.

* Consider a range of internal energy $[e_{\min}, e_{\max}]$.

* Provide a discretization of $[e_{\min}, e_{\max}]$ with $e_j = e_{\min} + j\Delta e$, $\Delta e = (e_{\max} - e_{\min})/N$.

* Set $\delta_j = \frac{e_j - e_{j-1}}{d}$.

* For $j = N, N-1, ..., 2$, find the first index $j_0$ such as $\mathcal{P}(e_{j_0-1})\mathcal{P}(e_{j_0}) < 0$.

* Launch a bisection algorithm on $[e_{j_0-1}, e_{j_0}]$.

* If the bisection fails, there must exist a real $\tilde{e} \in [e_{j_0-1}, e_{j_0}]$ where $\mathcal{P}$ is undefined. Hence, two cases must be considered.

1. If the previous iterations indicate that the root $\bar{e}$ is in the neighbourhood of $e_{j_0-1}$ (resp. $e_{j_0}$), relaunch the bisection on the ranges $[e_{j_0-1}, e_{j_0-1}+\delta_j]$, $[e_{j_0-1}+\delta_j, e_{j_0-1}+2\delta_j]$, ..., until convergence is reached (resp. $[e_{j_0}, e_{j_0}-\delta_j]$, $[e_{j_0}-\delta_j, e_{j_0}-2\delta_j]$, ...).

2. Else the last iteration $e_\infty$ provides an approximation of $\bar{e}$. If $e_\infty < \bar{e}$ (resp. $e_\infty > \bar{e}$), relaunch the bisection on the range $[e_{j_0-1}, e_\infty]$ (resp. on $[e_\infty, e_{j_0}]$).

In practice, we choose $d = 10$ and $N = 10^3$. Higher values of $N$ can be necessary on simulations at high pressure, see Section 4, Figure 11.
**Remark 4.** In order to reduce costs of this algorithm, a reminder can be implemented to launch first a bisection method around the solution found on the previous cell if it detects two similar consecutive states. Let $e_{old}$ and $e_{new}$ be the values of the internal energies on consecutive cells. At the beginning of the next computation of the energy, if $|e_{old} - e_{new}| < 1$, we first launch a bisection method on the interval $[e_{new} - 10; e_{new} + 10]$. If the solution is negative or if the bisection algorithm does not converge, it goes back to previous algorithm.

3.3. **Source term approximation.** The second step of the splitting strategy corresponds to the source term approximation.

As BGK-like source term is concerned, trajectories can be explicitly defined, see Proposition 2, as soon as the equilibrium vector $Y_{eq}$ is determined. Two difficulties arise. First, according to Remark 1, one could encounter numerical difficulties to compute the equilibrium. Moreover, when complex EoS are considered, the maximization process of the mixture entropy (14) under the constraints (11)-(12) has to be performed with an appropriate optimization algorithm. We refer to [24] which proposes a Broyden method to compute the equilibrium for a mixture of two stiffened gases and a NASG-Chemkin law for the liquid phase.

We focus in this section on the gradient of entropy source term $\Gamma_2$.

In order to capture the long-time asymptotic equilibria, an implicit approximation is considered. However, the complexity of the source term $\Gamma_2$ makes a full implicit discretization out of reach, especially when considering NASG-CK EoS. We propose hereafter a semi-implicit Euler algorithm, based on a Broyden method [6] to evaluate the implicit part.

On the time range $[t^{n+1,*}, t^{n+1}]$, one considers the Cauchy problem (32) with initial condition $W^{n+1,*}$ computed from the convective step. The implicit Euler method boils down to find the root of the nonlinear functional $G: [0, 1]^4 \to \mathbb{R}$

$$ (44) \qquad G(Y) = Y - \Delta t Y(1 - Y) \nabla_Y \sigma(\tau_{n+1;*}, e_{n+1;*}, Y, y_g) - Y_{n+1;*}, $$

To do so, we make use of an iterative Broyden algorithm, which requires to compute the Jacobian matrix $J(Y)$ of $G(Y)$ with coefficients

$$ (45) \qquad (J(Y))_{i,j} = \partial_{Y_j} Y_i - \Delta t \lambda_i \partial_{Y_j} [Y_i(1 - Y_i) \partial_{Y_i} \sigma(\tau_{n+1;*}, e_{n+1;*}, Y, y_g)], \quad 1 \le i, j \le 4, $$

where $Y_i$ (resp. $\lambda_i$) stands for the $i$-th coefficient of $Y$ (resp. $\lambda$). One remarks that only the initialization step requires a Jacobian computation. The algorithm consists in computing the sequence $(Y_k)_{k \in \mathbb{N}}$ given by

$$ (46) \qquad \left\{ \begin{array}{l} Y_k = Y_{k-1} - J_{k-1}^{-1} G(Y_{k-1}), \\ J_k^{-1} = J_{k-1}^{-1} + \frac{\Delta Y_k - J_{k-1}^{-1} \Delta G_k}{\Delta Y_k^T J_{k-1}^{-1} \Delta G_k} (\Delta Y_k)^T J_{k-1}^{-1}, \\ \Delta Y_k = Y_k - Y_{k-1}, \\ \Delta G_k = G(Y_k) - G(Y_{k-1}), \end{array} \right. $$

until $\|Y_k - Y_{k-1}\| < \eta$, where $\eta$ is a tolerance ($\eta = 10^{-10}$ in practice). The initialization requires to choose $Y_0$ to $Y_n$ given by the previous time step computation.

## 4. NUMERICAL RESULTS

This Section provides two types of numerical experiments. First verification test cases assess the good approximation of the convective part of the model (19). They consist in the approximation of exact Riemann problem solutions at low and high
pressure for three-phase mixture with either SG or NASG-CK EoS (for the liquid phase).

A comparison of CPU costs is presented in Section 4.1.3, highlighting the performances of the relaxation scheme compared to the VFRoe-ncv scheme.

Finally, a global simulation is presented, accounting for the contribution of the gradient of entropy source term.

**4.1. Infinitely slow relaxation.** These test cases correspond to the verification of the approximation of the convective part (when the source terms are inactive). For two different sets of EoS, we compare the numerical results of the different numerical fluxes presented in Section 3.1 on Riemann problems at low and high pressure. In both cases, the Riemann problem solutions are built according to Section 2.2.2. Considering a left state $W_L$, the state $W_1$ is computed considering that the first nonlinear wave is a ghost wave. Then the contact wave allows to deduce the state $W_2$, the pressure and velocity being preserved through the contact contrary to $Y$. Finally a shock wave of velocity $\xi$ separates the state $W_2$ and $W_R$ while the fractions $Y$ are preserved.

**4.1.1. Mixture of three Stiffened Gases.** Tables 4 and 5 present the exact Riemann data in the low and high pressure case respectively when considering that the three phases are depicted by Stiffened Gas laws with parameters 1. The four constant states are depicted in terms of their primitive variables, namely $Z = (y_l, y_g, \alpha_l, z_l, z_g, \tau, u, p)$.

<table><caption>TABLE 4. Mixture of three SG. Data for the low-pressure Riemann problem. $u + c$-shock of speed $\xi = 708.41563589800251$.</caption><thead><tr><th>Z<sub>L</sub> = Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>R</sub></th></tr></thead><tbody><tr><td>0.91826194646564185</td><td>0.89784395822033747</td><td>0.89784395822033747</td></tr><tr><td>5.0000000000000003×10<sup>-2</sup></td><td>5.99999999999994613×10<sup>-2</sup></td><td>5.99999999999994613×10<sup>-2</sup></td></tr><tr><td>1.6767800121816667×10<sup>-2</sup></td><td>1.0000000000001008×10<sup>-2</sup></td><td>1.0000000000001008×10<sup>-2</sup></td></tr><tr><td>0.41969031490727782</td><td>0.30000000000002064</td><td>0.30000000000002064</td></tr><tr><td>8.8350466782145878×10<sup>-2</sup></td><td>9.8999999999997021×10<sup>-2</sup></td><td>9.8999999999997021×10<sup>-2</sup></td></tr><tr><td>5.4427805491776261×10<sup>-2</sup></td><td>6.8863503931918257×10<sup>-2</sup></td><td>7.000000000000007×10<sup>-2</sup></td></tr><tr><td>10.0</td><td>10.0</td><td>-1.5263757834575244</td></tr><tr><td>1.8668002620104561×10<sup>5</sup></td><td>1.8668002620104561×10<sup>5</sup></td><td>6.9779191756698478×10<sup>4</sup></td></tr></tbody></table>

<table><caption>TABLE 5. Mixture of three SG. Data for the low-pressure Riemann problem. $u + c$-shock of speed $\xi = 578.35007234359045$.</caption><thead><tr><th>Z<sub>L</sub> = Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>R</sub></th></tr></thead><tbody><tr><td>0.39999999999999997</td><td>0.3799999999999906</td><td>0.3799999999999906</td></tr><tr><td>0.29999999999999999</td><td>0.3300000000000052</td><td>0.3300000000000052</td></tr><tr><td>8.0000000000000071×10<sup>-2</sup></td><td>6.999999999999840×10<sup>-2</sup></td><td>6.999999999999840×10<sup>-2</sup></td></tr><tr><td>0.25000000000000011</td><td>0.2199999999999936</td><td>0.2199999999999936</td></tr><tr><td>9.9999999999999978×10<sup>-2</sup></td><td>0.12000000000000005</td><td>0.12000000000000005</td></tr><tr><td>6.0000000000000010×10<sup>-3</sup></td><td>5.0000000000000062×10<sup>-3</sup></td><td>5.3000000000000000×10<sup>-3</sup></td></tr><tr><td>100.0</td><td>100.0</td><td>71.298995659385199</td></tr><tr><td>1.4999999999999996×10<sup>7</sup></td><td>1.499999999999996×10<sup>7</sup></td><td>1.2254174499466641×10<sup>7</sup></td></tr></tbody></table>
The three numerical fluxes are compared with the same CFL parameter set to 0.45, on several test cases that are set on the domain $[0, 1]$. Figures 7-8-10-11 were obtained with 500 cells. The convergence results, obtained for a final time $t = 10^{-3}$s., see Figure 6, show that the VFRoe-ncv and the relaxation schemes achieve expected convergence rates of 1/2 that is optimal for typical shock tube computations involving contact discontinuities.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Mixture of three Stiffened Gases. $L^1$-convergence for the low and high-pressure Riemann problems.</figcaption></figure>

Figure 7 presents the profiles of density, velocity, pressure and fractions obtained by the relaxation scheme in the low-pressure configuration. One observes the good approximations of the fractions. Some oscillations are present in the pressure and velocity profiles and are probably due to the nonlinear EoS. This phenomenon is a well known problem that is not specific to the multiphase context. It has been identified in [1] and possible consequences have been investigated in the single phase context in [33]. These oscillations also appear in the high-pressure configuration, see Figure 8, while the fraction profiles are still well captured.

4.1.2. *Mixture of two Stiffened Gases with a NASG-CK liquid phase.* Tables 6 and 7 present the exact Riemann data in the low and high pressure cases respectively when considering that the liquid phase is depicted by a NASG-CK law with parameters 3.

The convergence results in Figure 9 are obtained with the same discretization parameters and CFL condition as in paragraph 4.1.1. One observes that the VFRoe-ncv and the relaxation schemes achieve similar convergence rates.

The profiles of density, velocity, pressure, temperature and fractions are given in Figure 10 (resp. Figure 11) for the low-pressure (resp. high) configuration with the relaxation scheme. Similar observations can be done. The fractions are well captured whereas oscillations, due to the nonlinearity of the pressure law (with respect to the internal energy), arise. Note that considering the implicit NASG-CK EoS increases the amplitude of the oscillations, at least for the low-pressure test cases.

4.1.3. *Comparison between VFRoe-ncv and the relaxation scheme.* Although the previous section shows that the convergence for the VFRoe-ncv and relaxation schemes are similar, some remarks can be made on their performances.
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Mixture of three Stiffened Gases. Low-pressure test case and comparison with the exact solution. From top left to bottom right, profiles of density, velocity, pressure, temperature, $y_l$, $\alpha_l$, $z_l$ and $z_g$.</figcaption></figure>

As CPU time cost is concerned, a comparison of the two scheme performances is provided in Figure 12 for the mixture of three SG and a mixture of two SG and a NASG-CK liquid. Obviously the VFRoe-ncv scheme is much more CPU time-consuming. This is due to the fact that it requires more internal energy computations. These are done using the algorithms provided in Sections 3.2.1 and 3.2.2. On the other hand, the relaxation scheme only requires two internal energy
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Mixture of three Stiffened Gases. High-pressure test case and comparison with the exact solution. From top left to bottom right, profiles of density, velocity, pressure, temperature, $y_l$, $\alpha_l$, $z_l$ and $z_g$.</figcaption></figure>

computations: one at the beginning of each time step to go from conservative to primitive variables, the second at the end, for the inverse change of variables.

One can observe on Figure 12 that for a given number of cells, the VFRoe-ncv scheme needs on average two times more CPU time than the relaxation scheme, in the two SG mixture case.
<table><caption>TABLE 6. Mixture of two Stiffened Gases with a NASG-CK liquid phase. Data for the low-pressure Riemann problem. $u + c$-shock of speed $\xi = 1027.3711820835272$.</caption><thead><tr><th>Z<sub>L</sub> = Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>R</sub></th></tr></thead><tbody><tr><td>0.91841879730033971</td><td>0.89999999102958528</td><td>0.89999999102958528</td></tr><tr><td>5.0000000000000003E-002</td><td>5.5000004933728031E-002</td><td>5.5000004933728031E-002</td></tr><tr><td>1.6376791114071865E-002</td><td>1.4999998527357117E-002</td><td>1.4999998527357117E-002</td></tr><tr><td>0.92743215703897408</td><td>0.90000000855046036</td><td>0.90000000855046036</td></tr><tr><td>1.1475280368993632E-002</td><td>9.9999991449539437E-003</td><td>9.9999991449539437E-003</td></tr><tr><td>0.19492558898696019</td><td>0.19500001720076998</td><td>0.1970000000000001</td></tr><tr><td>100.0</td><td>100.0</td><td>90.488583338128308</td></tr><tr><td>2.4530790395664933×10<sup>5</sup></td><td>2.4530790395664933×10<sup>5</sup></td><td>2.0007399146972221×10<sup>5</sup></td></tr></tbody></table>

<table><caption>TABLE 7. Mixture of two Stiffened Gases with a NASG-CK liquid phase. Data for the high-pressure Riemann problem. $u + c$-shock of speed $\xi = 1591.8798489453186$.</caption><thead><tr><th>Z<sub>L</sub> = Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>Z<sub>R</sub></th></tr></thead><tbody><tr><td>0.92020382010222246</td><td>0.89999998651866320</td><td>0.89999998651866320</td></tr><tr><td>5.0000000000000003×10<sup>-2</sup></td><td>5.3000007145108524×10<sup>-2</sup></td><td>5.3000007145108524×10<sup>-2</sup></td></tr><tr><td>0.68595071360084603</td><td>0.59999996404976996</td><td>0.59999996404976996</td></tr><tr><td>0.95997763551798010</td><td>0.90000000093700905</td><td>0.90000000093700905</td></tr><tr><td>6.5876209523389063×10<sup>-3</sup></td><td>7.9999999250392681×10<sup>-3</sup></td><td>7.9999999250392681×10<sup>-3</sup></td></tr><tr><td>4.8078440722228851×10<sup>-3</sup></td><td>5.0000002246889472×10<sup>-3</sup></td><td>5.069999999999999×10<sup>-3</sup></td></tr><tr><td>100.0</td><td>100.0</td><td>79.113750095130015</td></tr><tr><td>1.7275285041205704×10<sup>7</sup></td><td>1.7275285041205704×10<sup>7</sup></td><td>1.1043330250633653×10<sup>7</sup></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 9. Mixture of two Stiffened Gases with a NASG-CK liquid phase. $L^1$-convergence for the low and high-pressure Riemann problems.</figcaption></figure>

In the case of a NASG-CK liquid, the gap is obviously wider and the factor between the two schemes varies from 20 (for finer meshes) to 100 (for coarser meshes). These variations are certainly due to the internal energy algorithm that struggles
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Mixture of two Stiffened Gases and a NASG-CK liquid. Low-pressure test case and comparison with the exact solution. From top left to bottom right, profiles of density, velocity, pressure, temperature, $y_l$, $\alpha_l$, $z_l$ and $z_g$.</figcaption></figure>

more on coarse meshes because of the large variations between two consecutive states, see Section 3.2.2.

Considering that the EOS computation is the crucial point of the CPU costs, Figure 13 represents the CPU time with respect to the number of EOS evaluations. For a given number of evaluations, the 2 SG + 1 NASG-CK law configuration requires approximatively twenty times more time than the 3 SG case.
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Mixture of two Stiffened Gases and a NASG-CK liquid. High-pressure test case and comparison with the exact solution. From top left to bottom right, profiles of density, velocity, pressure, temperature, $y_l$, $\alpha_l$, $z_l$ and $z_g$.</figcaption></figure>

4.2. **Global out-of-equilibrium simulation.** This final section provides numerical results of the fractional step method, using the relaxation scheme for the convective part and the Broyden algorithm for the source terms. The two source terms of Section 2.3 are compared and several (constant) relaxation time $\lambda$ as well. We focus on a mixture of three SG in the low-pressure configuration. The initial condition is $Z(0, x) = Z_L \mathbf{1}_{x<0.5} + Z_R \mathbf{1}_{x>0.5}$ with $Z_L$ and $Z_R$ given in Table 4. The
<figure><img src="image_14.png" /><figcaption>FIGURE 12. Comparison of CPU time costs with respect to the number of cells for the VFRoe-ncv and the relaxation schemes on the low pressure Riemann problem.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 13. Comparison of CPU time costs with respect to the number of equation of state evaluations for the low pressure Riemann problems solved with the relaxation scheme.</figcaption></figure>

computational domain $[0, 1]$ contains 1000 cells and the CFL is again 0.45. The final time is $t = 10^{-4}$. The source term $\Gamma_1$ requires relaxation time $\lambda$ lower than $5.10^3$ to avoid the appearance of non-admissible states.

One observes on Figure 14 an expected behaviour of the solutions with respect to the time scales $\lambda$. Indeed, the monotonicity property of the solutions of the dynamical system (22) is verified. The more $\lambda$ grows, the more the solution goes away from the out-of-equilibrium solution, which corresponds to the absence of
source terms ($\lambda = 0$). This is mathematically and physically relevant. However, the upper bound on $\lambda$ prevents us from obtaining solutions closer to the equilibrium. It is due to non-admissible states that can appear around the contact discontinuity, as one can see on the velocity plot, which seems to be the most sensitive variable on this test case. Finally, the left and right constant states on the fractions allow us to estimate the "distance" to the equilibrium. At final time $T = 0.001$ and for the value $\lambda = 5.10^3$, the system has done approximatively 40% of the path to the equilibrium. It is consistent with the previous observations done in Section 2.3.1.

As $\Gamma_2$ is concerned, Figure 15 illustrates the non-monotonic solutions for this source term, given in Section 2.3.2. For instance, one can observe intersections between the velocity profiles for $\lambda = 1$ and $\lambda = 1.1$. Because the source term is highly nonlinear, it is difficult, and even out of reach, to give quantitative informations about the possible trajectories. By comparing these profiles to the BGK source terms $\Gamma_1$, it clearly appears that the entropy gradient source term provides very complex intermediates states. The behaviors of the density, velocity and pressure profiles around the contact discontinuity are disorganized. In particular, it seems to have consequences on the right constant states of the pressure (for any values of $\lambda$), that are far from those of the BGK source term case. Finally, the left and right constant states of $y_l$ illustrate that time taken to reach equilibrium depends on the initial state (by supposing that the evolution of these latter are monotone, that seems true on this case).

## 5. CONCLUSION

This study is focused on the numerical simulation of the so-called homogeneous relaxation model for a three-phase flow with miscible constraints. As the numerical approximation is concerned, the work on the convective part leads to a preference for the relaxation scheme for CPU time and robustness reasons, although convergence results are similar for the relaxation and the VFRoe-ncv schemes. The computation of the internal energy plays a major role in the CPU time cost and robustness, that is why algorithms are detailed since they involve tricky computations. Secondly, as the thermodynamical modelling is concerned, both the thermodynamical equilibrium of the three-phase mixture and the relaxation towards this equilibrium are analyzed. The good definition of the equilibrium fractions and its constraints are fundamental features. We notably have addressed the question of disappearance of phases and the associated equilibria. When considering the time relaxation towards the thermodynamical equilibrium, two types of source terms have been compared. The BGK type one provides an efficient simulation, as soon as the computation of the equilibrium is reachable, which is a strong constraint. Moreover, it requires that relaxation time scales must be equal, which is again a strong assumption. On the other hand, the gradient type source term is more versatile, despite a far more complex analysis, that leads to numerical difficulties. These pros and cons are highlighted on a global simulation in the last part.

## Acknowledgments
<figure><img src="image_16.png" /><figcaption>FIGURE 14. Global simulation for a mixture of three SG in the low-pressure configuration with the BGK-like source term $\Gamma_1$. Comparison of the impact of the relaxation time. From top left to bottom right, profiles of density, velocity, pressure, $y_l$, $\alpha_l$. The plot for $\lambda = 1.10^4$ is not represented for the density and the volume fraction for the sake of clarity.</figcaption></figure>

This work has received the financial support from the ANR MoHyCon and the CNRS project Needs M2SIR. Part of the work has benefited from fruitful discussions with Mr Hérard from EDF Lab Chatou.

## APPENDIX A. EQUILIBRIUM COMPUTATION FOR 3 STIFFENED GAS

Let $(\tau, e, y_g)$ be a given state. We assume that the equilibrium $Y_{eq}$ is reached in the interior of $[0, 1]^4$, such that it is characterized by (18). Following [2], one characterizes the equilibrium considering first the thermal equilibrium, then the mechanical one and, to finish, the chemical potential equalities. According to the definition (5) of the Stiffened Gas temperature, the sum of the phase temperatures, weighted by the mass fractions $y_k$, reads

$$ (47) \qquad \sum_{k=l,g,v} y_k C_k T_k = e - \sum_k y_k Q_k - \sum_k \alpha_k \Pi_k \tau. $$
<figure><img src="image_17.png" /><figcaption>FIGURE 15. Global simulation for a mixture of three SG in the low-pressure configuration with the gradient entropy source term $\Gamma_2$. Comparison of the impact of the relaxation time. From top left to bottom right, profiles of density, velocity, pressure, $y_l$, $\alpha_l$.</figcaption></figure>

The thermal equilibrium $T_l = T_g = T_v = \bar{T}$ gives then

$$ (48) \qquad \bar{T}(y_v, \alpha_v) = \frac{e - Q(y_v) - \Pi(\alpha_v)\tau}{C(y_v)}, $$

where $C(y_v) = \sum_k y_k C_k$, $Q(y_v) = \sum_k y_k Q_k$ and $\Pi(\alpha_v) = \sum_k \alpha_k \Pi_k$. Using the intensive constraints (11)-(12), one observes that the coefficients $C(y_v)$ and $Q(y_v)$ are solely functions of the vapor mass fraction $y_v$, while $\Pi(\alpha_v)$ depends only on the volume fraction $\alpha_v$. Hence, the equilibrium temperature $\bar{T}(y_v, \alpha_v)$ does not depend on the energy fractions. This allows to split the computations in two steps: first determine the pair $(y_v, \alpha_v)$ and then, determine the energy fractions $(z_l, z_g)$.

We now turn to the mechanical equilibrium to express $\alpha_v$ as a function of $y_v$.

Indeed, the mechanical equilibrium corresponds to $p_l - (p_g + p_v) = 0$ or equivalently to $\alpha_v(1 - \alpha_v)[p_l - (p_g + p_v)] = 0$. It turns out that $F_{y_v}(\alpha_v) = \alpha_v(1 - \alpha_v)[p_l - (p_g + p_v)]$ is a second order polynomial in $\alpha_v$ which reads

$$ (49) \qquad F_{y_v}(\alpha_v) = A_F \alpha_v^2 + B_F \alpha_v + C_F, $$
with coefficients

$$
\begin{align}
A_F &= (\beta + 1)\Pi_0, \\
B_F &= -\Pi_0(\beta_v + \beta_g + 1) - \beta(\rho(e - Q) - \Pi_l), \\
C_F &= (\beta_v + \beta_g)(\rho(e - Q) - \Pi_l), \\
(50) \qquad \Pi_0 &= \Pi_g + \Pi_v - \Pi_l, \\
\beta_k &= (\gamma_k - 1)C_k y_k, \\
\beta &= \sum_k \beta_k.
\end{align}
$$

It remains to determine a root $\alpha_v$ of $F_{y_v}$, which only depends on $y_v$, according to the coefficient definitions (50). In order to assess the nonnegativity of the mixture temperature $\bar{T}(y_v, \alpha_v)$, this volume fraction $\alpha_v$ should comply with some constraints. Considering that $\Pi_g = 0$, $\Pi_v < 0$ and $\Pi_v - \Pi_l < 0$ (which is the case for considered data 1), the temperature positivity is equivalent to $\rho(e - Q(y_v) - \Pi(y_v)) \ge 0$, that is

$$
\alpha_v \geq \frac{\rho(e - Q) - \Pi_l}{\Pi_v - \Pi_l} =: \delta(y_v).
$$

Hence one must seek for $\alpha_v \in [\alpha_{v,inf}; 1]$, with $\alpha_{v,inf} = \max(0, \min(1, \delta(y_v))) \in [0; 1]$.

Generally, the equation $F_{y_v}(x) = 0$ does not necessarily admit a solution on this interval. If this happens, we must determine which limit case $\alpha_v = \alpha_{v,inf}$ or $\alpha_v = 1$ maximises the entropy.

Thanks to the concavity of $\alpha_v \mapsto \sigma$, one only has to study the sign of $\frac{\partial \sigma}{\partial \alpha_v}$ on the
border, and that sign is the same as the sign of $\Delta p = p_v + p_g - p_l$. Thus we have:

* If $\lim_{\alpha_v \to \alpha_{v,inf}} \Delta p(\alpha_v) < 0$, $\alpha_v = \alpha_{v,inf}$ corresponds to a stable state,

* If $\lim_{\alpha_v \to 1} \Delta p(\alpha_v) > 0$, $\alpha_v = 1$ corresponds to a stable state with only gases.

In practical simulations, we choose to solve first $F_{y_v}(x) = 0$, then verify if we
have an unique solution $\alpha_v$ in $[0; 1]$ such as $\bar{T}(y_v, \alpha_v) \ge 0$, and if not we consider
the limit cases.

At last, the chemical equilibrium $\mu_l = \mu_v$ allows us to determine $(\bar{y}_v, \bar{\alpha}_v)$ in
an implicit way, thanks to $F_{y_v}$. We can already define $y_{v,max} \in [0; 1 - y_g]$ as the
greatest mass fraction of vapor such as the temperature is positive.

Then we must solve the following equation on $[0; y_{v,max}]$:

$$
(51) \quad \Delta\mu(y_v) = \Delta\mu(y_v, \alpha_v(y_v)) = \mu_l(y_v, \alpha_v(y_v)) - \mu_v(y_v, \alpha_v(y_v)) = 0,
$$

where

$$
\left\{
\begin{array}{l}
\mu_l(y_v, \alpha_v(y_v)) = Q_l + C_l + \bar{T}(y_v, \alpha_v) \left( \gamma_l - \ln(C_l T) - \frac{s_l^0}{C_l} - (\gamma_l - 1) \ln \left( \frac{1-\alpha_v \tau}{1-y_v - y_g} \right) \right) \\
\mu_v(y_v, \alpha_v(y_v)) = Q_v + C_v + \bar{T}(y_v, \alpha_v) \left( \gamma_v - \ln(C_v T) - \frac{s_v^0}{C_v} - (\gamma_v - 1) \ln \left( \frac{\alpha_v \tau}{y_v} \right) \right)
\end{array}
\right.
$$

and $\alpha_v(y_v)$ is the solution of $F_{y_v}(x) = 0$. As previously, it is possible not to find
$y_v \in [0; y_{v,max}]$ and $\alpha_v$ solutions of this equation. This time we use the concavity
of $y_v \to \sigma$ and the fact that the sign of $\frac{\partial \sigma}{\partial y_v}$ is the same one as the sign of $\Delta \mu$ to
study the limit cases:

* If $\lim_{y_v \to 0} \Delta\mu(y_v) < 0$, $y_v = 0$ corresponds to a stable state, with only liquid and inert gas, i.e. $y_l = 1 - y_g$.
• If $\lim_{y_v \to y_{v,max}} \Delta\mu(y_v) > 0$, $y_v = y_{v,max}$ corresponds to a stable state with only gases, i.e. $y_l = 0$.

Finally, once $(\bar{y}_v, \bar{\alpha}_v)$ is determined, we can compute $(\bar{z}_v, \bar{z}_g)$ with the following relations:

$$ (52) \qquad \bar{z}_k = \frac{(\bar{T}(\bar{y}_v, \bar{\alpha}_v)C_k + Q_k)\bar{y}_k + \Pi_k \bar{\alpha}_k \tau}{e}. $$

REFERENCES

[1] Rémi Abgrall. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach. *Journal of Computational Physics*, 125(1):150–160, 1996.

[2] M. Bachmann, S. Müller, P. Helluy, and H. Mathis. A simple model for cavitation with non-condensable gases. In *Hyperbolic problems—theory, numerics and applications. Volume 1*, volume 17 of Ser. Contemp. Appl. Math. CAM, pages 289–296. World Sci. Publishing, Singapore, 2012.

[3] M. R. Baer and J. W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation transition (ddt) in reactive granular materials. *International journal of multiphase flow*, 12(6):861–889, 1986.

[4] T. Barberon and P. Helluy. Finite volume simulation of cavitating flows. *Computers & fluids*, 34(7):832–858, 2005.

[5] P. Boivin, M.A. Cannac, and O. Le Métayer. A thermodynamic closure for the simulation of multiphase reactive flows. *International Journal of Thermal Sciences*, 137:640–649, 2019. URL: https://hal.archives-ouvertes.fr/hal-01981954/document.

[6] C. G. Broyden. A class of methods for solving nonlinear simultaneous equations. *Mathematics of computation*, 19(92):577–593, 1965.

[7] T. Buffard, T. Gallouët, and J.-M. Hérard. A sequel to a rough godunov scheme: application to real gases. *Computers & fluids*, 29(7):813–847, 2000.

[8] J. Bussac. Study of relaxation processes in a two-phase flow model. Preprint, November 2021. URL: https://hal.archives-ouvertes.fr/hal-03418754.

[9] C. Chalons and J.-F. Coulombel. Relaxation approximation of the euler equations. *Journal of Mathematical Analysis and Applications*, 348(2):872–893, 2008.

[10] F. Coquel and B. Perthame. Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics. *SIAM J. Numer. Anal.*, 35(6):2223–2249, 1998. URL: http://dx.doi.org/10.1137/S0036142997318528.

[11] S. Gavrilyuk. The structure of pressure relaxation terms: the one-velocity case. Technical report, EDF, H-183-2014-0276-EN, 2014.

[12] U. S. NRC: Glossary. Loss of coolant accident (loca). URL: https://www.irsn.fr/FR/connaissances/Installations_nucleaires/Les-centrales-nucleaires/criteres_surete_ria_aprp/Pages/2-accident-perte-refrigerant-primaire.aspx?Id=6f765313-55e7-4581-8f5c-09b08a1ac859&dwId=3163662f-e01c-4895-af0f-0b82c2f2a3af# .Yd7ryC1h1qs.

[13] U. S. NRC: Glossary. Reactivity initiated accident (ria). URL: https://www.irsn.fr/FR/connaissances/Installations_nucleaires/Les-centrales-nucleaires/criteres_surete_ria_aprp/Pages/1-accident-reactivite-RIA.aspx#.XAkgQq17SE8.

[14] M. Hantke and S. Müller. Analysis and simulation of a new multi-component two-phase flow model with phase transitions and chemical reactions. *Quart. Appl. Math.*, 76(2):253–287, 2018. URL: https://doi.org/10.1090/qam/1498.

[15] P. Helluy, J.-M. Hérard, H. Mathis, and S. Müller. A simple parameter-free entropy correction for approximate riemann solvers. *Comptes rendus Mécanique*, 338(9):493–498, 2010.

[16] P. Helluy and N. Seguin. Relaxation models of phase transition flows. *M2AN Math. Model. Numer. Anal.*, 40(2):331–352, 2006. URL: http://dx.doi.org/10.1051/m2an:2006015, doi: 10.1051/m2an:2006015.

[17] J.-M. Hérard. A class of compressible multiphase flow models. *C. R. Math. Acad. Sci. Paris*, 354(9):954–959, 2016. URL: http://dx.doi.org/10.1016/j.crma.2016.07.004, doi: 10.1016/j.crma.2016.07.004.
[18] J.-M. Hérard, O. Hurisse, and L. Quibel. A four-field three-phase flow model with both miscible and immiscible components. *ESAIM Math. Model. Numer. Anal.*, 55(suppl.):S251-S278, 2021. doi:10.1051/m2an/2020037.

[19] J.-M. Hérard and H. Mathis. A three-phase flow model with two miscible phases. *ESAIM Math. Model. Numer. Anal.*, 53(4):1373-1389, 2019. doi:10.1051/m2an/2019028.

[20] O. Hurisse. Application of an homogeneous model to simulate the heating of two-phase flows. *International Journal on Finite Volumes*, 11, 2014. URL: https://hal.archives-ouvertes.fr/hal-01114808.

[21] O. Hurisse. Bgk source terms for out-of-equilibrium two-phase flow models. *Continuum Mechanics and Thermodynamics*, 2022. URL: https://doi.org/10.1007/s00161-022-01085-9.

[22] O. Hurisse and L. Quibel. A homogeneous model for compressible three-phase flows involving heat and mass transfer. *ESAIM: Proceedings and Surveys*, 66:84-108, 2019. URL: https://hal.archives-ouvertes.fr/hal-01976903.

[23] O. Hurisse and L. Quibel. Simulations of liquid-vapor water flows with non-condensable gases on the basis of a two-fluid model. *Appl. Math. Model.*, 99:514-537, 2021. doi:10.1016/j.apm.2021.06.020.

[24] O. Hurisse and L. Quibel. Simulations of water-vapor two-phase flows with non-condensable gas using a noble-able-chemkin stiffened gas equation of state. *Computers & Fluids*, 239:105399, 2022. URL: https://hal.archives-ouvertes.fr/hal-02963324.

[25] S. Jaouen. *Etude mathématique et numérique de stabilité pour des modèles hydrodynamiques avec transition de phase*. PhD thesis, Paris 6, 2001.

[26] J. Jung. *Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bi-fluides*. PhD thesis, Université de Strasbourg, 2013. URL: https://tel.archives-ouvertes.fr/tel-00876159v2.

[27] O. Le Métayer, J. Massoni, and R. Saurel. Élaboration des lois d'état d'un liquide et de sa vapeur pour les modèles d'écoulements diphasiques. *International Journal of Thermal Sciences*, 43(3):265-276, 2004.

[28] H. Mathis. A thermodynamically consistent model of a liquid-vapor fluid with a gas. *ESAIM: Mathematical Modelling and Numerical Analysis*, 53(1):63-84, 2019. URL: DOI:https://doi.org/10.1051/m2an/2018044, doi:10.1051/m2an/2018044.

[29] M. Pelanti and K.-M. Shyue. A numerical model for multiphase liquid-vapor-gas flows with interfaces and cavitation. *Int. J. Multiph. Flow*, 113:208-230, 2019. doi:10.1016/j.ijmultiphaseflow.2019.01.010.

[30] L. Quibel. *Simulation of water-vapor two-phase flows with non-condensable gas*. Theories, Université de Strasbourg, September 2020. URL: https://tel.archives-ouvertes.fr/tel-03505530.

[31] V. V. Rusanov. The calculation of the interaction of non-stationary shock waves with barriers. *Ž. Vyčisl. Mat i Mat. Fiz.*, 1:267-279, 1961.

[32] K. Saleh. A relaxation scheme for a hyperbolic multiphase flow model—Part I: Barotropic EOS. *ESAIM Math. Model. Numer. Anal.*, 53(5):1763-1795, 2019. doi:10.1051/m2an/2019034.

[33] Richard Saurel, Erwin Franquet, Eric Daniel, and Olivier Le Metayer. A relaxation-projection method for compressible flows. part i: The numerical equation of state for the euler equations. *Journal of Computational Physics*, 223(2):822-845, 2007.

[34] N. N. Yanenko. *The method of fractional steps. The solution of problems of mathematical physics in several variables*. Springer-Verlag, New York-Heidelberg, 1971. Translated from the Russian by T. Cheron. English translation edited by M. Holt.

NANTES UNIVERSITÉ, CNRS, LABORATOIRE DE MATHÉMATIQUES JEAN LERAY, LMJL, F-44000
NANTES, FRANCE

Email address: jean.bussac@univ-nantes.fr

NANTES UNIVERSITÉ, CNRS, LABORATOIRE DE MATHÉMATIQUES JEAN LERAY, LMJL, F-44000
NANTES, FRANCE

Email address: helene.mathis@univ-nantes.fr