Gloria Faccanoni
email: gloria.faccanoni@univ-tln.fr

Bérénice Grec
email: berenice.grec@u-paris.fr

COMPARISON OF SEVERAL COMPLETE CUBIC LAWS FOR TWO-PHASE FLOW MODELS *

In the present paper, we investigate several cubic equations of state widely used in the literature, for which we are able to construct analytically the complete law. In order to describe twophase flows, we use Maxwell's construction, which amounts to consider pure phases and a mixture zone at saturation. The parameters appearing in the different equations of state are fitted in order to be precise in the saturation zone at high pressures. The different laws are then compared in a large range of pressures, showing the best accuracy of Clausius equation of state.

Introduction

In order to describe the dynamical evolution of a fluid, the physical models stem from the principles of conservation of mass, momentum and energy. To obtain a complete description, these conservation laws must be supplemented with constitutive relations characterizing the material properties of the fluid. The thermodynamic properties of a material are given in a relation called Equation of State (EoS). Thermodynamic considerations imposes mathematical constraints on the admissible EoS.

Many attempts have been made to describe the thermodynamic behavior of fluids, predicting their physical properties under given conditions. Therefore, many forms of equations of state have been presented in the literature. There are two main approaches: either find an analytical EoS, or use tabulated values to construct an EoS (for example the IAPWS relations [START_REF] Wagner | Tables of the Properties of Water and Steam[END_REF] for water). Of course, the main advantage of the second approach is that such EoS are capable of expressing fluid thermodynamic properties very precisely, but the non-analytically defined EoS require time-consuming computations. As far as the first approach is concerned, a huge literature exists, trying to find a good compromise between the simplicity of the EoS (related to the number of parameters it involves) and the precision of such EoS.

An (incomplete) EoS is a thermodynamic expression that relates pressure (p), temperature (T ), and specific volume (τ ). The easiest EoS is the ideal gas law pτ = rT , which is satisfactory for gases at low pressures. However, the application of the ideal gas law at higher pressures may lead to very big errors. Indeed, real gases behave differently than ideal gases, since the ideal gas law was derived under the assumption that the molecular volume tends to zero and neither molecular attraction nor repulsion exists between them.

Later, cubic equations of state (CEoS) have been introduced, taking into account both attraction and repulsion effects. These CEoS happen to give reasonable results for the thermodynamic behavior of real fluids. In 1873 Van der Waals [START_REF] Van Der Waals | Continuity of the gaseous and liquid state of matter[END_REF] made one of the earliest attempts to describe the behavior of real gases by an analytical CEoS p = rT /(τb) -a/τ 2 . The repulsion pressure is described by the term rT /(τb), and the attraction pressure by a/τ 2 . In this equation, a is a measure of the intermolecular attractive forces between the molecules, and b is known as the co-volume (linked to the volume of molecules). The values of both parameters can be obtained from the critical properties of the fluid. Despite its simplicity, the Van der Waals EoS correctly describe the behavior of fluids in both liquid and gaseous phases. However, it is not accurate enough to be useful for design purposes. Other researchers began attempts to improve Van der Waals EoS. A usual approach consists in adapting the molecular attraction term a/τ 2 . In 1880, Clausius [START_REF] Clausius | Über das Verhalten der Kohlensäure in Bezug auf Druck, Volumen und Temperatur[END_REF] suggested that the molecular attraction term was inversely proportional to temperature, leading to the relation p = rT /(τb) -a/(T (τ + δ) 2 ). The addition of a new constant δ enabled better agreement with data. However, the mathematical manipulations required for thermodynamic calculations were more difficult. So Berthelot, in 1899 [START_REF] Berthelot | Sur une méthode purement physique pour la détermination des poids moléculaires des gaz et des poids atomiques de leurs éléments[END_REF], removed this additional constant, assuming δ = 0. Another important modification of the Van der Waals EoS is the Redlich-Kwong one in 1949 [START_REF] Redlich | On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions[END_REF], where the attraction pressure term was replaced with a generalized temperature-dependent term p = rT /(τb) -a/( √ T τ (τ + b)). This idea led to many improvements. A particular important one in the literature is the one by Soave [START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF], who replaced the term (a/ √ T ) by a more general temperature-dependent term, denoted by aα(T ), with α(T ) = (1 + m(1 -√ T )) 2 . In this relation, T is the reduced temperature and the parameter m is correlated to the centric factor.

In this work, we shall only focus on some CEoS, presenting the main features of the improvements we mentioned with very few parameters. However, note that following these founding papers, thousands of variations of CEoS have been proposed in the literature, involving many other parameters, see e.g. overview papers [START_REF] Economou | Cubic and generalized Van der Waals equations of state[END_REF][START_REF] Ghoderao | A four parameter cubic equation of state with temperature dependent covolume parameter[END_REF][START_REF] Lopez-Echeverry | Peng-Robinson equation of state: 40 years through cubics[END_REF][START_REF] Wilczek-Vera | Understanding cubic equations of state: A search for the hidden clues of their success[END_REF].

As indicated before, equations of state are used to describe the behavior of fluids coupled with PDE models. From the thermodynamic considerations used to derive an EoS, it is classical to write it under the form (τ, T ) → p(τ, T ). However, in order to describe all thermodynamic variables, one needs to construct another relation, for instance the one giving the specific internal energy e under the form (τ, T ) → e(τ, T ) [START_REF] Faccanoni | Admissible Equations of State for Immiscible and Miscible Mixtures[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF]. Such a relation has to be derived following standard thermodynamic principles, and might involve some new modeling choices. This allows to define a so-called complete EoS. Note that, depending on the nature of the PDE, several choices of thermodynamic variables can be made. In Computational Fluid Dynamics (CFD), it is usual to use other thermodynamic functions than e(τ, T ), which requires to invert the complete EoS. In compressible models, the conservation laws lead to express p(τ, e), whereas in low Mach number models, where the pressure is almost constant, the EoS is expressed as T (τ, p) or in terms of enthalpy h(τ, p). In practice, despite the wide variety of existing CEoS, only few are used in CFD, due to the difficulty of the EoS completion. For more references on the mathematical analysis of CEoS and their use in PDE, we refer e.g. to [START_REF] Dahmen | Riemann Problem for the Euler Equation with Non-Convex Equation of State including Phase Transitions[END_REF][START_REF] Giovangigli | Supercritical fluid thermodynamics from equations of state[END_REF][START_REF] Müller | The Riemann Problem for the Euler Equations with Nonconvex and Nonsmooth Equation of State: Construction of Wave Curves[END_REF][START_REF] Quartapelle | Solution of the Riemann problem of classical gasdynamics[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF].

CEoS aim at describing correctly a fluid behavior in both liquid and gaseous phases, and could therefore be used as a single EoS in order to take into account two-phase flows and phase transition. However, this approach leads to non physical regions, since in the coexistence zone of the two phases, the pressure is increasing along an isotherme curve. In this work, in order to account for phase change phenomena, the EoS is modified thanks to Maxwell's construction [START_REF] Mortimer | Physical Chemistry[END_REF]. Note that, for the specific construction needed in the mixture zone (where both phases are present), other approaches based on one EoS per phase exist [START_REF] Faccanoni | Admissible Equations of State for Immiscible and Miscible Mixtures[END_REF]. Maxwell's construction allows, from an equation of state representing both the liquid and gaseous states of a pure body, to compute the so-called dome at saturation, where both phases are present, and the (constant) pressure at which the phase transition occurs at a given temperature. Observe that, although the CEoS are given under an analytical form, Maxwell's construction is not analytical, and the boundaries of the coexistence zone and the saturation pressure (or temperature) have to be approximated numerically.

As any analytical law, (incomplete) CEoS involve some parameters, which have to be fitted. As we already mentioned, we focus on simple CEoS, for which the number of parameters remain small (up to 4). Since no choice of parameters allows to be precise for any values of the thermodynamic variables, the fitting of the parameters is done in order for the EoS to be precise in some pressure range. Since a CEoS describe both liquid and gaseous phases, only one set of parameters is used for both phases, which allows Maxwell's construction of the coexistence zone. To this end, it is natural to use the critical point (at which both phases coexist) as a reference experimental value for the fitting [START_REF] Economou | Cubic and generalized Van der Waals equations of state[END_REF][START_REF] Mortimer | Physical Chemistry[END_REF].

By computing analytically the critical point of the CEoS, three relations between the parameters and the experimental critical values are obtained. Depending on the number of free parameters, the system can be under-or over-determined. In the literature, in the case of an over-determined system, the usual approach consists of dropping one of the relations (and thus not satisfying one of the three critical values). In this work, we chose to compare these three possible choices and the approach consisting in relaxing the value of the parameter r, which is usually considered fixed at the ideal gas constant. In the case of an under-determined system, one degree of freedom remains for the choice of the parameters. In both of these cases, we decide in this paper to be precise on Maxwell's construction at a fixed pressure, which allows to compute optimal values for the parameters. This work is done for water. Further, the quality of the parameters fitting is assessed by comparing the different CEoS with experimental data on a large range of pressures.

The paper is organized as follows. In section 2, we first introduce the different families of (incomplete) EoS on which we shall focus in the following. Section 3 is devoted to the construction of a complete EoS for all considered CEoS. Further, we tackle in section 4 the description of the mixture zone by means of Maxwell's construction. Section 5 focuses on the best possible choice of parameters for applications with a large range of pressures, and the comparison of the performances of our different CEoS. Last, in appendix, we give additional useful thermodynamic quantities for the different EoS: the speed of sound for compressible CFD, and the inversion h(τ, p) as well as its derivative for use in low Mach number models.

Cubic equations of state

Thermodynamic considerations

In this work, in order to describe accurately both liquid water and vapor, we need to take into account phase transition. To this end, let us start with some thermodynamic considerations.

Experimentally, it is proved that the temperature at which two phases can coexist at equilibrium is a smooth function of the pressure. In a phase diagram (T, p), such a function is represented by a so-called coexistence curve (see Figure 1a). The three coexistence curves (solid-gas, liquid-gas, solid-liquid) separate the (T, p) plane into regions where the fluid is a gas, a liquid or a solid (see Figure 1a). Those curves meet at a point called the triple point, which corresponds to the unique value of the pressure and temperature at which the three phases can coexist. The liquid-vapor coexistence curve terminates at a point which is called the critical point. The temperature, specific volume, and pressure at this point are called the critical temperature, denoted by T c , the critical specific volume, denoted by τ c , and the critical pressure, denoted by p c . For water, the values of these critical constants are

p c,exp = 22.064 × 10 6 Pa, T c,exp = 647.096 K, τ c,exp = 1/322 m 3 kg -1 . (1)
At temperatures higher than the critical temperature and pressures higher than the critical pressure there is no transition between liquid and gas phases.

In the plane (τ, p) (see Figure 1b), the saturation curve is defined as the one separating the regions where the fluid is a pure phase (liquid or gas) from the region where the two phases coexist. The set of states delimited by the saturation curve is called "dome" of saturation and these states are called "saturated states" or "at saturation" or "equilibrium mixing state". Again, for temperatures higher than the critical temperature, the pure component is under a single phase for all pressures: there is no longer phase transition.

Incomplete Cubic Equations of State

In order to describe the thermodynamic behavior of a fluid, we use an Equation of State (EoS) which links the three variables p, τ and T , usually devised in the litterature under the form p(τ, T ).

p(τ, T ) = rT τ -b - aα(T ) τ (τ + d) + c(τ -d) , (2)
where r, a, b, c and d are constants and α is a temperature-dependent function. The parameter b is usually interpreted as the volume occupied by the molecules, and r as the specific ideal gas constant. The second term in (2) can be thought as representing the effect of intermolecular forces, and is non-positive. Such EoS are called cubic since the pressure definition leads to a polynomial equation of order 3 on τ in the following form

τ 3 + c + d -b - rT p τ 2 -cd + (c + d) b + rT p - aα(T ) p τ + cbd + cd rT p -b aα(T ) p = 0. (3)
Among the many CEoS, some well-known families are included in the general form [START_REF] Berthelot | Sur une méthode purement physique pour la détermination des poids moléculaires des gaz et des poids atomiques de leurs éléments[END_REF]. The Patel-Teja (PT) family of EoS results from setting d = b, the Peng-Robinson (PR) one from setting c = d = b (see [START_REF] Lopez-Echeverry | Peng-Robinson equation of state: 40 years through cubics[END_REF] and the references therein).

The temperature-dependent function α must satisfy the following properties [START_REF] Le Guennec | A consistency test for α-functions of cubic equations of state[END_REF][START_REF] Le Guennec | On the imperative need to use a consistent α-function for the prediction of pure-compound supercritical properties with a cubic equation of state[END_REF][START_REF] Lopez-Echeverry | Peng-Robinson equation of state: 40 years through cubics[END_REF]:

α ∈ C 2 (R + ) and α(T ) ≥ 0, α ′ (T ) ≤ 0, α ′′ (T ) ≥ 0, α ′′′ (T ) ≤ 0. (4)
Again, many choices can be made for the function α, leading to different well-known families of CEoS.

In this paper, we do not intend to review the numerous cubic equations of state available in the literature. We focus here on some CEoS that are widely used in the literature and simple enough to give analytical formulas. These CEoS, though simple, capture the essential physics of liquid-vapor phase transition. The EoS that we consider are the following • the Van der Waals (VdW) EoS [START_REF] Mortimer | Physical Chemistry[END_REF], obtained by setting c = d = 0 and α(T ) = 1:

p(τ, T ) = rT τ -b - a τ 2 ; (5)
• the Berthelot (Ber) EoS [START_REF] Berthelot | Sur une méthode purement physique pour la détermination des poids moléculaires des gaz et des poids atomiques de leurs éléments[END_REF][START_REF] Mortimer | Physical Chemistry[END_REF], 1 obtained by setting c = d = 0 and α(T ) = 1 T :

p(τ, T ) = rT τ -b - a T τ 2 ; (6)
• the Clausius (Cl) EoS [START_REF] Clausius | Über das Verhalten der Kohlensäure in Bezug auf Druck, Volumen und Temperatur[END_REF], 1 obtained by setting c and d such that c 2 + d 2 + 6cd = 0 and α(T ) = 1 T :

p(τ, T ) = rT τ -b - a T (τ + δ) 2 , (7)
where we denoted δ def = c+d 2 ; • the Redlich-Kwong (RK) EoS [START_REF] Mortimer | Physical Chemistry[END_REF], 2 obtained by setting c = 0, d = b and α(T

) = 1 √ T : p(τ, T ) = rT τ -b - a √ T τ (τ + b) ; (8)
• the Soave-Redlich-Kwong (SRK) EoS [START_REF] Le Guennec | A consistency test for α-functions of cubic equations of state[END_REF][START_REF] Mansour | Equation of State[END_REF][START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF], 3 obtained by setting c = 0, d = b, and

α(T ) = (σ - √ T ) 2 : p(τ, T ) = rT τ -b - a(σ - √ T ) 2 τ (τ + b) . (9)
Note that, in order to satisfy (4), the parameter σ has to be larger than √ T . In our case, since we are interested in temperatures less than the critical one, it is sufficient that σ > √ T c , which is the case by definition of σ.

The critical point in a cubic equation of state

On Figure 2, we plot schematically the pressure of a fluid as a function of its specific volume for several temperatures (isotherm curves) with a Van der Waals EoS. Since, at fixed pressure and temperature, the equation of state (3) is a cubic equation in τ , we have the following behavior:

• for each isotherm curve above the critical temperature, the pressure is a monotone decreasing function of the volume; • for isotherm curves under the critical temperature, the pressure is not monotone. This defines the so-called spinodal zone where the pressure is increasing. On the left of the spinodal zone, the fluid is liquid, whereas it is gaseous on the right of this zone. For any fixed pressure smaller than p c , there is a temperature range for which this pressure is associated with three volumes; a gas phase and a liquid phase can coexist. For a CEoS, we can compute analytically the critical point. Indeed, the critical isotherm curve is the one possessing one inflection point (horizontal tangent and change of convexity). The critical point is this inflection point of the critical isotherm curve, it thus satisfies

p(τ c , T c ) = p c , ∂p ∂τ T (τ c , T c ) = 0 and ∂ 2 p ∂τ 2 T (τ c , T c ) = 0. ( 10
)
1 In the litterature, the function α is usually given of the form 1/(T /Tc). This α can be rewritten under the form we propose by correcting the coefficient a by a multiplicative factor Tc. 2 In the litterature, the function α is usually given of the form 1/ T /Tc. This α can be rewritten under the form we propose by correcting the coefficient a by a multiplicative factor √ Tc. 3 In the litterature, the function α is usually given of the form

[1 + m(ω)(1 -T /Tc)] 2 = m 2 Tc √ Tc m + √ Tc - √ T 2
, where ω denotes a physical constant (acentric factor), and the function m is a polynomial of order two with given coefficients. This α can be rewritten under the form we propose by setting σ = √ Tc(1 + 1/m(ω)) and correcting the coefficient a by a multiplicative factor m(ω) 2 /Tc. For the general CEoS (2), this leads to the following system

                   rT c τ c -b - aα(T c ) c(τ c -d) + τ c (τ c + d) = p c , - rT c (τ c -b) 2 + a (2τ c + c + d) α(T c ) (c(τ c -d) + τ c (τ c + d)) 2 = 0, 2rT c (τ c -b) 3 + 2aα(T c ) (c(τ c -d) + τ c (τ c + d)) 2 - 2a (2τ c + c + d) 2 α(T c ) (c(τ c -d) + τ c (τ c + d)) 3 = 0. (11)
At this point, a fluid exhibits some unusual properties such as infinite heat capacity or infinite compressibility. An alternative method to solve [START_REF] Gavrilyuk | Theoretical model of the Leidenfrost temperature[END_REF] has been introduced in [START_REF] Martin | Development of an equation of state for gases[END_REF]. Instead of defining p(τ, T ), one could also define implicitly τ (p, T ) from (3). Then, system [START_REF] Gavrilyuk | Theoretical model of the Leidenfrost temperature[END_REF] defining the critical point becomes

τ (p c , T c ) = τ c , ∂τ ∂p T (p c , T c ) = 0 and ∂τ 2 ∂ 2 p T (p c , T c ) = 0. ( 12
)
From this and the fact that (3) is a unitary polynomial, it follows that the function τ (p, T ) should be equal at the critical point to (τ -

τ c ) 3 = τ 3 -3τ c τ 2 + 3τ 3 c τ -τ 3 c
. Thus, it is enough to identify the coefficients of this polynomial with the ones in (3) evaluated at (p c , T c ) to obtain the system satisfied by the critical point

               c + d -b - rT c p c = -3τ c , cd + (c + d) b + rT c p c - aα(T c ) p c = -3τ 2 c , cbd + cd rT c p c -b aα(T c ) p c = -τ 3 c . (13)
3. Construction of a complete cubic EoS using the variables (τ, T )

The goal of a CEoS is to predict thermodynamic properties of a fluid, such as pressure, temperature and density, but also the enthalpy of vaporization, the isochoric heat capacity and the speed of sound. To do this, we have to construct a complete EoS from the incomplete EoS given as p = p(τ, T ).

In this section we focus on determining the specific internal energy e = e(τ, T ), which is called the caloric equation of state. Other quantities (e.g. speed of sound) are computed in Appendix.

Main idea

We search for the specific internal energy e as a function of temperature T and specific volume τ . Denoting s the entropy, we differentiate the Gibbs relation de = T dspdτ with respect to τ to obtain ∂e ∂τ T = T ∂s ∂τ T -p. Using one of the Maxwell relations [START_REF] Faccanoni | Admissible Equations of State for Immiscible and Miscible Mixtures[END_REF] leads to ∂s ∂τ T = ∂p ∂T τ

. Thus the partial derivative of e with respect to τ can be expressed as

∂e ∂τ T = T ∂p ∂T τ -p. (14)
The other partial derivative defines the specific heat capacity at constant volume c v :

∂e ∂T τ def = c v . (15)
Since e(τ, T ) must be an exact differential form, the two mixed second partial derivatives have to be equal to each other:

∂ ∂τ ∂e ∂T τ T = ∂ ∂T ∂e ∂τ T τ .
We can write this as

∂c v ∂τ T = ∂ ∂T T ∂p ∂T τ -p τ = ∂p ∂T τ + T ∂ 2 p ∂T 2 τ - ∂p ∂T τ ,
which leads to the following compatibility condition on c v

∂c v ∂τ T = T ∂ 2 p ∂T 2 τ . (16)
Thus, c v can be computed by integrating this compatibility relation ( 16)

c v (τ, T ) = c v,0 (T ) + τ τ0 T ∂ 2 p(ς, T ) ∂T 2 ς dς. (17)
Observe that, since we are working on an incomplete equation of state, we only have a constraint on the dependence of c v with respect to τ , and we have some freedom on the dependence of the integration constant c v,0 with respect to T . From the partial derivatives of e(τ, T ), we deduce that e(τ, T ) =

The point (τ 0 , T 0 ) with respect to which the integration takes place can be chosen freely. A standard approach in the literature [START_REF] Balmer | Modern engineering thermodynamics-textbook with tables booklet[END_REF] is to choose some arbitrary reference temperature and essentially zero pressure. In our work, we choose to work with the critical point (τ c , T c ).

Assumption. For the sake of simplicity, we choose c v,0 (T ) to be independent on T ,4 so that

c v,0 (T ) = c 0 . ( 19
)
The choice of the value of c 0 is discussed in Remark 5.2.

For each of the considered EoS, the previous construction allows to define e(τ, T ) explicitly, as showed in the next subsections.

Van der Waals caloric equation of state

For the Van der Waals equation of state [START_REF] Dellacherie | Accurate steam-water equation of state for two-phase flow LMNC model with phase transition[END_REF], we compute

∂ 2 p ∂T 2 τ (τ, T ) = 0. (20)
Therefore, the compatibility condition [START_REF] Lopez-Echeverry | Peng-Robinson equation of state: 40 years through cubics[END_REF] implies that c v (τ, T ) does not depend on τ , and is thus a function of T only. With our previous assumption [START_REF] Mortimer | Physical Chemistry[END_REF], this implies that

c v (τ, T ) = c 0 . (21)
Further, we can obtain the energy e(τ, T ) by integrating

e(τ, T ) = e c + c 0 (T -T c ) -a 1 τ - 1 τ c . (22)

Berthelot and Clausius caloric equation of state

Observe that Berthelot EoS can be obtained from the Clausius one by setting δ = 0. We thus proceed here with the general Clausius EoS. For the Clausius equation of state [START_REF] Faccanoni | A homogeneous relaxation low Mach number model[END_REF], we compute

∂ 2 p ∂T 2 τ (τ, T ) = - 2a T 3 (τ + δ) 2 . ( 23
)
The compatibility condition gives the dependence of c v (τ, T ) as a function of τ . We use again [START_REF] Mortimer | Physical Chemistry[END_REF], which leads to

c v (τ, T ) = c 0 + 2a T 2 1 τ + δ - 1 τ c + δ , (24)
since we choose to fit the constants using the critical point. Further, we obtain the energy e(τ, T ) by integrating

e(τ, T ) = e c + c 0 (T -T c ) - 2a T 1 τ + δ - 1 τ c + δ . (25)

Redlich-Kwong caloric equation of state

For the Redlich-Kwong equation of state [START_REF] Faccanoni | Modelling and Simulation of Liquid-Vapor Phase Transition in Compressible Flows Based on Thermodynamical Equilibrium[END_REF], we compute

∂ 2 p ∂T 2 τ (τ, T ) = - 3a 4T 2 √ T τ (τ + b) . ( 26
)
As before, we compute c v (τ, T ) from the compatibility condition and ( 19)

c v (τ, T ) = c 0 + 3a 4bT √ T ln τ + b τ τ c τ c + b . ( 27
)
This leads to the energy

e(τ, T ) = e c + c 0 (T -T c ) - 3a 2b 1 √ T ln τ + b τ τ c τ c + b . ( 28
)

Soave-Redlich-Kwong caloric equation of state

For the Soave-Redlich-Kwong equation of state ( 9), we compute

∂ 2 p ∂T 2 τ (τ, T ) = - aσ 2T √ T τ (τ + b) . (29)
Again, we first obtain c v (τ, T ) as

c v (τ, T ) = c 0 + aσ 2b √ T ln τ + b τ τ c τ c + b , (30)
and the energy

e(τ, T ) = e c + c 0 (T -T c ) - aσ b √ T -σ ln τ τ + b τ c + b τ c . (31)

General caloric equation of state

In this subsection, we perform the formal computations for any generic CEoS, without using the assumption [START_REF] Mortimer | Physical Chemistry[END_REF]. For the general CEoS (2), we compute

∂ 2 p ∂T 2 τ (τ, T ) = -a τ (τ + d) + c(τ -d) α ′′ (T ). ( 32
)
Let us define ∆ def = c 2 + d 2 + 6cd, and introduce 5 Therefore, the compatibility condition implies that

C(τ ) def = τ - a ς(ς + d) + c(ς -d) dς =          - a √ ∆ ln 2τ + c + d - √ ∆ 2τ + c + d + √ ∆ if ∆ > 0, a τ + c+d 2 if ∆ = 0.
c v (τ, T ) = c v,0 (T ) + τ τc T ∂ 2 p(ς, T ) ∂T 2 ς dς = c v,0 (T ) + T α ′′ (T ) [C(τ )] τ τc . (34)
Further, we can compute the energy

e(τ, T ) = e c + T Tc c v,0 (θ)dθ + τ τc T ∂p(ς, T ) ∂T τ -p(ς, T ) dς = e c + T Tc c v,0 (θ)dθ + (T α ′ (T ) -α(T )) τ τc -a ς(ς + d) + c(ς -d) dς = e c + T Tc c v,0 (θ)dθ + (T α ′ (T ) -α(T )) [C(τ )] τ τc . ( 35
)
4. Phase Equilibrium

Maxwell's construction in incomplete EoS

Under the critical temperature, the spinodal zone (see Figure 2), situated inside the saturation dome, corresponds to a non physical part, since the pressure increases for increasing volume. Moreover, the part of the saturation dome outside the spinodal zone corresponds to the metastable states. Thus, cubic EoS can only describe in a satisfactory manner the properties of the fluid at the transition between liquid and vapor above the critical point or outside the saturation dome. Inside this saturation dome, the equation of state has to be modified in order to take into account phase change phenomena. To this end, we introduce a saturation plateau at a pressure depending only on the temperature. In thermodynamics, it is classical to use the so-called Maxwell's Area Rule, or Maxwell's construction, which allows to compute the saturation pressure at which the phase transition occurs at a given temperature. From a graphical point of view, this rule is equivalent to saying that the saturating vapor pressure of the isotherm is located on the level for which the two (positive) areas of the two blue surfaces of Figure 3a are equal.

Let us now give more details on the computations of Maxwell's construction. For fixed p 0 < p c , we have to compute τ s l , τ s g and T s such that the three following relations hold

                 p(τ s l , T s ) = p 0 , p(τ s g , T s ) = p 0 , τ s g τ s l (p(τ, T s ) -p 0 ) dτ = 0. (36)
Solving this system for any p < p c allows to define the following functions

p → τ s l (p), τ s g (p), T s (p) . (37)
Equivalently, for any fixed T 0 < T c ,we can look for τ s l , τ s g and p s such that p(τ s l , T 0 ) = p s , p(τ s g , T 0 ) = p s and τ s g τ s l (p(τ, T 0 ) -p s ) dτ = 0 (this is what is commonly done in the standard thermodynamics literature). Solving this system for any T < T c allows to define the following functions For the different EoS introduced in Section 2.2, we can compute explicitly the integral to rewrite system (36) in the following way: find (τ s l , τ s g , T s ) such that

T → τ s l (T ), τ s g (T ), p s (T ) . (38)
                   rT s τ s l -b - aα(T s ) τ s l (τ s l + d) + c(τ s l -d) = p 0 , rT s τ s g -b - aα(T s ) τ s g (τ s g + d) + cτ s g (τ s g -d) = p 0 , rT s ln τ s g -b τ s l -b + α(T s )(C(τ s g ) -C(τ s l )) -p 0 (τ s g -τ s l ) = 0.
where C(τ ) has been defined in (33). The functions C and α are (re)given here for the sake of completeness.

• Let us now explain how the complete EoS is constructed in the mixture zone. For (τ, T ) given such that τ s l (T ) < τ < τ gs (T ), we define the mass fraction of the mixture

φ(τ, T ) def = τ -τ s l (T ) τ s g (T ) -τ s l (T ) ∈ (0, 1).
The energy of the mixture at saturation is then given by e(τ, T

) def = φ(τ, T )e s l (T ) + (1 -φ(τ, T ))e s g (T ),
where we define the functions at saturation e s l/g (T ) def = e(τ s l/g (T ), T ) using ( 38) in (35). The complete EoS e Maxwell (τ, T ) is then defined piecewise using the saturation values:

e Maxwell (τ, T ) = φ(τ, T )e s l (T ) + (1 -φ(τ, T
))e s g (T ) in the mixture, i.e. τ s l (T ) < τ < τ s g (T ), Eq. ( 35) otherwise.

This construction ensures that the speed of sound remains positive [START_REF] Faccanoni | Modelling and Simulation of Liquid-Vapor Phase Transition in Compressible Flows Based on Thermodynamical Equilibrium[END_REF]. For the computation of the speed of sound for the different EoS, see Appendix A.

Remark 4.1 (Specific isochoric heat capacity). By definition of c v (15), the heat capacity of the mixture at saturation is given by

c v (τ, T ) = ∂φ ∂T τ (e s l (T ) -e s g (T )) + φ(τ, T )(e s l ) ′ (T ) + (1 -φ(τ, T ))(e s g ) ′ (T )
where

∂φ ∂T τ (τ, T ) = φ(τ, T )(τ s g ) ′ (T ) + (1 -φ(τ, T ))(τ s l ) ′ (T ) τ s g (T ) -τ s l (T ) , (e s κ ) ′ (T ) = ∂e ∂T τ (τ s κ (T ), T ) + (τ s κ ) ′ (T ) × ∂e ∂τ T (τ s κ (T ), T ) = c v (τ s κ (T ), T ) + (τ s κ ) ′ (T ) × (T × (p s ) ′ (T ) -p s (T )), (p s ) ′ (T ) = 1 T e s g (T ) -e s l (T ) τ s g (T ) -τ s l (T ) + p s (T ) T (Clapeyron relation).
Analytical expressions for the phase boundaries τ s l/g (T ) or their derivatives are not available, but they can be computed numerically.

Methodology to determine the parameters

Let us now describe shortly the methodology to determine the parameters of a CEoS for a given fluid. We proceed in two steps:

(1) First, by solving system (13), we establish a relation between the parameters (a, b, and possibly δ or σ), the experimental critical values (T c,exp , τ c,exp , p c,exp ) and r. A specific choice has to be made for the value of r: either treating it as a parameter, or using the experimental value r exp def = R/M , where R is the ideal gas constant and M the molar mass of the fluid. For water, r exp = 461.526 J K -1 kg -1 .

Case EoS

T c p c τ c r a b ① Van der Waals Berthelot T c,exp p c,exp 3rT c 8p c r exp a VdW = 27T 2 c r 2 64p c , a Ber = a VdW T c rT c 8p c Redlich-Kwong rT c 3p c 1 + 3 √ 2 + 3 √ 4 T 2 c √ T c r 2 9p c 3 √ 2 -1 rT c 3p c ② Van der Waals Berthelot T c,exp 3rT c 8τ c τ c,exp r exp a VdW = 9rT c τ c 8 , a Ber = a VdW T c τ c 3 Redlich-Kwong rT c 3τ c 1 + 3 √ 2 + 3 √ 4 rT 2 c √ T c τ c 3 3 √ 2 -1 τ c ③ Van der Waals Berthelot 8p c τ c 3r p c,exp τ c,exp r exp a VdW = 3p c τ 2 c , a Ber = a VdW T c τ c 3 Redlich-Kwong 3p c τ c r p c τ 2 c √ 3p c τ c 3 √ 2 -1 r 3 √ 2 -1 τ c ④ Van der Waals Berthelot T c,exp p c,exp τ c,exp 8p c τ c 3T c a VdW = 3p c τ 2 c , a Ber = a VdW T c τ c 3 Redlich-Kwong 3p c τ c T c 1 + 3 √ 2 + 3 √ 4 p c √ T c τ 2 c 3 √ 2 -1 τ c
Table 1. All possible choices of fixed parameters for Van der Waals, Berthelot and Redlich-Kwong EoS.

(2) When the problem is under-or over-determined, we optimize the choices of parameters of the parameter values in order to best fit the values at saturation at a given pressure p = p * with IAPWS data [START_REF] Wagner | Tables of the Properties of Water and Steam[END_REF]. After fixing these parameters, we compare the saturation temperature T s (p) (coexistence curves) and the phase boundaries τ s l/g (p) with respect to experimental data for different pressures to assess the quality of the EoS.

First step: expression of the parameters from the critical values

The critical point (T c , τ c , p c ) is obtained by solving the system (13).

Van der Waals, Berthelot and Redlich-Kwong EoS

For Van der Waals, Berthelot and Redlich-Kwong EoS, the problem is overdetermined. We explore all possible choices of fixed parameters.

① Usually, in the literature, T c = T c,exp , p c = p c,exp and r = r exp are fixed, to deduce a, b, τ c ̸ = τ c,exp .

Historically, this approach was chosen since the experimental value τ c,exp was not known [START_REF] Wilczek-Vera | Understanding cubic equations of state: A search for the hidden clues of their success[END_REF]. The obtained values are summarized in Table 1.

②
Remark 5.1. We treat the parameter r in the equation of state as an empirical one, different in general from the experimental value. We shall observe in section 5.2 that the experimental value of r does not lead to the best results when evaluating the saturation constants, and we shall therefore choose another optimized value.

In our setting, for high pressures, the value of r differs from the ideal gas constant. Thus, the CEoS does not converge to the ideal gas law when the pressure converges to zero. Note that, if one wants to be precise in another pressure regime, the parameters have to be fitted again. In particular, for low pressure values, the optimized value of r is close to the ideal gas constant. Indeed, this same approach has also been used independently from our work in [START_REF] Gavrilyuk | Theoretical model of the Leidenfrost temperature[END_REF] in the case of the Van der Waals EoS. In their work, the pressure is set at the atmospheric value, which leads to a value of r close to the experimental one.

Clausius EoS

For the Clausius EoS, we now have three parameters (a, b and δ). Thus, all three critical values can be fixed to the experimental values. We can also fix the value of r (see section 5.2.2), and a, b and δ are then recovered by the following relations

a(r) = 27T 3 c,exp r 2 64p c,exp , b(r) = τ c,exp - rT c,exp 4p c,exp , δ(r) = 3rT c,exp 8p c,exp -τ c,exp .

Soave-Redlich-Kwong EoS

For the Soave-Redlich-Kwong (SRK) EoS, we again have three parameters (a, b and σ). However, we observe that the first equation of ( 13) does not depend on the choice of the function α. Moreover, in this case, since c = 0 and d = b, no parameter (a, b and σ) appears in this equation. Thus, the system is overdetermined when fixing all three critical quantities and r to their experimental values. Guided by the best choice for Van der Waals, Berthelot and Redlich-Kwong EoS (cf. Section 5.2.1), we again allow r to vary, fixing the three critical quantities to their experimental values. When solving the system, we obtain

a(σ) = 1 + 3 √ 2 + 3 √ 4 p c,exp τ 2 c,exp (T c,exp -σ) 2 , b = 3 √ 2 -1 τ c,exp , r = 3p c,exp τ c,exp T c,exp ̸ = r exp .
Observe that the parameter σ remains free, and can thus be chosen in order to optimize some quantities of interest. We will explain in the following the choice we made for σ.

Reduced EoS

A usual approach to avoid computing the values of the parameters for different fluids is to consider the reduced EoS, for which the variables are the reduced temperature T = T /T c , reduced pressure p = p/p c and reduced specific volume τ = τ /τ c .

For the Van der Waals, the Berthelot and the Redlich-Kwong EoS, let us emphasize that all the different choices of fixed parameters as in Table 1 lead to the same reduced EoS:

p(τ , T ) =                  8 3 T τ -1/3 - 3 τ 2
for the Van der Waals EoS,

8 3 T τ -1/3 - 3 T τ 2 for the Berthelot EoS, 3T τ -3 √ 2 -1 - 1 + 3 √ 2 + 3 √ 4 √ T τ (τ + 3 √ 2 -1 )
for the Redlich-Kwong EoS, which do not involve any parameter. More precisely, this reduced form suggests to treat the parameter r as a variable one, in the same way as for a and b. Indeed, if the experimental values of the three critical quantities are known, the reduced EoS can be used as such. Of course, if one of these values is not known (corresponding to the cases ①-③), it is recovered by the corresponding formula in Table 1, involving then the (known experimental) value of r.

For Clausius EoS, we observe that, if we define the quantity ξ = 8p c,exp τ c,exp /T c,exp , the reduced Clausius EoS only involves this quantity ξ:

p = 8T 3ξτ + 2 -3ξ - 3 T (ξτ + 1 -ξ) 2 .
For Soave-Redlich-Kwong EoS, observe that, if we define the "reduced parameter" σ = σ/ √ T c , the reduced SRK EoS only involves this parameter σ:

p = 3T τ -3 √ 2 -1 - 1 + 3 √ 2 + 3 √ 4 σ -1 (σ - √ T ) 2 τ (τ + 3 √ 2 -1 )
.

Observe that the quantity ξ and the "reduced" parameter σ appearing in the Clausius and Soave-Redlich-Kwong reduced EoS do not play the same role. Indeed, for a fixed fluid, the critical values are given, and the constant ξ has a fixed value, whereas the reduced parameter σ involves the constant σ, and can therefore be chosen independently.

Second step: saturation at fixed pressure

Is it widely admitted in the literature [START_REF] Kontogeorgis | Taking Another Look at the Van der Waals Equation of State-Almost 150 Years Later[END_REF][START_REF] Wilczek-Vera | Understanding cubic equations of state: A search for the hidden clues of their success[END_REF] that the parameters have to be chosen in view of the application of the EoS. The novelty of our approach is to allow r to be a free parameter (which can be interpreted as a fraction of the experimental value of r). In our case, in order to handle diphasic flows, we want to be precise near to the saturation state around a fixed reference pressure. We here present the full procedure of parameter fitting, with the comparisons of all choices of fixed experimental values, for water at a fixed pressure p * = 155 bar. Of course, the same procedure can be applied for a different fluid, at a different pressure, or even using other reference values than the saturation values.

Van der Waals, Berthelot and Redlich-Kwong EoS

In Table 2 we present, for the Van der Waals, the Berthelot and the Redlich-Kwong EoS, the four possible choices of fixed experimental values described in the previous subsection (indicated in the first column), and show the values of the obtained parameters. Moreover, in the three last columns, we compute, using these values of the parameters, the saturation values by solving system (36). Note that, for the sake of readibility, we chose to display the density ϱ = 1/τ instead of the specific volume τ . We assess the quality of the EoS by the error made on these saturation values with respect to experimental IAPWS values [START_REF] Wagner | Tables of the Properties of Water and Steam[END_REF]. We observe that, for all three EoS, the best choice is to fix all three critical values (T c,exp , τ c,exp , p c,exp ), and to compute a, b and r. This choice leads to more precise saturation values, as can be seen on Figure 4, which depicts the temperature for these EoS in the plane (ϱ, p). Moreover, we observe that the Berthelot EoS leads to a far better description of the fluid at this pressure p * than the Van der Waals one: the saturation temperature is very precise, and the phase boundaries are closer to the experimental ones. As for the Redlich-Kwong EoS, it gives an even better approximation for the phase boundaries than the Berthelot EoS, but the saturation temperature slightly looses precision.

Fixed

param. a b r p c T c ϱ c ϱ s l (p * ) ϱ s g (p * ) T s (p * ) [m 3 J kg -2 ] [m 3 kg -1 ] [J K -1 kg -1 ] [Pa] [K] [kg m -3 ] [kg m -3 ] [kg m -3 ] [K] IAPWS 461

.2. Clausius EoS

As far as the Clausius EoS is concerned, we already mentioned that the three critical values can be fixed, as well as the value of r. A first choice for r is the experimental one given by the IAPWS. However, we observe in Table 3 and on Figure 5a that this choice leads to phase boundaries which are quite far from the experimental ones. To improve the parameter fitting, similarly to the Berthelot case for which it led to good results, we consider the possibility to relax the experimental value of r. In this setting, we have to choose a value for r. When choosing the same value of r as for the Berthelot EoS, we see in Table 3 that there still is some discrepancy on the phase boundaries. To improve the choice of r, we compare on Figure 5b the phase boundaries ϱ s l and ϱ s g for all values of r between the experimental value and the one obtained in the Berthelot EoS (we do not compare the saturation temperature since it does not vary for different values of r). We observe that, when trying to decrease the error on ϱ s g , the error on ϱ s l increases rapidly. Therefore, we can define an optimal value of r, corresponding to minimizing simultaneously both the normalized L 1 error on ϱ s g and the one on ϱ s l . Finally, we observe that the precision obtained with this choice of parameteres in the Clausius EoS is better than all other choices.

a b δ r ϱ s l (p * ) ϱ s g (p * ) T s (p * ) [m 3 J K kg -2 ] [m 3 kg -1 ] [m 3 kg -1 ] [J K -1 kg -1 ] [kg m -3 ] [kg m -3 ] [K] IAPWS 461

Soave-Redlich-Kwong EoS

For the Soave-Redlich-Kwong EoS, we already mentioned that the parameter σ remains free. As for the Clausius case, the idea is to optimize it in order to improve the saturation values (in this case, since the parameter σ is related to the temperature, it affects all three saturation values). We observe on Figure 6 that the optimal values of σ allowing to be precise on the phase boundaries or on the saturation temperature are very different. Thus, a choice has to be made on which saturation quantity needs the most precision. In our case, with applications in CFD in mind, it has been shown [START_REF] Dellacherie | Accurate steam-water equation of state for two-phase flow LMNC model with phase transition[END_REF] that the precise evaluation of the saturation temperature is of utmost importance. We therefore chose the optimal value σ opt = 48. For this choice, the parameters and saturation values are given in Table 4.

a b r ϱ s l (p * ) ϱ s g (p * ) T s (p * ) [m 3 J K kg -2 ] [m 3 kg -1 ] [J K -1 kg -1 ] [kg m -3 ] [kg m -3

All EoS

To summarize the performance of all EoS, we compare the temperature around saturation at pressure p * for the best choice of parameters in each EoS on Figure 7. We note that whereas the Van der Walls, and to a lesser extent, the Redlich-Kwong EoS are not very precise, the three other EoS are very accurate on the saturation temperature. The phase boundaries are best described by the Clausius and Soave-Redlich-Kwong EoS, with slightly more precision with the Clausius one.

Precision of CEoS with fitted parameters at all pressures

Fixing now the parameters at the optimal values obtained in the previous subsection for p = p * , we want to assess the quality of the different CEoS by comparing them with experimental data [START_REF] Wagner | Tables of the Properties of Water and Steam[END_REF] at all pressures.

Coexistence curve

More precisely, we plot on Figure 8 the coexistence curves p → T s (p) for variable pressures, starting from the smallest pressure in tabulated values up to the critical pressure p c , i.e. p ∈ [7 × 10 2 Pa, 22.06 × 106 Pa]. This is done for all four considered CEoS together with experimental values 6 . We observe a large discrepancy of the Van der Waals EoS with respect to IAPWS for pressures far from p c . On the other hand, the Berthelot, Clausius and Soave-Redlich-Kwong EoS show a very good accuracy with respect to IAPWS, even for pressures far from p * (at which the parameters were fitted). We also observe that there is no visible influence of the parameter δ, and both Berthelot and Clausius EoS are very satisfying. For the Soave-Redlich-Kwong EoS, the choice we made of optimal σ in order to lead maximal precision on the saturation temperature at p * leads to a very good accuracy on the whole pressure interval. As for the Redlich-Kwong EoS, it does not have the precision of SRK/Clausius EoS for the coexistence curve, but it still gives a satisfying approximation of experimental data on the whole pressure interval.

Phase boundaries

We further investigate the quality of the parameter fitting for all p ∈ [7×10 2 Pa, 22.06×10 6 Pa] by comparing the saturation zones. As before, for the sake of readibility, we chose to work in the plane (ϱ, p). In this plane, the mixture is at saturation in the set (ϱ, p) ϱ s g (p) ≤ ϱ ≤ ϱ s l (p) , the fluid is in the liquid phase for ϱ > ϱ s l (p) and in the gas phase for ϱ < ϱ s g (p). We can thus compare the phase boundaries ϱ s l/g (p) for each EoS with respect to experimental data as plotted on Figure 9. More precisely, for each pressure, we determine the saturation values, which allow to construct the solid lines describing the so-called saturation dome. Again, we observe a better accuracy for Berthelot, Redlich-Kwong, Clausius and SRK EoS than for the Van der Waals one. We also notice that in the lower range of pressures, the Redlich-Kwong/SRK EoS seem to be more adequate that the Berthelot/Clausius ones, when the parameters have been fitted at a higher pressure. We also plot on Figure 9 in dashed lines the isotherm curve ϱ → p(ϱ, T s (p * )) for each EoS. Since the value of p * = 155 bar is fixed, the height of the plateau of the curve is the same for any EoS. However, the value of T s (p * ) is different for each , by solving (36) with a simple fixed point procedure to obtain the saturation values. For Van der Waals and Berthelot EoS, a more refined procedure would be needed to compute the saturation, since the system is very stiff.

Internal energy

We can also observe the performance of the cubic laws by plotting the internal energy with respect to τ for different temperatures. To this end, we chose for the integration constants the experimental values given in [START_REF] Wagner | Tables of the Properties of Water and Steam[END_REF], namely e c = 2.019 × 10 6 J kg -1 , and c 0 ≈ 4.538 × 10 3 J kg -1 K -1 . We observed on Figure 7 that the best cubic laws were the Clausius and the Soave-Redlich-Kwong ones for evaluating the saturation temperature, but that the Clausius law was far better for evaluating the boundaries of the saturation zone. We thus chose to plot the internal energy for the Clausius EoS for different isotherm curves (see Figure 10).

Since the parameters have been fitted for the critical temperature, the energy is of course closer to experimental values for temperatures close to the critical one.

We observe that the error made on the saturation zone is inherent to the EoS. Since the evaluation of τ s l is quite precise, the liquid zone is almost exact, whereas the error on τ s g leads to a worse evaluation of the vapor zone.

It is worth noting that, on an isotherm curve, the error on the energy could be corrected by dropping assumption [START_REF] Mortimer | Physical Chemistry[END_REF], and choosing c v,0 (T ) depending on T . In particular, by defining c v,0 (T ) piecewise in the liquid and in the vapor phase, the two parts of the isotherm curves could be independently translated vertically, in order to improve the accuracy of the internal energy.

Remark 5.2. For any CEoS, when considering the monotony of c v with respect to τ on an isotherm curve, one notices that it is related to the sign of

- a τ (τ + d) + c(τ -d) ,
which is always non positive (on the whole physical domain for values of τ , and for any physical values of the parameters), leading to c v being non increasing. However, this is not the behavior observed for some fluids (for example, for water, c v is increasing in the liquid phase whereas it is decreasing in the vapor phase). Nevertheless, this problem does not occur on e, since it has been constructed to verify the correct properties, and the internal energy is thus always increasing. Note also that, with our assumption [START_REF] Mortimer | Physical Chemistry[END_REF], c v (τ, T ) given by [START_REF] Mansour | Equation of State[END_REF] does not diverge when approaching the critical point (τ c , T c , p c ), in contrast with experimental observations. Again, this problem does not appear on the internal energy, since c 0 is integrated between T c and T , which converges to zero at the critical point.

Conclusion

In this paper, we have studied several CEoS of different types (Van der Waals, Berthelot and Clausius, Redlich-Kwong, Soave-Redlich-Kwong), for which we are able to construct analytically the complete EoS. Moreover, we optimized the parameters appearing in the EoS in order to be precise on the coexistence zone of two phases (i.e. on Maxwell's construction). We investigated the performances of these EoS, and showed that, for a large range of pressures, Clausius (and to a lesser extent Soave-Redlich-Kwong) EoS lead to very satisfying results, improving greatly the accuracy with respect to Van der Waals EoS.

Following this work, some improvements could be investigated. In particular, it could be of interest to investigate possible corrections of the considered CEoS when changing the function α(T ), e.g. when choosing α(T ) = T ν , ν < 0. Moreover, we mentioned that the work could be extended by allowing c v,0 to depend on T . In this case, a choice has to be made for this function c v,0 (T ), which could be optimized to improve the performance of the considered complete EoS. Further, in the context of use in compressible CFD, the remaining step would be to develop an efficient strategy for inverting the complete EoS.

A. Speed of sound

In compressible flow models, knowledge of the speed of sound is fundamental. It is defined as

(c * ) 2 (τ, s) def = -τ 2 ∂p ∂τ s .
Let p(τ, s) = p(τ, T (τ, s)). By the composition rule, we have

(c * ) 2 (τ, T ) = -τ 2 ∂p ∂τ T + ∂p ∂T τ ∂T ∂τ s .
Using a Maxwell relation [START_REF] Faccanoni | Admissible Equations of State for Immiscible and Miscible Mixtures[END_REF] and the chain rule, we compute

(c * ) 2 (τ, T ) = -τ 2 ∂p ∂τ T - ∂p ∂T τ ∂p ∂s τ = -τ 2 ∂p ∂τ T - ∂p ∂T τ ∂p ∂T τ ∂T ∂e τ ∂e ∂s τ = -τ 2 ∂p ∂τ T - ∂p ∂T τ 2 ∂T ∂e τ T ,
where we used the definition of T def = ∂e ∂s τ in the last equality. We can rewrite this expression of (c * ) 2 using

c v (τ, T ) def = ∂e ∂T τ as (c * ) 2 (τ, T ) = -τ 2      ∂p ∂τ T -T ∂p ∂T τ 2 ∂e ∂T τ      = -τ 2 ∂p ∂τ T - ∂p ∂T τ 2 T c v (τ, T ) . ( 39
)
Note that, with Maxwell's construction, we never handle the spinodal zone, and thus, we only consider the EoS in the zone where ∂p ∂τ T is nonpositive. This guarantees that the speed of sound is positive.

In pure phases

In formula (39), all quantities are given by the EoS. More precisely, c v (τ, T ) is given for Van der Waals, the Berthelot/Clausius, the Redlich-Kwong and the Soave-Redlich-Kwong EoS respectively by ( 21), ( 24), ( 27) and (30), and the partial derivatives are given in Table 5 for each EoS.

In the mixture at saturation

In the mixture, ∂p ∂τ T = 0, since p(τ, T ) = p s (T ) in the mixture, and ∂p ∂T τ = (p s ) ′ (T ). The functions (p s ) ′ (T ) and c v (τ, T ) in the mixture are given in section 4.2, which defines all quantities in (39). When we are interested in simulating low Mach number flows, for which the thermodynamical pressure is constant, it is natural to work with the variables (τ, p) instead of (τ, T ). For some equations of state, it is possible to invert the formula of the pressure and express the temperature as a function of (τ, p). In this case, we can insert this expression of T in e(τ, T ) to gain ẽ(τ, p) def = e(τ, T (τ, p)),

which further leads to h(τ, p) = pτ + ẽ(τ, p). In the LMNC model [START_REF] Faccanoni | A homogeneous relaxation low Mach number model[END_REF], an important variable is the derivative of the enthalpy with respect to τ . We can derive h with respect to τ to obtain:

ζ(

B.1. Generic

We now consider the general cubic equation of state [START_REF] Berthelot | Sur une méthode purement physique pour la détermination des poids moléculaires des gaz et des poids atomiques de leurs éléments[END_REF]. Provided that the inversion of the temperature is possible, giving thus T (τ, p), we can do the same construction as before. In order to obtain ζ(τ, p) by (41), we need to compute the two partial derivatives of e given by (35):

B.2. Van der Waals

From the expression (5) of the pressure, we can invert the temperature as follows

T (τ, p) = (pτ 2 + a) (τ -b) rτ 2 .
Inserting this expression in relation ( 22) for e, we deduce ẽ(τ, p) and thus the enthalpy and its derivative

B.5. Redlich-Kwong

From the expression (8) of the pressure, we observe that the inversion of the temperature reduces to a third-order equation on

√ T : r τ -b √ T 3 -p √ T - a τ (τ + b) = 0.
Selecting the real root in Cardano's formula, we obtain

T (τ, p) =   3 2 3 p F + F 3 √ 18r   2 ,
where we defined We can thus compute the partial derivative

F def = 3 √ 3
∂T ∂τ p = 2 T (τ, p)   - 3 2 3 pF ′ F 2 + rF ′ -r ′ F 3 √ 18r 2   ,
with the following derivatives with respect to τ

F ′ = F 2 √ 3(9aa ′

Figure 1 .

 1 Figure 1. Schematic coexistence and saturation curves for water

Figure 2 .

 2 Figure 2. Isotherm curves for a cubic law in a reduced pressure-volume diagram (i.e. (τ /τ c , p/p c ) plane). The critical point is marked in red. The spinodal zone where the pressure is increasing (along an isotherm) is delimited by the gray curve.

  τ τ0 ∂e ∂τ T (ς, T )dς + T T0 ∂e ∂T τ (τ 0 , θ)dθ + e(τ 0 , T 0 ) and thus, using the expression of c v (τ, T ), we obtain e(τ, T ) = τ τ0 T ∂p(ς, T ) ∂T τ -p(ς, T ) dς + T T0 c v,0 (θ)dθ + e 0 .

5

  The Clausius EoS corresponds to ∆ = c 2 + d 2 + 6cd = 0 with δ = c+d 2 , and the (Soave)-Redlich-Kwong families correspond to ∆ = b 2 > 0.

Figure 3 .

 3 Figure 3. Isotherm curves for a cubic law in a pressure-volume plane (τ, P ) for the Maxwell construction

) 2 4. 2 .

 22 Van der Waals EoS: C(τ ) = 1/τ and α(T ) = 1 • Berthelot EoS: C(τ ) = 1/τ and α(T ) = 1/T • Clausius EoS: C(τ ) = 1/(τ + δ) and α(T ) = 1/T • Redlich-Kwong EoS: C(τ ) = (ln(τ + b) -ln(τ ))/b and α(T ) = 1/ √ T • Soave-Redlich-Kwong EoS: C(τ ) = (ln(τ + b) -ln(τ ))/b and α(T ) = (σ -√ T Phase transition in complete EoS in (τ, T )

  We fix τ c = τ c,exp , T c = T c,exp , r = r exp and deduce a, b, and p c ̸ = p c,exp . ③ We fix τ c = τ c,exp , p c = p c,exp , r = r exp and deduce a, b, and T c ̸ = T c,exp . ④ We fix T c = T c,exp , p c = p c,exp , τ = τ c,exp and deduce a, b, and r ̸ = r exp .

Figure 4 .

 4 Figure 4. Temperature T (ϱ, p * ) as a function of ϱ for all different choices of fixed parameters for Van der Waals, Berthelot and Redlich-Kwong EoS

( a )Figure 5 .

 a5 Figure 5. Comparison of different values of r for the Clausius EoS

Figure 6 .

 6 Figure 6. Saturation values σ → ϱ s l/g (p * ) and T s (p * ) as functions of σ for the SRK EoS

Figure 7 .

 7 Figure 7. T (ϱ, p * ) with the best parameters for the different EoS

Figure 8 .

 8 Figure 8. Coexistence curves p → T s (p) for the different EoS

Figure 9 .

 9 Figure 9. Saturation domes in the (ϱ, p) plane for the different EoS

Figure 10 .

 10 Figure 10. Internal energy as a function of τ at three constant temperatures for the Clausius EoS.

Table 5 .

 5 Partial derivatives of the pressure for the different EoS B. Use of the complete cubic EoS in a low Mach number model

  α ′ (T ) -α(T )) C ′ (τ ) = -a (T α ′ (T ) -α(T )) τ (τ + d) + c(τd) and ∂e ∂T τ = c v (τ, T ) = c v,0 (T ) + T α ′′ (T )(C(τ ) -C(τ c )).Remark B.1 (Phase transition in complete EoS in (τ, p)). Since the pressure is the same in both fluids in the mixture at saturation, the complete EoS h Maxwell (τ, p) is defined piecewise using the saturation values:h Maxwell (τ, p) = φ(τ, p)h(τ s l (p), T s (p)) + (1 -φ(τ, p))h(τ s g (p), T s (p)) in the mixture, i.e. τ s l (p) < τ < τ s g (p), pτ + ẽ(τ, p) s g (p),T s (p))-h(τ s l (p),T s (p)) τ s g (p)-τ s l (p),in the mixture, i.e. τ s l (p) < τ < τ s g (p) Eq. (41) otherwise.

  ζ(τ, p) = p + a τ 2 + c 0 pτ 3 + a(2b -τ ) rτ 3 .(42)B.3. Clausius (and Berthelot)From the expression[START_REF] Faccanoni | A homogeneous relaxation low Mach number model[END_REF] of the pressure, we can invert the temperature by solving a second-order equation on T and selecting the positive root, which is given byT (τ, p) = p(τb) 2r + D 2r(τ + δ) ,where we definedD def = (τb)(p 2 (τ + δ) 2 (τb) + 4ar).We can thus compute the partial derivative∂T ∂τ p = p 2r -D 2r(τ + δ) 2 + p 2 (τ + δ)(τb)(2τ + δb) + 2ar 2r(τ + δ)D ,which allows to compute ζ(τ, T ) by (41). Note that the Berthelot equation is recovered when δ = 0.B.4. Soave-Redlich-KwongFrom the expression (9) of the pressure, we can invert the temperature by solving a second-order equation on T and selecting the positive root, which is given byT (τ, p) = -aσ(τ -b) + D rτ (τ + b) -a(τb) ,where we definedD def = τ (τ 2 -b 2 )[a (rσ 2 -p(τb)) + prτ (τ + b)].We can thus compute the partial derivative∂T ∂τ p = -aσ + ( √ D) ′ -T (τ, p)(2rτ + rba) rτ (τ + b) -a(τb) ,where the derivative of √ D with respect to τ is given by( √ D) ′ = 1 2 √ D (3τ2 -b 2 ) a(rσ 2 -p(τb)) + prτ (τ + b) + pτ (τ 2 -b 2 )(2rτ + rba) which allows to compute ζ(τ, T ) by (41).

27a 2 r 4 -

 4 4r 3 p 3 -9ar 2 , r def = r τb , a def = -a τ (τ + b).

r 4 +

 4 18a 2 r ′ r 3 -2p 3 r ′ r 2 ) 27a 2 r 4 -4r 3 p 3 -3a ′ r 2 -6ar ′ r , r ′ = -r (τb) 2 , a ′ = a(2τ + b) τ 2 (τ + b) 2which allows to compute ζ(τ, T ) by (41).

Table 2 .

 2 Van

	.526	2.2064×10 7	647.096 322	594.358 101.925 617.942
		Van der Waals

der Waals, Berthelot and Redlich-Kwong EoS: comparison of the different choices of fixed parameters on computed critical and saturation values at p *

Table 3 .

 3 Comparison of the different choices of r on computed saturation values at p * for the Clausius EoS

	.526	594.358 101.925 617.942

Table 4 .

 4 Computed saturation values at p * for the Soave-Redlich-Kwong EoS

	] [K]

This assumption is referred to as the polytropic case in[START_REF] Quartapelle | Solution of the Riemann problem of classical gasdynamics[END_REF], but this denomination is more frequently used in the literature with another definition.

In fact, the smallest values of the pressure are only computed for Clausius, Redlich-Kwong and Soave-Redlich-Kwong EoS on

This work has been achieved within the framework of the M2SIR project, financially supported by CNRS AAP "NEEDS 2020", and its partners Andra, BRGM,