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 two-phase ﬂows, we use Maxwell’s construction, which amounts to consider pure phases and a mixture zone at saturation. The parameters appearing in the diﬀerent equations of state are ﬁtted in order to be precise in the saturation zone at high pressures. The diﬀerent 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 [26] 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 [25] 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 [3] 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 [2], removed this additional constant, assuming δ = 0. Another important modification of the Van der Waals EoS is the Redlich-Kwong one in 1949 [23], 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 [24], 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 [6,11,16,27].
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 ) [9,22].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 [4,12,[20][21][22].
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 [19].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 [9].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 [6,19].
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.

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 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 ).
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 Among the many CEoS, some well-known families are included in the general form (2). The Patel-Teja (PT) family of EoS results from setting d = b, the Peng-Robinson (PR) one from setting c = d = b (see [16] and the references therein).
The temperature-dependent function α must satisfy the following properties [14][15][16]: 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 [19], obtained by setting c = d = 0 and α(T ) = 1: • the Berthelot (Ber) EoS [2,19],1 obtained by setting c = d = 0 and α(T ) = 1 T : • the Clausius (Cl) EoS [3], 1 obtained by setting c and d such that c2 + d 2 + 6cd = 0 and α(T ) = 1 T : where we denoted δ def = c+d 2 ; • the Redlich-Kwong (RK) EoS [19], 2 obtained by setting c = 0, d = b and α(T • the Soave-Redlich-Kwong (SRK) EoS [14,17,24], 3 obtained by setting c = 0, d = b, and 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 by correcting the coefficient a by a multiplicative factor √ Tc. 3 In the litterature, the function α is usually given of the form , 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 At this point, a fluid exhibits some unusual properties such as infinite heat capacity or infinite compressibility.An alternative method to solve (10) has been introduced in [18].Instead of defining p(τ, T ), one could also define implicitly τ (p, T ) from (3).Then, system (10) defining the critical point becomes 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 . 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 ESAIM: PROCEEDINGS AND SURVEYS 123 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 ds − pdτ with respect to τ to obtain ∂e ∂τ T = T ∂s ∂τ T − p.Using one of the Maxwell relations [9] leads to ∂s ∂τ T = ∂p ∂T τ . Thus the partial derivative of e with respect to τ can be expressed as The other partial derivative defines the specific heat capacity at constant volume c v : Since e(τ, T ) must be an exact differential form, the two mixed second partial derivatives have to be equal to each other: We can write this as which leads to the following compatibility condition on c v Thus, c v can be computed by integrating this compatibility relation ( 16) 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 [1] 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 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 (5), we compute Therefore, the compatibility condition (16) implies that c v (τ, T ) does not depend on τ , and is thus a function of T only.With our previous assumption (19), this implies that Further, we can obtain the energy e(τ, T ) by integrating

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 (7), we compute The compatibility condition gives the dependence of c v (τ, T ) as a function of τ .We use again (19), which leads to since we choose to fit the constants using the critical point.Further, we obtain the energy e(τ, T ) by integrating

Redlich-Kwong caloric equation of state
For the Redlich-Kwong equation of state (8), we compute As before, we compute c v (τ, T ) from the compatibility condition and ( 19) This leads to the energy

Soave-Redlich-Kwong caloric equation of state
For the Soave-Redlich-Kwong equation of state ( 9), we compute Again, we first obtain c v (τ, T ) as and the energy

General caloric equation of state
In this subsection, we perform the formal computations for any generic CEoS, without using the assumption (19).For the general CEoS (2), we compute Let us define ∆ def = c 2 + d 2 + 6cd, and introduce5 Therefore, the compatibility condition implies that Further, we can compute the energy 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 Solving this system for any p < p c allows to define the following functions 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 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 The energy of the mixture at saturation is then given by e(τ, 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 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 [8].For the computation of the speed of sound for the different EoS, see Appendix A. where 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.526J K −1 kg −1 .

Case EoS
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 [26].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 [27].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 [10] 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

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 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: for the Van der Waals EoS, for the Berthelot EoS, 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 ξ: For Soave-Redlich-Kwong EoS, observe that, if we define the "reduced parameter" σ = σ/ √ T c , the reduced SRK EoS only involves this parameter σ: .
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 [13,27] 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 * = 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 [26].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.

.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.

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 [5] 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.

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 [26] 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 Figures 8-9, 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 [26], 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 (19), 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 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 (19), c v (τ, T ) given by (17) 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 Let p(τ, s) = p(τ, T (τ, s)).By the composition rule, we have Using a Maxwell relation [9] and the chain rule, we compute where we used the definition of T def = ∂e ∂s τ in the last equality.We can rewrite this expression of (c * ) 2 using 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  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 [7], 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 (2).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): Eq. ( 41) otherwise.

B.2. Van der Waals
From the expression (5) of the pressure, we can invert the temperature as follows Inserting this expression in relation ( 22) for e, we deduce ẽ(τ, p) and thus the enthalpy and its derivative

Figure 1 .
Figure 1.Schematic coexistence and saturation curves for water

Figure 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.

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

Figure 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 .
Figure 5.Comparison of different values of r for the Clausius EoS

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

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

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

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

Figure 10 .
Figure 10.Internal energy as a function of τ at three constant temperatures for the Clausius EoS.
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).

Table 2 .
Van 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 .
Comparison of the different choices of r on computed saturation values at p * for the Clausius EoS

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