Contents lists available atScienceDirect

## Journal of Computational Physics

www.elsevier.com/locate/jcp

## A mixture-energy-consistent six-equation two-phase numerical model for ﬂuids with interfaces,

## cavitation and evaporation waves

### Marica Pelanti

^{a}

^{,}^{∗}

### , Keh-Ming Shyue

^{b}

a*Department of Mechanical Engineering, École Nationale Supérieure de Techniques Avancées – ENSTA ParisTech,*
*828, Boulevard des Maréchaux, 91762 Palaiseau Cedex, France*

b*Department of Mathematics, National Taiwan University, Taipei 106, Taiwan*

a r t i c l e i n f o a b s t r a c t

*Article history:*

Received 31 March 2013

Received in revised form 30 October 2013 Accepted 1 December 2013

Available online 10 December 2013

*Keywords:*

Multiphase compressible ﬂow models Mechanical relaxation

Thermo-chemical relaxation Cavitation

Phase transition Finite volume schemes Wave propagation algorithms Riemann solvers

We model liquid–gas ﬂows with cavitation by a variant of the six-equation single-velocity two-phase model with stiff mechanical relaxation of Saurel–Petitpas–Berry (Saurel et al., 2009)[9]. In our approach we employ phasic total energy equations instead of the phasic internal energy equations of the classical six-equation system. This alternative formulation allows us to easily design a simple numerical method that ensures consistency with mixture total energy conservation at the discrete level and agreement of the relaxed pressure at equilibrium with the correct mixture equation of state. Temperature and Gibbs free energy exchange terms are included in the equations as relaxation terms to model heat and mass transfer and hence liquid–vapor transition. The algorithm uses a high-resolution wave propagation method for the numerical approximation of the homogeneous hyperbolic portion of the model. In two dimensions a fully-discretized scheme based on a hybrid HLLC/Roe Riemann solver is employed. Thermo-chemical terms are handled numerically via a stiff relaxation solver that forces thermodynamic equilibrium at liquid–vapor interfaces under metastable conditions. We present numerical results of sample tests in one and two space dimensions that show the ability of the proposed model to describe cavitation mechanisms and evaporation wave dynamics.

©2013 Elsevier Inc. All rights reserved.

**1. Introduction**

The modeling of cavitating ﬂows is relevant in numerous areas of engineering, from naval and submarine systems design to aerospace and nuclear power plants technologies. Cavitating ﬂuids are multiphase mixtures that often involve complex hydrodynamic and thermodynamic processes: liquid–vapor phase transition, dynamical creation of interfaces, vapor struc- tures collapse, and associated shock wave formation and interaction (cf.[1–3]). As a further reason of complexity, in many industrial applications these ﬂows occur in irregular geometries and they have a multi-dimensional character.

Extensive work has been dedicated in the past decades to the simulation of cavitating ﬂows and liquid–vapor ﬂows with phase change, see for instance[4–15]and the references therein. Among the different modeling approaches, the class of hyperbolic compressible multiphase models stemming from the original model of Baer–Nunziato[16] has shown great capabilities in describing the complex wave patterns and thermodynamic mechanisms of cavitation. A ﬁrst essential feature

### *

Corresponding author. Tel.: +33 1 69 31 98 19; fax: +33 1 69 31 99 97.*E-mail addresses:*marica.pelanti@ensta-paristech.fr(M. Pelanti),shyue@ntu.edu.tw(K.-M. Shyue).

0021-9991/$ – see front matter ©2013 Elsevier Inc. All rights reserved.

http://dx.doi.org/10.1016/j.jcp.2013.12.003

of these models is that compressibility is taken into account for all phases, vapor as well as liquid. This is fundamental to correctly capture wave propagation phenomena and acoustic perturbations, and it is particularly crucial when liquid–

vapor transition occurs[8]. Another important property is that these models can retain temperature and Gibbs free energy non-equilibrium effects, thus they are able to capture metastable states as well as evaporation fronts, when heat and mass transfer processes are included in the physical description through thermal and chemical relaxation source terms.

There exist various formulations of compressible temperature non-equilibrium multiphase ﬂow models, depending on the assumptions on mechanical and kinetic phase equilibrium. In choosing a particular model, one has to ﬁnd a good compromise between the accuracy of the description of the physical phenomena and the ability of conceiving robust and eﬃcient numerical methods. In the present work, we are interested in the hyperbolic single-velocity six-equation model proposed by Saurel, Petitpas and Berry in[9]for compressible two-phase ﬂows, see also Zein et al.[17]. This model consists of an advection equation for the volume fraction of one phase, mass and internal energy equations for each phase, and a mixture momentum equation. The six-equation model assumes instantaneous velocity equilibrium between the two phases, but it retains mechanical, thermal and chemical non-equilibrium effects. In the limit of instantaneous pressure relaxation the model reduces to the well known compressible two-phase ﬂow model of Kapila et al.[18]. Nonetheless, as emphasized in[9], and as we brieﬂy recall in Section2, numerically it is more advantageous to solve the six-equation system with stiff mechanical relaxation rather than the Kapila et al.[18]pressure-equilibrium ﬁve-equation model system.

The single-velocity six-equation two-phase model with stiff pressure relaxation was employed in [9]for applications to interface problems and mechanical cavitation processes (that is cavitation with no phase transition). It was later used by Zein et al. in [17] to simulate liquid–vapor transition in metastable liquids. One diﬃculty of the numerical algorithm illustrated in the latter work, as noted by the authors, is that it may require a very small time step for stability for some expansion problems with phase transition, due to the stiffness of the chemical relaxation terms. Only one-dimensional numerical results are presented by the authors in[17].

The aim of the present paper is to conceive a new multiphase ﬂow computational model on the basis of the six-equation system of[9]that could deal eﬃciently with interfaces, cavitation and evaporation waves, while retaining simplicity and time-affordability. The key idea of our approach is to employ an alternative mathematical formulation of the standard six-equation model system [9]in the numerical discretization. Rather than using the two phasic internal energy equations of the classical model, in our algorithm we employ two equations for the phasic total energies. Mathematically, these two model systems are equivalent. The present model, however, is numerically advantageous with respect to the standard one, since it allows us to easily design a simple numerical method that ensures important consistency properties with mixture total energy conservation and with the mixture thermodynamic state. More speciﬁcally, ﬁrst, we are able to automatically recover a conservative discrete form of the mixture total energy equation, whereas the classical six-equation model system needs to be augmented with an additional conservation law for the mixture total energy to correct the thermodynamic state [9,17]. Secondly, as a consequence of the mixture total energy conservation consistency property, we are able to easily ensure agreement of the relaxed pressure at equilibrium with the correct mixture equation of state for the full six-equation two-phase model that includes mechanical and thermo-chemical stiff relaxation effects. Relaxation terms are therefore eﬃciently handled.

To numerically solve the proposed two-phase model with pressure, temperature, and Gibbs free energy relaxation, we employ a simple fractional step approach that consists of the homogeneous hyperbolic system solution step, and a sequence of steps thereafter to solve systems of ordinary differential equations containing the relaxation source terms. A high- resolution wave propagation method based on Riemann solvers (HLLC and Roe) (cf. [19]) is employed for the numerical solution of the homogeneous hyperbolic system. The algorithm is easily implemented in the framework of the CLAWPACK software package[20]. For solving the ordinary differential equations with stiff relaxation sources, we have devised robust solvers that drive the mixture to the desired equilibrium conditions in a sequence of relaxation processes (cf. [21,10,8,17, 14]). In this procedure, similar to[8,22], thermodynamic equilibrium is forced at liquid–vapor interfaces under metastable conditions. Numerically for this task we employ an idea similar to[14,23]that uses the thermodynamic equilibrium condi- tions to reduce the solution of the ODEs relaxation problem to the solution of a simple system of algebraic equations for the equilibrium state variables.

This paper is organized as follows. In Section2.1, we begin by recalling the six-equation single-velocity model with stiff mechanical relaxation of Saurel–Petitpas–Berry[9]for compressible two-phase ﬂows. We then propose in Section2.2a vari- ant of this model system, by employing phasic total energy equations in the mathematical formulation instead of the phasic internal energy equations of the classical approach. The extended model that includes thermal and chemical relaxation terms to model heat and mass transfer is described in Section2.3. In Section3we illustrate the numerical method to solve the basic model system with mechanical relaxation only. In this section we also discuss the mixture-energy-consistency property of the algorithm. The numerical treatment of temperature and Gibbs free energy relaxation source terms is de- scribed in Section 4. Finally, in Section5we present a selection of numerical results obtained by employing the proposed method with and without activation of heat and mass transfer.

**2. The six-equation single-velocity two-phase ﬂow model**
*2.1. Phasic-internal-energy-based formulation*

The six-equation single-velocity compressible two-phase ﬂow model with stiff mechanical relaxation proposed by Saurel et al.[9]has the form

*∂*

*t*

*α*

1### +

*u*

### · ∇ *α*

1### = *μ* *(*

*p*1

### −

*p*2

*),*

(1a)
*∂*

*t*

*(* *α*

1*ρ*

1*)* *+ ∇ · (* *α*

1*ρ*

1

^{u}### *)* =

^{0}

*,*

(1b)
*∂*

*t*

*(* *α*

2*ρ*

2*)* *+ ∇ · (* *α*

2*ρ*

2

^{u}### *)* =

^{0}

*,*

(1c)
*∂*

*t*

*(* *ρ*

*u*

*)* *+ ∇ · (* *ρ*

*u*

### ⊗

^{u}*)* *+ ∇(* *α*

1*p*

_{1}

### + *α*

2*p*

_{2}

*)* =

^{0}

*,*

(1d)
*∂*

*t*

*(* *α*

1*E*1

*)* *+ ∇ · (* *α*

1*E*1

*u*

### *)* + *α*

1*p*

_{1}

### ∇ ·

^{u}### = − *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*),*

(1e)
*∂*

_{t}*(* *α*

_{2}

*E*2

*)* *+ ∇ · (* *α*

_{2}

*E*2

*u*

### *)* + *α*

_{2}

*p*

_{2}

### ∇ ·

^{u}### = *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*).*

(1f)
Here *α*

*k*

*is the volume fraction of phase k, k*=

^{1}

*,*2 (

*α*

1+*α*

2=^{1),}

*ρ*

*k*

*is the phasic density, p*

*is the phasic pressure, and*

_{k}*E*

*k*=

*ρ*

*k*

*ε*

*k*is the phasic internal energy, with

*ε*

*k*denoting the phasic speciﬁc internal energy. We have

*ρ*

=*α*

1*ρ*

1+*α*

2*ρ*

2
denoting the mixture density, and*u representing the ﬂow velocity vector. The source terms appearing in*(1a),(1e), and(1f) model mechanical relaxation. In these terms

*μ*

*>0 represents the pressure relaxation parameter and p*

_{I}is the interface

*pressure, p*I=

^{Z}^{2}

^{p}

_{Z}^{1}

_{1}

^{+}

_{+}

^{Z}

_{Z}^{1}

_{2}

^{p}^{2}

^{, where Z}*k*=

*ρ*

_{k}*c*

_{k}^{2}

*is the acoustic impedance of phase k, and c*

_{k}*is the sound speed of phase k. We*assume an inﬁnite-rate pressure relaxation with

*μ*

→ ∞, therefore mechanical equilibrium is reached instantaneously.
It is known (cf. [9]) that in this instantaneous pressure relaxation limit the six-equation model above reduces to the single-velocity single-pressure ﬁve-equation model of Kapila et al.[18](see also e.g.[24]). The ﬁve-equation model system is composed of two phasic mass balance equations, the mixture momentum equation, the mixture total energy equation, and an evolution equation for the volume fraction of one phase with a source term that results from the asymptotic limit of instantaneous velocity and pressure equilibrium of the non-equilibrium compressible two-phase ﬂow model of Baer and Nunziato[16]. Despite this equivalence, the six-equation model(1)offers signiﬁcant advantages in the numerical approxi- mation, compared to the ﬁve-equation model[18]. The main numerical issues in the solution of the ﬁve-equation system come from the non-conservative contribution in the volume fraction equation that depends on the divergence of the ﬂow velocity and on the phasic impedances [9,59]. The variation of the volume fraction across acoustic waves associated to this term makes the construction of approximate Riemann solvers more challenging. In particular, the presence of this non-conservative contribution makes it diﬃcult to preserve volume fraction positivity, especially when shocks and strong rarefaction waves are involved[9].

*2.1.1. System’s closure*

The closure of system(1)is obtained through the speciﬁcation of an equation of state for each phase, which we choose
to express in terms of *E**k* and

*ρ*

*k*

*, p*

*=*

_{k}

^{p}*k*

*(E*

*k*

*,*

*ρ*

*k*

*), k*=

^{1}

*,2. The phasic sound speed can be written as c*

*=*

_{k}*κ*

*k*

*h*

*+*

_{k}*χ*

*k*,

*where h*

*=*

_{k}

^{E}

^{k}

_{ρ}^{+}

_{k}

^{p}*is the phasic speciﬁc enthalpy,*

^{k}*κ*

*k*=

^{∂}

^{p}

^{k}

^{(}_{∂}^{E}_{E}

^{k}

_{k}

^{,}^{ρ}

^{k}

^{)}^{and}

*χ*

*k*=

^{∂}

^{p}

^{k}

^{(}_{∂}^{E}_{ρ}

^{k}

_{k}

^{,}^{ρ}

^{k}*. The mixture sound speed of this model is*

^{)}*c*

### =

*Y*_{1}*c*^{2}_{1}

### +

*2*

^{Y}*c*

^{2}

_{2}(2)

*where Y** _{k}*=

^{α}

^{k}

_{ρ}^{ρ}

^{k}*is the mass fraction of phase k (Y*1+

*2=1). Notice the monotonic character of(2)with respect to the volume fraction in contrast to the non-monotonic behavior of the Wood’s sound speed[25]of the Kapila et al. model[18].*

^{Y}This feature of the 6-equation model also represents an advantage over the 5-equation model in the numerical approxima- tion. As explained in [26], the non-monotonic behavior of the 5-equation model’s sound speed in the numerical diffusion zone of an interface may result in the presence of two sonic points in this region even when the ﬂow is subsonic in the two pure ﬂuids. The construction of approximate Riemann solvers able to handle robustly and eﬃciently these sonic points does not appear a simple task. In [9]the authors also explain that the sound speed non-monotonicity might affect the propagation of acoustic waves interacting with the interfacial zone, and result in a temporal delay in the wave transmission.

*Let us remark that the pressure law p** _{k}*=

^{p}*k*

*(E*

*k*

*,*

*ρ*

*k*

*), by using a more rigorous terminology, represents an incomplete*equation of state. This law suﬃces to determine the ﬂuid dynamics when thermal and chemical phenomena are neglected,

*as in the model above. However, a caloric law T*

*k*=

^{T}*k*

*(p*

*k*

*,*

*ρ*

*k*

*)for each phasic temperature T*

*k*is also needed in order to completely describe the thermodynamic state of the ﬂuid (see e.g.[27]). The full thermodynamic characterization is required when heat and mass transfer terms are included in the ﬂow model for problems with phase change, see Section2.3.

Here we will restrict our study to the case of species governed by the stiffened gas equation of state (SG EOS):

*p*_{k}

*(*

*E*

*k*

*,* *ρ*

_{k}*)* *= (* *γ*

_{k}_{−}

1*)*

*E*

*k*

### − *γ*

_{k}*π*

_{k}_{− (} *γ*

_{− (}

_{k}_{−}

1*)* *η*

_{k}*ρ*

_{k}*,*

(3a)
*T*

_{k}*(*

*p*

_{k}*,* *ρ*

*k*

*)* =

^{p}

^{k}### + *π*

*k*

*C*_{vk}

*ρ*

*k*

*(* *γ*

*k*

### −

^{1}

*)*

^{(3b)}

*for k*=^{1}*,*2, where

*γ*

*k*,

*π*

*k*,

*η*

*k*

*, and C*

*vk*

*are material-dependent parameters. The associated entropy s*

*k*

*(p*

*k*

*,T*

*k*

*)*and the Gibbs

*free energy (chemical potential) g*

_{k}*(p*

_{k}*,T*

_{k}*)*=

^{h}*k*−

^{T}*k*

*s*

*are given by*

_{k}*s*_{k}

*(*

*p*

_{k}*,*

*T*

_{k}*)* =

^{C}*vk*log

*T*

*k*

*γ*

*k*

*(*

*p*

_{k}### + *π*

*k*

*)*

^{γ}

^{k}^{−}

^{1}

### + *η*

^{}

_{k}*,*

(3c)
*g*_{k}

*(*

*p*

_{k}*,*

*T*

_{k}*)* =

*γ*

_{k}*C*

_{vk}### − *η*

^{}

_{k}^{}

*T*

_{k}### −

^{C}*vk*

*T*

*log*

_{k}*T*

_{k}

^{γ}

^{k}*(*

*p*

_{k}### + *π*

_{k}*)*

^{γ}

^{k}^{−}

^{1}

### + *η*

_{k}*,*

(3d)
where

*η*

^{}

*is a constant. With the SG EOS assumption(3), we have*

_{k}*κ*

*k*

*= (*

*γ*

*k*−

^{1}

*)*,

*χ*

*k*

*= −(*

*γ*

*k*−

^{1}

*)*

*η*

*k*, and the phasic sound

*speed can be expressed as c*

*=*

_{k}*γ*

*k*

*p*

*+*

_{k}*π*

*k*

*ρ**k* .

The mixture speciﬁc internal energy for the model considered here is deﬁned as

*ε*

_{=}

*Y*1

*ε*

1+*2*

^{Y}*ε*

2, and, equivalently, the
mixture internal energy per unit volume is *E =*

*ρε*

_{=}

*α*

1*E*1+

*α*

2*E*2. The latter relation, by using the isobaric assumption

*p*1=

*2=*

^{p}*p in the energy laws*

*E*

*k*

*(p*

*k*

*,*

*ρ*

*k*

*), k*=

^{1}

*,*2, gives the mixture equation of state, which determines implicitly the

*mixture pressure law p*=

^{p}(E,*ρ*

1*,*

*ρ*

2*,*

*α*

1*)*:

*E*

### = *α*

1*E*1

*(*

*p*

*,* *ρ*

1*)* + *α*

2*E*2

*(*

*p*

*,* *ρ*

2*).*

(4)
In the case with the SG EOS(3)*we ﬁnd an explicit expression for the mixture pressure p* [28]:

*p*

*(*

*E*

*,* *ρ*

1*,* *ρ*

2*,* *α*

1*)* =

^{E}*− (* *α*

_{1}

*ρ*

_{1}

*η*

_{1}

_{+} *α*

_{2}

*ρ*

_{2}

*η*

_{2}

*)* *− (*

^{α}_{γ}^{1}

_{1}

^{γ}_{−}

^{1}

^{π}_{1}

^{1}

### +

^{α}_{γ}^{2}

_{2}

^{γ}_{−}

^{2}

^{π}_{1}

^{2}

*)*

*α*1

*γ*1−1

### +

_{γ}_{2}

^{α}_{−}

^{2}

_{1}

*.*

(5)
When physically relevant values of the ﬂow state variables are deﬁned in the region of thermodynamic stability, the single-
velocity six-equation model (1) is hyperbolic, that is it has real eigenvalues and a complete set of eigenvectors, see for
instance [17]andAppendix A. Let us also recall that(1)represents the asymptotic limit of the hyperbolic seven-equation
two-velocity two-phase ﬂow model of Saurel and Abgrall [21] for instantaneous kinetic relaxation (see e.g. [17,18]). One
*advantage of the six-equation model with respect to the two-velocity model is that the order of the eigenvalues is a priori*
known, yielding an easy decomposition of waves in computing approximate solutions to Riemann problems (cf.[29–31]).

*2.2. Phasic-total-energy-based formulation*

For numerical reasons that we will discuss in the following, we propose to consider a mathematically equivalent for-
mulation of the model (1), obtained by replacing the two phasic internal energy equations with two phasic total energy
*equations. We denote with E** _{k}*=

*E*

*k*+

^{1}

_{2}

*ρ*

*k*

*·*

^{u}*u the total energy of phase k. The alternative form of the six-equation model*reads

*∂*

*t*

*α*

1### +

^{u}### · ∇ *α*

1### = *μ* *(*

*p*

_{1}

### −

*2*

^{p}*),*

(6a)
*∂*

*t*

*(* *α*

1*ρ*

1*)* *+ ∇ · (* *α*

1*ρ*

1*u*

*)* =

0*,*

(6b)
*∂*

_{t}*(* *α*

_{2}

*ρ*

_{2}

*)* *+ ∇ · (* *α*

_{2}

*ρ*

_{2}

*u*

*)* =

^{0}

*,*

(6c)
*∂*

*t*

*(* *ρ*

^{u}### *)* *+ ∇ · (* *ρ*

^{u}### ⊗

^{u}*)* *+ ∇(* *α*

1*p*

_{1}

### + *α*

2*p*

_{2}

*)* =

^{0}

*,*

(6d)
*∂*

*t*

*(* *α*

1*E*1

*)* *+ ∇ · (* *α*

1*E*1

*u*

### + *α*

1*p*1

*u*

### *)* *+ Σ(*

^{q}*,* ∇

^{q}*)* = − *μ*

*p*I

*(*

*p*1

### −

*2*

^{p}*),*

(6e)
*∂*

*t*

*(* *α*

2*E*

_{2}

*)* *+ ∇ · (* *α*

2*E*

_{2}

*u*

### + *α*

2*p*

_{2}

*u*

### *)* *− Σ(*

^{q}*,∇*

^{q}*)* = *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*),*

(6f)
where the non-conservative term*Σ*appearing in the phasic total energy equations is

*Σ (*

*q*

*,* ∇

^{q}*)* = −

^{u}### ·

*Y*_{2}

*∇(* *α*

_{1}

*p*

_{1}

*)* −

*1*

^{Y}*∇(* *α*

_{2}

*p*

_{2}

*)*

### = −

^{u}### ·

*(*

*Y*

_{2}

*p*

_{1}

### +

*1*

^{Y}*p*

_{2}

*)* ∇ *α*

_{1}

_{+} *α*

_{1}

*Y*

_{2}

### ∇

*1*

^{p}### − *α*

_{2}

*Y*

_{1}

### ∇

*2*

^{p}### *,*

(6g)
*with q denoting the vector of the system unknowns, see* (6). Note that unlike the previous model (1) with the phasic
internal energy equations, here the simple sum of the two non-conservative phasic total energy equations (6e) and(6f)
*recovers the equation expressing conservation of the mixture total energy E*=*E +*^{1}_{2}

*ρ*

*·*

^{u}*=*

^{u}*α*

1*E*

_{1}+

*α*

2*E*

_{2}:

*∂*

*t*

*E*

*+ ∇ · (*

^{E}^{u}### + *α*

1*p*

_{1}

^{u}### + *α*

2*p*

_{2}

^{u}*)* =

^{0}

*.*

(7)
*This feature is beneﬁcial in the numerical approximation of the model to ensure consistency with the conservation of E, see*Section3.1.

*2.3. Model with heat and mass transfer*

A classical way to model thermal and chemical inter-phase phenomena mathematically consists in introducing additional
heat and mass transfer source terms*Q*and*m, respectively, into the original model system (cf.*˙ [14,32,21,8,11,17,22]). In the
present work we follow in particular the modeling approach of Saurel and co-workers[8,22], employing the six-equation
model with phasic total energy equations(6)as the basic system. The ﬂow model augmented with heat and mass transfer
terms takes the form:

*∂*

*t*

*α*

1### +

^{u}### · ∇ *α*

1### = *μ* *(*

*p*1

### −

*2*

^{p}*)* +

*m*

### ˙ *ρ*

I
*,*

(8a)
*∂*

*t*

*(* *α*

1*ρ*

1*)* *+ ∇ · (* *α*

1*ρ*

1*u*

### *)* = ˙

*m*

*,*

(8b)
*∂*

*t*

*(* *α*

2*ρ*

2*)* *+ ∇ · (* *α*

2*ρ*

2

^{u}### *)* = − ˙

^{m}*,*

(8c)
*∂*

*t*

*(* *ρ*

^{u}*)* *+ ∇ · (* *ρ*

^{u}### ⊗

^{u}*)* *+ ∇(* *α*

1*p*

_{1}

### + *α*

2*p*

_{2}

*)* =

^{0}

*,*

(8d)
*∂*

*t*

*(* *α*

1*E*

_{1}

*)* *+ ∇ · (* *α*

1*E*

_{1}

*u*

### + *α*

1*p*

_{1}

*u*

### *)* *+ Σ(*

^{q}*,* ∇

^{q}*)* = − *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*)* +

*Q*

### +

*I*

^{e}*m*

### ˙ *,*

(8e)
*∂*

*t*

*(* *α*

2*E*

_{2}

*)* *+ ∇ · (* *α*

2*E*

_{2}

*u*

### + *α*

2*p*

_{2}

*u*

### *)* *− Σ(*

^{q}*,* ∇

^{q}*)* = *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*)* −

*Q*

### −

*I*

^{e}*m*

### ˙ *,*

(8f)
where the terms*Q*

^{and}

*m can be written as*˙

*Q*

*= θ(*

*T*2

### −

*T*1

*),*

(9a)
### ˙

*m*

### = *ν* *(*

*g*

_{2}

### −

*1*

^{g}*).*

(9b)
Here*θ*and

*ν*

are the thermal and chemical relaxation parameters, respectively. These parameters are assumed to be inﬁnite
at selected locations, while they are set to zero elsewhere. More speciﬁcally, thermal transfer is activated at liquid–vapor
interfaces, and thermo-chemical transfer is activated at liquid–vapor interfaces under metastable thermodynamic conditions
*(liquid temperature T*liq

*higher than the saturation temperature T*sat at the given pressure):

*θ* =

### ∞

^{if}

### I

### *α*

1^{1}

### −

I*,*

0 otherwise

*,*

^{(10a)}

*ν* =

### ∞

^{if}

### I

### *α*

1^{1}

### −

I*and T*

_{liq}

*>*

*T*

_{sat}

*,*

0 otherwise

*.*

^{(10b)}

Here the parameter

### Iidentiﬁes liquid–vapor interface locations (e.g.

### I=

^{10}

^{−}

^{4}

^{). See}

^{[8]}for further discussion on the choice of tolerances.

We refer the reader to[33,34]*for a rigorous and systematic analytical study of Liu’s subcharacteristic condition*[35]for the
six-equation two-phase model with the three levels of relaxation considered here, namely, pressure relaxation, simultaneous
pressure and temperature relaxation, simultaneous pressure, temperature and chemical potential relaxation. In particular,
the works[33,34]show that the subcharacteristic condition (characteristic speeds of the relaxation system at least as large
as the characteristic speeds of the relaxed equilibrium system) is satisﬁed for all cases as long as physical fundamental
positivity conditions on a set of physically positive thermodynamic variables are fulﬁlled.

*Let us notice that a theoretical pressure–temperature (p–T ) saturation curve can be obtained by imposing the equilib-*
*rium condition on the Gibbs free energy g*_{1}=* ^{g}*2 for the liquid and vapor phases (cf.[36,37,8]). With the SG EOS(3a)–(3d),

*the p–T saturation curve is deﬁned by the equation*

*A*_{S}

### +

^{B}^{S}

*T*

### +

*S*

^{C}*log T*

### +

*Slog*

^{D}*(*

*p*

### + *π*

1*)* −

^{log}

*(*

*p*

### + *π*

2*)* =

^{0}

*,*

(11a)
with

*A*_{S}

### =

^{C}

^{p1}### −

^{C}*p2*

### + *η*

_{2}

^{}

_{−} *η*

^{}

_{1}

*C*_{p2}

### −

^{C}*v2*

*,*

*B*

_{S}

### = *η*

1### − *η*

2
*C*_{p2}

### −

^{C}*v2*

*,*

*C*

_{S}

### =

^{C}

^{p2}### −

^{C}*p1*

*C*_{p2}

### −

^{C}*v2*

*,*

*D*

_{S}

### =

^{C}

^{p1}### −

^{C}*v1*

*C*_{p2}

### −

^{C}*v2*

*.*

(11b)
*Here C**pk* *denotes the heat capacity at constant pressure, C**pk*=^{C}*vk*

*γ*

*k*

*, k*=

^{1}

*,*2. The parameters in the SG EOS (3a)–(3d)

*can be chosen to ﬁt the above theoretical p–T curve with the experimental curve for the considered material*[36]. The

*theoretical p–T curve is used in the numerical algorithm, see Section*4.2.

Concerning the interfacial density

*ρ*

I*and the interfacial speciﬁc total energy e*Ithat appear in the mass and heat transfer source terms, their expression can be obtained by imposing appropriate thermodynamic constraints on the thermal and chemical processes. For instance, in [17], Zein et al. determine these quantities in a model similar to (8) by assuming pressure and temperature equilibrium during the relaxation processes. However, thanks to the assumption (10) on the relaxation parameters and to the particular numerical algorithm used for the treatment of thermal and chemical source terms, there is no need to specify the expressions of

*ρ*

I*and e*

_{I}(see Section4). Let us also note that, due to the symmetry of the source terms in(8), the mixture mass and total energy equations are recovered as for the model in(6).

To end this section, for the ease of later reference, we write(8)in a compact form as

*∂*

*t*

*q*

### + ∇ ·

^{f}*(*

*q*

*)* + *σ* *(*

*q*

*,* ∇

^{q}*)* *= ψ*

*μ*

*(*

*q*

*)* *+ ψ*

*θ*

*(*

*q*

*)* *+ ψ*

*ν*

*(*

*q*

*),*

(12a)
where
*q*

### =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝ *α*

1
*α*

1*ρ*

1
*α*

2*ρ*

2
*ρ*

*u*

### *α*

1*E*1

*α*

2*E*2

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

*,*

*f*

*(*

*q*

*)* =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝

0

*α*

1*ρ*

1*u*

### *α*

2*ρ*

2

^{u}*ρ*

*u*

### ⊗

*u*

*+ (* *α*

1*p*1

### + *α*

2*p*2

*)I* *α*

1*(*

*E*1

### +

*p*1

*)*

*u*

*α*

2*(*

*E*2

### +

*2*

^{p}*)*

*u*

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

*,* *σ* *(*

*q*

*,* ∇

^{q}*)* =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝

*u*

### · ∇ *α*

1
0 0 0

*Σ (*

*q*

*,* ∇

*q*

*)*

*−Σ(*

^{q}*,* ∇

^{q}*)*

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

*,*

(12b)
*ψ*

*μ*

*(*

*q*

*)* =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝

*μ* *(*

*p*1

### −

*p*2

*)*

0
0
0
### − *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*)* *μ*

*p*

_{I}

*(*

*p*

_{1}

### −

*2*

^{p}*)*

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

*,* *ψ*

*θ*

*(*

*q*

*)* =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝

0 0 0 0*θ (*

*T*

_{2}

### −

*1*

^{T}*)*

*−θ(*

*2*

^{T}### −

*1*

^{T}*)*

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

*,* *ψ*

*ν*

*(*

*q*

*)* =

### ⎛

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎜ ⎜

### ⎝

*ν*

^{g}^{2}

^{−}

_{ρ}

^{g}^{1}

I

*ν* *(*

*g*2

### −

*g*1

*)*

### − *ν* *(*

*g*2

### −

*g*1

*)*

0
*ν*

*e*I

*(*

*g*2

### −

*1*

^{g}*)*

### − *ν*

*e*

_{I}

*(*

*g*

_{2}

### −

*1*

^{g}*)*

### ⎞

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎟ ⎟

### ⎠

(12c)

with*Σ* as in(6g). Above we have put into evidence the conservative portion of the spatial derivative contributions in the
system as∇ · * ^{f}(q)*, and we have indicated the non-conservative term as

*σ*

*(q,*∇

*. The source terms*

^{q})*ψμ(q)*,

*ψθ(q)*,

*ψν(q)*contain mechanical, thermal and chemical relaxation terms, respectively.

**3. Numerical solution of the model system with mechanical relaxation**

We begin by considering the solution of the model system(12)with only mechanical relaxation:

*∂*

*t*

*q*

### + ∇ ·

^{f}*(*

*q*

*)* + *σ* *(*

*q*

*,* ∇

^{q}*)* *= ψ*

*μ*

*(*

*q*

*).*

(13)
The treatment of thermal and chemical relaxation terms will be described in Section 4. To numerically solve the system
above we use a fractional step technique, similar to [9,17], where we alternate between the solution of the homogeneous
hyperbolic system and the solution of a system of ordinary differential equations that takes into account pressure relaxation
source terms. That is, the algorithm consists of the following steps:
*1. Homogeneous hyperbolic system. We solve over a time intervalt the homogeneous hyperbolic portion of*(13):

*∂*

*t*

*q*

### + ∇ ·

^{f}*(*

*q*

*)* + *σ* *(*

*q*

*,* ∇

^{q}*)* =

^{0}

*.*

(14)
*2. Stiff mechanical relaxation. We solve in the limit*

*μ*

_{→ ∞}the system of ordinary differential equations (ODEs)

*∂*

*t*

*q*

*= ψ*

*μ*

*(*

*q*

*).*

(15)
*This step drives the two-phase ﬂow to mechanical equilibrium with an equilibrium relaxed pressure p*1=* ^{p}*2=

*step the partial densities, the mixture momentum, the mixture total energy, and the mixture internal energy remain constant. The volume fraction*

^{p. In this}*α*

1*, the mixture pressure p, and the phasic internal energies*

*α*

*k*

*E*

*k*,

*E*

*k*=

*E*

*k*

*(p, (*

*α*

*k*

*ρ*

*k*

*)/*

*α*

*k*

*)*

*for k*=

^{1}

*,*2, are updated before returning to Step 1.

*3.1. Mixture-energy-consistent discretization*

Before illustrating each step of the algorithm in more details, we emphasize that in the design of a numerical method for computing approximate solutions of the two-phase model (12) it is important to ensure at all times conservation of the quantities that are physically conserved, namely the partial densities

*α*

_{k}*ρ*

_{k}*, k*=

^{1}

*,*2 (together with the mixture density

*ρ*

_{=}

*α*

1*ρ*

1+*α*

2*ρ*

2), the mixture momentum*ρ*

*u, and the mixture total energy E*=

*α*

1*E*1+

*α*

2*E*2. Note in particular that the values of the equilibrium pressure to be used at the beginning of the homogeneous system solution step must satisfy the mixture EOS (4)for values of

*E =*−

^{E}^{1}

_{2}

*ρ*

*·*

^{u}*u that correspond to conservation-consistent discrete values of E, so as to*approximate correctly the ﬂow thermodynamic state. Godunov-type schemes can be relatively easily formulated to preserve conservation at the discrete level of quantities that are governed by conservation laws (cf.[38,19]). However, the two-phase mathematical model (either in the form (1) or(6)) does not contain the conservation law for the total energy, but two phasic energy equations from which the total energy is recovered. The diﬃculty then is to discretize the phasic energy

equations in such a way that total energy conservation is fulﬁlled at the discrete level, and that consistency with the correct thermodynamic state is ensured in the sense made precise hereafter.

Let us denote with superscript 0 the quantities computed by solving the homogeneous system in Step 1 of the algorithm above, and with superscript∗the quantities at mechanical equilibrium computed in Step 2. As mentioned above (see also Section3.3), we have

*ρ*

_{k}^{∗}=

*ρ*

_{k}^{0}

*, k*=

^{1}

*,*2,

*(*

*ρ*

^{u}*)*

^{∗}

*= (*

*ρ*

^{u}*)*

^{0},

*E*

^{∗}=

*E*

^{0}

^{, and E}^{∗}=

^{E}^{0}

*. Let us also denote with E*

^{0}

^{,}^{C}discrete values

*of the mixture total energy that come from a conservative approximation of the conservation law for E in*(7).

* Deﬁnition 3.1. We say that the numerical scheme based on the fractional step algorithm above is mixture-energy-consistent*
if the following two properties are satisﬁed

(i) Mixture total energy conservation consistency:

*E*^{0}

### =

*E*

^{0}

^{,}^{C}

*,*

(16a)
*where E*^{0}*= (*

*α*

1*E*1

*)*

^{0}

*+ (*

*α*

2*E*2

*)*

^{0}. (ii) Relaxed pressure consistency:

*E*^{0}^{,}^{C}

### = *α*

^{∗}

_{1}

*E*1

*p*^{∗}

*,* *(* *α*

1*ρ*

1*)*

^{0}

*α*

_{1}

^{∗}

### + *α*

_{2}

^{∗}

*E*2

*p*^{∗}

*,* *(* *α*

2*ρ*

2*)*

^{0}

*α*

^{∗}

_{2}

*,*

(16b)
where*E*^{0}^{,}^{C}=^{E}^{0}^{,}^{C}−^{(}^{ρ}^{u}^{}^{)}_{2}^{0}* _{ρ}^{·(}*0

^{ρ}

^{u}^{}

^{)}^{0}.

The ﬁrst property (i) means that the sum of the discrete values of the phasic total energies given by the solution of
the homogeneous system must recover discrete values of the mixture total energy that are consistent with a conservative
discrete form of(7). The second property (ii) means that the value of the relaxed (equilibrium) pressure p^{∗} predicted in
*the relaxation step must be equal to the pressure as computed through the mixture equation of state p(E*^{0}^{,}^{C}*,*

*α*

^{∗}

_{1}

*,*

*ρ*

^{∗}

_{1}

*,*

*ρ*

_{2}

^{∗}

*)*, deﬁned by(4). That is, with the SG EOS, this consistency condition reads (cf.(5))

*p*^{∗}

### =

^{E}0*,*C

*− ((* *α*

1*ρ*

1*)*

^{0}

*η*

1*+ (* *α*

2*ρ*

2*)*

^{0}

*η*

2*)* *− (*

^{α}_{γ}^{1}

^{∗}

_{1}

^{γ}_{−}

^{1}

^{π}_{1}

^{1}

### +

^{α}_{γ}^{∗}

^{2}

_{2}

^{γ}_{−}

^{2}

^{π}_{1}

^{2}

*)*

*α*_{1}^{∗}

*γ*1−^{1}

### +

_{γ}^{α}_{2}

_{−}

^{2}

^{∗}

_{1}

*.*

(17)
The mathematical formulation of the two-phase model with the phasic total energy equations (6)easily allows us to
satisfy both properties (i) and (ii). To ensure the property (i), it suﬃces in Step 1 to apply a standard conservative scheme
to the conservative portion of the energy equations(6e)and(6f), that is*∂**t**(*

*α*

*k*

*E*

_{k}*)+ ∇ · (*

*α*

*k*

*E*

_{k}*+*

^{u}*α*

*k*

*p*

_{k}

^{u}*), k*=

^{1}

*,*2, and to discretize symmetrically the non-conservative contribution

*Σ*appearing there. In such a way, the sum of the discrete non- conservative energy equations recovers a conservative discrete form of the mixture energy equation, as a consequence of the cancellation of non-conservative discrete contributions. The fulﬁllment of mixture total energy conservation consistency then easily enables us to ensure also the property (ii), the agreement of the relaxed equilibrium pressure with the correct mixture equation of state. See the simple pressure relaxation procedure for Step 2 in Section3.3.

On the other hand, it appears diﬃcult to obtain a mixture-energy-consistent scheme if we apply an analogous frac-
tional step algorithm to the classical six-equation two-phase model (1). Although clearly both formulations (1) and (6)
mathematically recover the conservation law for the mixture total energy, it seems hard to discretize the phasic internal
energies equations(1e)and(1f)in a way that recovers a conservative discrete form of(7). Indeed, numerical models such
as[9,17] built on the formulation (1)*need to augment the six-equation system with the equation for E. The additional*
*conservation law for E is solved through a standard conservative scheme to obtain consistent discrete values E*^{0}^{,}^{C}. These
values are then used to correct the thermodynamic state predicted by the non-conservative internal energy equations via
the mixture equation of state. Note that this approach in general does not guarantee the consistency property (ii), that is
*p*^{∗}=^{p}(E^{0}^{,}^{C}*,*

*α*

^{∗}

_{1}

*,*

*ρ*

^{∗}

_{1}

*,*

*ρ*

_{2}

^{∗}

*)*.

*3.2. Homogeneous hyperbolic system solution step*

In Step 1 of the algorithm we employ the wave propagation method of[19,39]to compute approximate solutions of the homogeneous system(14). This method belongs to a class of Godunov-type ﬁnite volume schemes[40,41,38,19]for solving hyperbolic systems of partial differential equations. We describe hereafter the basic ideas of the method in one dimension.

The two-dimensional scheme will be brieﬂy recalled in Section3.2.4.

*3.2.1. One-dimensional wave-propagation scheme*

We consider the solution of the one dimensional system*∂**t**q+ ∂**x**f(q)*+

*σ*

*(q, ∂*

*x*

*q)*=0 (as obtained by setting

*u*=

^{u and}*∇ = ∂**x* in(12)). We assume a uniform grid with cells of width*x, and we denote with Q*_{i}* ^{n}* the approximate solution at the

*ith cell at time t*

^{n}*, i*∈ Z

*∈ N*

^{, n}^{. Setting}

*t*=

^{t}

^{n}^{+}

^{1}−

^{t}*, the second-order one-dimensional wave-propagation scheme has the form*

^{n}*Q*_{i}^{n}^{+}^{1}

### =

^{Q}

_{i}

^{n}### −

*t*

*x*

*A*

^{+}

*Q*

_{i}_{−}

_{1}

_{/}_{2}

### +

*A*

^{−}

*Q*

_{i}_{+}

_{1}

_{/}_{2}

### −

*t*

*x*

*F*

_{i}^{h}

_{+}

_{1}

_{/}_{2}

### −

^{F}^{h}

_{i}_{−}

_{1}

_{/}_{2}

*,*

(18)
where*A*^{∓}*Q*_{i}_{+}_{1}_{/}_{2}*are the so-called ﬂuctuations at interfaces x*_{i}_{+}_{1}_{/}_{2}*between cells i and(i*+^{1}*), and F*_{i}^{h}_{+}_{1}_{/}_{2}are second-order
correction ﬂuxes for higher resolution.

Here the ﬂuctuations*A*^{∓}*Q*_{i}_{+}_{1}_{/}_{2}*are computed by solving local Riemann problems at cell interfaces x*_{i}_{+}_{1}_{/}_{2} for each pair
*of data Q*_{i}^{n}*, Q*_{i}^{n}_{+}_{1}corresponding to adjacent cells. A Riemann solver (cf.[41,38,19]) must be provided to perform this task.

Let us specify the jump relations that need to be satisﬁed by an approximate Riemann solver for a Riemann problem
*with left and right initial data q _{}and q*

*. The Riemann solution structure deﬁned by the solver can be expressed in general*

_{r}*by a set of M waves*

*W*

^{l}*and corresponding speeds s*

^{l}*, M*

*6. For example, see the HLLC-type solver (M*=

^{3) and the}

*Roe-type solver (M*=6) presented below. By using the formalism introduced in[42], we also deﬁne the f -waves

*Z*

*=*

^{l}

^{s}

^{l}*W*

*,*

^{l}*l*=

^{1}

*, . . . ,M, which have the dimension of a ﬂux. The sum of the waves must be equal to the initial jump in the vector q of*the system variables:

*q*

### ≡

^{q}*r*

### −

^{q}### =

*M*

*l*=

^{1}

*W*^{l}

*.*

(19)
Moreover, for any variable of the model system governed by a conservative equation the initial jump in the associated
*ﬂux function must be recovered by the sum of the associated f -wave components. In the considered model the conserved*
quantities are

*α*

*k*

*ρ*

*k*

*, k*=

^{1}

*,*2, and

*ρ*

*u, therefore in order to guarantee conservation we need:*

*f*

^{(ξ )}### ≡

^{f}^{(ξ )}*(*

*q*

*r*

*)* −

^{f}^{(ξ )}*(*

*q*

_{}*)* =

*M*

*l*=

^{1}

*s*^{l}*W _{ξ}*

^{l}### =

*M*

*l*=

^{1}

*Z*^{l}* _{ξ}* (20)

for *ξ*=^{2}*,*3*,4, where f ^{(ξ )}* is the

*ξth component of the ﬂux vector f , and*

*W*

^{l}

_{ξ}^{and}

*Z*

_{ξ}

^{l}^{denote the}

*ξ*th component of the

*lth wave and of the lth f -wave, l*=

^{1}

*, . . . ,M, respectively. It is clear that conservation of the partial densities ensures*conservation of the mixture density

*ρ*

=*α*

1*ρ*

1+*α*

2*ρ*

2. In addition, we must ensure conservation of the mixture total energy,
that is the consistency condition (i) in(16a):
*f*

_{E}### ≡

^{f}*E*

*(*

*q*

_{r}*)* −

^{f}*E*

*(*

*q*

_{}*)* =

*M*

*l*=

^{1}

*s*^{l}

*W*^{l}_{5}

### +

*W*

_{6}

^{l}### =

*M*

*l*=

^{1}

*Z*^{l}_{5}

### +

*Z*

_{6}

^{l}*,*

(21)
*where f**E*=*u(E*+

*α*

1*p*1+

*α*

2*p*2

*)is the ﬂux function associated to the mixture total energy E.*

Once the Riemann solution structure {W_{i}^{l}_{+}_{1}_{/}_{2}*,s*^{l}_{i}_{+}_{1}_{/}_{2}}*l*=1*,...,**M* *arising at each cell edge x*_{i}_{+}_{1}_{/}_{2} is deﬁned through a Rie-
mann solver, the ﬂuctuations*A*^{∓}*Q*_{i}_{+}_{1}_{/}_{2} *and the high-resolution correction ﬂuxes F*_{i}^{h}_{+}_{1}_{/}_{2} in(18)are computed as

*A*^{±}

*Q*

_{i}_{+}

_{1}

_{/}_{2}

### =

*M*

*l*=1

*s*

^{l}*i*+1*/*2

_{±}

*W*_{i}^{l}_{+}_{1}_{/}_{2}

*,*

(22)
*s*^{+}=^{max}*(s,*0*), s*^{−}=^{min}*(s,*0*)*, and
*F*^{h}_{i}_{+}_{1}_{/}_{2}

### =

^{1}

2

*M*

*l*=

^{1}

*s*

^{l}

_{i}_{+}

_{1}

_{/}_{2}

1

### −

*t*

*x*

*s*

^{l}

_{i}_{+}

_{1}

_{/}_{2}

*W*_{i}^{h}_{+}^{,}^{l}_{1}_{/}_{2}

*,*

(23)
where*W*_{i}^{h}_{+}^{,}^{l}_{1}_{/}_{2}are a modiﬁed version of*W*^{l}_{i}_{+}_{1}_{/}_{2}obtained by applying to*W*_{i}^{l}_{+}_{1}_{/}_{2}a limiter function (cf.[19]).

We present in the next sections two approximate Riemann solvers that we have developed for the model system(12):

a HLLC-type solver[43,44,38], similar to the solvers described in[9,17], and a new Roe-type solver[45]. The conservation consistency condition(21)(equivalently(16a)) is easily fulﬁlled in both cases.

*3.2.2. HLLC-type solver*

To begin with, we deﬁne an approximate solver for the two-phase system(14)by applying the idea of the HLLC solver
of Toro et al. [44] (see also[38]). One diﬃculty in designing a HLLC-type solver for the present model is related to the
non-conservative character of the phasic energy equations, for which we lack a notion of weak solution in the distributional
framework (see for instance on this subject[46–49]). Nonetheless, note that to correctly set the initial thermodynamic state
for the solution of the homogeneous system we only require that the sum of the phasic total energies computed at the pre-
vious time level fulﬁll the consistency condition*(*

*α*

1*E*

_{1}

*)*

^{0}

*+ (*

*α*

2*E*

_{2}

*)*

^{0}=

^{E}^{0}

^{,}^{C}. The individual phasic energies are re-set at the

*beginning of Step 1 by using the relaxed variables p*

^{∗}and

*α*

^{∗}

_{1}, obtained in Step 2 through a procedure that by construction ensures that the energies’ sum recovers the correct mixture energy state (Section3.3). The individual phasic energies values that come from the solution of the homogeneous system have only a role in the initial condition for the mechanical relax- ation step. Aiming at designing the simplest method that could provide reasonable mixture-energy-consistent estimates for