# A 2 x 2 hyperbolic system modelling incompressible two phase flows: theory and numerics

M. Ndjinga, Tp. Nguyen, Christophe Chalons

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

## To cite this version:

M. Ndjinga, Tp. Nguyen, Christophe Chalons. A 2 x 2 hyperbolic system modelling incompressible two phase flows: theory and numerics. Nonlinear Differential Equations and Applications, 2017, 24, 10.1007/s00030-017-0458-6. hal-02339649

HAL Id: hal-02339649

https://hal.science/hal-02339649v1

Submitted on 16 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A 2 × 2 hyperbolic system modelling incompressible two phase flows: theory and numerics

M. Ndjinga, T. P. K. Nguyen and C. Chalons

**Abstract.** We propose a 2 × 2 hyperbolic system of conservation laws to model the dynamics of two incompressible fluids in mechanical disequilibrium. In the theoretical part of the paper we show that this 1D system is not strictly hyperbolic, that the characteristic speed can not a priori be ordered and that the characteristic fields are neither genuinely nonlinear, nor linearly degenerate. We nevertheless prove the existence and uniqueness of an admissible solution to the Riemann problem. This solution remains bounded with positive volume fractions even when one the phases vanishes. We conclude that the multiphase/single phase transition does not imply mechanical equilibrium but displays a non classical wave structure. In the numerical part of the paper we propose some approximate Riemann solvers to simulate the model, especially the multiphase/single phase transition. The classical Riemann solvers have been considered as Godunov scheme, Roe scheme with or without entropy fix. We also propose an in-cell discontinuous reconstruction method which proves to be successful, whereas the other schemes may show some spurious oscillations in some Riemann problem. Finally, as an application we study and simulate the problem of phase separation by gravity.

**Mathematics Subject Classification.** 35L65, 35L67, 35Q35, 35D30, 65M08, 76T10.

**Keywords.** Riemann problem, Incompressible fluid, Hyperbolic system, Non genuinely non linear, Two-fluid model.

## 1. Introduction

The flow regime involved in a nuclear reactor core is purely liquid in normal operating conditions but may become a liquid-gas mixture in incidental conditions or even purely gaseous in the case of a serious accident involving a core dewatering. The simulation of the single phase/multiphase transition is numerically challenging and has been a major difficulty in the design of new

<figure><img src="image_3.png" /></figure>
simulation platforms based on advanced two-fluid models, see for instance [1,2] and references therein. An important issue is to guarantee the positivity of the volume fraction of each phase. There is an open debate as to whether this positivity is intrinsic to the conservation laws or requires some adequate source terms such as inter-phase friction. The thermal hydraulics platform CATHARE [3] assumes that when a phase disappears, its velocity is equal to the velocity of the other phase. In order to strongly couple the two phase velocities, they use a very high interfacial friction term to deal numerically with these transitions. This paper intends to prove in the case of incompressible phases that the Riemann problem admits a positive solution without any frictional term and that the two velocities are not necessarily equal (Sect. 2), to propose some numerical methods able to deal with vanishing phases (Sect. 3) and then presents some numerical results (Sect. 4). Allowing phase disappearance with mechanical disequilibrium enables for example the modelling of bubbles ascending in a liquid as a consequence of Archimedes' principle.

## 1.1. The compressible model

We consider a one dimensional isentropic two phase flow involving two fluids with densities $ρ_1 < ρ_2$, pressures $P_1$ and $P_2$, volume fractions $α_1 \in [0, 1]$, $α_2 = 1 - α_1$, and velocities $u_1$ and $u_2$. After averaging the mass and momentum balance equations for each phase (see [4–6]), and neglecting mass and momentum transfer terms, the two-fluid model consists in the following four equations

$$ \partial_t α_1 ρ_1 + \partial_x (α_1 ρ_1 u_1) = 0, \qquad (1.1a) $$

$$ \partial_t α_2 ρ_2 + \partial_x (α_2 ρ_2 u_2) = 0, \qquad (1.1b) $$

$$ \partial_t (α_1 ρ_1 u_1) + \partial_x (α_1 ρ_1 u_1^2) + α_1 \partial_x P_1 = α_1 ρ_1 g, \qquad (1.1c) $$

$$ \partial_t (α_2 ρ_2 u_2) + \partial_x (α_2 ρ_2 u_2^2) + α_2 \partial_x P_2 = α_2 ρ_2 g, \qquad (1.1d) $$

where $g$ is the gravitational acceleration. Unlike [1,2,7,8] we do not assume pressure equality $P_1 = P_2$ nor do we introduce an interfacial pressure default term $Δρ∂_xα_k$ in the governing equations (1.1c) and (1.1d). Instead we consider a non zero pressure difference of the form $P_1 - P_2 = ρ_1 ρ_2 (u_1 - u_2)^2$ which is the minimal pressure gap yielding a hyperbolic system.

It is known since the 70s [9,10] that the single pressure two-fluid model may suffers instabilities due to the non real nature of its characteristic waves. The compressible six equations two-fluid model display six characteristic waves in 1D. Two of those are acoustic connected to pressure variations, and two are entropy waves connected to temperature and/or internal energy variations. These first four waves are necessarily real and similar to those encountered in the Euler equations for single phase flows. There are however two possibly complex waves that are specific to two phase flow and are strongly connected to dynamics of the void fraction $α$, and called void waves (see [4,11] for more details). The possibly unstable modes of the two-fluid model have been connected by authors such as W. Fullmer to the shear flow instability known as Kelvin–Helmholtz [12]. This instability is thus related to the surface instabilities that are captured by the model through the void fraction parameter (see
for example [13,14]). Some authors such as Keyfitz et al. [7,15] have studied measure type solution to the non hyperbolic system, whilst other authors such as [9,10] have considered relaxing the single pressure assumption through the introduction of an interfacial pressure. In this paper we will follow the second approach and allow a small pressure difference between the two phases that will yield mathematically an hyperbolic system with no unstable mode. This pressure difference originates physically from the surface tension stabilizing effect connected to the curvature of the interface and the surface tension coefficient $\gamma$ via the Young-Laplace equation [12]

$$P_1 - P_2 = \gamma \frac{1}{R_1} + \frac{1}{R_2} \qquad (1.2)$$

where $R_1$ and $R_2$ are the principal radii of curvature. In a non viscous two-phase flow, the Kelvin-Helmholtz instability arises at the interface when surface tension effects are dominated by shear flow instability. This usually generates a change of shape that distorts the bubbles or droplets and can lead to them breaking up. The consequence of those distortions and breaking up is an increase of the mean curvature $\frac{1}{R_T} + \frac{1}{R_2}$ and therefore of the surface tension stress. We are not interested in the modeling of bubbles/droplets distortion and break up. We will assume that surface tension dynamically adapts to oppose instantly the shear mechanism that triggers the instability. In the absence of surface tension, the Kelvin-Helmholtz instability may appear whenever there is a velocity jump across an interface and its amplitude is proportional to the square of the velocity jump. We therefore chose a dynamic surface tension model of the form

$$P_1 - P_2 = \rho^*(u_1 - u_2)^2 \qquad (1.3)$$

where $\rho^*$, should the smallest constant that stabilizes the system. Looking for the smallest perturbation of B. Keyfitz model that yields a well posed PDE problem we find the critical value

$$\rho^* = \frac{\rho_1 \rho_2}{2(\rho_1 - \rho_2)} (u_1 - u_2)^2 \qquad (1.4)$$

Any correlation $p_1 - p_2 = \rho(u_1 - u_2)^2$ will yield complex characteristic waves if $\rho < \rho^*$ and real characteristic waves if $\rho > \rho^*$.

The system (1.1a-1.1d) has four main unknowns: $\alpha_1, P_1, u_1, u_2$, the other unknowns can be obtained using the equations of state $\rho_k(P_k)$ and the pressure gap law

$$P_1 - P_2 = \frac{\rho_1 \rho_2}{2(\rho_1 - \rho_2)} (u_1 - u_2)^2. \qquad (1.5)$$

Defining the mixture sound wave $c_m = \frac{(\alpha_1 \rho_2 + \alpha_2 \rho_1) c_1^2 c_2^2}{\alpha_1 \rho_2 c_2^2 + \alpha_2 \rho_1 c_1^2}$ where $c_k = c_k(P_k)$,

$k = 1, 2$ are the sound speeds of each phase, the characteristic waves of the system are given by the roots of a fourth degree polynomial involving $c_m$, $u_1, u_2, \rho_1, \rho_2$ and $\alpha$. It is not possible to obtain simple closed expressions for these eigenvalues unless $u_1 = u_2$. However, following the work in [16], the Taylor expansion of the four eigenvalues of the system (1.1a-1.1d) can be derived near the mechanical equilibrium $u_1 = u_2$. If $u_1 - u_2 = c_m$, $\frac{u_1 - u_2}{c_m}$
can be considered as a small perturbation and the Taylor expansion can be written as

$$ \lambda_1 = \frac{\rho_1 u_1 - \rho_2 u_2}{\rho_1 - \rho_2} - \frac{\rho_1 \rho_2}{(\alpha_1 \rho_2 + \alpha_2 \rho_1)^2} + O \frac{u_1 - u_2}{c_m} \quad \text{void fraction (or transport) wave,} $$

$$ \lambda = \frac{\rho_1 u_1 - \rho_2 u_2}{\rho_1^2 - \rho_2^2} + O \frac{u_1 - u_2}{c_m} \quad \text{void fraction (or transport) wave,} $$

$$ \lambda_3 = \frac{\alpha_1 \rho_2 u_1 + \alpha_2 \rho_1 u_2}{\alpha_1 \rho_2 + \alpha_2 \rho_1} + c_m + O \frac{u_1 - u_2}{c_m} \quad \text{acoustic wave,} $$

$$ \lambda_4 = \frac{\alpha_1 \rho_2 u_1 + \alpha_2 \rho_1 u_2}{\alpha_1 \rho_2 + \alpha_2 \rho_1} - c_m + O \frac{u_1 - u_2}{c_m} \quad \text{acoustic wave.} $$

Thus for small relative velocities $u_1 - u_2$ $c_m$, the system (1.1a-1.1d) is hyperbolic with 2 acoustic waves involving the mixture sound speed $\lambda_1$ and $\lambda_2$ that are similar to those encountered in the Euler system for single phase flows, and two void waves that are specific to the two phase dynamics. Since we are interested by the void wave dynamics and flows at low Mach numbers, we devote more attention to the incompressible limit of the system (1.1a-1.1d).

## 1.2. The incompressible model

In order to study more precisely the volume fraction waves involved in applications where the fluid densities are almost constant, we follow [7] and assume that both phases are incompressible with constant densities $\rho_1$ and $\rho_2$. We will consider the incompressible model which is obtained by applying the same method as in [7] but using the pressure law (1.5) instead of considering $P_1 = P_2$ as was done in [7]. Recalling the closure laws, $\alpha_1 + \alpha_2 = 1$, the four equation system (1.1a-1.1d) should be solved for the four unknowns ($\alpha_1, u_1, u_2$, and $P_1$).

Inspired by the method in [7], we derive a system of two equations, which allows for the study of the void waves and avoids singularities when one phase disappears.

The sum of $\frac{1}{\rho_1}$ (1.1a) and $\frac{1}{\rho_2}$ (1.1b) yields

$$ \partial_x(\alpha_1 u_1 + \alpha_2 u_2) = 0. \quad (1.6) $$

Therefore the quantity

$$ K = \alpha_1 u_1 + \alpha_2 u_2 \quad (1.7) $$

is constant in space and can be determined from the boundary conditions. Assuming that the boundary condition is independent of time, we obtain that $K$ is constant both in space and time.

From the equation of momentum conservation (1.1c) and (1.1d), we derive

$$ u_1 \rho_1 (\partial_t \alpha_1 + \partial_x \alpha_1 u_1) + \alpha_1 \rho_1 (\partial_t u_1 + u_1 \partial_x u_1) + \alpha_1 \partial_x P_1 = \alpha_1 \rho_1 g, \quad (1.8a) $$

$$ u_2 \rho_2 (\partial_t \alpha_2 + \partial_x \alpha_2 u_2) + \alpha_2 \rho_2 (\partial_t u_2 + u_2 \partial_x u_2) + \alpha_2 \partial_x P_2 = \alpha_2 \rho_2 g. \quad (1.8b) $$
Thank to (1.1a) and (1.1b), these equations can be simplified:

$$ \alpha_1 \partial_t (\rho_1 u_1) + \partial_x \frac{\rho_1 u_1^2}{2} + \alpha_1 \partial_x P_1 = \alpha_1 \rho_1 g, \quad (1.9a) $$

$$ \alpha_2 \partial_t (\rho_2 u_2) + \partial_x \frac{\rho_2 u_2^2}{2} + \alpha_2 \partial_x P_2 = \alpha_2 \rho_2 g. \quad (1.9b) $$

Assuming that initially $\alpha_1 \alpha_2 \neq 0$, we can simplify by $\alpha_1$ in (1.9a) and by $\alpha_2$ in (1.9b), then subtracting the two equations yields

$$ \partial_t (\rho_1 u_1 - \rho_2 u_2) + \partial_x \frac{1}{2} (\rho_1 u_1^2 - \rho_2 u_2^2) + P_1 - P_2 = (\rho_1 - \rho_2)g. \quad (1.10) $$

The space differential in (1.10) can be simplified thanks to
$$ \begin{aligned} \frac{1}{2} (\rho u_1^2 - \rho u_2^2) + P_1 - P_2 \frac{\rho u_1^2 - \rho u_2^2}{2} + \frac{\rho_1 \rho_2}{2} (u_1 - u_2)^2 &= \frac{1}{2} \begin{pmatrix} 1 & 1 & 2 & 2 \\ 2 & 1 & 1 & 2 \\ 2 & 2 & 1 & 1 \\ 2 & 2 & 2 & 1 \end{pmatrix} = \frac{1}{2(\rho_1 - \rho_2)} \begin{pmatrix} \rho u_1 - \rho u_2 \\ 1 & 1 & 2 & 2 \end{pmatrix}^2. \end{aligned} \quad (1.11) $$

We set the new unknowns $(\alpha, \omega)$ as

$$ \begin{aligned} \alpha &= \alpha_1, \\ \omega &= \rho_1 u_1 - \rho_2 u_2. \end{aligned} \quad (1.12) $$

The original unknowns $u_1$ and $u_2$ can be recovered from (1.7) and (1.12):

$$ \begin{aligned} u_1 &= \frac{(1 - \alpha)\omega}{\alpha(\rho_2 - \rho_1) + \rho_1} + \frac{K\rho_2}{\alpha(\rho_2 - \rho_1) + \rho_1}, \\ u_2 &= \frac{-\alpha\omega}{\alpha(\rho_2 - \rho_1) + \rho_1} + \frac{K\rho_1}{\alpha(\rho_2 - \rho_1) + \rho_1}. \end{aligned} $$

Finally, from (1.1a), (1.10) and (1.11) we obtain the 2 × 2 system. More precisely, the incompressible two-fluid model can be written in closed form as

$$ \left\{
\begin{aligned}
\partial_t \alpha + \partial_x \frac{\alpha(1-\alpha)\omega}{\alpha(\rho_2-\rho_1)+\rho_1} + (\rho_1-\rho_2)\frac{K\rho_1\rho_2}{\alpha(\rho_2-\rho_1)+\rho_1} &= 0, && (1.13) \\
\partial_t \omega + \partial_x \frac{\omega^2}{2(\rho_1-\rho_2)} &= (\rho_1-\rho_2)g,
\end{aligned}
\right. $$

since $K$ is constant in space and time.

For simplification, we will assume that $K = 0$, which is true for example if there is a wall boundary condition. If $K \neq 0$ is constant both in time and space, a Galilean change of reference frame $u'_k = u_k - K$ yields a new system with $K' = \alpha_1 u'_1 + \alpha_2 u'_2 = 0$ (see [7]). With the assumption $K = 0$, the system (1.13) can be therefore rewritten as

$$ \left\{
\begin{aligned}
\partial_t \alpha + \partial_x \frac{\alpha(1-\alpha)\omega}{\alpha(\rho_2-\rho_1)+\rho_1} &= 0, \\
\partial_t \omega + \partial_x \frac{\omega^2}{2(\rho_1-\rho_2)} &= (\rho_1-\rho_2)g,
\end{aligned}
\right. $$

or in the following compact form

$$ \partial_t U + \partial_x F(U) = S, \quad (1.14) $$
with

$$ U = \begin{pmatrix} \alpha \\ \omega \end{pmatrix}, \quad F(U) = \begin{pmatrix} \frac{\alpha(1-\alpha)\omega}{\alpha(\rho_2-\rho_1)+\rho_1} \\ \frac{\omega^2}{2(\rho_1-\rho_2)} \end{pmatrix}, \quad S = \begin{pmatrix} 0 \\ (\rho_1-\rho_2)g \end{pmatrix}, \quad (1.15) $$

$$ \nabla F(U) = \frac{\omega}{\rho_1 - \rho_2} \left( 1 - \frac{\frac{\rho_1 \rho_2}{\alpha(\rho_2 - \rho_1) + \rho_1}}{0} \right)^2 \frac{\alpha(1-\alpha)}{\alpha(\rho_2 - \rho_1) + \rho_1} \frac{\omega}{\rho_1 - \rho_2} \quad (1.16) $$

## 2. Theoretical study

### 2.1. Hyperbolicity and characteristic fields

The following theorem is a direct consequence of formula (1.16).

**Theorem 1.** (Hyperbolicity of system (1.14)) The jacobian matrix (1.16) of the system (1.14) always admits real eigenvalues

$$ \lambda_1 = \frac{\omega}{\rho_1 - \rho_2} \left( 1 - \frac{\rho_1 \rho_2}{(\alpha(\rho_2 - \rho_1) + \rho_1)^2} \right), \quad \lambda_2 = \frac{\omega}{\rho_1 - \rho_2}, \quad (2.17) $$

provided $U \in \mathbb{R}$ where

$$ 7 = \{(\alpha, \omega), \alpha \in [0, 1], \omega \in \mathbb{R}\}. $$

Moreover, such a matrix is diagonalisable provided $(\alpha, \omega) \in \mathbb{R}_*$ where

$$ \begin{gathered} 7_* = 7_+ \cup 7_- \cup \{(0, \omega), \omega \in \mathbb{R}\} \cup \{(1, \omega), \omega \in \mathbb{R}\}, \\ \text{and } 7_{\pm} = \{(\alpha, \omega), \omega \in \mathbb{R}_{\pm} \setminus \{0\}, \text{ and } \alpha \in (0, 1)\}. \end{gathered} $$

In general, the system (1.14) is weakly hyperbolic on the domain $7 = 7_* \cup \{(\alpha, 0), \alpha \in (0, 1)\}$.

The states $(\omega = 0, \alpha = 0)$ and $(\omega = 0, \alpha = 1)$ will play an important role in connecting states in $7_+$ to states in $7_-$ and will be called **critical states**.

$7$ is neither an open nor a simply connected subset of $\mathbb{R}^2$, see Fig. 1. The eigenvectors of $\nabla F$ are

$$ \rightarrow r_1 = {}^t(1, 0), \quad \rightarrow r_2 = {}^t(\alpha(1 - \alpha)(\rho_1 - \rho_2)(\alpha(\rho_2 - \rho_1) + \rho_1), \rho_1\rho_2\omega). \quad (2.18) $$

<figure><img src="image_4.png" /><figcaption>FIGURE 1. Strict hyperbolicity domain $7_* = 7_+ \cup 7_- \cup \{(0, \omega), \omega \in \mathbb{R}\} \cup \{(1, \omega), \omega \in \mathbb{R}\}$</figcaption></figure>
Assuming that $\rho_1 < \rho_2$, we notice that the two eigenvalues are not a priori ordered since

$$ \left\{ \begin{array}{ll} \lambda_1 < \lambda_2 & \text{if } \omega < 0, \\ \lambda_1 > \lambda_2 & \text{if } \omega > 0, \\ \lambda_1 = \lambda_2 & \text{if } \omega = 0. \end{array} \right. \qquad (2.19) $$

Moreover, the signs of $\nabla\lambda_1 \cdot \nabla r_1 = \frac{-2\rho_1\rho_2\omega}{(\alpha(\rho_2-\rho_1)+\rho_1)}$ and $\nabla\lambda_2 \cdot \nabla r_2 = \frac{\rho_1\rho_2\omega}{(\rho_1-\rho_2)}$ are not known a priori since

$$ \left\{ \begin{array}{ll} \nabla\lambda_k \cdot \nabla r_k > 0 & \text{if } \omega < 0, \\ \nabla\lambda_k \cdot \nabla r_k < 0 & \text{if } \omega > 0, \\ \nabla\lambda_k \cdot \nabla r_k = 0 & \text{if } \omega = 0. \end{array} \right. \qquad (2.20) $$

where $k = 1, 2$. Therefore the characteristic fields associated to $\lambda_1$ and $\lambda_2$ are genuinely nonlinear in each domain $7_+$ and $7_-$, but are neither genuinely non linear, nor linearly degenerate in general.

## 2.2. Triangular systems of conservation laws

The system (1.14) without source terms is a particular case of triangular systems of conservation laws due to the fact that the second equation is independent of the first one. Let us consider the simplest triangular system of conservation laws

$$ \partial_t \alpha + \partial_x g(\alpha, \omega) = 0, \qquad (2.21a) $$

$$ \partial_t \omega + \partial_x f(\omega) = 0, \qquad (2.21b) $$

which has been studied largely in the literature. For example, the reader is referred to [17–21] and references therein for more details. Inspired by the theory of scalar conservation laws, a simple way to solve a triangular system of 2 × 2 equations is to compute the solution of the independent equation then replace it to the remaining one. However for the general results of existence and uniqueness to the Cauchy problem (or even the Riemann problem) for the system (2.21a–2.21b) with regular functions of $f(\omega)$ and $g(\alpha, \omega)$ is still open. The difficulties encountered include the non strict hyperbolicity and the resonance.

Let us introduce well-known results for the triangular system as well as the approaches so that we have a general point of view and can figure out our contribution in this interest. The first approach is to generalize the weak solutions, i.e. to extend the space of an admissible solution in the sense that they are not necessarily bounded. Following this approach, some authors have been studying non strictly hyperbolic triangular systems. More precisely, considering a particular system of (2.21a–2.21b) where $g(\alpha, \omega) = \alpha\omega$ and $f(\omega) = \omega^2$, the authors in [17,21] introduced the class of admissible solutions to the Riemann problems admitting delta-shocks. Such study has been extended by many researchers, for example in [18,22] and references therein. The general results obtained for the Cauchy problem in such references require the function $g(\alpha, \omega)$ to be linear with respect to $\alpha$ for each $\omega$, i.e. $g(\alpha, \omega)$ can be rewritten as
$g(\alpha, \omega) = h(\omega)\alpha$. The key idea in this case is to use the study of the linear transport equation where the velocity admits the discontinuities in space and in time. The resulting admissible solutions include the delta-shock waves, which appear around the configuration of non strict hyperbolicity.

Another approach can be found in [19,20], where the authors used the theory of compensated compactness to prove the existence and uniqueness of the Cauchy problem for a triangular hyperbolic system. The main idea is to use numerical schemes, as Godunov-type in [19] or the relaxation scheme in [20], in order to prove the convergence of the schemes which implies the existence of the solution. However, the proof of the convergence strongly depends on the assumption of the function $g(\alpha, \omega)$. More precisely, the function $g(\alpha, \omega)$ must be genuinely nonlinear in $\alpha$ for any $\omega$ in the interesting domain, i.e. $6\omega$, $g_{\alpha\alpha}(\alpha, \omega) \neq 0$.

We cannot apply directly the results from the references presented above to the system (1.14) although there are some similar properties between our system and the one in [17,21], such as not strict hyperbolicity when $\omega = 0$, two eigenvalues that are not ordered, etc. The first reason is that we were not able to exhibit entropy pairs for our system, which is consistent with the fact that incompressible Euler equations do not admit obvious entropies, and the system (1.14) is therefore not symmetrizable. The second reason is that our flux function in Eq. (2.21a) is not linear with respect to $\alpha$, therefore we are not able to use directly the theory of transport equations. Moreover, the function $g(\alpha, \omega)$ is not genuinely nonlinear with respect to $\alpha$. We therefore present in this paper a new approach yielding the existence and uniqueness of the admissible solution to the Riemann problem for the system (1.14) without source terms.

## 2.3. Admissible solutions of the Riemann problem

The fact that the domain $\mathcal{J}$ is not open, that the system is not strictly hyperbolic and that the characteristic fields are neither genuinely nonlinear neither linearly degenerate raises many theoretical as well as numerical difficulties. We cannot use the classical Lax theorem (see [23]) to obtain solutions to the Riemann problem but will however build solutions to the Riemann problem having a non classical wave structure for any pair of left and right states data in $\mathcal{J}$.

**Definition 2.1.** (Hugoniot locus and Riemann invariants) Given a state $U \in \mathcal{J}$, the Hugoniot locus $\mathcal{F}(U)$ associated to (2.28) and $U$ is the set of states that can be connected to $U$ via a shock wave:

$$ \mathcal{F}(U) = \{V \in \mathcal{J}, \exists \sigma(U, V) \in \mathbb{R}, F(U) - F(V) = \sigma(U, V)(U - V)\}. $$

For any $k \in \{1, 2\}$, a $k$-Riemann invariant associated to (2.28) is a function $R_k$ defined on $\mathcal{J}$ such that

$$ 6U \in \mathcal{J}, \quad \nabla R_k \cdot \rightarrow r_k = 0. $$

and the $k$-rarefaction wave associated to the Riemann invariant $R_k$ at a left state $U$ is

$$ k\text{-}R(U) = \{V \in \mathcal{J}, \quad R_k(V) = R_k(U), \quad \lambda_k(V) \geq \lambda_k(U)\}. $$
We did not find an entropy pair to our system and instead of an entropy based selection criterion we use the Liu criterion to define admissible solutions.

**Definition 2.2.** (Admissible solution of the Riemann problem) An admissible solution to (2.28) is a weak solution that is composed of a finite number of constant states connected by a rarefaction wave or a shock wave, connecting left and right states $U_L$ and $U_R$ and propagating at a speed $\sigma$, that satisfies the Liu criterion:

$$ \sigma(U_L, U_R) \leq \sigma(U_L U), \quad 6U \in \mathcal{F}(U_L) \text{ between } U_L \text{ and } U_R. $$

Liu's shock admissibility criterion is more general than the classical Lax's shock admissibility criterion which is often used to select the entropy solution, see [24,25]. These two criteria coincide in the cases where the system is genuinely nonlinear [26]. We will use Liu's criterion since in our case, the incompressible system (1.14) is not genuinely nonlinear.

## 2.4. Shock and rarefaction waves

In the theory of 2 × 2 strictly hyperbolic systems with genuinely nonlinear characteristic fields, each eigenvector family $\#_{k}$, $k = 1, 2$ defines at any state $U_0$ one single shock and one single rarefaction curve, both starting and having tangent vector $\rightarrow r_k$ at $U = U_0$. These two curves are defined in an open neighborhood of $U_0$ and do not necessarily extend to the entire domain, see [23,27] for example.

Our system is neither strictly hyperbolic neither genuinely nonlinear, our strict hyperbolicity domain is not open. However we prove below that there are two families of shock “curves” (1-shocks and 2-shocks) and rarefaction “curves” (1-rarefactions and 2-rarefactions) and that they extend to the entire domain.

The originality in our system is the fact that the shock “curve” associated to a state is not necessarily a connected set. Indeed the 1-shock “curve” associated to the states $(\alpha \in [0, 1], \omega = 0)$ is a contact discontinuity whose speed is $\sigma = 0$, whereas there is no 1-rarefaction associated to such states. Moreover the 2-shock family of “curves” can be composed of the two branches of a hyperbola and there are more than one 2-rarefaction curves passing through the critical state $(\alpha = 1, \omega = 0)$.

We now characterise the shock and rarefaction curves and illustrate them on Fig. 2. We note that the admissibility criterion is not taken into account in the following theorem.

**Theorem 2.** (Existence of shock curves) Any state $U_0 = (\alpha_0, \omega_0) \in \mathbb{R}\,(\alpha = 1, \omega = 0)$ belongs to two shock curves in 7. The state $(\alpha = 1, \omega = 0)$ belong to three shock curves in 7.

The equation of the 1-shock family (1-S) of states $(\alpha, \omega_1(\alpha))$ connected to the state $U_0 = (\alpha_0, \omega_0)$ is

* If $\omega_0 \neq 0$, $\omega_1(\alpha) = \omega_0$

* If $\omega_0 = 0$, $\omega_1(\alpha) = 0$, the 1-shock “curve” degenerates a contact discontinuity whose the speed is $\sigma = 0$.
The equation of the 2-shock family (2-S) of states $(\alpha, \omega_2(\alpha))$ connected to the state $U_0 = (\alpha_0, \omega_0)$ is

• If $\omega_0 \neq 0$ and $\alpha_0 \in (0, 1)$

$$ \omega_2(\alpha) = \frac{\alpha(\rho_2 - \rho_1) + \rho_1}{\alpha_0(\rho_2 - \rho_1) + \rho_1} \times \frac{\alpha^2(\rho_2 - \rho_1) + \alpha_0(\rho_1 - 2\rho_2 + \alpha\rho_2 - \alpha\rho_1) + \alpha\rho_1}{\alpha^2(\rho_2 - \rho_1) + \alpha(\rho_1 - 2\rho_2 + \alpha_0\rho_2 - \alpha_0\rho_1) + \alpha_0\rho_1} \omega_0 $$

• If $\omega_0 \neq 0$ and $\alpha_0 = 0$, then either $\alpha = 0$ or $\omega_2 = \frac{(\alpha(\rho_2 - \rho_1) + \rho_1)\omega_0}{\alpha(\rho_2 - \rho_1) + \rho_1 - 2\rho_2}$.

• If $\omega_0 \neq 0$ and $\alpha_0 = 1$, then either $\alpha = 1$ or $\omega_2 = \frac{-(\alpha(\rho_2 - \rho_1) + \rho_1)\omega_0}{-\alpha(\rho_2 - \rho_1) + \rho_1}$.

• If $\omega_0 = 0$ and $\alpha_0 \in [0, 1)$, the 2-shock curve is straight line defined by the constant $\alpha$ which is the unique solution in $[0, 1)$ of the following equation

$$ \alpha^2(\rho_2 - \rho_1) + \alpha(\rho_1 - 2\rho_2 + \alpha_0\rho_2 - \alpha_0\rho_1) + \alpha_0\rho_1 = 0. $$

• If $\omega_0 = 0$ and $\alpha_0 = 1$, then either $\alpha = 1$ or $\alpha = \frac{\rho_1}{\rho_2 - \rho_1}$ if $\rho_2 > 2\rho_1$.

*Proof.* A state $U$ in $\mathcal{F}(U_0)$ must fulfil the Rankine-Hugoniot condition

$$ F(U) - F(U_0) = \sigma(U, U_0)(U - U_0), \quad (2.22) $$

which for the system (2.28) takes the form

$$ \left\{ \begin{array}{l} \frac{\alpha(1-\alpha)\omega}{\omega^2} - \frac{\alpha_0(1-\alpha_0)\omega_0}{\alpha(\rho_2-\rho_1)+\rho_1} = \sigma(U, U_0)(\alpha - \alpha_0), \\ \frac{\omega_0^2}{2(\rho_1-\rho_2)} - \frac{\omega_0^2}{2(\rho_1-\rho_2)} = \sigma(U, U_0)(\omega - \omega_0). \end{array} \right. \quad (2.23) $$

• If $\omega = \omega_0$, there always exists $\sigma(U, U_0)$ such that (2.23). This comes from the fact that the second equation is independent from $\alpha$. The 1-shock family is made of states sharing the same values of $\omega$ which is consistent with the expression of the eigenvector $\to r_1$ (Eq. 2.18).

• If $\omega \neq \omega_0$, the second equation of (2.23) implies $\sigma(U, U_0) = \frac{\omega_0 + \omega_1}{2(\rho_1 - \rho_2)}$, which yields in the first equation

$$ \frac{\alpha^2(\rho_2 - \rho_1) + \alpha(\rho_1 - 2\rho_2 + \alpha_0\rho_2 - \alpha_0\rho_1) + \alpha_0\rho_1}{\alpha(\rho_2 - \rho_1) + \rho_1} \omega - \frac{\alpha_0^2(\rho_2 - \rho_1) + \alpha_0(\rho_1 - 2\rho_2 + \alpha\rho_2 - \alpha\rho_1) + \alpha\rho_1}{\alpha_0(\rho_2 - \rho_1) + \rho_1} \omega_0 = 0. \quad (2.24) $$

— If $\omega_0 = 0$, (2.24) implies

$$ \alpha^2(\rho_2 - \rho_1) + \alpha(\rho_1 - 2\rho_2 + \alpha_0\rho_2 - \alpha_0\rho_1) + \alpha_0\rho_1 = 0. \quad (2.25) $$

It is easy to see that $\delta\alpha \in [0, 1)$, the Eq. (2.25) has a unique solution $\alpha \in [0, 1)$, and such an $\alpha$ defines the straight line 2-shock curve connecting to $(\alpha_0, \omega_0 = 0)$. Moreover, the state $(\alpha_0 = 1, \omega_0 = 0)$ connects to the straight line $\alpha = 1$ and can connect to all states $(\alpha = \frac{\rho_1}{\rho_2 - \rho_1}, \omega \in \mathbb{R}^*)$ under the assumption that $\rho_2 > 2\rho_1$.

— If $\omega_0 \neq 0$ and $\alpha_0 = 0$, then (2.24) implies either $\alpha = 0$ or $\omega = \frac{(\alpha(\rho_2 - \rho_1) + \rho_1)\omega_0}{\alpha(\rho_2 - \rho_1) + \rho_1 - 2\rho_2}$.

— If $\omega_0 \neq 0$ and $\alpha_0 = 1$, then (2.24) implies either $\alpha = 1$ or $\omega = \frac{-(\alpha(\rho_2 - \rho_1) + \rho_1)\omega_0}{-\alpha(\rho_2 - \rho_1) + \rho_1}$.
— If $\omega_0 \neq 0$ and $\alpha_0 \in (0, 1)$, then then (2.24) implies

$$ \omega(\alpha) = \frac{\alpha(\rho_2 - \rho_1) + \rho_1}{\alpha_0(\rho_2 - \rho_1) + \rho_1} \times \frac{\alpha^2(\rho_2 - \rho_1) + \alpha_0(\rho_1 - 2\rho_2 + \alpha\rho_2 - \alpha\rho_1) + \alpha\rho_1}{\alpha^2(\rho_2 - \rho_1) + \alpha(\rho_1 - 2\rho_2 + \alpha_0\rho_2 - \alpha_0\rho_1) + \alpha_0\rho_1} \omega_0. $$

The shape of the shock curves can be seen on Fig. 2.

**Theorem 3.** (Existence of rarefaction curves)

* Any state $U_0 = (\alpha_0, \omega_0) \in \mathbb{R}\{(\alpha, \omega = 0)\}$ belongs to two rarefaction curves, the 1-rarefaction curve (1-R) is described by the equation $\omega_1(\alpha) = \omega_0$ and the 2-rarefaction curve (2-R) is described by the equation

$$ \omega_2(\alpha) = \omega_0 \frac{\alpha(\rho_2 - \rho_1) + \rho_1}{\alpha_0(\rho_2 - \rho_1) + \rho_1} \frac{\alpha}{\alpha_0} \frac{\frac{\rho_2}{\rho_1 - \rho_2}}{1 - \alpha_0} \frac{1 - \alpha}{1 - \alpha_0} \quad (2.26) $$

* The state $(\alpha_0 = 0, \omega_0 = 0)$ belongs to a single rarefaction curve, the 2-rarefaction curve described by $\alpha = 0$.

* The state $(\alpha_0 = 1, \omega_0 = 0)$ belongs to all 2-rarefaction curves going through $U_0 = (\alpha_0, \omega_0) \in \mathbb{R}\{(\alpha = 0, \omega = 0)\}$ and described by (2.26) if $\alpha_0 \neq 1$ or $\alpha = 1$, otherwise.

* There is no rarefaction going through the state $(\alpha_0 \in (0, 1), \omega = 0)$.

*Proof.* From the definition of Riemann invariants, $R_1$ and $R_2$ must satisfy

$$ \frac{\partial R_1}{\partial \alpha} = 0; \quad \alpha(1 - \alpha)(\rho_1 - \rho_2)(\alpha(\rho_2 - \rho_1) + \rho_1)\frac{\partial R_2}{\partial \alpha} + \rho_1 \rho_2 \omega \frac{\partial R_2}{\partial \omega} = 0. $$

Since $R_1$ is function of only $\omega$, it is easy to obtain the equation of the 1-rarefaction curves.

Considering the Riemann invariant $R_2$,

* if $\alpha \in (0, 1)$, we have

$$ \frac{\partial \omega}{\partial \alpha} = \frac{\rho_1 \rho_2 \omega}{\alpha(1 - \alpha)(\rho_1 - \rho_2)(\alpha(\rho_2 - \rho_1) + \rho_1)}. \quad (2.27) $$

Solving the linear ordinary differential equation (2.27), we obtain explicitly the equation of the rarefaction curves (2.26) in the theorem.

* if $\alpha \in \{0, 1\}$, then

- either $\frac{\partial R_2}{\partial \omega} = 0$, the rarefaction curves are $\alpha(\omega) = 0$ (or $\alpha(\omega) = 1$) corresponding to $\alpha_0 = 0$ (or $\alpha_0 = 1$),

- or $\omega = 0$. However, only the state $(\alpha = 1, \omega = 0)$ belongs to all $R_2$ described by (2.26) due to continuing property $(\lim_{\alpha \to 1} \omega_2(\alpha) = 0, \lim_{\alpha \to 0} \omega_2(\alpha) = \pm\infty)$.

The structure of the rarefaction curves is illustrated on Fig. 2.

## 2.5. Solution of the Riemann problem

We consider the Riemann problem for the conservative system (1.14) in the case $g = 0$ with a piecewise constant initial data:
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Shock curves and Rarefaction curves. $U_0 = (0, 40):$ green. $U_0 = (0.5, -10):$ red. $U_0 = (0.5, 10):$ violet. $U_0 = (1, -20):$ blue. a 1-Shock curves, b 1-rarefaction curves, c 2-shock curves, d 2-rarefaction curves (color figure online)</figcaption></figure>

$$ \begin{gathered} \partial_t U + \partial_x F(U) = 0 \\ U_0(x) = \begin{cases} U_L(\alpha_L, \omega_L) & \text{if } x \le 0, \\ U_R(\alpha_R, \omega_R) & \text{if } x > 0, \end{cases} \end{gathered} \quad (2.28) $$

where $U =^t (\alpha, \omega) \in \mathbb{R}$, $F(U)$ are defined in (1.15), and $U_L, U_R \in \mathbb{R}$. We start by stating and proving in Sect. 2.5.1 existence and uniqueness of an admissible solution to (2.28) in Theorem 4. We then give in Sect. 2.5.2 four important examples of non classical solutions to the Riemann problem in order to illustrate the proof of Theorem 4 and the unusual behaviour of the system (2.28) as well.

**2.5.1. Main result.** We intend to prove existence and uniqueness of an admissible solution. The weak solutions are found using the classical characteristic method with the shock and the rarefaction curves found in Theorems 2 and 3 correspondingly. Due to the complex structure of the shock curves, there are numerous different cases to be considered. We first describe the structure of admissible solutions in the following lemma.

**Theorem 1.** Assume that an admissible solution of the Riemann problem contains two adjacent waves which connect the left state $U_L = (\alpha_L, \omega_L) \in \mathbb{R}$ to the
right state $U_R = (\alpha_R, \omega_R) \in \mathbb{N}_0$ through the intermediate state $U_I = (\alpha_I, \omega_I)$. Then, these two adjacent waves must be in different family (i.e. one wave from a 1-family and the other from a 2-family) at the exception that a 2-rarefaction may be followed by a 2-rarefaction.

Moreover, if the two adjacent waves are a 1-wave followed by a 2-wave (resp. a 2-wave followed by a 1-wave), then the intermediate state satisfies $\omega_I = \omega_L \le 0$ (resp. $\omega_I = \omega_R \ge 0$). Therefore, the value of $\omega$ satisfies the maximum principle.

**Proof.** This lemma is a direct consequence of the characterisation of k-shock curves and k-rarefaction curves (Theorems 2 and 3) and of the speed order criterion. Assume first that two adjacent waves are two 1-shock waves connecting the left state $U_L = (\alpha_L, \omega_L)$ to the right state $U_R = (\alpha_R, \omega_R)$ through the intermediate state $U_I = (\alpha_I, \omega_I)$. A consequence of Theorem 2 is that $\omega_L = \omega_R = \omega_I$. From the first equation in (2.23), the speed of the first 1-shock and the second 1-shock is

$$ \sigma_{1_1} = \frac{\omega_L}{\rho_1 - \rho_2} - \frac{\rho_1 \rho_2}{(\alpha_I (\rho_2 - \rho_1) + \rho_1)(\alpha_L (\rho_2 - \rho_1) + \rho_1)} \quad \text{and} \quad (2.29) $$

$$ \sigma_{1_2} = \frac{\omega_R}{\rho_1 - \rho_2} - \frac{\rho_1 \rho_2}{(\alpha_I (\rho_2 - \rho_1) + \rho_1)(\alpha_R (\rho_2 - \rho_1) + \rho_1)} \quad \text{and} \quad (2.30) $$

Liu's criterion for both shocks gives $\alpha_R > \alpha_I > \alpha_L$ if $\omega_L = \omega_R = \omega_I > 0$ and $\alpha_R < \alpha_I < \alpha_L$ if $\omega_L = \omega_R = \omega_I < 0$, while the speed order criterion $\sigma_{1_1} < \sigma_{1_2}$ gives $\alpha_R < \alpha_L$ if $\omega_L = \omega_R = \omega_I > 0$ and $\alpha_R > \alpha_L$ if $\omega_L = \omega_R = \omega_I < 0$. Therefore, the admissible solutions of the Riemann problem do not admit two adjacent waves in the same 1-shock family.

The conclusion for the 1-rarefaction family and 2-shock family are completely the same while the conclusion for the 2-rarefaction family is an exception.

Since the two 2-rarefaction curves join at unique point $(\alpha = 1, \omega = 0)$, if two adjacent waves are two 2-rarefaction, then $\omega_L > 0 > \omega_R$. Such two 2-rarefaction satisfy the speed criterion although there exists only an intermediate point $(\alpha = 1, \omega = 0)$ (no intermediate constant state). It is the same conclusion of a system in [17].

As for the second statement of the lemma, assume that the intermediate state $U_I = (\alpha_I, \omega_I)$ connects to the left state (resp. right state) by a 1-wave and connects to the right state (resp. the left state) by a 2-wave, due to the fact that the equation of 1-wave is that $\omega$ is constant, see theorem of existence of shock curves Theorem 2 and rarefaction curves Theorem 3, then $\omega_I = \omega_L$ (or $\omega_I = \omega_R$). In order to prove the sign property of $\omega_I = \omega_L$ in the case of a 1-wave followed by a 2-wave, let us notice that for the solution to be admissible, the speed of propagation of the 1-wave must be smaller or equal to the one of the 2-wave. (2.19) thus imposes $\omega_I = \omega_L \le 0$. Similarly, we also obtain the sign property $\omega_I = \omega_R \ge 0$ in the case of a 2-wave followed by a 1-wave.

**Theorem 4.** (Existence and uniqueness for the Riemann problem) Let us consider two states $U_L(\alpha_L, \omega_L)$ and $U_R(\alpha_R, \omega_R)$ in 7. The Riemann problem
(2.28) admits a unique admissible solution (in the sense of Sect. 2.3) $U(x, t) \in \mathcal{J}$ which depends continuously on the initial condition.

*Proof.* Let $U_L = (\alpha_L, \omega_L) \neq U_R = (\alpha_R, \omega_R)$ be in $\mathcal{J}$. Assuming that $\rho_1 < \rho_2$, we find admissible solutions satisfying the initial data (2.28).

**Case 1** $\omega_L = \omega_R$. The solution is a single 1-wave. In particular, if $\omega_L = \omega_R = 0$ and $\alpha_L \neq \alpha_R$, such a 1-wave is the degenerate 1-contact discontinuity.

**Case 2** $\omega_L > \omega_R$. An admissible solution will contain at least one 2-rarefaction (see Theorem 3 for the existence of the rarefaction curves which are illustrated in Fig. 2). In detail, we consider two possibilities, $\omega_L\omega_R \ge 0$ or $\omega_L\omega_R < 0$.

* **Case 2.1** $\omega_L\omega_R \ge 0$. Without loss of generality, we assume that $\omega_L$ and $\omega_R$ are non-positive. From (2.19), an admissible solution is a 2-rarefaction followed by a 1-wave. In addition, the 2-rarefaction monotonic curve $\omega_2(\alpha)$ intersects the 1-wave curve, $\omega(\alpha) = \omega_R$, at a unique point, see also Fig. 2. The uniqueness of the admissible solution is thus obtained.

* **Case 2.2** $\omega_L > 0 > \omega_R$. Due to Lemma 1, the left state must connect to a 2-rarefaction and the right state does the same (otherwise it violates Lemma 1). This property then implies the uniqueness of the admissible solution whose structure depends on the values of $\alpha_L$ and $\alpha_R$, there are three possibilities

- If $\alpha_L \neq 0$ and $\alpha_R \neq 0$, the admissible solution is a 2-rarefaction followed by another 2-rarefaction. See Example 2.3 and Fig. 4a, c. In particular, if $(\alpha_L = 1, \omega_L = 0)$ or $(\alpha_R = 1, \omega_R = 0)$, the solution is a single 2-rarefaction.

- If $\alpha_L = \alpha_R = 0$, the admissible solution is a single 2-rarefaction due to existence of a 2-rarefaction which goes through $\alpha = 0$, see Theorem 3.

- If $(\alpha_L = 0$ and $\alpha_R \neq 0)$ or $(\alpha_L \neq 0$ and $\alpha_R = 0)$, due to Theorem 3 the 2-rarefaction connecting the left state and the one connecting the right state are neither coincide nor intersecting, the admissible therefore contains more than two waves. The unique result solution which satisfies the speed criterion is a 2-rarefaction attached to a 1-contact discontinuity and then followed by another 2-rarefaction. See Example 2.4 and Fig. 4b, d.

**Case 3** $\omega_L < \omega_R$. This case is more technical and the uniqueness of the admissible solution has to be carefully studied. An admissible solution in this case contains a 2-shock, since $\omega$ must increase from $\omega_L$ to $\omega_R$, see Fig. 2 for the shock and rarefaction curves.

We introduce a new variable

$$ \theta = \alpha(\rho_2 - \rho_1) + \rho_1. \qquad (2.31) $$

Since $\alpha \in [0, 1]$, we have $\theta \in [\rho_2, \rho_2]$.

First of all, we look for admissible solutions whose structure is a 2-shock followed by a 1-wave with $U^*$ as intermediate state. The necessary condition $\omega^* = \omega_R \ge 0$ follows from Lemma 1. In order to select an admissible solution, we introduce the speed order criterion in this specific case, which is
$\sigma_2 \le \lambda_1(U^*)$ if the solution is a 2-shock followed by a 1-rarefaction, (2.32)

or $\sigma_2 < \sigma_1$ if the solution is a 2-shock followed by a 1-shock, (2.33)

where $\sigma_2 = \frac{\omega_L + \omega_R}{2(\rho_1 - \rho_2)}$, $\sigma_1 = \frac{\omega_R}{\rho_1 - \rho_2}$, $1 - \frac{\rho_1}{\sigma_1^2} = \frac{\rho_1}{\sigma_2^2}$, $U^* = (\alpha^*, \omega^*)$ such that $\alpha^* = \frac{\rho_1 - \beta^*}{\rho_1 - \rho_2^*}$. For simplicity, we can rewrite the inequalities (2.32) as $\beta^* \le \frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L}$

and (2.33) as $\beta^* \le \frac{2\rho_1\rho_2\omega_R}{\beta R(\omega_R - \omega_L)}$.

Replacing $\alpha$ by $\frac{\beta - \rho_1}{\rho_2 - \rho_1}$ in the 2-shock curve equation in Theorem 2 yields after some calculations that all states $(\alpha, \omega_R)$ which are connected to $U_L$ by a 2-shock satisfy the following quadratic equation

$$ \beta_L(\omega_L - \omega_R)\beta^2 - (\beta_L(\omega_L - \omega_R)(2(\rho_1 + \rho_2) - \beta_L) - 2\rho_1\rho_2\omega_L)\beta - 2\rho_1\rho_2\beta_L\omega_R = 0. \quad (2.34) $$

Let $G(\beta_L, \beta)$ be the left hand side of the Eq. (2.34), then

$$ G(\beta_L, \rho_1) = -\rho_1(\beta_L - \rho_1)(\beta_L(\omega_R - \omega_L) + 2\rho_2\omega_L) \quad \text{and} \quad (2.35) $$

$$ G(\beta_L, \rho_2) = \rho_2(\rho_2 - \beta_L)(\beta_L(\omega_R - \omega_L) + 2\rho_1\omega_L). \quad (2.36) $$

Recall that $\omega_L < \omega_R$, the concave quadratic function $G(\beta_L, \beta)$ may have no solution or more than one solution. We look for a condition on $\omega_R$ such that the Eq. (2.34) has a non-negative solution $\alpha^* \in [0, 1]$, namely $\beta^* \in [\rho_1, \rho_2]$, or equivalently $G(\beta_L, \rho_1)G(\beta_L, \rho_2) \le 0$ because $G(\beta_L, \rho_1)$ and $G(\beta_L, \rho_2)$ can not be negative at the same time. By considering a variation of $\alpha_L$, we get different cases.

* **Case 3.1** If $\alpha_L \in \{0, 1\}$, then $G(\beta_L, \rho_1)G(\beta_L, \rho_2) = 0$ and the Eq. (2.34) may have two solutions with $\beta^* \in [\rho_1, \rho_2]$. Using the speed order criterion, we obtain:

- **Case 3.1.1** If $\alpha_L = 0$, i.e. $\beta_L = \rho_1$ the two potential solutions of (2.34) are $\beta^* = \rho_1$ and $\beta^* = \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$. However $\frac{2\rho_2\omega_R}{\omega_R - \omega_L}$ violates the criteria (2.32) and (2.33). Therefore, only $\beta^* = \rho_1$ is acceptable and the admissible solution is a 2-shock followed by a 1-shock (not a 1-rarefaction) since $\omega_L = \omega_R > 0$ (see 1-shock curves and 1-rarefaction curves on Fig. 2). Such an admissible solution satisfies the criterion (2.33) if and only if $\beta_R \le \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$.

- **Case 3.1.2** If $\alpha_L = 1$ and $\beta_L \ge \frac{2\rho_2\omega_L}{\omega_R - \omega_L}$ (which is equivalent to $\omega_R \ge -\omega_L$ and also to $\frac{2\rho_1\omega_R}{\omega_R - \omega_L} \ge \rho_1$), then there exists a unique solution $\beta^* = \min\{\rho_2, \frac{2\rho_1\omega_R}{\omega_R - \omega_L}\}$ satisfying the criteria (2.32) and (2.33).

* **Case 3.2** If $\alpha_L \in (0, 1)$, or equivalently $\beta_L \in (\rho_1, \rho_2)$, then it is obvious that $G(\beta_L, \rho_1)$ and $G(\beta_L, \rho_2)$ can not be non-positive at the same time, the Eq. (2.34) therefore has at most a solution $\beta^* \in [\rho_1, \rho_2]$.

The first possibility is $G(\beta_L, \rho_1) \le 0$ and $G(\beta_L, \rho_2) \ge 0$, or equivalently

$$ \beta_L \ge \frac{-2\rho_2\omega_L}{\omega_R - \omega_L}, \quad (2.37) $$

i.e. $G(\beta_L, \beta) \ge 0$ implies $\beta \in [\beta^*, +\infty)$.
The second possibility is $G(\mathcal{B}_L, \rho_1) \ge 0$ and $G(\mathcal{B}_L, \rho_2) \le 0$, or equivalently

$$ \mathcal{B}_L \le \frac{-2\rho_1\rho_L\omega_L}{\omega_R - \omega_L}. \qquad (2.38) $$

i.e. $G(\mathcal{B}_L, \mathcal{B}) \ge 0$ implies $\mathcal{B} \in (-\infty, \mathcal{B}^*]$.

As long as the condition (2.37) or (2.38) is satisfied, the speed order criterion (2.32) and (2.33) will help us to select an admissible solution.

- **Case 3.2.1** Assume first that the admissible solution is a 2-shock followed by a 1-rarefaction and use the speed order criterion (2.32). We denote $H(\mathcal{B}_L) = G \mathcal{B}_L \frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L}$ and first prove that $H(\mathcal{B}_L) > 0$ for all $\mathcal{B}_L \in (\rho_1, \rho_2)$. Writing $H(\mathcal{B}_L)$ explicitly,

$$ H(\mathcal{B}_L) = -4\rho_1\rho_2\mathcal{B}_L\omega_R + (\mathcal{B}_L(\omega_R - \omega_L)(2(\rho_1 + \rho_2) - \mathcal{B}_L) + 2\rho_1\rho_2\omega_L) \frac{\overline{2\rho_1\rho_2\omega_R}}{\omega_R - \omega_L} $$

Calculating and evaluating $H(\rho_1)$,

$$ \begin{aligned} H(\rho_1) &= -4\rho_1^2\rho_2\omega_R + \rho_1(\rho_1(\omega_R - \omega_L) + 2\rho_2\omega_R) \frac{\overline{2\rho_1\rho_2\omega_R}}{\omega_R - \omega_L} \\ &\geq -4\rho_1^2\rho_2\omega_L + 2\rho_1\rho_2\rho_1(\omega_R - \omega_L) \frac{\overline{2\rho_1\rho_2\omega_R}}{\omega_R - \omega_L} \\ &= 0. \end{aligned} \qquad (2.39) $$

The inequality (2.39) is obtained by the Cauchy's inequality for two non-negative numbers $\rho_1(\omega_R - \omega_L)$ and $2\rho_2\omega_R$. The result $H(\rho_2) \ge 0$ is obtained similarly. Moreover, considering $\mathcal{B}_L$ as a variable of the quadratic function $H(\mathcal{B}_L)$ whose highest order's coefficient is negative and both $H(\rho_1)$ and $H(\rho_2)$ are non-negative, we achieve $H(\mathcal{B}_L) > 0$ for all $\mathcal{B}_L \in (\rho_1, \rho_2)$. This result shows that the condition (2.37) satisfies the criterion (2.32) (since $G \mathcal{B}_L \frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L} = H(\mathcal{B}_L)$ and $G(\mathcal{B}_L, \mathcal{B}) \ge 0$ implies $\mathcal{B} \in [\mathcal{B}^*, +\infty)$, so that $\mathcal{B}^* \le \frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L}$ while the condition (2.38) is impossible to satisfy (since it implies $\mathcal{B}^* \ge \frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L}$ because $G(\mathcal{B}_L, \mathcal{B}) \ge 0$ implies $\mathcal{B} \in (-\infty, \mathcal{B}^*]$ in this case).

- **Case 3.2.2** Assume that the admissible solution is a 2-shock followed by a 1-shock and use the speed order criterion (2.33). If $\frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L} > \rho_2$, the criterion (2.33) is always satisfied. So we will merely consider

$\frac{2\rho_1\rho_2\omega_R}{\omega_R - \omega_L} \le \rho_2$, where we define $M(\mathcal{B}_L) = G \mathcal{B}_L \frac{2\omega_R\rho_1\rho_2}{\omega_R - \omega_L}$. Rewrite $M(\mathcal{B}_L)$ as the following

$$ M(\mathcal{B}_L) = \frac{2\rho_1\rho_2\omega_R}{\omega_R^2(\omega_R - \omega_L)} \left( (\mathcal{B}_R\omega_L - \mathcal{B}_L\omega_R) + \mathcal{B}_L\mathcal{B}_R(\omega_R - \omega_L)(2\rho_1 + 2\rho_2 - \mathcal{B}_L - \mathcal{B}_R) \right). $$
We compute

$$
\begin{align*}
M(\rho_1) &= \frac{2\rho_1\rho_2\omega_R}{\beta_R} (\rho_1 - \beta_R) \quad \rho_1 - \frac{2\rho_1\rho_2\omega_R}{\beta_R(\omega_R - \omega_L)}, \\
M(\rho_2) &= \frac{2\rho_1\rho_2\omega_R}{\beta_R} (\rho_2 - \beta_R) \quad \rho_2 - \frac{2\rho_1\rho_2\omega_R}{\beta_R(\omega_R - \omega_L)}.
\end{align*}
$$

The condition (2.37) implies $M(\rho_1) \ge 0$ and $M(\rho_2) \ge 0$, and then obviously $M(\boldsymbol{\beta}_L) > 0$, $6\boldsymbol{\beta}_L \in (\rho_1, \rho_2)$ (since the second order polynomial $\boldsymbol{\beta}_L \to M(\boldsymbol{\beta}_L)$ is concave), this result satisfies the criterion (2.33).

On the other hand, the condition (2.38) violates the criterion (2.33)
as long as we assume the solution has more than one wave. By the
continuity, this solution is not admitted.

We summarize the structure of the admissible solution according to the initial
data.

• If $\alpha_L \in (0, 1]$ and $\boldsymbol{\beta}_L \ge \frac{2\rho_2\omega_L}{\omega_R - \omega_L}$, the solution is a 2-shock followed by a 1-wave.

• If $\alpha_L = 0$ and $\boldsymbol{\beta}_R \le \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$ the solution is a 2-shock followed by a 1-shock ($\alpha^* = 0$).

Similarly, we obtain the following results

• If $\alpha_R \in (0, 1]$ and $\boldsymbol{\beta}_R \ge \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$ the solution is a 1-wave followed by a 2-shock.

• If $\alpha_R = 0$ and $\boldsymbol{\beta}_L \le \frac{2\rho_2\omega_L}{\omega_R - \omega_L}$ the solution is a 1-shock followed by a 2-shock ($\alpha^* = 0$).

Before continuing the proof, we can conclude that if $\alpha_L \ne 0$ and $\alpha_R \ne 0$, the solution consisting of a 2-shock (resp. 1-wave) followed by a 1-wave (resp. a 2-shock) is admissible if $\boldsymbol{\beta}_L \ge \frac{-2\rho_2\omega_L}{\omega_R - \omega_L}$ (resp. $\boldsymbol{\beta}_R \ge \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$).

The rest of our proof considers the initial data which are not studied above,
i.e. $\alpha_L \alpha_R \ne 0$ and $\boldsymbol{\beta}_L < \frac{-2\rho_2\omega_L}{\omega_R - \omega_L}$ and $\boldsymbol{\beta}_R < \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$. According to Lemma 1,
an admissible solution must be a 1-wave followed by a 2-shock connected to
another 1-wave. Let us denote the two intermediate states ordered from the
left to the right by $U^*(\alpha^*, \omega_L)$ and $U^{**}(\alpha^{**}, \omega_R)$.

• If $\alpha^* \ne 0$ and $\alpha^{**} \ne 0$, due to the previous results a 1-wave followed by a 2-shock is admissible if $\boldsymbol{\beta}^{**} \ge \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$ and this 2-shock followed by another 1-wave is admissible if $\boldsymbol{\beta}^* \ge \frac{-2\rho_2\omega_L}{\omega_R - \omega_L}$. Both $\boldsymbol{\beta}^{**} \ge \frac{2\rho_2\omega_R}{\omega_R - \omega_L}$ and $\boldsymbol{\beta}^* \ge \frac{-2\rho_2\omega_L}{\omega_R - \omega_L}$ are satisfied if and only if $\omega_L = -\omega_R$ and $\boldsymbol{\beta}^* = \boldsymbol{\beta}^{**} = \rho_2$, this condition however violates the speed order criterion.

• If $(\alpha^* \ne 0$ and $\alpha^{**} = 0$) or $(\alpha^* = 0$ and $\alpha^{**} \ne 0)$, Theorem 2 shows that it is impossible.

$\alpha^* = \alpha^{**} = 0$ is admissible, this solution satisfies all criteria of speed order. See Fig. 3b for the construction of such an admissible solution. An example is performed in Example 2.5, see Fig. 5a, c.
<figure><img src="image_6.png" /><figcaption>FIGURE 3. The admissible solution 2-shock 1-shock (a), 1-shock 2-shock 1-shock (b)</figcaption></figure>

The continuity to the initial condition of the admissible solution is obvious
due to the construction of such an admissible solution.

We would like to remark that an admissible solution of a Riemann problem
may admit a 2-shock wave which connects a left state which is on a branch
of a hyperbola to a right state located on the other branch, see Fig. 3a. In
the following section, this 2-shock will be called a **non classical shock wave**. It
turns out that such a non classical shock wave is not easily captured by some
classical numerical methods.

**2.5.2. Key examples and important comments.** The first three Riemann problems (Examples 2.3, 2.4, 2.5) we present do not give rise to classical weak solution made of two waves of different families while the last one (Example 2.6) produces a non classical shock wave. The first Riemann problem gives rise two rarefactions of the 2-family, the second one leads to three waves (a 2-wave followed by and attached to a 1-wave followed by a 2-wave) and the third one produces a pure phase ($\alpha = 0$) starting from a mixture and contains three shocks. It is interesting to notice that in the third Riemann problem (Example 2.5) the value of the velocity of the vanishing phase does not necessarily equal the one of the non vanishing phase. Finally, the last Riemann problem illustrates an admissible solution consisting in two waves of different families whose one is a non classical wave.

*Example 2.3.* The configuration $\alpha_L = \alpha_R = 0.5$, $\omega_L = -\omega_R = 3$ generates a pure gas intermediate value $U^* = (1, 0)$. The admissible solution consists of two rarefaction waves. See Fig. 4a, c.

*Example 2.4.* The configuration $\alpha_L = 0, \alpha_R > 0, \omega_L > 0 > \omega_R$ is an example where the solution is a 2-rarefaction touching the degenerate 1-shock followed by another 2-rarefaction such that $\lambda_2(U^*) = \lambda_2(U^{**}) = 0$ and the speed of degenerate 1-shock is also zero, where $U^* = (0, 0)$, $U^{**} = (1, 0)$ are the intermediate values. See Fig. 4b, d.
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Example 2.3: 2-rarefaction 2-rarefaction (a, c). Example 2.4: 2-rarefaction 1-shock 2-rarefaction (b, d)</figcaption></figure>

*Example 2.5.* The configuration $\alpha_L = \alpha_R = 0.5$, $\omega_L = -\omega_R = -5$ generates a pure liquid and the solution consists in three shocks (a 1-shock connects to 2-shock followed by another 1-shock). See Fig. 5a, c.

*Example 2.6.* The configuration $\alpha_L = 0.8$, $\alpha_R = 0.5$, $\omega_L = -3$, $\omega_R = 5$ generates a non classical shock and the solution consists in two shocks (a non classical 2-shock connects to a 1-shock). See Fig. 5b, d.

### 3. Numerical study

We now investigate the numerical simulation of the system (1.14) and show that the basic Roe scheme fails to capture the expected dynamics whereas the Godunov scheme and the Roe scheme with a Harten type correction capture the analytic solution. However in the non classical shock wave (corresponding to a passage through the domain $\mathcal{J}_+$ and $\mathcal{J}_-$, both of these schemes show oscillations. We then propose a reconstruction method, see in [28–32] and
<figure><img src="image_8.png" /><figcaption>FIGURE 5. Example 2.5: 1-shock 2-shock 1-shock (a, c); Example 2.6: 2-rarefaction 1-shock 2-rarefaction (b, d)</figcaption></figure>

references therein, which significantly improves the numerical result in this case.

We consider a uniform mesh of the computational domain $[0, 1]$ whose $N$ cells are centered at $x_i, i = 1, \dots, N$. The space step $\Delta x = x_i - x_{i-1}$ is constant whereas the time step $\Delta t(U^n) > 0$ depends on the discrete field $U^n = (U_i^n)_{i=1,\dots,N}$ which approximates the exact solution $U(x, t)$ at cells $i$ and time $t^n = \sum_{k=0}^{n-1} \Delta t(U^k)$. The time step should satisfy the following CFL condition in order to ensure the stability of the explicit schemes: $\Delta t \le \max_i \{\lambda_1(U_i^n) \lambda_2(U_{i+1})\}$, where $\lambda_k(U_{i\nu}, U_{i+1})$, $k = 1, 2$ is the largest value of $|\lambda_k|$ on the path connecting $U_i$ to $U_{i+1}$ using the rarefactions and admissible shock waves computed in Theorems 2 and 3. We point out that $\lambda_k(U_{i\nu}, U_{i+1})$ may be different from $|\lambda_k(U_i)|$ and $|\lambda_k(U_{i+1})|$ because the characteristic fields are non genuinely nonlinear. Denote $U^*$ intermediate states, then

$$ \lambda_k(U_{i\nu}, U_{i+1}) = \max_{U^*} \{|\lambda_k(U_i)|, |\lambda_k(U_{i+1})|, |\lambda_k(U^*)|\}, \quad k = 1, 2. \quad (3.40) $$
We consider conservative finite volume schemes in the following explicit form:

$$U_i^{n+1} = U_i^n - \frac{\Delta t}{\Delta x} \Phi_{i+1/2}^n - \Phi_{i-1/2}^n, \quad (3.41)$$

where $\Phi_{i+1/2}^n$ is the numerical flux function at the interface between cells $i$ and $i+1$, and at time $t^n$. We compute the numerical flux $\Phi_{i+1/2}^n$ using one of the following strategy.

### 3.1. Godunov scheme

$$\Phi_{i+1/2}^n = \mathbf{F}_{j+1/2}^{\text{God}} = F(U_i^*(U_i^n, U_{i+1}^n)),$$

where $U_i^*(U_i^n, U_{i+1}^n)$ is the value taken by the solution of the Riemann problem between the left state $U_i^n$ and the right state $U_{i+1}^n$ at the interface.

### 3.2. Roe scheme with a Harten type correction

$$\Phi_{i+1/2}^n = \mathbf{F}_{j+1/2}^{\text{Har}} = \frac{F(U_i^{(n)} + F(U^n))}{2} \Big|_{\Delta U_{i+1/2} = \frac{2}{\Delta x}} - |A^{\text{Roe}}(U_i^n, U_{i+1}^n)| + \text{har}_{i,i+1}^n \text{Id}$$

where $\Delta U_{i+1/2} = U_{i+1}^n - U_i^n$, $A^{\text{Roe}}(U_i^n, U_{i+1}^n)$ is the Roe matrix, (see (3.42) for the derivation of a Roe matrix), and $\text{har}_{i,i+1}^n = C \max(|\lambda_1(U_i^n) - \lambda_1(U_{i+1}^n)|, |\lambda_2(U_i^n) - \lambda_2(U_{i+1}^n)|)$. If $C=0$ we recover the standard Roe scheme. However it is well-known that the Roe scheme may capture non admissible solutions (see [33]). Hence we used a constant value $C = \frac{1}{5}$ to include a Harten type entropic correction in the Roe scheme.

A Roe matrix $A^{\text{Roe}}(U_L, U_R)$ for the system (1.14) and two states $U_L, U_R \in \mathbb{R}$ is a diagonalisable matrix such that

$$F(U_L) - F(U_R) = A^{\text{Roe}}(U_L, U_R)(U_L - U_R)$$

$$A^{\text{Roe}}(U, U) = \nabla F(U)$$

After some calculations, we obtained and used the following Roe matrix

$$A^{\text{Roe}}(U_L, U_R) = \begin{pmatrix} A_{1,1} & A_{1,2} \\ A_{2,1} & A_{2,2} \end{pmatrix}, \quad (3.42)$$

where

$$A_{1,1} = \frac{w_L + w_R}{2(\rho_1 - \rho_2)} \left( 1 - \frac{\rho_1 \rho_2}{\beta_L \beta_R} \right),$$

$$A_{1,2} = \frac{1}{2(\rho_1 - \rho_2)} \frac{(\beta_L - \rho_1)(\beta_L - \rho_2)}{\beta_L} + \frac{(\beta_R - \rho_1)(\beta_R - \rho_2)}{\beta_R},$$

$$A_{2,1} = 0,$$

$$A_{2,2} = \frac{\omega_L + \omega_R}{2(\rho_1 - \rho_2)},$$

and $\beta_k = \alpha_k(\rho_2 - \rho_1) + \rho_1, k = L, R$.
### 3.3. Reconstruction scheme

As we will see in the next section devoted to the numerical experiments, the Roe scheme (with Harten entropy correction) and the Godunov scheme provide good numerical results when the solution stays in one of the two domains $7_-$ and $7_+$, or when the solution belongs to $7_+$ but crosses the line $w = 0$ with a pure phase for which $\alpha$ is a constant (more precisely, $\alpha = 0$ or $\alpha = 1$). In these situations, the model is strictly hyperbolic and the characteristic fields are genuinely non linear. In other words, the model satisfies the usual properties which ensure that the solution is classical and do not develop non classical shocks. On the contrary, we will see that when the solution goes from $7_-$ to $7_+$, both methods fail in producing good numerical results. More precisely, the methods are not able to properly compute the non classical discontinuity joining the two sets $7_-$ and $7_+$. Recall that the model is hyperbolic in $7_- \cup 7_+$ but the characteristic fields are not genuinely non linear (nor linearly degenerate) in this domain, which gives rise to a non classical shock.

From a numerical point of view, it is now very well-known that approximating non classical shocks is a challenging issue because of the dependence on the underlying diffusion mechanisms. Standard techniques (like Godunov's method and Roe's method) are useless and a deeper analysis shows that the failure of these techniques can be related to the (un)control of the underlying numerical diffusion. We refer for instance the reader to [27,34-39]... and the references therein.

With this in mind and in order to obtain numerical solutions in agreement with the exact ones, we suggest to develop a new method able to control the numerical diffusion across the non classical discontinuity. The strategy is based on in-cell discontinuous reconstructions and was first proposed in [40] to solve transport equations with no numerical diffusion (see also [28,41] for non linear scalar equations). More precisely, the aim of this non-dissipative finite volume approach is to follow exactly the propagation of a given discontinuity by means of in-cell reconstructions. Then, the usual projection onto the piecewise constant solutions of the time evolution of the reconstructed discontinuities makes the capture of isolated discontinuities exact. In particular, the numerical diffusion is perfectly controlled and made of one point only.

Extending this approach to our model is not immediate since we first face a *system* of conservation law (the above papers consider scalar equations), and secondly we aim at following the non classical discontinuity of a *non linear* model and not a contact discontinuity of a (linear) transport equation as above. However, the in-cell discontinuous reconstruction strategy has already been extended to the computation of non-classical shocks for non linear scalar equations, see for instance [30,31], but also to the system case, see for instance [29,32,42] and references therein for more details.

In what follows, we therefore follow the same approach and propose an in-cell discontinuous reconstruction strategy of the non classical discontinuities of our model. Again, the strategy allows for a sharp computation of these non classical discontinuities. By sharp, we mean that such isolated discontinuities are exactly computed with only one point of numerical diffusion, the value of
which corresponds to the mean value of the exact solution in the corresponding cell. It is clearly the best result achievable by a conservative scheme.

Let us now present the main ideas of the method and summarize the computation of the numerical flux function, which is what is really needed to implement the method. For more details, we refer the reader to the references above.

Considering our system, the unknown variables is $U = (\alpha, \omega)$ but only $\alpha$ varies in a 1-shock, whereas the 2-shock corresponds to the case where both $\alpha$ and $\omega$ vary at the jump. Numerical methods in general capture well the shock when only one of the two variables varies, and show difficulties in the cases where the two variables vary at the jump. This is especially the case of the non classical 2-shock i.e. when the two states located on different branches of a hyperbola. Therefore, we will develop the in-cell reconstruction corresponding to this specific configuration.

**In-cell reconstruction criteria.** First of all, one has to decide when it is necessary to introduce a discontinuous reconstruction in a given cell $j$. In other words, it amounts to wonder when a non classical discontinuity is expected to be present in the cell $j$. With this in mind, we suggest to consider the Riemann solution associated with the left and right initial states $U_{j-1}$ and $U_{j+1}$. We denote $RP(U_{j-1}, U_{j+1})$ this Riemann solution and our detection criteria of non classical discontinuity in the cell $j$ is the following. If the solution of $RP(U_{j-1}, U_{j+1})$ contains an admissible discontinuity between the left state $U_j^-= (\alpha_j^-(U_{j-1}, U_{j+1}), \omega_j^-(U_{j-1}, U_{j+1}))$ and right state $U_j^+ = (\alpha_j^+(U_{j-1}, U_{j+1}), \omega_j^+(U_{j-1}, U_{j+1}))$ such that $\alpha_j^- \neq \alpha_j^+$ and $\omega_j^- \neq \omega_j^+$, we propose a discontinuous in-cell reconstruction between $\dot{U}_j^-$ and $\dot{U}_j^+$ in cell $j$ as shown on Fig. 6. Otherwise, $U_j^- = U_j^+ = U_j$.

**Location of the reconstructed discontinuity and conservativity.** It is important to note that the in-cell reconstructed discontinuity will not be necessarily located at the same place for both variables $\alpha$ and $\omega$. Indeed, our objective is to preserve conservation for both variables $\alpha$ and $\omega$, which means that we require the average value of the reconstructed discontinuity to be equal for each variable to the initial constant value on the cell $j$. We thus define the coefficient $\vartheta_j^\alpha$ (resp. $\vartheta_j^\omega$) such that the distance from $x_{j-1/2}$ to the discontinuity of the variable $\alpha$ (resp. variable $\omega$) is $\vartheta_j^\alpha \Delta x$ (resp. $\vartheta_j^\omega \Delta x$), see Fig. 6. We would like to locate the discontinuities of $\alpha$ and $\omega$ in a way that yields a conservative scheme, so that $\vartheta_j^\alpha$ and $\vartheta_j^\omega$ must satisfy

$$ \left\{ \begin{array}{l} \vartheta_j^\alpha \alpha_j^- + 1 - \vartheta_j^\alpha \alpha_j^+ = \alpha_j^n, \\ \\ \vartheta_j^\omega \omega_j^- + 1 - \vartheta_j^\omega \omega_j^+ = \omega_j^n. \end{array} \right. $$

At this stage it may happen that the values of $\vartheta$ given by these relations lie outside of $[0, 1]$, meaning that the reconstructed discontinuity should be located outside of the cell $j$. In this case, we simply suggest to give up the reconstruction process in the cell $j$, more precisely we set
<figure><img src="image_9.png" /><figcaption>FIGURE 6. Reconstruct $U_j$ by $U^- = (\alpha_j^-, \omega_j^-)$, $U_j^+ = (\alpha_j^+, \omega_j^+)$ and $d_j^\alpha = \vartheta_j^\alpha \Delta x$, $d_j^\omega = \vartheta_j^\omega \Delta x$</figcaption></figure>

• If $\vartheta_j^\alpha \notin [0, 1]$ : no reconstruction for $\alpha$ i.e.

$$ \alpha_j^+ = \alpha_j^- = \alpha_j^n. $$

• If $\vartheta_j^\omega \notin [0, 1]$ : no reconstruction for $\omega$ i.e.

$$ \omega_j^+ = \omega_j^- = \omega_j^n. $$

• If $\vartheta_j^\alpha \in [0, 1]$ (and/or $\vartheta_j^\omega \in [0, 1]$), reconstructing for $\alpha$ (and/or $\omega$).

**Definition of the numerical flux using the reconstructed discontinuity.** In this paragraph, we give the definition of the numerical fluxes taking into account the presence of in-cell reconstructed discontinuities. The idea is to evaluate the numerical fluxes by considering the propagation of the reconstructed discontinuities and by calculating the value of the exact flux either on the left or the right state of the discontinuity depending on its (varying in time) position inside or outside the cell $j$. With this in mind, let us denote $\sigma_j$ the exact value of the speed of propagation of the discontinuity $U_j^-$, $U_j^+$. We then compute the numerical flux function between $t^n$ and $t^n + \Delta t$ by using the reconstructed discontinuities rather than the average values. More precisely, if $\sigma_j > 0$ (resp. $\sigma_j < 0$), we are going to calculate the flux at interface $j + 1/2$ (resp. $j - 1/2$) by considering that the numerical flux equals the exact flux evaluated on the right value $U_j^+$, until the corresponding discontinuity reaches the interface $j + 1/2$ (resp. $j - 1/2$), and the exact flux evaluated on the left value $U_j^-$ afterwards. Therefore, such a flux function will be computed relying on the speed of shock propagation $\sigma_j$ of the reconstructed discontinuity and on the times $\Delta t_{j+1/2}^\alpha$, $\Delta t_{j+1/2}^\omega$ needed by this discontinuity to reach the interface $j \pm 1/2$ depending on the sign of $\sigma_j$. More explicitly:

• If $\sigma_j > 0$. Denote $\Delta t^\omega = \frac{(1-\vartheta_j^\omega)\Delta x}{t^{n+1/2}}$, $\Delta t^\alpha = \frac{(1-\vartheta_j^\alpha)\Delta x}{j+1/2}$, $\sigma_j^\alpha$ is computed by using $U_j^-$ and $U_j^+$.
— If $\vartheta_j^\alpha \le \vartheta_j^\omega$, then

$$ \Delta t \mathbf{F}^{\text{Rec}}_{j+1/2} = \min_j \Delta t^\omega_{j+1/2}, \Delta t F(U_j^+) + \max_j \Delta t - \Delta t^\alpha_{j+1/2}, \begin{array}{l} 0, F(U_j^-) \\ \Delta t^\omega_{j+1/2}, \Delta t - \Delta t^\omega_{j+1/2}, \end{array} $$

— If $\vartheta_j^\alpha > \vartheta_j^\omega$, then

$$ \Delta t \mathbf{F}^{\text{Rec}}_{j+1/2} = \min_j \Delta t^\alpha_{j+1/2}, \Delta t F(U_j^+) + \max_j \Delta t - \Delta t^\omega_{j+1/2}, \begin{array}{l} 0, F(U_j^-) \\ \Delta t^\omega_{j+1/2}, \Delta t - \Delta t^\omega_{j+1/2}, \end{array} $$

$$ \cdot \text{ If } \sigma_j < 0. \text{ Denote } \Delta t^\omega_{j-1/2} = \frac{\vartheta_j^\omega \Delta x}{-\sigma_j}, \quad \Delta t^\alpha_{j-1/2} = \frac{\vartheta_j^\alpha \Delta x}{-\sigma_j}. $$

The numerical flux function $\mathbf{F}^{\text{Rec}}_{j-1/2}$ is computed by using $U_j^-$ and

$U_j^+$. If $\vartheta_j^\alpha \le \vartheta_j^\omega$, then

$$ \Delta t \mathbf{F}^{\text{Rec}}_{j-1/2} = \min_j \Delta t^\alpha_{j-1/2}, \Delta t F(U_j^-) + \max_j \Delta t - \Delta t^\omega_{j-1/2}, \begin{array}{l} 0, F(U_j^+) \\ \Delta t^\omega_{j-1/2}, \Delta t - \Delta t^\omega_{j-1/2}, \end{array} $$

— If $\vartheta_j^\alpha > \vartheta_j^\omega$, then

$$ \Delta t \mathbf{F}^{\text{Rec}}_{j-1/2} = \min_j \Delta t^\omega_{j-1/2}, \Delta t F(U_j^-) + \max_j \Delta t - \Delta t^\alpha_{j-1/2}, \begin{array}{l} 0, F(U_j^+) \\ \Delta t^\omega_{j-1/2}, \Delta t - \Delta t^\omega_{j-1/2}, \end{array} $$

## 4. Numerical results

We present some numerical results obtained with the constant densities $\rho_1 = 1$, $\rho_2 = 3$, which give a good overview of the wave structure. Moreover, the simulation is implemented on a spacial domain $[0, 1]$, uniform mesh with space step $\Delta x$ and CFL number is less than or equal to 1. The time step is defined by

$$ \Delta t = \text{CFL} \times \frac{\Delta x}{\max\{\lambda_1(U_{i\uparrow}, U_{i\downarrow+1}), \lambda_2(U_{i\uparrow}, U_{i\downarrow+1})\}} \quad (4.43) $$

where $\lambda_k(U_{i\uparrow}, U_{i\downarrow+1})$, $k = 1, 2$ are defined by (3.40).

We first show in Sect. 4.1 that the Godunov scheme and the Roe scheme with Harten type correction are able to capture the non classical wave structure joining two states in different domains $\mathbf{7}_-$ and $\mathbf{7}_+$ in the Riemann problem involving a pure phase intermediate state (Examples 2.3 and 2.5). These schemes however show strong oscillation in capturing the non classical 2-shock wave in Example 2.6, see this configuration in Fig. 3a, whereas the reconstructing method show very good results, Fig. 9.

Then in Sect. 4.2 we simulate the classical problem of phase separation under gravity.
<figure><img src="image_10.png" /><figcaption>FIGURE 7. Solution of the Riemann problem at time $t = 0.15$ for the initial data $\alpha_1 = \alpha_2 = 0.5$ and $\omega_L = -\omega_R = 5$; 100 cells and CFL = 0.9</figcaption></figure>

## 4.1. The Riemann problem

The Riemann problem consists in solving the system (1.14) with $K = 0, S = 0$ and the initial data

$$ U(x, 0) = \begin{cases} (\beta_L, \omega_L) & \text{if } x \le 0, \\ (\beta_R, \omega_R) & \text{if } x > 0. \end{cases} \qquad (4.44) $$

From Theorem 4, this problem admits a unique admissible solution satisfying Liu's criterion with $\alpha_1, \alpha_2 \in [0, 1]$. In the special case where $\omega_L = -\omega_R$, the solution involves a pure phase: the lighter if $\omega_L > 0$, and the heavier if $\omega_L < 0$. It consists of two transonic rarefactions in the former case and three shocks waves in the latter. We present in Figs. 7 and 8, the numerical results obtained using the Godunov scheme, the Roe scheme, and the Roe scheme with the Harten type entropy fix presented at Sect. 3. In the first case $\omega_L > 0$, Fig. 7, the original Roe scheme is unable to capture the admissible solution and captures instead an inadmissible shock, i.e. does not satisfy the Liu criterion. We remark that the velocity of the liquid in the pure gas region is smooth. In the second case $\omega_L < 0$, Fig. 8, the original Roe scheme and others schemes capture well the pure liquid state. In this case, the gas velocity in the pure liquid region includes of three shocks and is bounded. The velocity of the vanishing phase in both of cases is not necessarily equal to the one of the pure phase.
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Solution of the Riemann problem at time $t = 0.15$ for the initial data $\alpha_1 = \alpha_2 = 0.5$ and $\omega_L = -\omega_R = -5$; 100 cells and CFL = 0.9</figcaption></figure>

The third numerical simulation of the Riemann problem is the non classical 2-shock wave as in Fig. 3a. We recall that this 2-shock goes through the domain $7_-$ and $7_+$, connects the left state $U_L$ to the intermediate state $U_{int}$ such that each component of $U_L$ and $U_{int}$ is different and the speed propagation of the 2-shock is not equal to zero. These challenges lead to oscillations given by both the Godunov scheme and the Roe scheme with Harten entropy fix while the Roe scheme without entropy fix yields strong oscillations. The admissible solution is well captured only by the reconstruction method, see Fig. 9 (a uniform mesh with 100 cells), Fig. 11 (a uniform mesh with 500 cells) and Fig. 10 for the convergence of these schemes.

## 4.2. The phase separation under gravity

This is a classical test case in the assessment of numerical methods in the modeling of counter-current two phase flows with steep transition (see [2]). We consider the model (1.1) with $g = -10m/s^2$, $K = 0$ and $x \in [0, 1]$ with the initial data $u_1(x, 0) = 0$, $u_2(x, 0) = 0$, $\alpha_1(x, 0) = 0.5$, $\alpha_2(x, 0) = 0.5$ and boundary data $u_1(0, t) = u_2(0, t) = u_1(1, t) = u_2(1, t) = 0$. The transient result in Fig. 12 (left) shows that the Roe scheme captures an inadmissible shock departing from $x = 1$. This is consistent with the results shown in the previous section since the Riemann problems at the walls yield pure phases
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Solution of the Riemann problem at time $t = 0.12$ for the initial data $\alpha_L = 0.9$, $\alpha_R = 0.3$ and $\omega_L = -3$, $\omega_R = 5$; 100 cells and CFL = 0.5</figcaption></figure>

intermediate states and a transonic rarefaction fan for the lighter phase. However, the Roe scheme with Harten entropic correction gives a similar result to the Godunov scheme, both of them being consistent with the analysis of the Riemann problem (Figs. 13, 14).

Both the physical and mathematical analysis agree that the expected stationary state for the volume fraction and velocities should satisfy $\alpha_1 = 0$ on $[0, 0.5]$ and $\alpha_1 = 1$ on $[0.5, 1]$ velocity $u_2 = 0$ on $[0, 0.5]$ and $u_1 = 0$ on $[0.5, 1]$. However there is a debate as to what should be the value of $u_2$ (resp. $u_1$) on $[0.5, 1]$ (resp. $[0, 0.5]$) since in that region the liquid (resp. the gas) is absent). In our model we can compute the stationary velocity of the liquid which is not zero hence there is no mechanical equilibrium. However there does not exist a stationary value for the gas velocity on the whole of domain, we refer the reader to Appendix 5 for details. During the numerical simulation, all schemes except the original Roe scheme captured well the vanishing velocity of liquid in the pure gas domain as well as the (non stationary) vanishing velocity of gas in the pure liquid domain.
<figure><img src="image_13.png" /><figcaption>FIGURE 10. Solution of the Riemann problem at time $t = 0.12$ for the initial data $\alpha_L = 0.9$, $\alpha_R = 0.3$ and $\omega_L = -3$, $\omega_R = 5$; 500 cells and CFL = 0.5</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 11. Convergence curves (mesh refinement for the Riemann problem with initial data $\alpha_L = 0.9$, $\alpha_R = 0.3$ and $\omega_L = -3$, $\omega_R = 5$; CFL = 0.5)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 12. Volume fraction $\alpha_1$ for the sedimentation problem. a Transient, b stationary</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13. The velocity $u_1$ for the sedimentation problem. a Transient, b stationary</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 14. The velocity $u_2$ for the sedimentation problem. a Transient, b stationary</figcaption></figure>
## 5. Appendix: Vanishing velocity of the phase separation under gravity

This appendix demonstrates the derivations of the vanishing velocity, i.e the stationary state of the phase separation under gravity in Sect. 4.2.

We consider the stationary state of the 2×2 system (1.13) assuming that $g < 0$ and $\rho_1 < \rho_2$. The first equation in (1.13) at stationary state and the boundary condition together with (1.7) yield

$$ \alpha_1 u_1 = 0, $$

$$ \alpha_2 u_2 = 0. $$

We seek a solution consisting of two zones. A bottom zone with pure phase 2: $\alpha_2 = 1$ and constant velocity $u_2 = 0$ in the region $x \in [0, 0.5]$ and a top zone with pure phase 1: $\alpha_1 = 1$ and constant velocity $u_1 = 0$ in the region $x \in [0.5, 1]$. We are going to use the second equation of the system (1.13) to determine the vanishing velocity of phase 2 in the region $x \in [0.5, 1]$. The second equation of the system (1.13), using the physical variables $u_1$ and $u_2$ (or equivalently equation 1.10) is

$$ \partial_x \rho_1 \frac{u_1^2}{2} - \rho_2 \frac{u_2^2}{2} + \frac{\rho_1 \rho_2}{2(\rho_1 - \rho_2)} (u_1 - u_2)^2 = (\rho_1 - \rho_2)g. \quad (5.45) $$

Integrating (5.45) we obtain the **two phase Bernoulli's principle**:

$$ \rho_1 \frac{u_1^2}{2} - \rho_2 \frac{u_2^2}{2} + \frac{\rho_1 \rho_2}{2(\rho_1 - \rho_2)} (u_1 - u_2)^2 - (\rho_1 - \rho_2)gx = \text{constant}. \quad (5.46) $$

In order to compute the vanishing phase velocity of phase 2, we remark that the velocities at the walls $x = 1$ are $u_1 = u_2 = 0$. Hence the constant in (5.46) equals $-(\rho_1 - \rho_2)g$, and since $u_1 = 0$ for $x \in [0.5, 1]$ the two phase Bernoulli's principle becomes

$$ \frac{\rho_2^2}{\rho_1 - \rho_2} \frac{u_2^2(x)}{2} = (\rho_1 - \rho_2)g(x - 1) \quad \text{for } x \in [0.5, 1]. $$

Hence

$$ u_2(x) = -2 \left(1 - \frac{\rho_1}{\rho_2} g(x - 1)\right) \quad \text{for } x \in [0.5, 1]. \quad (5.47) $$

The velocity profile is therefore not constant and furthermore shows a discontinuity at the interface $x = 0.5$.

We remark that we cannot determine a stationary vanishing velocities for phase 1 in the region $x \in [0, 0.5]$. Indeed since $u_2 = 0$ for $x \in [0, 0.5]$, the two-phase Bernoulli's principle takes the form

$$ \frac{\rho_1^2}{\rho_1 - \rho_2} \frac{u_1^2(x)}{2} = (\rho_1 - \rho_2)gx \quad \text{for } x \in [0, 0.5]. \quad (5.48) $$

Since $g = -10m/s < 0$, the Eq. (5.48) has no solution as it yields $u_1^2 < 0$. In the stationary state of the sedimentation problem, the vanishing liquid velocity in the pure gas region is given by the formula (5.47), however the
stationary vanishing gas velocity in the pure liquid phase need not exist as explained above. This result is consistent with the other similar results in the simulation of compressible two-phase flows, for example in [43–45], where the authors can capture and study the vanishing liquid velocity but no information is given about the gas vanishing velocity. It is in fact a difficult problem encountered in the transition between the two-phase flow model to the single one.

## References

[1] Cordier, F., Degond, P., Kumbaro, A.: Phase appearance or disappearance in two-phase flows. J. Sci. Comput. **58**(1), 115–148 (2014)

[2] Jeong, J.J., Yoon, H.Y., Cho, H.K., Jim, J.: A semi-implicit numerical scheme for transient two-phase flows on unstructured grids. Nucl. Eng. Des. **238**, 3403–3412 (2008)

[3] Bestion, D.: The appearance and disappearance in the CATHARE code. Trends in Numerical and Physical Modeling for Industrial Multiphase Flows, Cargese. 27th–29th Sept (2000)

[4] Drew, D.A., Passman, S.L.: Theory of Multicomponents Fluids. Springer, New York (1999)

[5] Ishii, M.: Thermo-Fluid Dynamic Theory of Two-Phase Flow. Eyrolles, Paris (1975)

[6] Stewart, H.B., Wendroff, B.: Two-phase flow: models and methods. J. Comput. Phys. **56**, 3 (1984)

[7] Keyfitz, B.L., Sanders, R., Sever, M.: Lack of Hyperbolicity in the two-fluid model for two-phase incompressible flow. Discret. Contin. Dyn. Syst. Ser. B **3**(4), 541–563 (2003)

[8] Ndjinga, M., Kumbaro, A., De Vuyst, F., Laurent-Gengoux, P.: Numerical simulation of hyperbolic two-phase flow models using a Roe-type solver. Nucl. Eng. Des. **238**(8), 2075–2083 (2008)

[9] Stuhmiller, J.H.: The influence of interfacial pressure forces on the character of two-phase flow model equations. Int. J. Multiph. Flow **3**, 551–560 (1977)

[10] Stewart, H.B.: Stability of two-phase flow calculation using two-fluid models. J. Comput. Phys. **33**(2), 259–270 (1979)

[11] Park, J.-W., Drew, D.A., Lahey Jr., R.T.: The analysis of void wave propagation in adiabatic monodispersed bubbly two-phase flows using an ensemble-averaged two-fluid model. Int. J. Multiphase Flow **24**, 1205–1244 (1999)

[12] Batchelor, G.K.: An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge (1967)
[13] de Bertodano, M.L., Fullmer, W., Clausse, A., Ransom, V.: Two-Fluid Model Stability, Simulation and Chaos. Springer (2017). doi:10.1007/978-3-319-44968-5

[14] Fullmer, W.D., de Bertodano, M.L., Zhang, X.: Verification of a higher-order finite difference scheme for the one-dimensional two-fluid model. J. Comput. Multiphase Flows 5(2), 139–155 (2013)

[15] Keyfitz, B.L., Sever, M., Zhang, F.: Viscous singular shock structure for a non-hyperbolic two-fluid model. Nonlinearity 17(5), 1731 (2004)

[16] Ndjinga, M.: Quelques espects de modélisation et d'analyse des systèmes issus des 'ecoulements diphasiques. PhD Thesis (2007)

[17] Tan, D., Chang, T., Zheng, Y.: Delta-shock waves as limits of vanishing viscosity for hyperbolic systems of conservation laws. J. Differ. Equ. 112, 1–32 (1994)

[18] Danilov, V.G., Mitrovic, D.: Delta shock wave formation in the case of triangular system of conservation laws. J. Differ. Equ. 245, 3707–3734 (2008)

[19] Karlsen, K.H., Mishra, S., Risebro, N.H.: Convergence of finite volume schemes for triangular systems of conservation laws. Numer. Math. 111, 559–589 (2009)

[20] Agrell, C., Risebro, N.H.: Convergence of a relaxation scheme for a 2 × 2 triangular system of conservation laws. Int. J. Numer. Anal. Model. 11, 148–171 (2014)

[21] Tan, D.C., Zhang, T.: Two-dimensional riemann problem for a hyperbolic system of nonlinear conservation laws: II. Initial data involving some rarefaction waves. J. Differ. Equ. 111(2), 255–282 (1994)

[22] Selkovich, V.M.: The Riemann problem admitting $\delta$-, $\delta'$-shocks and vacuum states (the vanishing viscosity approach). J. Differ. Equ. 231, 459–500 (2006)

[23] Godlewski, E., Raviart, P.-A.: Numerical Approximation of Hyperbolic Systems of Conservation Laws. Applied Mathematical Sciences, vol. 118. Springer, New York (1996)

[24] Dafermos, C.M.: The entropy rate admissibility criterion for solutions of hyperbolic conservation laws. J. Differ. Equ. 14, 202–212 (1973)

[25] Lax, P.D.: Shock waves and entropy. In: Zarantonello, E.H. (ed.) Contributions to Nonlinear Functional Analysis, pp. 603–634. Academic Press, New York (1971)

[26] Dafermos, C.M.: Hyperbolic Conservation Laws in Continuum Physics, 3rd edn. Springer, Berlin (2010)

[27] LeFloch, P.G.: Hyperbolic Systems of Conservation Laws: The Theory of Classical and Nonclassical Shock Waves. Birkhauser, Basel (2002)

[28] Lagouti'ere, F.: A non-dissipative entropic scheme for convex scalar equations via discontinuous cell-reconstruction. C. R. Acad. Sci. Paris Ser. I 338, 549–554 (2004)
[29] Aguillon, N., Chalons, C.: Nondiffusive conservative schemes based on approximate Riemann solvers for Lagrangian gas dynamics (submitted)

[30] Boutin, B., Chalons, C., Lagoutire, F., LeFloch, P.G.: Convergent and conservative schemes for non-classical solutions based on kinetic relations. Interfaces Free Bound. **10**, 399–421 (2008)

[31] Chalons, C., Delle Monache, M.-L., Goatin, P.: A numerical scheme for moving bottlenecks in traffic flow. In: Proceedings of HYP2014 International Conference (submitted)

[32] Aguillon, N.: Problèmes d’interfaces et couplages singuliers dans les systèmes hyperboliques: analyse et analyse numérique. Ph.D. Thesis (2015)

[33] Harten, A.: High resolution schemes for hyperbolic conservation laws. J. Comput. Phys. **49**, 357–393 (1983)

[34] Hayes, B.T., LeFloch, P.G.: Nonclassical shocks and kinetic relations: finite difference schemes. SIAM J. Numer. Anal. **35**, 2169–2194 (1998)

[35] Kissling, F., Rohde, C.: The computation of nonclassical shock waves with a heterogeneous multiscale method. Netw. Heterog. Media **5**, 661–674 (2010)

[36] LeFloch, P.G., Mercier, J.M., Rohde, C.: Fully discrete, entropy conservative schemes of arbitrary order. SIAM J. Numer. Anal. **40**, 1968–1992 (2002)

[37] LeFloch, P.G., Rohde, C.: High-order schemes, entropy inequalities, and non-classical shocks. SIAM J. Numer. Anal. **37**, 2023–2060 (2000)

[38] Merkle, C., Rohde, C.: Computation of dynamical phase transitions in solids. Appl. Numer. Math. **56**, 1450–1463 (2006)

[39] Chalons, C., Engel, P., Rohde, C.: A conservative and convergent scheme for undercompressive shock waves. SIAM J. Numer. Anal. **52**(1), 554–579 (2014)

[40] Desprès, B., Lagoutière, F.: Contact discontinuity capturing schemes for linear advection and compressible gas dynamics. J. Sci. Comput. **16**(4), 479–524 (2001)

[41] Lagoutière, F.: Stability of reconstruction schemes for scalar hyperbolic conservation laws. Commun. Math. Sci. **6**(1), 57–70 (2008)

[42] Aguillon, N., Chalons, C.: Non diffusive conservative schemes based on approximate Riemann solvers for Lagrangian gas dynamics. Math. Model. Numer. Anal. (M2AN) **50**(6), 1887–1916 (2016)

[43] Coquel, F., El Amine, K., Godlewski, E., Perthame, B., Rascle, P.: A numerical method using upwind schemes for the resolution of two-phase flows. J. Comput. Phys. **136**(2), 272–288 (1997)

[44] Munkejord, S.T., Evje, S., Flatten, T.: A MUSTA scheme for a nonconservative two-fluid model. J. Sci. Comput. **31**(4), 2587–2622 (2009)

[45] Shekari, Y., Hajidavalloo, E.: Application of Osher and PRICE-C schemes to solve compressible isothermal two-fluid models of two-phase flow. J. Comput. Fluids **86**, 363–379 (2013)
M. Ndjinga and T. P. K. Nguyen

Den-Service de Thermo-Hydraulique et de M'ecanique des Fluides (STMF)

CEA

Universit'e Paris-Saclay

91191 Gif-sur-Yvette

France

e-mail: tpkieu.nguyen@gmail.com

M. Ndjinga

e-mail: michael.ndjinga@cea.fr

T. P. K. Nguyen and C. Chalons

LMV, UMR 8100

UVSQ

45 Avenue des États-Unis

78035 Versailles Cedex

France

e-mail: christophe.chalons@uvsq.fr

Received: 19 April 2016.

Accepted: 29 May 2017.