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

HAL
open science

# Comparison of several complete cubic laws for two-phase flow models

Gloria Faccanoni, Bérénice Grec

► To cite this version:

Gloria Faccanoni, Bérénice Grec. Comparison of several complete cubic laws for two-phase flow models. ESAIM: Proceedings and Surveys, 2023, Workshop on Compressible Multiphase Flows: Derivation, closure laws, thermodynamics, 72, pp.117-142. 10.1051/proc/202372117. hal-03863134v2

HAL Id: hal-03863134

https://hal.science/hal-03863134v2

Submitted on 6 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
COMPARISON OF SEVERAL COMPLETE CUBIC LAWS FOR TWO-PHASE
FLOW MODELS *

GLORIA FACCANONI<fnref n="1" /> AND BÉRÉNICE GREC<fnref n="2" />

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

# 1. 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 repul-
sion 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

<footnote>* This work has been achieved within the framework of the M2SIR project, financially supported by CNRS AAP "NEEDS 2020", and its partners Andra, BRGM, CEA, EDF, Framatome, IRSN et Orano.</footnote>

<footnote n="1">¹ Université de Toulon - IMATH, EA 2134, avenue de l'Université, 83957 La Garde, France – e-mail: gloria.faccanoni@univ-tln.fr</footnote>

<footnote n="2">² Université Paris Cité, CNRS, MAP5, F-75006 Paris, France – e-mail: berenice.grec@u-paris.fr</footnote>
CEoS $p = rT/(\tau - b) - a/\tau^2$. The repulsion pressure is described by the term $rT/(\tau - b)$, and the attraction pressure by $a/\tau^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/\tau^2$. In 1880, Clausius [3] suggested that the molecular attraction term was inversely proportional to temperature, leading to the relation $p = rT/(\tau - b) - a/(T(\tau + \delta)^2)$. The addition of a new constant $\delta$ 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 $\delta = 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/(\tau - b) - a/(\sqrt{T}\tau(\tau + 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/\sqrt{T})$ by a more general temperature-dependent term, denoted by $a\alpha(T)$, with $\alpha(T) = (1 + m(1 - \sqrt{T}))^2$. In this relation, $\bar{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 $(\tau, T) \mapsto p(\tau, 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 $(\tau, T) \mapsto e(\tau, 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(\tau, T)$, which requires to invert the complete EoS. In compressible models, the conservation laws lead to express $p(\tau, e)$, whereas in low Mach number models, where the pressure is almost constant, the EoS is expressed as $T(\tau, p)$ or in terms of enthalpy $h(\tau, 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–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(\tau, p)$ as well as its derivative for use in low Mach number models.

## 2. CUBIC EQUATIONS OF STATE

### 2.1. 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 $\tau_c$, and the critical pressure, denoted by $p_c$. For water, the values of these critical constants are

$$p_{c,\text{exp}} = 22.064 \times 10^6 \text{ Pa}, \quad T_{c,\text{exp}} = 647.096 \text{ K}, \quad \tau_{c,\text{exp}} = 1/322 \text{ m}^3 \text{ kg}^{-1}. \quad (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 $(\tau, 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.

### 2.2. 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$, $\tau$ and $T$, usually devised in the literature under the form $p(\tau, T)$.
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Schematic coexistence and saturation curves for water</figcaption></figure>

Cubic EoS are very widely used due to their simplicity. A four parameters Cubic Equation Of State (CEoS) [11] can be written as

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a\alpha(T)}{\tau(\tau + d) + c(\tau - d)}, \quad (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

$$\tau^3 + \left(c+d-b-\frac{rT}{p}\right)\tau^2 - \left dah + (c+d)\left(b+\frac{rT}{p}\right) - \frac{a\alpha(T)}{p} \tau + \left(cbd+cd\frac{rT}{p}-b\frac{a\alpha(T)}{p}\right) = 0. \quad (3)$$

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–16]: α ∈ ℂ²(ℝ⁺) and

$$\alpha(T) \ge 0, \qquad \alpha'(T) \le 0, \qquad \alpha''(T) \ge 0, \qquad \alpha'''(T) \le 0. \quad (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 [19], obtained by setting c = d = 0 and α(T) = 1:

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a}{\tau^2}; \quad (5)$$
• the Berthelot (Ber) EoS [2, 19],<fnref n="1" /> obtained by setting $c = d = 0$ and $\alpha(T) = \frac{1}{T}$:

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a}{T\tau^2}; \qquad (6)$$

• the Clausius (Cl) EoS [3],<fnref n="1" /> obtained by setting $c$ and $d$ such that $c^2 + d^2 + 6cd = 0$ and $\alpha(T) = \frac{1}{T}$:

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a}{T(\tau + \delta)^2}, \qquad (7)$$

where we denoted $\delta \stackrel{\text{def}}{=} \frac{c+d}{2}$;

• the Redlich-Kwong (RK) EoS [19],<fnref n="2" /> obtained by setting $c = 0$, $d = b$ and $\alpha(T) = \frac{1}{\sqrt{T}}$:

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a}{\sqrt{T}\tau(\tau + b)}; \qquad (8)$$

• the Soave-Redlich-Kwong (SRK) EoS [14, 17, 24],<fnref n="3" /> obtained by setting $c = 0$, $d = b$, and $\alpha(T) = (\sigma - \sqrt{T})^2$:

$$p(\tau, T) = \frac{rT}{\tau - b} - \frac{a(\sigma - \sqrt{T})^2}{\tau(\tau + b)}. \qquad (9)$$

Note that, in order to satisfy (4), the parameter $\sigma$ has to be larger than $\sqrt{T}$. In our case, since we are interested in temperatures less than the critical one, it is sufficient that $\sigma > \sqrt{T_c}$, which is the case by definition of $\sigma$.

## 2.3. 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 $\tau$, 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(\tau_c, T_c) = p_c, \quad \left. \frac{\partial p}{\partial \tau} \right|_T (\tau_c, T_c) = 0 \quad \text{and} \quad \left. \frac{\partial^2 p}{\partial \tau^2} \right|_T (\tau_c, T_c) = 0. \qquad (10)$$

<footnote n="1">¹In the literature, the function $\alpha$ is usually given of the form $1/(T/T_c)$. This $\alpha$ can be rewritten under the form we propose by correcting the coefficient $a$ by a multiplicative factor $T_c$.</footnote>

<footnote n="2">²In the literature, the function $\alpha$ is usually given of the form $1/\sqrt{T/T_c}$. This $\alpha$ can be rewritten under the form we propose by correcting the coefficient $a$ by a multiplicative factor $\sqrt{T_c}$.</footnote>

<footnote n="3">³In the literature, the function $\alpha$ is usually given of the form $[1 + m(\omega)(1 - \sqrt{T/T_c})]^2 = \frac{m^2}{T_c^2} [\frac{\sqrt{T_c}}{m} + \sqrt{T_c} - \sqrt{T}]^2$, where $\omega$ denotes a physical constant (acentric factor), and the function $m$ is a polynomial of order two with given coefficients. This $\alpha$ can be rewritten under the form we propose by setting $\sigma = \sqrt{T_c}(1 + 1/m(\omega))$ and correcting the coefficient $a$ by a multiplicative factor $m(\omega)^2/T_c$.</footnote>
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Isotherm curves for a cubic law in a reduced pressure-volume diagram (i.e. ($\tau/\tau_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.</figcaption></figure>

For the general CEoS (2), this leads to the following system

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{rT_c}{\tau_c - b} - \frac{a\alpha(T_c)}{c(\tau_c - d) + \tau_c(\tau_c + d)} = p_c, \\[1.5em]
\displaystyle -\frac{rT_c}{(\tau_c - b)^2} + \frac{a(2\tau_c + c + d)\alpha(T_c)}{(c(\tau_c - d) + \tau_c(\tau_c + d))^2} = 0, \\[1.5em]
\displaystyle \frac{2rT_c}{(\tau_c - b)^3} + \frac{2a\alpha(T_c)}{(c(\tau_c - d) + \tau_c(\tau_c + d))^2} - \frac{2a(2\tau_c + c + d)^2\alpha(T_c)}{(c(\tau_c - d) + \tau_c(\tau_c + d))^3} = 0.
\end{array}
\right.
\qquad (11) $$

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(\tau, T)$, one could also define implicitly $\tau(p, T)$ from (3). Then, system (10) defining the critical point becomes

$$ \tau(p_c, T_c) = \tau_c, \quad \left. \frac{\partial \tau}{\partial p} \right|_T (p_c, T_c) = 0 \quad \text{and} \quad \left. \frac{\partial^2 \tau}{\partial^2 p} \right|_T (p_c, T_c) = 0. \qquad (12) $$

From this and the fact that (3) is a unitary polynomial, it follows that the function $\tau(p, T)$ should be equal at the critical point to $(\tau - \tau_c)^3 = \tau^3 - 3\tau_c\tau^2 + 3\tau_c^3\tau - \tau_c^3$. 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

$$ 
\left\{ 
\begin{array}{l}
c + d - b - \displaystyle\frac{rT_c}{p_c} = -3\tau_c, \\[1.5em]
cd + (c+d)\left(b + \displaystyle\frac{rT_c}{p_c}\right) - \displaystyle\frac{a\alpha(T_c)}{p_c} = -3\tau_c^2, \\[1.5em]
cbd + cd\displaystyle\frac{rT_c}{p_c} - b\displaystyle\frac{a\alpha(T_c)}{p_c} = -\tau_c^3.
\end{array}
\right.
\qquad (13) $$
### 3. CONSTRUCTION OF A COMPLETE CUBIC EOS USING THE VARIABLES $(\tau, 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(\tau, T)$.

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

#### 3.1. Main idea

We search for the specific internal energy $e$ as a function of temperature $T$ and specific volume $\tau$. Denoting $s$ the entropy, we differentiate the Gibbs relation $de = Tds - pd\tau$ with respect to $\tau$ to obtain $\left.\frac{\partial e}{\partial \tau}\right|_T = T \left.\frac{\partial s}{\partial \tau}\right|_T - p$. Using one of the Maxwell relations [9] leads to $\left.\frac{\partial p}{\partial \tau}\right|_T = \left.\frac{\partial s}{\partial T}\right|_T$. Thus the partial derivative of $e$ with respect to $\tau$ can be expressed as

$$ \left.\frac{\partial e}{\partial \tau}\right|_T = T \left.\frac{\partial p}{\partial T}\right|_{\tau} - p. \qquad (14) $$

The other partial derivative defines the specific heat capacity at constant volume $c_v$:

$$ \left.\frac{\partial e}{\partial T}\right|_{\tau} \stackrel{\text{def}}{=} c_v. \qquad (15) $$

Since $e(\tau, T)$ must be an exact differential form, the two mixed second partial derivatives have to be equal to each other:

$$ \left.\frac{\partial}{\partial \tau}\left(\left.\frac{\partial e}{\partial T}\right|_{\tau}\right)\right|_{T} = \left.\frac{\partial}{\partial T}\left(\left.\frac{\partial e}{\partial \tau}\right|_{T}\right)\right|_{\tau}. $$

We can write this as

$$ \left.\frac{\partial c_v}{\partial \tau}\right|_T = \left.\frac{\partial}{\partial T}\left(T \left.\frac{\partial p}{\partial T}\right|_{\tau} - p\right)\right|_{\tau} = \left.\frac{\partial p}{\partial T}\right|_{\tau} + T \left.\frac{\partial^2 p}{\partial T^2}\right|_{\tau} - \left.\frac{\partial p}{\partial T}\right|_{\tau}, $$

which leads to the following compatibility condition on $c_v$

$$ \left.\frac{\partial c_v}{\partial \tau}\right|_T = T \left.\frac{\partial^2 p}{\partial T^2}\right|_{\tau}. \qquad (16) $$

Thus, $c_v$ can be computed by integrating this compatibility relation (16)

$$ c_v(\tau, T) = c_{v,0}(T) + \int_{\tau_0}^{\tau} T \left.\frac{\partial^2 p(\zeta, T)}{\partial T^2}\right|_{\zeta} d\zeta. \qquad (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 $\tau$, 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(\tau, T)$, we deduce that

$$ e(\tau, T) = \int_{\tau_0}^{\tau} \left.\frac{\partial e}{\partial \tau}\right|_{T}(\zeta, T)d\zeta + \int_{T_0}^{T} \left.\frac{\partial e}{\partial T}\right|_{\tau}(\tau_0, \theta)d\theta + e(\tau_0, T_0) $$
and thus, using the expression of $c_v(\tau, T)$, we obtain

$$e(\tau, T) = \int_{\tau_0}^{\tau} \left( T \left. \frac{\partial p(\zeta, T)}{\partial T} \right|_{\tau} - p(\zeta, T) \right) d\zeta + \int_{T_0}^{T} c_{v,0}(\theta) d\theta + e_0. \quad (18)$$

The point $(\tau_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 $(\tau_c, T_c)$.

**Assumption.** For the sake of simplicity, we choose $c_{v,0}(T)$ to be independent on $T$,<fnref n="4" /> so that

$$c_{v,0}(T) = c_0. \quad (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(\tau, T)$ explicitly, as showed in the next subsections.

## 3.2. Van der Waals caloric equation of state

For the Van der Waals equation of state (5), we compute

$$\left. \frac{\partial^2 p}{\partial T^2} \right|_{\tau} (\tau, T) = 0. \quad (20)$$

Therefore, the compatibility condition (16) implies that $c_v(\tau, T)$ does not depend on $\tau$, and is thus a function of $T$ only. With our previous assumption (19), this implies that

$$c_v(\tau, T) = c_0. \quad (21)$$

Further, we can obtain the energy $e(\tau, T)$ by integrating

$$e(\tau, T) = e_c + c_0(T - T_c) - a \left( \frac{1}{\tau} - \frac{1}{\tau_c} \right). \quad (22)$$

## 3.3. Berthelot and Clausius caloric equation of state

Observe that Berthelot EoS can be obtained from the Clausius one by setting $\delta = 0$. We thus proceed here with the general Clausius EoS. For the Clausius equation of state (7), we compute

$$\left. \frac{\partial^2 p}{\partial T^2} \right|_{\tau} (\tau, T) = - \frac{2a}{T^3(\tau + \delta)^2}. \quad (23)$$

The compatibility condition gives the dependence of $c_v(\tau, T)$ as a function of $\tau$. We use again (19), which leads to

$$c_v(\tau, T) = c_0 + \frac{2a}{T^2} \left( \frac{1}{\tau + \delta} - \frac{1}{\tau_c + \delta} \right), \quad (24)$$

since we choose to fit the constants using the critical point. Further, we obtain the energy $e(\tau, T)$ by integrating

$$e(\tau, T) = e_c + c_0(T - T_c) - \frac{2a}{T} \left( \frac{1}{\tau + \delta} - \frac{1}{\tau_c + \delta} \right). \quad (25)$$

<footnote n="4">⁴This assumption is referred to as the polytropic case in [21], but this denomination is more frequently used in the literature with another definition.</footnote>
### 3.4. Redlich-Kwong caloric equation of state

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

$$ \left. \frac{\partial^2 p}{\partial T^2} \right|_{\tau} (\tau, T) = - \frac{3a}{4T^2 \sqrt{T} \tau (\tau + b)}. \qquad (26) $$

As before, we compute $c_v(\tau, T)$ from the compatibility condition and (19)

$$ c_v(\tau, T) = c_0 + \frac{3a}{4bT\sqrt{T}} \ln \left( \frac{\tau + b}{\tau} \frac{\tau_c}{\tau_c + b} \right). \qquad (27) $$

This leads to the energy

$$ e(\tau, T) = e_c + c_0(T - T_c) - \frac{3a}{2b} \frac{1}{\sqrt{T}} \ln \left( \frac{\tau + b}{\tau} \frac{\tau_c}{\tau_c + b} \right). \qquad (28) $$

### 3.5. Soave-Redlich-Kwong caloric equation of state

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

$$ \left. \frac{\partial^2 p}{\partial T^2} \right|_{\tau} (\tau, T) = - \frac{a\sigma}{2T\sqrt{T}\tau(\tau + b)}. \qquad (29) $$

Again, we first obtain $c_v(\tau, T)$ as

$$ c_v(\tau, T) = c_0 + \frac{a\sigma}{2b\sqrt{T}} \ln \left( \frac{\tau + b}{\tau} \frac{\tau_c}{\tau_c + b} \right), \qquad (30) $$

and the energy

$$ e(\tau, T) = e_c + c_0(T - T_c) - \frac{a\sigma}{b} (\sqrt{T} - \sigma) \ln \left( \frac{\tau}{\tau + b} \frac{\tau_c + b}{\tau_c} \right). \qquad (31) $$

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

$$ \left. \frac{\partial^2 p}{\partial T^2} \right|_{\tau} (\tau, T) = \frac{-a}{\tau(\tau + d) + c(\tau - d)} \alpha''(T). \qquad (32) $$

Let us define $\Delta \stackrel{\text{def}}{=} c^2 + d^2 + 6cd$, and introduce$^5$

$$ \mathfrak{E}(\tau) \stackrel{\text{def}}{=} \int_{\tau} - \frac{a}{\zeta(\zeta + d) + c(\zeta - d)} d\zeta = \begin{cases} - \frac{a}{\sqrt{\Delta}} \ln \left( \frac{2\tau + c + d - \sqrt{\Delta}}{2\tau + c + d + \sqrt{\Delta}} \right) & \text{if } \Delta > 0, \\ \frac{a}{\tau + \frac{c+d}{2}} & \text{if } \Delta = 0. \end{cases} \qquad (33) $$

<footnote>$^5$The Clausius EoS corresponds to $\Delta = c^2 + d^2 + 6cd = 0$ with $\delta = \frac{c+d}{2}$, and the (Soave)-Redlich-Kwong families correspond to $\Delta = b^2 > 0$.</footnote>
Therefore, the compatibility condition implies that

$$c_v(\tau, T) = c_{v,0}(T) + \int_{\tau_c}^{\tau} T \left. \frac{\partial^2 p(\zeta, T)}{\partial T^2} \right|_{\zeta} d\zeta = c_{v,0}(T) + T\alpha''(T) [\mathfrak{E}(\tau)]_{\tau_c}^{\tau}. \quad (34)$$

Further, we can compute the energy

$$\begin{aligned}
e(\tau, T) &= e_c + \int_{T_c}^{T} c_{v,0}(\theta) d\theta + \int_{\tau_c}^{\tau} \left( T \left. \frac{\partial p(\zeta, T)}{\partial T} \right|_{\tau} - p(\zeta, T) \right) d\zeta \\
&= e_c + \int_{T_c}^{T} c_{v,0}(\theta) d\theta + (T\alpha'(T) - \alpha(T)) \int_{\tau_c}^{\tau} \frac{-a}{\zeta(\zeta+d) + c(\zeta-d)} d\zeta \\
&= e_c + \int_{T_c}^{T} c_{v,0}(\theta) d\theta + (T\alpha'(T) - \alpha(T)) [\mathfrak{E}(\tau)]_{\tau_c}^{\tau}.
\end{aligned}
\quad (35)$$

# 4. PHASE EQUILIBRIUM

## 4.1. 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 $\tau_l^s$, $\tau_g^s$ and $T^s$ such that the three following relations hold

$$\left\{
\begin{array}{l}
p(\tau_l^s, T^s) = p_0, \\
p(\tau_g^s, T^s) = p_0, \\
\displaystyle\int_{\tau_l^s}^{\tau_g^s} (p(\tau, T^s) - p_0) d\tau = 0.
\end{array}
\right.
\quad (36)$$

Solving this system for any $p < p_c$ allows to define the following functions

$$p \mapsto (\tau_l^s(p), \tau_g^s(p), T^s(p)). \quad (37)$$

Equivalently, for any fixed $T_0 < T_c$, we can look for $\tau_l^s$, $\tau_g^s$ and $p^s$ such that $p(\tau_l^s, T_0) = p^s$, $p(\tau_g^s, T_0) = p^s$ and $\int_{\tau_l^s}^{\tau_g^s} (p(\tau, T_0) - p^s) d\tau = 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

$$T \mapsto (\tau_l^s(T), \tau_g^s(T), p^s(T)). \quad (38)$$
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Isotherm curves for a cubic law in a pressure-volume plane (τ, P) for the Maxwell construction</figcaption></figure>

After such a construction, we define the EoS in a piecewise manner (see Figure 3b):

$$p_{\text{Maxwell}}(\tau, T) = \begin{cases} p(\tau, T) & \text{in pure phases, i.e. } \tau \le \tau_{\text{l}}^s(T) \text{ or } \tau \ge \tau_{\text{g}}^s(T), \\ p^s(T) & \text{in the mixture, i.e. } \tau_{\text{l}}^s(T) \le \tau \le \tau_{\text{g}}^s(T). \end{cases}$$

For the different EoS introduced in Section 2.2, we can compute explicitly the integral to rewrite system (36) in the following way: find $(\tau_{\text{l}}^s, \tau_{\text{g}}^s, T^s)$ such that

$$\left\{ \begin{aligned} & \frac{rT^s}{\tau_{\text{l}}^s - b} - \frac{a\alpha(T^s)}{\tau_{\text{l}}^s(\tau_{\text{l}}^s + d) + c(\tau_{\text{l}}^s - d)} = p_0, \\ & \frac{rT^s}{\tau_{\text{g}}^s - b} - \frac{a\alpha(T^s)}{\tau_{\text{g}}^s(\tau_{\text{g}}^s + d) + c\tau_{\text{g}}^s(\tau_{\text{g}}^s - d)} = p_0, \\ & rT^s \ln\left(\frac{\tau_{\text{g}}^s - b}{\tau_{\text{l}}^s - b}\right) + \alpha(T^s)(\mathcal{C}(\tau_{\text{g}}^s) - \mathcal{C}(\tau_{\text{l}}^s)) - p_0(\tau_{\text{g}}^s - \tau_{\text{l}}^s) = 0. \end{aligned} \right.$$

where $\mathcal{C}(\tau)$ has been defined in (33). The functions $\mathcal{C}$ and $\alpha$ are (re)given here for the sake of completeness.

* Van der Waals EoS: $\mathcal{C}(\tau) = 1/\tau$ and $\alpha(T) = 1$

* Berthelot EoS: $\mathcal{C}(\tau) = 1/\tau$ and $\alpha(T) = 1/T$

* Clausius EoS: $\mathcal{C}(\tau) = 1/(\tau + \delta)$ and $\alpha(T) = 1/T$

* Redlich-Kwong EoS: $\mathcal{C}(\tau) = (\ln(\tau + b) - \ln(\tau))/b$ and $\alpha(T) = 1/\sqrt{T}$

* Soave-Redlich-Kwong EoS: $\mathcal{C}(\tau) = (\ln(\tau + b) - \ln(\tau))/b$ and $\alpha(T) = (\sigma - \sqrt{T})^2$
## 4.2. Phase transition in complete EoS in $(\tau, T)$

Let us now explain how the complete EoS is constructed in the mixture zone. For $(\tau, T)$ given such that $\tau_{\|}^s(T) < \tau < \tau_{\|}^s(T)$, we define the mass fraction of the mixture

$$ \varphi(\tau, T) \stackrel{\text{def}}{=} \frac{\tau - \tau_{\|}^s(T)}{\tau_{\|}^s(T) - \tau_{\|}^s(T)} \in (0, 1). $$

The energy of the mixture at saturation is then given by

$$ e(\tau, T) \stackrel{\text{def}}{=} \varphi(\tau, T)e_{\|}^s(T) + (1 - \varphi(\tau, T))e_{\|}^s(T), $$

where we define the functions at saturation $e_{\|/g}^s(T) \stackrel{\text{def}}{=} e(\tau_{\|/g}^s(T), T)$ using (38) in (35).

The complete EoS $e_{\text{Maxwell}}(\tau, T)$ is then defined piecewise using the saturation values:

$$ e_{\text{Maxwell}}(\tau, T) = \begin{cases} \varphi(\tau, T)e_{\|}^s(T) + (1 - \varphi(\tau, T))e_{\|}^s(T) & \text{in the mixture, i.e. } \tau_{\|}^s(T) < \tau < \tau_{\|}^s(T), \\ \text{Eq. (35)} & \text{otherwise.} \end{cases} $$

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.

**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(\tau, T) = \left. \frac{\partial \varphi}{\partial T} \right|_{\tau} (e_{\|}^s(T) - e_{\|}^s(T)) + \varphi(\tau, T)(e_{\|}^s)'(T) + (1 - \varphi(\tau, T))(e_{\|}^s)'(T) $$

where

$$ \begin{aligned} \left. \frac{\partial \varphi}{\partial T} \right|_{\tau} (\tau, T) &= \frac{\varphi(\tau, T)(\tau_{\|}^s)'(T) + (1 - \varphi(\tau, T))(\tau_{\|}^s)'(T)}{\tau_{\|}^s(T) - \tau_{\|}^s(T)}, \\ (e_{\kappa}^s)'(T) &= \left. \frac{\partial e}{\partial T} \right|_{\tau} (\tau_{\kappa}^s(T), T) + (\tau_{\kappa}^s)'(T) \times \left. \frac{\partial e}{\partial \tau} \right|_{T} (\tau_{\kappa}^s(T), T) = c_v(\tau_{\kappa}^s(T), T) + (\tau_{\kappa}^s)'(T) \times (T \times (p^s)'(T) - p^s(T)), \\ (p^s)'(T) &= \frac{1}{T} \frac{e_{\|}^s(T) - e_{\|}^s(T)}{\tau_{\|}^s(T) - \tau_{\|}^s(T)} + \frac{p^s(T)}{T} \quad (\text{Clapeyron relation}). \end{aligned} $$

Analytical expressions for the phase boundaries $\tau_{\|/g}^s(T)$ or their derivatives are not available, but they can be computed numerically.

# 5. 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 $\delta$ or $\sigma$), the experimental critical values $(T_{c,\text{exp}}, \tau_{c,\text{exp}}, p_{c,\text{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_{\text{exp}} \stackrel{\text{def}}{=} R/M$, where $R$ is the ideal gas constant and $M$ the molar mass of the fluid. For water, $r_{\text{exp}} = 461.526 \text{ J K}^{-1} \text{ kg}^{-1}$.
<table><caption>TABLE 1. All possible choices of fixed parameters for Van der Waals, Berthelot and Redlich-Kwong EoS.</caption><thead><tr><th>Case</th><th>EoS</th><th>T<sub>c</sub></th><th>p<sub>c</sub></th><th>τ<sub>c</sub></th><th>r</th><th>a</th><th>b</th></tr></thead><tbody><tr><td rowspan="2">①</td><td>Van der Waals Berthelot</td><td rowspan="2">T<sub>c,exp</sub></td><td rowspan="2">p<sub>c,exp</sub></td><td><sup>3rT<sub>c</sub></sup>&frasl;<sub>8p<sub>c</sub></sub></td><td rowspan="2">r<sub>exp</sub></td><td>a<sub>VdW</sub> = <sup>27T<sub>c</sub><sup>2</sup>r<sup>2</sup></sup>&frasl;<sub>64p<sub>c</sub></sub>, a<sub>Ber</sub> = a<sub>VdW</sub>T<sub>c</sub></td><td><sup>rT<sub>c</sub></sup>&frasl;<sub>8p<sub>c</sub></sub></td></tr><tr><td>Redlich-Kwong</td><td><sup>rT<sub>c</sub></sup>&frasl;<sub>3p<sub>c</sub></sub></td><td><sup>(1 + <sup>3</sup>&radic;<span style="text-decoration: overline">2</span> + <sup>3</sup>&radic;<span style="text-decoration: overline">4</span>) T<sub>c</sub><sup>2</sup> &radic;<span style="text-decoration: overline">T<sub>c</sub></span>r<sup>2</sup></sup></sup>&frasl;<sub>9p<sub>c</sub></sub></td><td><sup>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) rT<sub>c</sub></sup>&frasl;<sub>3p<sub>c</sub></sub></td></tr><tr><td rowspan="2">②</td><td>Van der Waals Berthelot</td><td rowspan="2">T<sub>c,exp</sub></td><td><sup>3rT<sub>c</sub></sup>&frasl;<sub>8τ<sub>c</sub></sub></td><td rowspan="2">τ<sub>c,exp</sub></td><td rowspan="2">r<sub>exp</sub></td><td>a<sub>VdW</sub> = <sup>9rT<sub>c</sub>τ<sub>c</sub></sup>&frasl;<sub>8</sub>, a<sub>Ber</sub> = a<sub>VdW</sub>T<sub>c</sub></td><td><sup>τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr><tr><td>Redlich-Kwong</td><td><sup>rT<sub>c</sub></sup>&frasl;<sub>3τ<sub>c</sub></sub></td><td><sup>(1 + <sup>3</sup>&radic;<span style="text-decoration: overline">2</span> + <sup>3</sup>&radic;<span style="text-decoration: overline">4</span>) rT<sub>c</sub><sup>2</sup> &radic;<span style="text-decoration: overline">T<sub>c</sub></span>τ<sub>c</sub></sup></sup>&frasl;<sub>3</sub></td><td><sup>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr><tr><td rowspan="2">③</td><td>Van der Waals Berthelot</td><td><sup>8p<sub>c</sub>τ<sub>c</sub></sup>&frasl;<sub>3r</sub></td><td rowspan="2">p<sub>c,exp</sub></td><td rowspan="2">τ<sub>c,exp</sub></td><td rowspan="2">r<sub>exp</sub></td><td>a<sub>VdW</sub> = 3p<sub>c</sub>τ<sub>c</sub><sup>2</sup>, a<sub>Ber</sub> = a<sub>VdW</sub>T<sub>c</sub></td><td><sup>τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr><tr><td>Redlich-Kwong</td><td><sup>3p<sub>c</sub>τ<sub>c</sub></sup>&frasl;<sub>r</sub></td><td><sup>p<sub>c</sub>τ<sub>c</sub><sup>2</sup> &radic;<span style="text-decoration: overline">3p<sub>c</sub>τ<sub>c</sub></span></sup>&frasl;<sub>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) r</sub></td><td><sup>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr><tr><td rowspan="2">④</td><td>Van der Waals Berthelot</td><td rowspan="2">T<sub>c,exp</sub></td><td rowspan="2">p<sub>c,exp</sub></td><td rowspan="2">τ<sub>c,exp</sub></td><td><sup>8p<sub>c</sub>τ<sub>c</sub></sup>&frasl;<sub>3T<sub>c</sub></sub></td><td>a<sub>VdW</sub> = 3p<sub>c</sub>τ<sub>c</sub><sup>2</sup>, a<sub>Ber</sub> = a<sub>VdW</sub>T<sub>c</sub></td><td><sup>τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr><tr><td>Redlich-Kwong</td><td><sup>3p<sub>c</sub>τ<sub>c</sub></sup>&frasl;<sub>T<sub>c</sub></sub></td><td><sup>(1 + <sup>3</sup>&radic;<span style="text-decoration: overline">2</span> + <sup>3</sup>&radic;<span style="text-decoration: overline">4</span>) p<sub>c</sub> &radic;<span style="text-decoration: overline">T<sub>c</sub></span>τ<sub>c</sub><sup>2</sup></sup>&frasl;<sub>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) τ<sub>c</sub></sub></td><td><sup>(<sup>3</sup>&radic;<span style="text-decoration: overline">2</span> - 1) τ<sub>c</sub></sup>&frasl;<sub>3</sub></td></tr></tbody></table>

(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 $\tau_{l/g}^s(p)$ with respect to experimental data for different pressures to assess the quality of the EoS.

## 5.1. First step: expression of the parameters from the critical values

The critical point $(T_c, \tau_c, p_c)$ is obtained by solving the system (13).

### 5.1.1. *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, \tau_c \neq \tau_{c,exp}$. Historically, this approach was chosen since the experimental value $\tau_{c,exp}$ was not known [27].

② We fix $\tau_c = \tau_{c,exp}$, $T_c = T_{c,exp}$, $r = r_{exp}$ and deduce $a, b$, and $p_c \neq p_{c,exp}$.

③ We fix $\tau_c = \tau_{c,exp}$, $p_c = p_{c,exp}$, $r = r_{exp}$ and deduce $a, b$, and $T_c \neq T_{c,exp}$.

④ We fix $T_c = T_{c,exp}$, $p_c = p_{c,exp}$, $\tau = \tau_{c,exp}$ and deduce $a, b$, and $r \neq r_{exp}$.

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.

### 5.1.2. Clausius EoS

For the Clausius EoS, we now have three parameters ($a$, $b$ and $\delta$). 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 $\delta$ are then recovered by the following relations

$$a(r) = \frac{27 T_{c, \exp}^3 r^2}{64 p_{c, \exp}}, \quad b(r) = \tau_{c, \exp} - \frac{r T_{c, \exp}}{4 p_{c, \exp}}, \quad \delta(r) = \frac{3 r T_{c, \exp}}{8 p_{c, \exp}} - \tau_{c, \exp}.$$

### 5.1.3. Soave-Redlich-Kwong EoS

For the Soave-Redlich-Kwong (SRK) EoS, we again have three parameters ($a$, $b$ and $\sigma$). However, we observe that the first equation of (13) does not depend on the choice of the function $\alpha$. Moreover, in this case, since $c = 0$ and $d = b$, no parameter ($a$, $b$ and $\sigma$) 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(\sigma) = \frac{(1 + \sqrt[3]{2} + \sqrt[3]{4}) p_{c, \exp} \tau_{c, \exp}^2}{(T_{c, \exp} - \sigma)^2}, \quad b = (\sqrt[3]{2} - 1) \tau_{c, \exp}, \quad r = \frac{3 p_{c, \exp} \tau_{c, \exp}}{T_{c, \exp}} \neq r_{\exp}.$$

Observe that the parameter $\sigma$ 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 $\sigma$.

### 5.1.4. 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 $\bar{T} = T/T_c$, reduced pressure $\bar{p} = p/p_c$ and reduced specific volume $\bar{\tau} = \tau/\tau_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:

$$\bar{p}(\bar{\tau}, \bar{T}) = \begin{cases} \frac{8}{3} \frac{\bar{T}}{\bar{\tau} - 1/3} - \frac{3}{\bar{\tau}^2} & \text{for the Van der Waals EoS,} \\ \frac{8}{3} \frac{\bar{T}}{\bar{\tau} - 1/3} - \frac{3}{\bar{T}\bar{\tau}^2} & \text{for the Berthelot EoS,} \\ \frac{3\bar{T}}{\bar{\tau} - (\sqrt[3]{2} - 1)} - \frac{(1 + \sqrt[3]{2} + \sqrt[3]{4})}{\sqrt{\bar{T}}\bar{\tau}(\bar{\tau} + (\sqrt[3]{2} - 1))} & \text{for the Redlich-Kwong EoS,} \end{cases}$$

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 $\bar{\xi} = 8p_{c,exp}\tau_{c,exp}/T_{c,exp}$, the reduced Clausius EoS only involves this quantity $\bar{\xi}$:

$$ \bar{p} = \frac{8\bar{T}}{3\bar{\xi}\bar{\tau} + 2 - 3\bar{\xi}} - \frac{3}{\bar{T}(\bar{\xi}\bar{\tau} + 1 - \bar{\xi})^2} $$

For Soave-Redlich-Kwong EoS, observe that, if we define the “reduced parameter” $\bar{\sigma} = \sigma/\sqrt{T_c}$, the reduced SRK EoS only involves this parameter $\bar{\sigma}$:

$$ \bar{p} = \frac{3\bar{T}}{\bar{\tau} - (\sqrt[3]{2} - 1)} - \frac{(1 + \sqrt[3]{2} + \sqrt[3]{4})}{\bar{\sigma} - 1} \frac{(\bar{\sigma} - \sqrt{\bar{T}})^2}{\bar{\tau}(\bar{\tau} + (\sqrt[3]{2} - 1))} $$

Observe that the quantity $\bar{\xi}$ and the “reduced” parameter $\bar{\sigma}$ 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 $\bar{\xi}$ has a fixed value, whereas the reduced parameter $\bar{\sigma}$ involves the constant $\sigma$, and can therefore be chosen independently.

## 5.2. 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_* = 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.

### 5.2.1. 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 readability, we chose to display the density $\rho = 1/\tau$ instead of the specific volume $\tau$. 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}$, $\tau_{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 $(\rho, 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 loses precision.
<table><caption>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_*$</caption><thead><tr><th>Fixed param.</th><th>a<br>[m<sup>3</sup> J kg<sup>-2</sup>]</th><th>b<br>[m<sup>3</sup> kg<sup>-1</sup>]</th><th>r<br>[J K<sup>-1</sup> kg<sup>-1</sup>]</th><th>p<sub>c</sub><br>[Pa]</th><th>T<sub>c</sub><br>[K]</th><th>ρ<sub>c</sub><br>[kg m<sup>-3</sup>]</th><th>ρ<sub>c</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>ρ<sub>c</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>T<sup>s</sup>(p<sub>*</sub>)<br>[K]</th></tr></thead><tbody><tr><td colspan="10" style="text-align: center; font-weight: bold;">IAPWS</td></tr><tr><td></td><td></td><td></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>322</td><td>594.358</td><td>101.925</td><td>617.942</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Van der Waals</td></tr><tr><td>pTτ</td><td>638.401</td><td>1.0352×10<sup>-3</sup></td><td>282.376</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>322</td><td>513.77</td><td>152.106</td><td>593.788</td></tr><tr><td>Tτr</td><td>1043.43</td><td>1.0352×10<sup>-3</sup></td><td>461.526</td><td>3.60622×10<sup>7</sup></td><td>647.096</td><td>322</td><td>606.717</td><td>86.8006</td><td>530.753</td></tr><tr><td>pTr</td><td>1705.41</td><td>1.69196×10<sup>-3</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>197.01</td><td>314.341</td><td>93.0631</td><td>593.788</td></tr><tr><td>pτr</td><td>638.401</td><td>1.0352×10<sup>-3</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>395.914</td><td>322</td><td>513.77</td><td>152.106</td><td>363.299</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Berthelot</td></tr><tr><td>pTτ</td><td>413107</td><td>1.0352×10<sup>-3</sup></td><td>282.376</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>322</td><td>525.844</td><td>142.917</td><td>616.312</td></tr><tr><td>Tτr</td><td>675197</td><td>1.0352×10<sup>-3</sup></td><td>461.526</td><td>3.60622×10<sup>7</sup></td><td>647.096</td><td>322</td><td>621.923</td><td>77.4192</td><td>578.977</td></tr><tr><td>pTr</td><td>1.10357×10<sup>6</sup></td><td>1.69196×10<sup>-3</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>197.01</td><td>321.728</td><td>87.4414</td><td>616.312</td></tr><tr><td>pτr</td><td>252752</td><td>1.0352×10<sup>-3</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>395.914</td><td>322</td><td>525.844</td><td>142.917</td><td>377.079</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">Redlich-Kwong</td></tr><tr><td>pTτ</td><td>20826.5</td><td>8.07208×10<sup>-4</sup></td><td>317.673</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>322</td><td>562.68</td><td>136.755</td><td>608.57</td></tr><tr><td>Tτr</td><td>30257.3</td><td>8.07208×10<sup>-4</sup></td><td>461.526</td><td>3.20553×10<sup>7</sup></td><td>647.096</td><td>322</td><td>664.945</td><td>84.3884</td><td>572.579</td></tr><tr><td>pTr</td><td>43958.8</td><td>1.17274×10<sup>-3</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>647.096</td><td>221.636</td><td>387.299</td><td>94.1299</td><td>608.57</td></tr><tr><td>pτr</td><td>17278.6</td><td>8.07208×10<sup>-4</sup></td><td>461.526</td><td>2.2064×10<sup>7</sup></td><td>445.403</td><td>322</td><td>562.68</td><td>136.755</td><td>418.886</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Temperature $T(\rho, p_*)$ as a function of $\rho$ for all different choices of fixed parameters for Van der Waals, Berthelot and Redlich-Kwong EoS</figcaption></figure>

### 5.2.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 $\rho_c^s$ and
<table><caption>TABLE 3. Comparison of the different choices of *r* on computed saturation values at *p*<sub>*</sub> for the Clausius EoS</caption><thead><tr><th></th><th>a<br>[m<sup>3</sup> J K kg<sup>-2</sup>]</th><th>b<br>[m<sup>3</sup> kg<sup>-1</sup>]</th><th>δ<br>[m<sup>3</sup> kg<sup>-1</sup>]</th><th>r<br>[J K<sup>-1</sup> kg<sup>-1</sup>]</th><th>ρ<sub>l</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>ρ<sub>g</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>T<sup>s</sup>(p<sub>*</sub>)<br>[K]</th></tr></thead><tbody><tr><td>IAPWS</td><td></td><td></td><td></td><td>461.526</td><td>594.358</td><td>101.925</td><td>617.942</td></tr><tr><td>r<sub>exp</sub></td><td>1.10357×10<sup>6</sup></td><td>-2.78334×10<sup>-4</sup></td><td>0.0019703</td><td>461.526</td><td>878.8</td><td>105.642</td><td>616.312</td></tr><tr><td>r<sub>Ber</sub></td><td>413107</td><td>1.0352×10<sup>-3</sup></td><td>0</td><td>282.376</td><td>525.844</td><td>142.917</td><td>616.312</td></tr><tr><td>r<sub>opt</sub></td><td>577587</td><td>6.57483×10<sup>-4</sup></td><td>0.00056657</td><td>333.892</td><td>594.505</td><td>129.752</td><td>616.312</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Comparison of different values of *r* for the Clausius EoS</figcaption></figure>

$\rho_g^s$ 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 $\rho_g^s$, the error on $\rho_l^s$ increases rapidly. Therefore, we can define an optimal value of *r*, corresponding to minimizing simultaneously both the normalized $L^1$ error on $\rho_g^s$ and the one on $\rho_l^s$. Finally, we observe that the precision obtained with this choice of parameteres in the Clausius EoS is better than all other choices.

### 5.2.3. Soave-Redlich-Kwong EoS

For the Soave-Redlich-Kwong EoS, we already mentioned that the parameter $\sigma$ 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 $\sigma$ is related to the temperature, it affects all three saturation values). We observe on Figure 6 that the optimal values of $\sigma$ 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 $\sigma_{opt} = 48$. For this choice, the parameters and saturation values are given in Table 4.
<table><caption>TABLE 4. Computed saturation values at $p_*$ for the Soave-Redlich-Kwong EoS</caption><thead><tr><th></th><th>a<br>[m<sup>3</sup> J K kg<sup>-2</sup>]</th><th>b<br>[m<sup>3</sup> kg<sup>-1</sup>]</th><th>r<br>[J K<sup>-1</sup> kg<sup>-1</sup>]</th><th>ρ<sub>l</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>ρ<sub>g</sub><sup>s</sup>(p<sub>*</sub>)<br>[kg m<sup>-3</sup>]</th><th>T<sup>s</sup>(p<sub>*</sub>)<br>[K]</th></tr></thead><tbody><tr><td>IAPWS</td><td></td><td></td><td>461.526</td><td>594.358</td><td>101.925</td><td>617.942</td></tr><tr><td>σ<sub>opt</sub> = 48</td><td>1.60834</td><td>0.0008.07208×10<sup>-3</sup></td><td>317.673</td><td>568.93</td><td>133.162</td><td>617.974</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Saturation values $\sigma \mapsto \rho_{l/g}^s(p_*)$ and $T^s(p_*)$ as functions of $\sigma$ for the SRK EoS</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>FIGURE 7. $T(\rho, p_*)$ with the best parameters for the different EoS</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Coexistence curves $p \mapsto T^s(p)$ for the different EoS</figcaption></figure>

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

### 5.3.1. Coexistence curve

More precisely, we plot on Figure 8 the coexistence curves $p \mapsto T^s(p)$ for variable pressures, starting from the smallest pressure in tabulated values up to the critical pressure $p_c$, i.e. $p \in [7 \times 10^2 \text{Pa}, 22.06 \times 10^6 \text{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 $\delta$, and both Berthelot and Clausius EoS are very satisfying. For the Soave-Redlich-Kwong EoS, the choice we made of optimal $\sigma$ 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.

### 5.3.2. Phase boundaries

We further investigate the quality of the parameter fitting for all $p \in [7 \times 10^2 \text{Pa}, 22.06 \times 10^6 \text{Pa}]$ by comparing the saturation zones. As before, for the sake of readibility, we chose to work in the plane $(\varrho, p)$. In this plane, the mixture is at saturation in the set $\{ (\varrho, p) \mid \varrho_g^s(p) \le \varrho \le \varrho_l^s(p) \}$, the fluid is in the liquid phase for $\varrho > \varrho_l^s(p)$ and in the gas phase for $\varrho < \varrho_g^s(p)$. We can thus compare the phase boundaries $\varrho_{l/g}^s(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 $\varrho \mapsto p(\varrho, T^s(p_*))$ for each EoS. Since the value of $p_* = 155 \text{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

<footnote>$^6$In fact, the smallest values of the pressure are only computed for Clausius, Redlich-Kwong and Soave-Redlich-Kwong EoS on 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.</footnote>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Saturation domes in the $(\rho, p)$ plane for the different EoS</figcaption></figure>

EoS, and thus these isotherm curves do not all correspond to the same temperature. Again, the discrepancy between experimental values and CEoS is bigger for Van der Waals and Berthelot EoS, and smaller for Clausius and (Soave)-Redlich-Kwong ones.

### 5.3.3. Internal energy

We can also observe the performance of the cubic laws by plotting the internal energy with respect to $\tau$ for different temperatures. To this end, we chose for the integration constants the experimental values given in [26], namely $e_c = 2.019 \times 10^6 \text{ J kg}^{-1}$, and $c_0 \approx 4.538 \times 10^3 \text{ J kg}^{-1} \text{ K}^{-1}$.
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Internal energy as a function of $\tau$ at three constant temperatures for the Clausius EoS.</figcaption></figure>

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 $\tau_c$ is quite precise, the liquid zone is almost exact, whereas the error on $\tau_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 $\tau$ on an isotherm curve, one notices that it is related to the sign of

$$-\frac{a}{\tau(\tau+d) + c(\tau-d)}$$

which is always non positive (on the whole physical domain for values of $\tau$, 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(\tau, T)$ given by (17) does not diverge when approaching the critical point $(\tau_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.

## 6. 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 $\alpha(T)$, e.g. when choosing $\alpha(T) = T^{\nu}$, $\nu < 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(\tau, s) \stackrel{\text{def}}{=} -\tau^2 \left. \frac{\partial p}{\partial \tau} \right|_s . $$

Let $p(\tau, s) = p(\tau, T(\tau, s))$. By the composition rule, we have

$$ (c^*)^2(\tau, T) = -\tau^2 \left( \left. \frac{\partial p}{\partial \tau} \right|_T + \left. \frac{\partial p}{\partial T} \right|_{\tau} \left. \frac{\partial T}{\partial \tau} \right|_s \right). $$

Using a Maxwell relation [9] and the chain rule, we compute

$$ \begin{aligned} (c^*)^2(\tau, T) &= -\tau^2 \left( \left. \frac{\partial p}{\partial \tau} \right|_T - \left. \frac{\partial p}{\partial T} \right|_{\tau} \left. \frac{\partial p}{\partial s} \right|_{\tau} \right) = -\tau^2 \left( \left. \frac{\partial p}{\partial \tau} \right|_T - \left. \frac{\partial p}{\partial T} \right|_{\tau} \left. \frac{\partial p}{\partial T} \right|_{\tau} \left. \frac{\partial T}{\partial e} \right|_{\tau} \left. \frac{\partial e}{\partial s} \right|_{\tau} \right) \\ &= -\tau^2 \left( \left. \frac{\partial p}{\partial \tau} \right|_T - \left( \left. \frac{\partial p}{\partial T} \right|_{\tau} \right)^2 \left. \frac{\partial T}{\partial e} \right|_{\tau} T \right), \end{aligned} $$

where we used the definition of $T \stackrel{\text{def}}{=} \left. \frac{\partial e}{\partial s} \right|_{\tau}$ in the last equality. We can rewrite this expression of $(c^*)^2$ using $c_v(\tau, T) \stackrel{\text{def}}{=} \left. \frac{\partial e}{\partial T} \right|_{\tau}$ as

$$ (c^*)^2(\tau, T) = -\tau^2 \left( \left. \frac{\partial p}{\partial \tau} \right|_T - T \left. \frac{\left( \left. \frac{\partial p}{\partial T} \right|_{\tau} \right)^2}{\left. \frac{\partial e}{\partial T} \right|_{\tau}} \right) = -\tau^2 \left[ \left. \frac{\partial p}{\partial \tau} \right|_T - \left( \left. \frac{\partial p}{\partial T} \right|_{\tau} \right)^2 \left. \frac{T}{c_v(\tau, T)} \right|_{\tau} \right]. \quad (39) $$

Note that, with Maxwell's construction, we never handle the spinodal zone, and thus, we only consider the EoS in the zone where $\left. \frac{\partial p}{\partial \tau} \right|_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(\tau, 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, $\left. \frac{\partial p}{\partial \tau} \right|_T = 0$, since $p(\tau, T) = p^s(T)$ in the mixture, and $\left. \frac{\partial p}{\partial T} \right|_\tau = (p^s)^\prime(T)$. The functions $(p^s)^\prime(T)$ and $c_v(\tau, T)$ in the mixture are given in section 4.2, which defines all quantities in (39).
<table><caption>TABLE 5. Partial derivatives of the pressure for the different EoS</caption><thead><tr><th></th><th>∂p/∂τ|<sub>T</sub></th><th>∂p/∂T|<sub>τ</sub></th></tr></thead><tbody><tr><td>Van der Waals</td><td>-rT/(τ-b)<sup>2</sup> + 2a/τ<sup>3</sup></td><td>r/τ-b</td></tr><tr><td>Berthelot and Clausius</td><td>-rT/(τ-b)<sup>2</sup> + 2a/T(τ+δ)<sup>3</sup></td><td>r/τ-b + a/T<sup>2</sup>(τ+δ)<sup>2</sup></td></tr><tr><td>Redlich-Kwong</td><td>-rT/(τ-b)<sup>2</sup> + a(2τ+b)/√Tτ<sup>2</sup>(τ+b)<sup>2</sup></td><td>r/τ-b + a/2T√Tτ(τ+b)</td></tr><tr><td>Soave-Redlich-Kwong</td><td>-rT/(τ-b)<sup>2</sup> + a(σ-√T)<sup>2</sup>(2τ+b)/τ<sup>2</sup>(τ+b)<sup>2</sup></td><td>r/τ-b + a(σ-√T)/√Tτ(τ+b)</td></tr></tbody></table>

## B. USE OF THE COMPLETE CUBIC EOS IN A LOW MACH NUMBER MODEL

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 $(\tau, p)$ instead of $(\tau, T)$. For some equations of state, it is possible to invert the formula of the pressure and express the temperature as a function of $(\tau, p)$. In this case, we can insert this expression of $T$ in $e(\tau, T)$ to gain

$$ \tilde{e}(\tau, p) \stackrel{\text{def}}{=} e(\tau, T(\tau, p)), \qquad (40) $$

which further leads to $h(\tau, p) = p\tau + \tilde{e}(\tau, p)$. In the LMNC model [7], an important variable is the derivative of the enthalpy with respect to $\tau$. We can derive $h$ with respect to $\tau$ to obtain:

$$ \zeta(\tau, p) = \left. \frac{\partial h}{\partial \tau} \right|_{p} = p + \left. \frac{\partial \tilde{e}}{\partial \tau} \right|_{p} = p + \left. \frac{\partial e}{\partial \tau} \right|_{T} + \left. \frac{\partial e}{\partial T} \right|_{\tau} \left. \frac{\partial T}{\partial \tau} \right|_{p}. \qquad (41) $$

### 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(\tau, p)$, we can do the same construction as before. In order to obtain $\zeta(\tau, p)$ by (41), we need to compute the two partial derivatives of $e$ given by (35):

$$ \left. \frac{\partial e}{\partial \tau} \right|_T = (T\alpha'(T) - \alpha(T)) \mathfrak{C}'(\tau) = - \frac{a(T\alpha'(T) - \alpha(T))}{\tau(\tau+d) + c(\tau-d)} $$

and

$$ \left. \frac{\partial e}{\partial T} \right|_{\tau} = c_v(\tau, T) = c_{v,0}(T) + T\alpha''(T)(\mathfrak{C}(\tau) - \mathfrak{C}(\tau_c)). $$

**Remark B.1** (Phase transition in complete EoS in $(\tau, p)$). Since the pressure is the same in both fluids in the mixture at saturation, the complete EoS $h_{\text{Maxwell}}(\tau, p)$ is defined piecewise using the saturation values:

$$ h_{\text{Maxwell}}(\tau, p) = \begin{cases} \varphi(\tau, p)h(\tau_{\text{l}}^s(p), T^s(p)) + (1 - \varphi(\tau, p))h(\tau_{\text{g}}^s(p), T^s(p)) & \text{in the mixture, i.e. } \tau_{\text{l}}^s(p) < \tau < \tau_{\text{g}}^s(p), \\ p\tau + \tilde{e}(\tau, p) & \text{otherwise,} \end{cases} $$
so that

$$ \zeta_{\text{Maxwell}}(\tau, p) = \begin{cases} \frac{h(\tau_g^s(p), T_g^s(p)) - h(\tau_l^s(p), T_l^s(p))}{\tau_g^s(p) - \tau_l^s(p)}, & \text{in the mixture, i.e. } \tau_l^s(p) < \tau < \tau_g^s(p) \\ \text{Eq. (41)} & \text{otherwise.} \end{cases} $$

## B.2. Van der Waals

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

$$ T(\tau, p) = \frac{(p\tau^2 + a)(\tau - b)}{r\tau^2}. $$

Inserting this expression in relation (22) for $e$, we deduce $\tilde{e}(\tau, p)$ and thus the enthalpy and its derivative

$$ \zeta(\tau, p) = p + \frac{a}{\tau^2} + c_0 \frac{p\tau^3 + a(2b - \tau)}{r\tau^3}. \quad (42) $$

## B.3. Clausius (and Berthelot)

From the expression (7) of the pressure, we can invert the temperature by solving a second-order equation on $T$ and selecting the positive root, which is given by

$$ T(\tau, p) = \frac{p(\tau - b)}{2r} + \frac{\mathfrak{D}}{2r(\tau + \delta)}, $$

where we defined $\mathfrak{D} \stackrel{\text{def}}{=} \sqrt{(\tau - b)(p^2(\tau + \delta)^2(\tau - b) + 4ar)}$. We can thus compute the partial derivative

$$ \left. \frac{\partial T}{\partial \tau} \right|_p = \frac{p}{2r} - \frac{\mathfrak{D}}{2r(\tau + \delta)^2} + \frac{p^2(\tau + \delta)(\tau - b)(2\tau + \delta - b) + 2ar}{2r(\tau + \delta)\mathfrak{D}}, $$

which allows to compute $\zeta(\tau, T)$ by (41). Note that the Berthelot equation is recovered when $\delta = 0$.

## B.4. Soave-Redlich-Kwong

From 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 by

$$ T(\tau, p) = \frac{-a\sigma(\tau - b) + \mathfrak{D}}{r\tau(\tau + b) - a(\tau - b)}, $$

where we defined $\mathfrak{D} \stackrel{\text{def}}{=} \sqrt{\tau(\tau^2 - b^2)[a(r\sigma^2 - p(\tau - b)) + pr\tau(\tau + b)]}$. We can thus compute the partial derivative

$$ \left. \frac{\partial T}{\partial \tau} \right|_p = \frac{-a\sigma + (\sqrt{\mathfrak{D}})' - T(\tau, p)(2r\tau + rb - a)}{r\tau(\tau + b) - a(\tau - b)}, $$

where the derivative of $\sqrt{\mathfrak{D}}$ with respect to $\tau$ is given by

$$ (\sqrt{\mathfrak{D}})' = \frac{1}{2\sqrt{\mathfrak{D}}} \left\{ (3\tau^2 - b^2) [a(r\sigma^2 - p(\tau - b)) + pr\tau(\tau + b)] + p\tau(\tau^2 - b^2)(2r\tau + rb - a) \right\} $$

which allows to compute $\zeta(\tau, T)$ by (41).
## 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 $\sqrt{T}$:

$$ \frac{r}{\tau - b} \sqrt{T^3} - p\sqrt{T} - \frac{a}{\tau(\tau + b)} = 0. $$

Selecting the real root in Cardano's formula, we obtain

$$ T(\tau, p) = \left( \frac{\sqrt[3]{\frac{2}{3}} p}{\mathfrak{F}} + \frac{\mathfrak{F}}{\sqrt[3]{18\tau}} \right)^2, $$

where we defined

$$ \mathfrak{F} \stackrel{\text{def}}{=} \sqrt[3]{\sqrt{3} \sqrt{27 a^2 r^4 - 4 r^3 p^3} - 9 a r^2}, $$

$$ r \stackrel{\text{def}}{=} \frac{r}{\tau - b}, $$

$$ a \stackrel{\text{def}}{=} \frac{-a}{\tau(\tau + b)}. $$

We can thus compute the partial derivative

$$ \left. \frac{\partial T}{\partial \tau} \right|_p = 2 \sqrt{T(\tau, p)} \left( - \frac{\sqrt[3]{\frac{2}{3}} p \mathfrak{F}'}{\mathfrak{F}^2} + \frac{r \mathfrak{F}' - r' \mathfrak{F}}{\sqrt[3]{18 r^2}} \right), $$

with the following derivatives with respect to $\tau$

$$ \mathfrak{F}' = \mathfrak{F}^2 \left[ \frac{\sqrt{3}(9a a' r^4 + 18a^2 r' r^3 - 2p^3 r' r^2)}{\sqrt{27a^2 r^4 - 4r^3 p^3}} - 3a' r^2 - 6a r' r \right], \quad r' = -\frac{r}{(\tau - b)^2}, \quad a' = \frac{a(2\tau + b)}{\tau^2(\tau + b)^2} $$

which allows to compute $\zeta(\tau, T)$ by (41).

## REFERENCES

[1] R. T. Balmer. *Modern engineering thermodynamics-textbook with tables booklet*. Academic Press, 2011.

[2] D. 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". In: *Journal de Physique Théorique et Appliquée* 8.1 (1899), pp. 263–274.

[3] R. Clausius. "Über das Verhalten der Kohlensäure in Bezug auf Druck, Volumen und Temperatur". In: *Annalen der Physik* 245.3 (1880), pp. 337–357.

[4] W. Dahmen, S. Müller, and A. Vöß. "Riemann Problem for the Euler Equation with Non-Convex Equation of State including Phase Transitions". In: *Analysis and Numerics for Conservation Laws*. Springer Berlin Heidelberg, pp. 137–162. DOI: 10.1007/3-540-27907-5_7.

[5] S. Dellacherie, G. Faccanoni, B. Grec, and Y. Penel. "Accurate steam-water equation of state for two-phase flow LMNC model with phase transition". In: *Applied Mathematical Modelling* 65 (2019), pp. 207–233. DOI: 10.1016/j.apm.2018.07.028.

[6] I. G. Economou. "Cubic and generalized Van der Waals equations of state". In: *Applied Thermodynamics of Fluids* 4.1 (2010), p. 53.

[7] G. Faccanoni, B. Grec, and Y. Penel. "A homogeneous relaxation low Mach number model". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 55 (2021), pp. 1569–1598. DOI: 10.1051/m2an/2021032.

[8] G. Faccanoni, S. Kokh, and G. Allaire. "Modelling and Simulation of Liquid-Vapor Phase Transition in Compressible Flows Based on Thermodynamical Equilibrium". In: *ESAIM Math. Model. Numer. Anal.* 46 (Sept. 2012), pp. 1029–1054. URL: https://hal.archives-ouvertes.fr/hal-00976983.
[9] G. Faccanoni and H. Mathis. "Admissible Equations of State for Immiscible and Miscible Mixtures". In: *ESAIM: Proceedings and Surveys*. Workshop on Compressible Multiphase Flows: Derivation, closure laws, thermodynamics 66 (2019), pp. 1–21. DOI: 10.1051/proc/201966001.

[10] S. Gavrilyuk and H. Gouin. "Theoretical model of the Leidenfrost temperature". In: *Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics* 106.5 (Nov. 2022). DOI: 10.1103/physRevE.106.055.102. URL: https://hal.archives-ouvertes.fr/hal-03617770.

[11] P. N. Ghoderao, V. H. Dalvi, and M. Narayan. "A four parameter cubic equation of state with temperature dependent covolume parameter". In: *Chinese Journal of Chemical Engineering* 27.5 (May 2019), pp. 1132–1148. DOI: 10.1016/j.cjche.2018.08.013.

[12] V. Giovangigli and L. Matuszewski. "Supercritical fluid thermodynamics from equations of state". In: *Physica D: Nonlinear Phenomena* 241.6 (2012), pp. 649–670.

[13] G. M. Kontogeorgis, R. Privat, and J.-N. Jaubert. "Taking Another Look at the Van der Waals Equation of State–Almost 150 Years Later". In: *Journal of Chemical & Engineering Data* 64.11 (Aug. 2019), pp. 4619–4637. DOI: 10.1021/acs.jced.9b00264.

[14] Y. Le Guennec, S. Lasala, R. Privat, and J.-N. Jaubert. "A consistency test for $\alpha$-functions of cubic equations of state". In: *Fluid Phase Equilibria* 427 (2016), pp. 513–538.

[15] Y. Le Guennec, R. Privat, S. Lasala, and J.-N. Jaubert. "On the imperative need to use a consistent $\alpha$-function for the prediction of pure-compound supercritical properties with a cubic equation of state". In: *Fluid Phase Equilibria* 445 (2017), pp. 45–53.

[16] J. S. Lopez-Echeverry, S. Reif-Acherman, and E. Araujo-Lopez. "Peng-Robinson equation of state: 40 years through cubics". In: *Fluid Phase Equilibria* 447 (Sept. 2017), pp. 39–71. DOI: 10.1016/j.fluid.2017.05.007.

[17] E. M. Mansour. "Equation of State". In: *Inverse Heat Conduction and Heat Exchangers*. IntechOpen, Dec. 2020. DOI: 10.5772/intechopen.89919.

[18] J. J. Martin and Y.-C. Hou. "Development of an equation of state for gases". In: *American Institute of Chemical Engineers* (1955).

[19] R. Mortimer. *Physical Chemistry*. Elsevier Science, 2008.

[20] S. Müller and A. Voss. "The Riemann Problem for the Euler Equations with Nonconvex and Nonsmooth Equation of State: Construction of Wave Curves". In: *SIAM Journal on Scientific Computing* 28.2 (Jan. 2006), pp. 651–681. DOI: 10.1137/040619909.

[21] L. Quartapelle, L. Castelletti, A. Guardone, and G. Quaranta. "Solution of the Riemann problem of classical gasdynamics". In: *Journal of Computational Physics* 190.1 (2003), pp. 118–140.

[22] L. Quibel. "Simulation of water-vapor two-phase flows with non-condensable gas." Theses. Université de Strasbourg, Sept. 2020. URL: https://theses.hal.science/tel-03505530v4.

[23] O. Redlich and J. N. Kwong. "On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions." In: *Chemical reviews* 44.1 (1949), pp. 233–244.

[24] G. Soave. "Equilibrium constants from a modified Redlich-Kwong equation of state". In: *Chemical engineering science* 27.6 (1972), pp. 1197–1203.

[25] J. Van der Waals. "Continuity of the gaseous and liquid state of matter". In: *University of Leiden* (1873).

[26] W. Wagner and H.-J. Kretzschmar. "Tables of the Properties of Water and Steam". In: *International Steam Tables: Properties of Water and Steam Based on the Industrial Formulation IAPWS-IF97* (2008), pp. 169–343.

[27] G. Wilczek-Vera and J. H. Vera. "Understanding cubic equations of state: A search for the hidden clues of their success". In: *AIChE Journal* 61.9 (Mar. 2015), pp. 2824–2831. DOI: 10.1002/aic.14741.