Exact and numerical solutions of the Riemann problem for a conservative model of compressible two-phase flows

In this work we study the solution of the Riemann problem for the barotropic version of the conservative symmetric hyperbolic and thermodynamically compatible (SHTC) two-phase flow model introduced in \cite{Romenski2007,Romenski2009}. All characteristic fields are carefully studied and explicit expressions are derived for the Riemann invariants and the Rankine-Hugoniot conditions. Due to the presence of multiple characteristics in the system under consideration, non-standard wave phenomena can occur. Therefore we briefly review admissibility conditions for discontinuities and then discuss possible wave interactions. In particular we will show that overlapping rarefaction waves are possible and moreover we may have shocks that lie inside a rarefaction wave. In contrast to nonconservative two phase flow models, such as the Baer-Nunziato system, we can use the advantage of the conservative form of the model under consideration. Furthermore, we show the relation between the considered conservative SHTC system and the corresponding barotropic version of the nonconservative Baer-Nunziato model. Additionally, we derive the reduced four equation Kapila system for the case of instantaneous relaxation, which is the common limit system of both, the conservative SHTC model and the non-conservative Baer-Nunziato model. Finally, we compare exact solutions of the Riemann problem with numerical results obtained for the conservative two-phase flow model under consideration, for the non-conservative Baer-Nunziato system and for the Kapila limit. The examples underline the previous analysis of the different wave phenomena, as well as differences and similarities of the three systems.


Introduction
The aim of the paper is to construct exact solutions for the Riemann problem of the one-dimensional barotropic version of the conservative Symmetric Hyperbolic and Thermodynamically Compatible (SHTC) model of compressible two-phase flows introduced in [63,58]. The results obtained can be useful for a qualitative understanding of the physical processes occurring in two-phase flows, for a comparative analysis of various models, and also for testing numerical methods.
Note that the existence of a well-developed theory of exact solutions of the Riemann problem ensured the success in the development of modern shock-capturing numerical methods [28] for solving the Euler equations of compressible gasdynamics, which are of fundamental importance in science and engineering, in particular aerospace engineering and astrophysics. See [74] for an exhaustive overview of shock capturing schemes based on the exact or approximate solution of the Riemann problem. And even at present, the known exact solutions of the Euler equations are successfully used to test new numerical methods for solving hyperbolic systems of equations. Therefore, in our opinion, the construction of exact solutions for the equations of two-phase flows will influence the formation of a common point of view in the field of multiphase flow modeling and the development of new numerical methods for solving problems related to this area.
In contrast to single-phase gas dynamics, there is still no universally accepted model of compressible multiphase flows, even for the case of only two phases, see, for example [68]. The generally accepted approach to design a two-phase flow model is based on the assumption that a mixture is a system of two interacting single phase continua, see, for example, [46]. The most widely used PDE system is the Baer-Nunziato model [2], various modifications of which have been applied by many authors to study different types of flows, including flows with phase transitions and chemical reactions, see, for example recent papers [11,26,69] and references therein.
In this paper, we will study a model developed on the basis of the theory of SHTC systems, which was developed in [31,29,30,56,66,53]. This theory connects the local well-posedness of the governing partial differential equations of continuum physics (symmetric hyperbolicity in the sense of Friedrichs) with the fulfillment of the laws of thermodynamics (the law of conservation of energy and the law of increasing entropy). The general master system of SHTC equations can be derived from an underlying variational principle [53], and the study of numerous models of continuum mechanics has shown that their governing differential equations belong to the SHTC class of PDEs. This theory can be used to formulate new, well-posed models of processes in complex media, including unified models for the description of viscous Newtonian and non-Newtonian fluids and nonlinear elasto-plastic solids at large deformations [54,17,18,52], its extension to general relativity [60], SHTC models of rupture dynamics [72,27] and flows in deformable porous media [62,61].
The object of study in this paper is the SHTC model of compressible two-phase flow, the general equations of which are formulated and discussed in [64,63,58,57]. In the model under consideration, the two-phase medium is assumed to be a single continuum, the properties of which take into account the features of the two-phase flow. This means that the element of the medium is described by the average field of phase velocities and pressures for the mixture, but the flow of phases through this element is allowed. We consider the simplest model of a one-dimensional flow of a barotropic mixture. In this case, the SHTC equations can be written in a fully conservative form and therefore allow a direct formulation of discontinuous solutions. Nevertheless, the construction of exact solutions for the Riemann problem turns out to be a rather difficult task.
It is necessary to note that the governing equations of the Baer-Nunziato model [2,1] in general can not be completely written in a conservative form, even in the onedimensional case. This creates difficulties in the definition of discontinuous solutions of the shock-wave type. But it turns out that the barotropic Baer-Nunziato system can be presented in a conservative form for one particular choice of interfacial velocity and interfacial pressure. And this conservative system is exactly the same as the barotropic case of the SHTC model studied in this paper. We discuss the similarities and differences between the different models by rewriting the equations of the SHTC model in the form of the Baer-Nunziato model for two interacting single phase continua. We also consider the reduced one-dimensional limit system obtained in the stiff pressure and velocity relaxation limit i.e. the resulting single-velocity single-pressure approximation of the SHTC and Baer-Nunziato models for the barotropic flow and show that they are identical and reduce to the well-known Kapila model [37].
The conservative form of the SHTC model of compressible two-phase flows has advantages not only for the construction of exact solutions of the Riemann problem, but also when using advanced shock capturing numerical methods. For the numerical simulations shown in this paper we will therefore rely on classical second order high resolution shock capturing TVD finite volume schemes as presented in [74].
Degenerate behaviour in hyperbolic systems is a topic of constant interest. However, one has to be careful identifying the source of degeneracy. Non-strictly hyperbolic systems, i.e. systems with multiple eigenvalues, naturally arise in multi dimensions [15]. Specific cases of non-strictly hyperbolic systems where studied in [38,3,44] and recently by Freistühler and Pellhammer [24]. There mostly systems of two equations where studied and the coincidence of eigenvalues often occurred in points where also the character of the related field changed from genuine nonlinear to linearly degenerated. In particular Freistühler showed in [22] that coinciding eigenvalues in the presence of a discontinuity must belong linearly degenerated fields. Considering the Euler equations or the present system such a situation is related to the vanishing of the fundamental derivative G. This may happen for certain equations of state, but in general does not necessarily imply multiple eigenvalues. A prominent example is the system of Euler equations and we highly recommend [75,45,49]. Another reason for degeneracies may be the loss of an eigenvector as described in [70,13]. Systems with missing eigenvectors are often called weakly hyperbolic and systems with coinciding eigenvalues are sometimes also called hyperbolic degenerate. For further reading and a more detailed survey we recommend [10]. Coinciding eigenvalues may lead to difficult situations, but as long as there is a full set of eigenvectors which span the complete space the system is still diagonisable. For the construction of complete Riemann solvers a full set of eigenvectors is of great importance, see e.g. [55,50,19]. However, we want to emphasize that the consequences depend crucially on the system under consideration. In some cases the numerical methods break down when applied to weakly hyperbolic systems, see e.g. [12], whereas in other situations a numerical treatment of weakly hyperbolic systems is still possible, see e.g. [33,32].
The rest of the paper is organized as follows: in Section 2 we present the SHTC system under consideration in this paper and study its eigenstructure. In Section 3 we briefly summarize degeneracies and admissibility conditions. The wave relations, i.e. the Riemann invariants and the Rankine-Hugoniot conditions of the model are presented in Section 4, while the possible wave configurations are shown in Section 5. The relation of the conservative SHTC system with the non-conservative Baer-Nunziato model and the common Kapila limit are shown in Section 6. Some examples of exact solutions and corresponding numerical results are presented in Section 7. The paper closes with some concluding remarks and an outlook to future work in Section 8.
2 Conservative barotropic SHTC model for compressible two-phase flows 2.1 Multi-dimensional case The PDE system for compressible barotropic two-phase flows was discussed in Romenski et al. [63,58]. Written in terms of the specific energy E = E(α 1 , c 1 , ρ, w k ) it reads Here, α 1 is the volume fraction of the first phase, which is connected with the volume fraction of the second phase α 2 by the saturation law α 1 + α 2 = 1, ρ is the mixture mass density, which is connected with the phase mass densities ρ 1 , ρ 2 by the relation ρ = α 1 ρ 1 + α 2 ρ 2 . The phase mass fractions are defined as c 1 = α 1 ρ 1 /ρ, c 2 = α 2 ρ 2 /ρ and it is easy to see that c 1 + c 2 = 1. The mixture velocity is given by u i = c 1 u i 1 + c 2 u i 2 and w i = u i 1 − u i 2 is the relative phase velocity. The equations describe the balance law for the volume fraction, the balance law for the mass fraction, the conservation of total mass, the total momentum balance law and the balance for the relative velocity. The phase interaction is present via algebraic source terms in (2.1a) and (2.1e), which are proportional to thermodynamic forces. These source terms are phase pressure relaxation to the common value and interfacial friction The coefficients θ 1 , θ 2 characterize the rate of pressure and velocity relaxation and can depend on parameters of state. Due to mass and momentum conservation throughout this work we assume ξ 2 = ξ 3 = ξ 4 = 0.

Discussion of the mixture equation of state
We now want to further specify the derivatives of the generalized energy E. Due to the relation for the mass fractions and the concentrations we write α ≡ α 1 , α 2 = 1 − α and c ≡ c 1 , c 2 = 1 − c, when appropriate. The mixture equation of state (EOS) is defined as the sum of the mass averaged phase equations of state and the kinematic energy of relative motion where e i (ρ i ) is the specific internal energy of the i-th phase and is assumed to be known. Special attention has to be paid when the derivatives of the internal energies are needed. However, if the temperature is held constant the derivative is given by Note that the relations for the barotropic case were not given in previous work and in particular the isothermal case is no straightforward simplification of the general case. The detailed calculations are given in the Appendix A and we will just summarize the results needed here. We now introduce the mixture pressure and the specific Gibbs energy of phase i The Gibbs energy, or sometimes free enthalpy, may be generalized to the chemical potential µ when more substances are involved, see [40]. For the mixture EOS we have the following derivatives The derivatives for the internal energy e(α, c, ρ) = c 1 e 1 (ρ 1 ) + c 2 e 2 (ρ 2 ) of the mixture are for the isentropic case and 13) for the isothermal case.

One-dimensional model
In this paper we focus on the one dimensional case and thus the equations simplify to Note, that the curl term in equation (2.1e) vanishes. Applying the obtained results and relations and introducing the notation the system can be rewritten in the following form: For the isentropic case the total energy inequality, which serves as mathematical entropy inequality, reads while for the isothermal case we have the inequality Inequality (2.17) is the mathematical formulation of the physical statement that the free energy of a system under consideration is minimized in an isothermal process. Moreover, this formulation is a straightforward generalization of the inequality obtained for the isothermal Euler equations, see [15,67,73]. We want to note that it is beneficial to derive the isothermal model from the more general model including the thermal impulse. More detailed information are given in the Appendix.

Conservative formulation
The system (2.15a) -(2.15e) can be written in the conservative form given by where the vector of conserved quantities reads Using the equations obtained so far we can write the conservative flux as follows (2.18)

Primitive formulation
We want to reformulate the barotropic system (2.15) in terms of the primitive variables α 1 , ρ 1 , ρ 2 , u 1 and u 2 . The other quantities are then obtained using the relations Here we assume α 1 ∈ (0, 1) and ρ 1 , ρ 2 > 0, i.e. we exclude vacuum states. We further introduce the speed of sound a i of phase i given by the following relation For a thermodynamically consistent equation of state the speed of sound is well defined and we do not have to differ between the two cases for the mathematical considerations. With W ≡ (α 1 , ρ 1 , ρ 2 , u 1 , u 2 ) the Jacobian of this system is given by 2.20) and the source terms can be transformed using the following matrix Now we can write the system in the following compact form The eigenvalues can be computed as and we have (up to scaling) the following right eigenvectors (2.23) Here we introduced the following abbreviations The subscripts 1, 2 refer to the corresponding phases and the pair (λ C , R C ) takes a special role as we will see in a moment. We further want to investigate the fields and see whether they are genuine nonlinear or linearly degenerated. The gradients of the eigenvalues with respect to the given variables W = (α 1 , ρ 1 , ρ 2 , u 1 , u 2 ) are given by We immediately obtain Here we have introduced the fundamental derivative G and the following relation holds For more details on G and its crucial influence on the flow we recommend Menikoff and Plohr [45] and Müller and Voss [49]. Throughout this work we assume G > 0 and thus the fields 1± and 2± are genuine nonlinear. Indeed for an ideal gas we have We further have the remarkable property that Thus the eigenvectors R 1± and R 2± span a four dimensional hyperplane in the state space where α is constant. It remains to discuss the field C. Therefore we use We obtain Thus this field is linearly degenerated and hence discontinuities associated to this field are contact waves. In view of the above results it is clear that each character of the present fields is independent of the flow. Note that up to now we cannot exclude situations where eigenvalues coincide, i.e. have a multiplicity larger than one. Different possible phenomena related to these special situations will be discussed in Section 5.

Degeneracies of the system and admissibility conditions
By construction, the system under consideration is symmetric hyperbolic in the sense of Friedrichs [25]. The initial value problem (Cauchy problem) for such a system is wellposed locally in time [15,4]. But the question of the solvability of the Cauchy problem in the large for symmetric hyperbolic systems is still a challenging problem and, for example, nontrivial interesting phenomena (such as the resonance effect) can arise. As mentioned in the previous Section we now will have to study under which conditions this system is hyperbolic in the sense that we have real eigenvalues and a full set of eigenvectors, or, more precisely, under which conditions certain degeneracies may occur. Furthermore another crucial point is, as mentioned before, that up to now the order of the waves is not clear, i.e. the order of the eigenvalues. Thus we will briefly review results concerning the admissibility conditions for discontinuities which will play an important role throughout this work.

Coincidence of Eigenvalues, Parabolic Degeneracy and Hyperbolic Resonance
We now want to investigate situations where the eigenvectors may become linearly dependent. Considering the eigenvectors (2.23) it is obvious that the eigenvectors R i± become linearly dependent in each phase iff a i = 0. Since we assume a strictly positive speed of sound in each phase this is excluded and the four eigenvectors R i± remain linearly independent. For the field (λ C , R C ) we can state for a fixed field (λ i± , R i± ) that This can be seen by direct calculations which are given in the Appendix A for convenience. We want to discuss the consequences and interpretation of this case later on. We note that in the extreme situation that three eigenvalues coincide, which is the maximum for reasonable EOS, the eigenvector R C becomes the null vector. This situation can only happen when λ C coincides with one eigenvalue of each phase and thus we have

Admissibility Conditions for Discontinuities Revisited
In the previous Section 2 we presented the eigenvalues (2.22) and eigenvectors (2.23) of the system under consideration (2.15). As noted before different situations may occur due to coinciding eigenvalues. Thus it is important to review suitable criteria to single out admissible solutions. In particular we are interested in the case when discontinuities occur. For a given and fixed state W we have for the eigenvalues with p ≤ n and n being the total number of the possible eigenvalues. In our case we have n = 5 and further With this we allow for situations where the order of the eigenvalues changes and that eigenvalues may coincide. We start with the original work by Lax [41] and assume p = n.
In particular we consider a strictly hyperbolic system. Further we consider the states left and right of a discontinuity denoted by W L and W R . According to Lax a k-shock with the speed S satisfies the condition From these inequalities we deduce that n − (k − 1) characteristics impinge from the left of the discontinuity and k from the right. Following the presentation given in Dafermos [15] the situation can be generalized as follows. We consider the eigenvalues for the left and right state with the agreement of λ (0) (W L ) = −∞ and λ (n+1) (W R ) = ∞. If the inequalities are satisfied with i = j we obtain the Lax condition (3.3) given above and the shock is called compressive. In the case of i < j the shock is called overcompressive and for i > j the shock is called undercompressive. Now we want to deal with the non strict situation, i.e. there is no well-defined ordering of the eigenvalues and they may coincide. In Keyfitz and Kranzer [38] a generalization to non-strictly hyperbolic systems is given as follows 1) n + 1 characteristics enter the shock and n − 1 leave it or 2) n − 1 characteristics enter and leave the shock whereas the remaining two are tangent to the shock and belong to linearly degenerated fields.
When an i-shock has n − 1 leaving characteristics and the corresponding eigenvectors fulfill then the shock is called evolutionary, see again [15] or the book of Kulikovskii et al. [39]. We want to end this brief review of admissibility conditions with the results presented in [1]. There the results obtained in the references given above are basically collected and generalized to incorporate different situations that are of interest for non-strictly hyperbolic systems. Let us consider the situation (3.2). Having p eigenvalues we conclude that we have p unknowns on each side of the discontinuity and additionally the shock speed S, i.e. N = 2p + 1 unknowns. These N unknowns can be determined as follows.
Let us assume we have m relations across the discontinuity, e.g. the jump conditions. Second, the unknowns should be determined by the flow using the characteristics. Let the number of incoming characteristics be i, the number of outgoing characteristics o and the number of coinciding characteristics be c. The incoming and coinciding characteristics are determined by the past and thus provide further information. Hence in order to determine the unknowns we demand In [1] this is called evolutionarity condition and a discontinuity is called evolutionary if the condition (3.5) holds (in agreement with the results above). This implies that a discontinuity is evolutionary iff o = m − 1, see [1]. A closely related concept was introduced by Freistühler in [23]. The idea is quite similar as again a linearised problem is investigated. The eigenvalues are distinguished, using there relation to the speed of the discontinuity, into slow and fast characteristics. In this sense outgoing characteristics are slow characteristics on the left side of the discontinuity and fast ones on the right side. For the incoming characteristics the situation is reversed. Coinciding characteristics move at the characteristic speed of the discontinuity.

Wave relations
We now want to obtain the relations that are valid across the different types of waves. In particular we derive the Riemann invariants which are constant across rarefaction waves and the Rankine Hugoniot jump conditions across discontinuities.

Rarefaction waves
Given a smooth solution we may apply a nonlinear transformation of the variables aiming to simplify the system using another appropriate choice of variables. Such a special set of variables is given by the Riemann invariants, cf. [15,21,71]. The Riemann invariants for R 1± and R 2± can be calculated using Here µ, ν ∈ {1, 2}, µ = ν where µ = 1 for R 1± and µ = 2 for R 2± . The relations for ρ µ and u µ can be combined to Thus we obtain the following invariants Furthermore, the slope inside a left rarefaction wave is given by and hence we obtain that the solution inside the rarefaction fan is given by Here, ρ is obtained as the root of F (ρ). Similar we obtain the results for right rarefaction waves corresponding to λ µ+ Although the Riemann invariants (4.1) state that ρ ν and u ν remain constant, we prefer the notion that these quantities are not affected by the rarefaction wave. As we see later on they might change due to other waves.

Shock waves
In the presence of discontinuities, such as shock waves, the situation is different. However, since the present system is conservative, corresponding Rankine Hugoniot jump conditions have to hold at a shock with speed S, see for example [15,42]. In the present case we obtain the following conditions across discontinuities The jump conditions (4.5a)-(4.5c) can be reformulated to Using the third jump condition (4.8) we obtain for the first equation (4.6) We introduce the abbreviations for the mixture mass flux and the mass fluxes of the phases, respectively. In particular we have Thus we can derive a jump condition for the partial mass flux of the second phase, using the continuity of the mixture mass flux (4.8) and the continuity mass flux of the first phase (4.7), i.e.
Using the jump conditions for the partial mass fluxes (4.7) and (4.10) we obtain for the fourth jump condition (4.5d) (4.11) The fifth jump condition (4.5e) can be reformulated as follows Summarizing we have the following system of jump conditions Note that (4.13c) may be replaced by (4.10). Further, equations (4.13d) and (4.13e) may also be reformulated to To give the complete picture we briefly summarize the jump conditions according to the system (2.14) in terms of the mixture EOS We also want to emphasize the analogy to the jump conditions of the Euler equations, cf. [15,74]. Equation (4.13d) is basically a volume fraction weighted combination of an individual momentum jump condition for each phase as it appears in the Euler equations. Moreover, a single jump term in (4.13e) agrees with the jump bracket of the energy equation for the Euler equations.

Exploiting the jump conditions for a Lax shock
In the following we assume that we have a Lax-shock, i.e. with no tangential eigenvalues.
In particular this implies u = S with S being the shock speed. Indeed, as we will see this cannot happen for shock waves. From equation (4.13a) we thus have that α 1 is continuous across a shock corresponding to the eigenvalues λ 1± and λ 2± . Up to now we have made no further assumption on α 1 . With the aim to further exploit and simplify the jump conditions (4.13a) -(4.13e) we assume α 1 ∈ (0, 1) from now on. First we have, due to the continuity of α 1 (4.13a), that also the mass fluxes Q 1 and Q 2 are continuous, see (4.13b) and (4.10). Hence the partial mass fluxes may be written outside the jump brackets in (4.15), i.e.
Further we can use the continuity of the partial mass fluxes and write with µ ∈ {1, 2} Using (4.18) together with the momentum jump condition (4.11) we get Thus equations (4.17) and (4.19) form a linear system for the partial mass fluxes The inverse is given by Note that across a shock the phase densities necessarily jump and thus the determinant is not equal to zero, see Appendix A. Hence we obtain for the partial mass fluxes Furthermore, these partial mass fluxes are also related to each other using Thus one may express one partial mass flux by the other. Up to now we only have the squares of the mass flux. The correct sign of the corresponding root is given by the Lax condition for the present shock. The Lax criterion states for a shock in phase µ that where W L and W R are the left and right states adjacent to the shock. For a left shock in phase µ ∈ {1, 2} we have Similar we obtain for a right shock Thus the partial mass flux has a strict sign. Hence across a λ 1± -shock the sign of Q 1 is determined and for a λ 2± -shock the sign of Q 2 is given, respectively. Therefore we choose the following The values for the densities and w may assumed to be given, i.e. the values on one side of the shock. Once we have obtained the partial mass flux we can eliminate a further unknown using the velocity jump condition (4.18)

Entropy inequality
Given the (mathematical) entropy inequalities (2.16) and (2.17) we have in the presence of a discontinuity for the isentropic case and for the isothermal case, respectively. With algebraic manipulations using the continuity of the partial mass fluxes (4.7), (4.10) and the jump condition for the momentum (4.13d) we obtain for both cases Note that the bracket terms correspond to the terms in the jump condition for the relative velocity (4.13e). Thus we can replace the term corresponding to one phase by the other if (4.13e) holds and obtain In particular the inequality for the isothermal case is perfectly analogous to the one used in [34].

Contact wave
The wave corresponding to the linearly degenerated field (λ C , R C ) is a contact wave and because of it is immediately clear that we have the Riemann invariant The characteristic condition for the contact wave is λ C (W L ) = S = λ C (W R ) which further gives for the velocity of the contact S = u. Thus we have Q = −ρ(u − S) = 0 and hence α 1 may be non-zero. Indeed, since α 1 is continuous across the other waves the jump of α 1 is given by the initial data. From the continuity of the partial mass fluxes we obtain For the mixture momentum at the contact we obtain using (4.14) and (4.31) From the jump condition for the relative velocity (4.13e) we get Thus we have the following equations at the contact Note that ρc 1 c 2 w 2 is a dynamic mixture pressure related to the relative velocity. By introducing the generalized total mixture pressurep = ρc 1 c 2 w 2 + p as the sum of the dynamic pressure and the static mixture pressure p we can rewrite (4.34c) as We again give the jump conditions in terms of the mixture EOS for completeness: Finally, it is easy to see from (4.29) and u = S that the entropy inequality is fulfilled with the right side being identically zero.

Wave configurations and relations
In this section we want to study particular wave configurations. The crucial point is, as mentioned before, that up to now the order of the waves is not clear, i.e. the order of the eigenvalues. In particular eigenvalues may coincide in certain points. Thus the reviewed results concerning the admissibility conditions for discontinuities will play a crucial role throughout this section.

Contact
We first want to discuss phenomena related to the linearly degenerates field (λ C , R C ), i.e. the contact wave. Obviously we can determine the position of the contact wave due to the fact that λ C = u is the convex combination of the individual phase velocities, i.e.
Let us now consider the situation that we have an isolated shock (w.l.o.g.) corresponding to λ µ− moving with speed S. A priori it is not obvious at all whether we may encounter the situation that u = S. In the following we will show that this is not possible. In this situation we have a discontinuity with the tangential eigenvalues λ − C = λ + C = S. Here a superscript − refers to the state W − left of the discontinuity and a superscript + refers to the state W + right of the discontinuity, respectively. According to the results obtained by Keyfitz et al. [38] cited above this corresponds to the second situation. Thus we would further need n − 1 incoming and n − 1 outgoing characteristics, with n = 5 in our case. If we consider the results given in [1] we obtain the same results with m = 5 jump conditions, c = 2 tangential eigenvalues, N = 2n + 1 = 2 · 5 + 1 = 11 unknowns and thus i = 4 incoming and o = m − 1 = 4 outgoing characteristics. Hence it is clear that for an evolutionary discontinuity we must not have further coinciding eigenvalues. Let us picture the situation more precisely and assume w.l.o.g. that we have a shock corresponding to λ 1− . Thus we would have We can conclude immediately that λ − 1+ > S and S > λ − 2− . Using the continuity of α 1 Q 1 it follows that λ + 1+ > S and hence S > λ + 2− . Since we now already have four ingoing characteristics it follows that λ − 2+ and λ + 2+ must be outgoing characteristics. The situation is depicted as an example in Figure 1. So clearly this situation is excluded, since the Figure 1: Shock with u = S (example).
information related to the eigenvalue λ 2+ is coming out of the discontinuity. Another short argument would be to say that we are in a strictly hyperbolic situation on both sides of the discontinuity and hence only a classical wave is allowed. Note that this is not by far that obvious for systems with multiple linearly degenerated fields. Thus we can state that in our system an evolutionary discontinuity with u = S is a contact and vice versa. Even more important is the statement that due to u = S for a shock we always have the continuity of α across the shock. This also gives another view on the results obtained above in Section 3 and in particular the situation of coinciding eigenvectors described by (3.1). In this particular situation α will not even jump across R C and thus remains constant in the complete fan. Hence we can exclude this situation by simply prescribing different values for α initially. Or in other words this situation may only occur when α L = α R holds for the initial states of the Riemann problem. We therefore could interpret the case (3.1) as a consequence of a redundant α equation. With α constant everywhere we can reformulate the system as The situation λ i± = λ C then corresponds to the case of hyperbolic resonance discussed by Isaacson and Temple [36]. In the literature you also find the phrasing parabolic degeneracy or weak hyperbolicity for missing eigenvectors, but we think resonance is the term best suited here.

Contact inside Rarefaction
Let us assume that a contact lies inside a rarefaction wave, see e.g. the sketch shown in Figure 2. It cannot be attached to one side (or both) of the rarefaction, because then we would have the situation of hyperbolic resonance discussed above (i.e. coincidence of eigenvalues and eigenvectors). Thus the contact will tear the rarefaction wave into two parts. Assume we have a rarefaction wave corresponding to λ µ− and thus we consider the following situation According to the results obtained above concerning the admissibility of discontinuities we have c = 2 coinciding characteristics and thus need i = o = 4 characteristics going in and out. From the given relation for the eigenvalues λ µ− and λ C we can directly conclude λ + µ+ > S. We can then conclude that that λ + ν− < S. Using the inequalities for λ µ− we see that Q µ < 0 across the contact. Since 0 = Q = α µ Q µ + α ν Q ν we thus have Q ν > 0. Using Q µ < 0 we yield λ − µ+ > S. From Q ν > 0 we yield λ − ν− < S and hence finally λ + ν+ < S. Altogether we therefore obtain Such a discontinuity violates the admissibility criteria given above due to the fact, that the characteristics λ − µ− and λ + µ− do not contribute any information to the discontinuity. Note that as mentioned before the phrasing outgoing characteristic seems to be not quite suited here for these two eigenvalues. However, the results remain valid and one could use the notation of slow and fast characteristics with respect to the side of the discontinuity, cf. [23]. In this sense the situation is analogue to that of a contact coinciding with a shock.

Overlapping Rarefaction Waves
A possible wave configuration might be two rarefaction waves that overlap, see Figure 3. In this situation one rarefaction wave belongs to phase one and the other to phase two, Figure 2: Contact inside a rarefaction (example). Figure 3: Overlapping rarefaction waves (example).
respectively. The rarefaction fans are cones given by In the case of an empty intersection we have the classical situation and the solution is obtained using the corresponding eigenvector. If however, C * = ∅ we have Due to the special structure of the eigenvectors the obtained invariants remain unchanged and the formulas (4.3) and (4.4) stay valid. Indeed if we write down the (homogeneous) system using the primitive variables W = (α 1 , ρ 1 , ρ 2 , u 1 , u 2 ) we have the Jacobian (2.20).
For a rarefaction wave α 1 remains constant and thus the system simplifies to Hence the system decouples into two barotropic Euler systems for each phase and the rarefaction waves can be obtained individually. The solution is then obtained as the superposition of the individual solutions.

Shock interacting with a Rarefaction Wave
Another situation that can occur is a shock which lies inside a rarefaction fan, see Figure  4. Let us assume for the moment that we have a left shock corresponding to phase Figure 4: Shock inside a rarefaction wave -Case (i).
µ ∈ {1, 2} and a left rarefaction corresponding to phase ν ∈ {1, 2}, µ = ν. Clearly, the phase µ is only affected by the shock wave due to the structure of the eigenvector corresponding to the rarefaction wave. For the phase ν the situation is more complicated.
We have for the eigenvalue λ ν− = ξ := x/t corresponding to the rarefaction wave where W L and W R denote the states left and right of the rarefaction wave. Let us denote the shock speed with S and quantities left of the shock are denoted with a superscript − and a + when they are on the right, respectively. There are four cases which are possible in this situation Case (i): The first case can be excluded since it implies linear degeneracy of the field (λ ν− , R ν− ), see [38,22]. This is obviously not the case as long as we have G ν = 0, which we may assume for our EOS. Further discussion of the fundamental derivative can be found in [45,49].

Case (ii):
Considering the characteristics in the second case we obviously have Thus we directly obtain Due to the continuity of the mass fluxes we further yield Since u is a convex combination of the phase velocities we also conclude λ − C > S and λ + C > S. Summarizing we have the following situation Therefore this situation is not admissible. According to the admissibility criteria we have N = 11 unknowns m = 5 equations and thus we would need i = 6 incoming and o = 4 outgoing characteristics. Case (iii): Now, for the third case we have Figure 5: Shock inside a rarefaction wave -Case (ii).
We have for the mass flux of phase µ that Q µ < 0. Thus we have u + µ > S and hence λ + µ+ > S. For the mass flux of phase ν we obtain Thus we also have Q < 0 and hence We therefore obtain that λ − C is an ingoing characteristic and λ + C an outgoing characteristic. In the situation under consideration we have as before N = 11 unknowns m = 5 equations. Additionally we have c = 1 coinciding wave and thus we need i = 5 incoming and o = 4 outgoing characteristics. This implies that λ − ν+ is also an ingoing characteristic. Summing up we have the following situation Case (iv): Finally, for the fourth case we have Using similar arguments as for the third case we obtain the following situation Figure 6: Shock inside a rarefaction wave -Case (iii). Figure 7: Shock inside a rarefaction wave -Case (iv).
Now we have two possible configuration which seem to be allowed. A priori it is not obvious whether one of these cases can be ruled out or if both may occur. Thus we will make use of the energy inequality (4.28) to investigate both cases.
Energy Inequality Case (iii): As in (4.29) we use the entropy inequality for the phase ν which defines the rarefaction, i.e.
From the characteristics λ − C , λ + C we directly conclude Q < 0 and hence we expect for this configuration in order to be admissible. We rewrite the kinetic energy in terms of the mass flux and assume the right state to be given. Thus it is possible to write the jump bracket as a function in the density ρ − ν of the left side of the discontinuity, i.e.
Investigating this function gives the following In order to get further insight we need to discuss the function g(ρ) = ρa ν (ρ) which is sometimes called Lagrangian wave speed [45]. More precisely, we already discussed it investigating the characteristic fields and we have for phase ν Thus we can conclude for all 0 Due to this monotonicity behaviour, ρ − ν < ρ + ν and since f (ρ + ν ) = 0 we can conclude Hence configuration (iii) respects the mathematical entropy inequality.
Energy Inequality Case (iv): We use a similar argumentation as before. Now we assume the left state to be given and write the jump bracket as a function in the density ρ + ν of the right side of the discontinuity, i.e.
Investigating this function gives the following Again we use the monotonicity of the Lagrangian wave speed (5.4) to conclude for all Due to this monotonicity behaviour, ρ − ν < ρ + ν and since f (ρ − ν ) = 0 we can conclude Hence configuration (iv) violates the energy inequality and finally we can say that only case (iii) is possible. The treatment for the case considering λ µ+ and λ ν+ is completely analogous showing that the shock and the rarefaction characteristic coincide now on the left side.
If we reconsider the arguments above it becomes clear that the side on which the rarefaction characteristic coincides with the discontinuity is basically defined by the mixture mass flux Q. Thus we can extend the above results to the cases of a λ µ+ shock interacting with a λ ν− rarefaction or the analogous case with λ µ− and λ ν+ . For Q > 0 in the first case we exemplary have the situation and for Q < 0 The other cases can be discussed as before. In particular the case Q = 0 is not admissible since it would contradict the genuine non-linearity of the field λ ν .

Shock Resonance
Due to the already mentioned result by Freistühler [22] and the results above we can exclude multiple eigenvalues near a shock discontinuity. However, we have to discuss the situation that we have a shock in each phase and both move at the same speed S. Hence we first consider the situation We directly verify u − µ > S, u − ν > S and hence u − > S. Due to the continuity of Q we also have u + > S. Furthermore we have λ − µ+ > S and λ − ν+ > S and this gives in total seven incoming characteristics, i.e.
There is no coinciding characteristic and we have three outgoing characteristics Thus such a configuration violates the admissibility conditions. Let us close with the case We directly verify u − µ > S, S > u + ν and hence we have From the continuity of the partial mass fluxes we conclude u + µ > S and S > u − ν and hence we already have the following six incoming and two outgoing characteristics Thus the situation of λ − C = λ + C = S can be excluded. But due to the continuity of the mass flux Q this implies that for Q < 0 or Q > 0 we have seven incoming and three outgoing characteristics and thus this situation is also not admissible. With this we have discussed most wave patterns that can appear in a Riemann problem for the model under consideration, including the interaction of two characteristic families.

Related Models
It is further interesting to discuss how other well established models are related to the studied system at hand. In particular we discuss Kapila's limit of the barotropic SHTC model and the relation to the Baer-Nunziato model in the following.

Instantaneous Relaxation Limit of the Barotropic SHTC Model
As it is noted above there are two relaxation terms in the system (2.1). One for the pressure relaxation (2.2) and one for the velocity relaxation (2.3), respectively. Since in real processes these relaxation processes can be quite fast, it is useful to obtain asymptotic limits of the solution for small values of the relaxation times. In the paper [47] the simplified equations for instantaneous pressure and relative velocity relaxations are accurately derived by the asymptotic analysis as a relaxation limit of the general Baer-Nunziato two-pressure two-velocity model for two-phase compressible fluid flows. In this section, based on the results obtained in [47] we derive a reduced SHTC system of PDEs for the instantaneous velocity and pressure relaxation. Note that it seems intuitive that the time scale of pressure relaxation is less than the time scale of velocity relaxation, because the physical mechanism of pressure relaxation is the pressure waves propagation. Thus, one can consider a possibility to study the instantaneous relaxation limit separately for the pressure and then for the velocity. If to consider only the pressure relaxation limit then we arrive at a single pressure two-velocity system of governing equations. It turns out that this PDE system is not hyperbolic. If we want to deal with hyperbolic equations we should consider instantaneous velocity relaxation together with the instantaneous pressure relaxation.
We will not repeat the rigorous asymptotic analysis as done in the [47], but will rely on the assumptions that follow from the rigorous theory. So, for the sake of simplicity we start with the instantaneous velocity relaxation. Let us consider equation (2.1e) An asymptotic expansion of the relative velocity w k = w k 0 + θ 2 w k 1 + ... for small θ 2 gives us For our purpose we only need to account for the zeroth order term of the expansion w k 0 = 0. This gives us a single velocity approximation of the model, whereas the second term of the expansion w k 1 = ∂E c 1 ∂x k gives us a phase diffusion Fick's law.
Now we substitute w k = 0 into the system (2.1) and remove the equation for the relative velocity. This gives a single velocity model for two-phase flows where u i is the single velocity of the flow.
Let us now consider the instantaneous pressure relaxation limit for (6.4) assuming θ 1 → 0. Note that θ 1 → 0 gives us p 1 = p 2 as a zeroth order approximation, but it is not correct to simply put p 1 = p 2 = P into the equations, because this would give us This equation for α 1 means that the volume fraction does not change along the trajectory, although the pressure can change. But if the phase pressures change, then α 1 must change, too, due to the different phase compressibility coefficients. The correct way to derive equations for the instantaneous pressure relaxation limit is to account for the following consequences of p 1 = p 2 = P : where K i = ρ i a 2 i is the phase bulk modulus, a i is the phase speed of sound. Now we use the above equation (6.6), the phase mass conservation equations (6.4b) and (6.4c) to derive the following equation for the volume fraction Thus, in case of instantaneous velocity and pressure relaxation we arrive at the following system The above equations are exactly the same as is in the Kapila model.

Comparison with the Baer-Nunziato Model
It is of further interest to compare the system (2.15) with the one dimensional barotropic Baer-Nunziato model given by the following equations It was already mentioned in [64] and [63] that for smooth solutions the systems can be reformulated into each other in one space dimension. However, we want to recall this transformation with a slightly different purpose. In the following we want to show the equivalence of these two systems (2.15) and (6.9) for smooth solutions and in particular that there is no freedom of choice for the interface quantities u I , p I and the sources in this case. This is of special interest since there are several possible choices for the interface quantities in the context of these Baer-Nunziato type models, see e.g. [1].
To show how the barotropic SHTC model (2.15) can be derived from the Baer-Nunziato model (6.9), we need to derive the balance for the relative velocity (2.15e). Therefore we want to reformulate the partial momentum balances (6.9d) and (6.9e) Subtracting equation (6.11) from equation (6.10) gives Using relation (2.19) and p I = (α 2 ρ 2 p 1 + α 1 ρ 1 p 2 )/ρ we obtain Thus we have shown that the SHTC model can be derived from the Baer-Nunziato model with the following choices It remains to show how the partial momentum balances can be derived using the SHTC model. We start with the balance for the relative velocity (2.15e) Now we multiply the mixture momentum balance (2.15d) with α 1 ρ 1 and add it to the previous equation (6.12) to obtain With the same choice for the interface pressure as before and the corresponding source term we have obtained the partial momentum balance for the first phase. The second balance can be obtained in a similar way (multiply with α 2 ρ 2 and subtract) Summarizing the relations we have One immediately verifies BC = I and thus the equivalence in the smooth case is proven. Following from this equivalence, we verify that the systems share the same Jacobian (2.20) and thus the same eigenvalues (2.22) and eigenvectors (2.23). However, there are crucial differences between the Baer-Nunziato system (6.9) and the SHTC system (2.15). Apart from the obvious distinction that system (2.15) can be written in conservative form, the most remarkable difference are the jump conditions for discontinuities. For system (6.9) the equation for α is not in conservative form and hence it is in general not possible to write down jump conditions, cf. [47]. However, by the argument that α stays constant across the shock the equations are decouple into Euler systems for each phase with corresponding jump conditions. In contrast to system (2.15) where the phases remain coupled across discontinuities and a shock in one phase also affects the other. Note that the particular choice for u I and p I made above is also mentioned as a consistent choice in a work by Hérard [35]. In the isentropic single temperature case it also coincides with the choice presented in Coquel et al. [14].

Numerical Results
For the SHTC model, all numerical results shown in this section have been obtained with a classical second order MUSCL-Hancock scheme, see [74] for details. For the nonconservative Baer-Nunziato system, which we solve for comparison at the end of this section, we employ a second-order path-conservative version of the MUSCL Hancock scheme, see e.g. [9,51,20,43,16]. In all cases we use the simple Rusanov-flux (local Lax-Friedrichs flux) as approximate Riemann solver.

Exact vs. numerical solution of the homogeneous system
The exact solution for the examples shown in this section is obtained as follows: 1) Assume that the eigenvalues λ i− are located left of λ C and the eigenvalues λ i+ on the right, respectively.
2) Prescribe the state left of the contact, α 1 on the right and then solve for remaining quantities on the right.
3) Choose eigenvalue next to the contact on the left side. Prescribe wave type (shock/rarefaction) and choose wave speed (head speed for a rarefaction).

4)
Repeat with remaining eigenvalue on the left side while the admissibility is constantly checked and thus certain waves might be excluded already.

5)
Repeat the afore mentioned steps for the right side.
6) Sample the complete solution.

Shock in rarefaction
The exact solution is obtained by using a contact centered inverse construction of the solution as described above. The states of this first Riemann problem (RP1 )are given in Table 1. Due to the shock inside the rarefaction we have an additional state behind the shock given byŪ . The following computation was performed using the ideal gas EOS  In Figures 8 and 9 the numerical results together with the exact solution are shown. The overall wave structure of the solution is depicted in Figure 9. It consists of a λ 2− -rarefaction (green) which is completely contained inside the λ 1− -rarefaction (blue). The contact (red) is (by construction) the middle wave. On the right we have a λ 1+rarefaction (blue) which contains a λ 2+ -shock (green) inside. The dashed blue line marks the intermediate tail of the rarefaction, i.e. the position when the rarefaction starts again after the shock. Looking at the densities and velocities in Figure 8 we verify that a rarefaction only affects the related phase. However, the interaction of the two rarefactions manifests itself in the mixture quantities ρ, u, w. For the volume fraction we observe that it only jumps at the contact as shown before and takes the initial values on each side. Special attention has to be paid to the density of the first phase. Right of the contact the rarefaction starts (not affecting phase two) until the shock occurs. There the density jumps according to the jump conditions. Continuing on the right we have a plateau of the right state of the shock and then the rarefaction continues again. The eigenvalues are shown in the upper right panel of Figure 9. One clearly sees the overlapping eigenvalues in the region where the rarefaction waves coincide. Moreover, one can see that in most of the states the system is strictly hyperbolic, but the ordering of the eigenvalues changes throughout the whole domain.

Solution without contact
The states of this second Riemann problem RP2 are given in Table 2 and we have α 1 = 0.5 for all states. The following computation was performed using the ideal gas EOS  Table 2: Primitive states of Riemann problem RP2. In Figures 10 and 11 the numerical results together with the exact solution are shown. It consists of a λ 2− -rarefaction (green) which overlaps with the λ 1− -rarefaction (blue). The contact is (by construction) not visible since all quantities stay constant across it. On the right we have an isolated λ 1+ -shock (blue) which is followed by an isolated λ 2+shock (green). Looking at the densities and velocities in Figure 10 we again verify that a rarefaction only affects the related phase. The interaction of the two rarefactions where they overlap can be observed by the changing slope in the mixture quantities ρ, u, w. The eigenvalues are shown at the right in Figure 11. One clearly sees the overlapping eigenvalues in the region where the rarefaction waves coincide. Moreover one can see that in most of the states the system is strictly hyperbolic, but the ordering of the eigenvalues changes throughout the whole domain three times.   Figure 11: Wave structure of RP2 (top): phase one (blue), contact (red) and phase 2 (green). Eigenvalues of RP2 (bottom).)

Symmetric double rarefaction
The following example (RP3) is taken from a paper of Romenski and Toro [64]. The exact solution is obtained by using a contact centered inverse construction of the solution as described previously. The states are given in Table 3 and we have α 1 = 0.9 for all states.
The following computation was performed using the following EOS   In Figures 12 and 13 the numerical results together with the exact solution are shown.
For the numerical solution we used the Force Flux together with Godunov's method as exemplary shown in [74]. This example again shows the behaviour of rarefaction waves quite nicely. Each rarefaction is seen individually in the corresponding phase but the interaction of the overlapping rarefaction waves is observed in the mixture quantities. In Figure 13 on can see that the rarefaction waves of phase one are contained in the rarefaction waves of phase two and that the corresponding eigenvalues coincide where they overlap.

Symmetric Double Shock
The following example is taken from a paper of Romenski and Toro [64]. The exact solution is obtained by using a contact centered inverse construction of the solution as described above. The states are given in Table 4 and we have α 1 = 0.9 for all states.
The following computation was performed using the following EOS   For the numerical solution we used the Force Flux together with Godunov's method as exemplary shown in [74]. This is another good test as it shows four separated shock waves one in each phase. As noted before, in contrast to the Baer-Nunziato type systems, the shocks affect every phase. In Figure 15 we can see that the shock waves related to phase one are slower than the shocks related to phase two. Moreover the order of the eigenvalues changes across every shock.

Comparison of the SHTC model with the Baer-Nunziato model
In this last section we present a comparison of the numerical solutions obtained for the conservative SHTC model discussed in this paper, and the well-known Baer-Nunziato model. We show the solution of two Riemann problems, without and with stiff relaxation source terms. The test problem RP5 has the following initial data: ρ 1,L = ρ 1,R = 2, ρ 2,L = ρ 2,R = 1, u 1,L = −2, u 1,R = +2, u 2,L = −1, u 2,R = +1, α L = 0.7, α R = 0.3. The computational domain is the interval [−1, 1], which is discretized with 10000 uniform grid cells and the final time of the simulation is t = 0.1. In Figure 16 we compare the exact solution of the Riemann problem for the SHTC system with the numerical solutions obtained for the SHTC system and the Baer-Nunziato model without any relaxation source terms. Since the wave structure consists of two rarefaction waves in both phases, the numerical results of the SHTC model and the numerical results obtained for the Baer-Nunziato system agree perfectly well with each other and with the exact solution of the Riemann problem, as expected. In Figure 17 we show the results obtained for the same initial data, but with stiff pressure and velocity relaxation, choosing the relaxation parameters as θ 1 = 10 −3 and θ 2 = 10 −8 , which corresponds to the Kapila limit of both systems. We find that the numerical solutions obtained for the conservative SHTC model and the Baer-Nunziato model agree perfectly well with each other.
The last test problem RP6 has the following initial data: ρ 1,L = 2, ρ 1,R = 1, ρ 2,L = 1, ρ 2,R = 2, u 1,L = u 1,R = 0, u 2,L = u 2,R = 0, α L = 0.7, α R = 0.3. The computational domain is again the interval [−1, 1], which is discretized with 10000 uniform grid cells and the final time of the simulation is t = 0.25. In Figure 18 we compare the exact solution of the Riemann problem for the SHTC system with the numerical solutions obtained for the SHTC system and the Baer-Nunziato model without any relaxation source terms. Since the wave structure consists of a shock and a rarefaction wave in both phases, the numerical results of the SHTC model and the numerical results obtained for the Baer-Nunziato system do not agree with each other any more, since the jump conditions of the non-conservative Baer-Nunziato system and the conservative SHTC model do not coincide. In particular, one can observe how the shock waves in one phase leave the other phase unchanged in the Baer-Nunziato model, while in the conservative SHTC system a shock in one phase always affects the other phase, which is a consequence of the jump conditions. In Figure 19 we show the results obtained for the same initial data, but with stiff pressure and velocity relaxation, choosing the relaxation parameters as θ 1 = 10 −3 and θ 2 = 10 −8 , which corresponds again to the Kapila limit of both systems. Despite the visible discrepancies observed in the homogeneous case, we find that the numerical solutions obtained in the stiff relaxation limit of the conservative SHTC model and the Baer-Nunziato model agree perfectly well with each other. This indicates that for numerical purposes it may be more beneficial to solve the SHTC system in the stiff relaxation limit, since standard numerical schemes for conservation laws can be applied, while in the Baer-Nunziato model special numerical techniques for the treatment of nonconservative terms are needed.   (a) Densities ρ1, ρ2, ρ and volume fraction α1.

Conclusion
In this paper we have presented the exact solution of the Riemann problem for the barotropic conservative two-phase model proposed in [64,63,58], which belongs to the SHTC class of symmetric hyperbolic and thermodynamically compatible systems. We have discussed the characteristic fields of the system and the admissibility criteria of shock waves. The Riemann invariants as well as the Rankine-Hugoniot jump conditions have been provided. A particular feature of the model under consideration is that some eigenvalues of the two phases may coincide, creating particular wave phenomena for the case where a shock in one phase interacts with a rarefaction wave in the other phase. We have discussed possible wave patterns and based on the mathematical entropy inequality associated with the system and the admissibility criteria of shock waves we have ruled out non-admissible wave configurations.
Independent of the potential physical applications of the model discussed in this paper, its mathematical structure and properties as well as the possible wave configurations that may appear in the solution of the Riemann problem make it a very interesting object of study from a mathematical point of view.
At the end of the paper we show exact solutions of some example Riemann problems, providing also a detailed comparison with numerical results obtained with a second order TVD scheme. For two problems, we also provide a comparison with numerical solutions obtained for the barotropic Baer-Nunziato model, both, for the homogeneous case without source terms and for the case of stiff pressure and velocity relaxation. In the case of smooth solutions, the results of the SHTC model and the Baer-Nunziato model agree perfectly well with each other, as expected, while in the presence of shock waves the numerical solutions disagree in the homogeneous case. On the contrary, in the presence of stiff velocity and pressure relaxation source terms, the numerical results obtained in the relaxation limit of the SHTC system and of the Baer-Nunziato model (Kapila limit) agree perfectly well with each other.
Future work will consider the extension of the Riemann solver presented in this paper to the non-barotropic case, solving the full model [63,58] and as also shown in the appendix. We furthermore plan to develop provably thermodynamically compatible finite volume schemes for the SHTC system, following the ideas outlined in [6,7,8].

A.1 Derivation of the Barotropic Submodel
A PDE system for compressible two-phase, two-temperature flow using a hyperbolic heat conduction model, was previously discussed in Romenski et al. [63,58]. Note, that there is another version of the SHTC model for hyperbolic heat transfer originally proposed in [65] and used in [17,18,59,5]. The mentioned version seems to be more suitable from the point of view of the general theory of SHTC models, since its PDE system can be derived from the variational principle [53]. However, since this heat transfer model is not implemented in the SHTC two-phase flow model, we use the equations formulated in [63]. Moreover, the conclusions drawn in this section should be the same for both models, since the difference in the equations for the thermal impulse does not affect the final result. Written in terms of the generalized energy the system under consideration reads as Here, α 1 is the volume fraction of the first phase which is connected with the volume fraction of the second phase α 2 by the saturation law α 1 + α 2 = 1, ρ is the mixture mass density which is connected with the phase mass densities ρ 1 , ρ 2 by the relation ρ = α 1 ρ 1 + α 2 ρ 2 . The phase mass fractions are defined as c 1 = α 1 ρ 1 /ρ, c 2 = α 2 ρ 2 /ρ and it is easy to see that c 1 + c 2 = 1. Eventually, u i = c 1 u i 1 + c 2 u i 2 is the mixture velocity, w i = u i 1 − u i 2 is the phase relative velocity. The quantities j k i is the thermal impulse of phase i and S i = c i s i is the partial entropy of phase i. The equations describe the balance law for the volume fraction, the balance law for the mass fraction, the conservation of total mass, the total momentum conservation law, the balance for the relative velocity, the two equations for the partial thermal impulses and the tow equations for the partial entropies. Moreover by construction the system is equipped with a conservation law for the total energy ∂ρ E + 1 2 u l u l ∂t + In this work the source terms are of no special interest. The details can be found in [63,58]. We want to present a derivation of the barotropic submodel in order to emphasize the differences in the conservation law for the total energy since this will affect the mathematical entropy for the system given above. By barotropic we understand that the internal energy of each phase solely depends on the corresponding phase density. More precisely we consider either an isentropic (constant entropy) or an isothermal (constant temperature) process. Let us first consider an isentropic process. This implies that the derivatives E S i vanish which leads to the following conservation law for the total energy ∂ρ E + 1 2 u l u l ∂t + ∂ ρu k E + 1 2 u l u l + p ρ + w l E w l + ρE c E w k ∂x k = 0. (A. 3) The isothermal case is a bit more involved. The reason is, that for an isothermal process we necessarily need the thermal impulse. We first rewrite the equation for the total energy (A.2) Since we have an isothermal process the partial space derivative of E S i = T i vanishes identically. We further reformulate the partial entropy balances (A.1h) and (A.1i), i.e.
and insert these two equations into the equation for the total energy. From now on we also assume that the flow has a single temperature T = T 1 = T 2 and thus the sources π i vanish. This gives ∂ρ E + 1 2 u l u l ∂t + ∂ ρu k E + 1 2 u l u l + p ρ + w l E w l + ρE c E w k ∂x k Before we further simplify the energy equations for the isentropic and the isothermal case we need to precisely discuss the mixture equation of state in each case. In the following we will discuss the isentropic and isothermal case separately.
well known Maxwell relation (omitting the phase index for the following considerations) ∂p ∂T ρ = −ρ 2 ∂s ∂ρ T (A. 11) which bases on the equality of second order derivatives of the internal energy depending on the canonical variables ρ and s, exemplary see [40]. This relation then can be used to derive a more complicated Maxwell relation which is also known but not that obvious. More precisely we want to show the following relation (A.12) This is also a well established relation, see [48,40], but for the convenience of the reader we show how it can be derived. Here g denotes the specific Gibbs energy which may be generalized to the chemical potential µ when more substances are involved, see [40]. The Gibbs energy (or sometimes free enthalpy) is given by The differentials for the internal energy and the Gibbs energy depending on the variables ρ and T are given by

A.1.4 Simplifying the Energy Equations
Using the results obtained above we can further simplify the energy equations. For the isentropic case we the relations (A.7) -(A.10) can be used to yield For the isothermal case we first rewrite the energy equality again using S i = c i s i and S = S 1 + S 2 , i.e.
∂ρ E + 1 2 u l u l ∂t + ∂ρ E − T S + 1 2 u l u l ∂t + ∂ ρu k E − T S + 1 2 u l u l + p ρ + w l E w l + ρE c E w k ∂x k .
Here F = E −T S is generalized mixture free energy and we can further write Π = Π 1 +Π 2 for the non-negative total entropy production. Now we use relations (A.18) -(A.21) to obtain −T Π = ∂ρ E − T S + 1 2 u l u l ∂t + ∂ ρu k E − T S + 1 2 u l u l + p ρ + w l E w l + ρE c E w k ∂x k In [58] a single temperature model is presented. The governing equations are the same as in the isothermal case, the energy equation reads However, after some manipulations this equation may also be written in the same form as it is used in this work. Thus the thermal impulses are again basically hidden in the energy equality.

A.2 Calculation for Coinciding Eigenvectors
We want to give the detailed calculation for the statement (3.1). Let us first assume that the vector R C is a multiple of a fixed vector R i± , see (2.23). We can directly conclude that Note that the specific correct sign ± is determined by the non-zero components of the eigenvector. To show the reverse direction we assume that λ C coincides with a fixed λ i± , see (2.22). From this we get that the corresponding ε i vanishes. To simplify the notation let us assume w.l.o.g. that i = 1 and thus ε 1 = 0. We further obtain Thus we get This shows the desired relation.

A.3 Regularity of the Mass Flux System Matrix
In Section 4.2 we derived a linear system for the squares of the partial mass fluxes (4.20). It remains to verify that the matrix is regular, i.e. the determinant and thus the determinant will only vanish iff at least one of densities does not jump across the shock. If the partial density would not jump across the shock we obtain from the continuity of the corresponding partial mass flux that the associated partial velocity also would not jump. From the Lax condition of the eigenvalues we hence can conclude that for a shock in phase µ ∈ {1, 2} the partial density of phase µ must jump. Let us now consider the other phase ν ∈ {1, 2}, ν = µ. If we assume ρ ν = 0 we obtain for (4.18) Thus we have for the momentum jump condition (4.13d) and for the relative velocity jump condition (4.13e) −Q µ u µ + p µ = 0, (A.22)