# Generalization of the Wong-Sandler mixing rule to a generic cubic equation of state: examples of use for systems of industrial interest (Hydrogen, CCUS, refrigeration)

Salaheddine Chabab, Christophe Coquelet, Fabien Rivollet

► To cite this version:

Salaheddine Chabab, Christophe Coquelet, Fabien Rivollet. Generalization of the Wong-Sandler mixing rule to a generic cubic equation of state: examples of use for systems of industrial interest (Hydrogen, CCUS, refrigeration). Journal of Supercritical Fluids, 2024, 212, pp.106336. 10.1016/j.supflu.2024.106336. hal-04624419

HAL Id: hal-04624419

https://imt-mines-albi.hal.science/hal-04624419v1

Submitted on 3 Jul 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Generalization of the Wong-Sandler mixing rule to a generic cubic equation of state: Examples of use for systems of industrial interest (Hydrogen, CCUS, refrigeration)

Salaheddine Chabab<sup>a,b,*</sup>, Christophe Coquelet<sup>c</sup>, Fabien Rivollet<sup>d</sup>

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

## HIGHLIGHTS

* Wong-Sandler mixing rules has been generalized (g-WS MR) and validated against existing versions (PR and PT EoSs).

* The implementation procedure of g-WS MR is presented.

* 3-parameter EoSs (NEOS, PT and mPT) were tested with g-WS MR for different applications (refrigeration, oil&gas, H₂, CCUS).

* Excellent results are obtained with g-WS MR, in particular the prediction of the critical locus and supercritical densities.

* The g-WS MR can be easily implemented in existing process simulators and thermodynamic calculation software.

## GRAPHICAL ABSTRACT

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

## ARTICLE INFO

**Keywords:**

Critical locus

Phase Equilibria

Density

Refrigerent

Hydrogen

CCUS

## ABSTRACT

Cubic equations of state (EoS) remain the most widely used models in chemical engineering, and their predictive capabilities are highly dependent on mixing rules (MR). Considered as one of the most accurate, Wong Sandler's (WS) MR, have been generalized in this work, to allow its use with all 2-parameter (PR, SRK, ...) and 3-parameter (PT, mPT, NEOS, ...) cubic EoSs. The procedure for implementing the generalized WS MR is presented and tested with PR, PT and its variants (mPT and NEOS) on vapor-liquid equilibrium (VLE) and density data on various classes of mixtures for different types of use (hydrogen, CCUS, refrigeration). The results clearly show the benefits of this type of MR, especially around the critical point and for density data without considering volume translation frequently used with PR and SRK EoS.

<footnote><sup>a</sup> Université de Pau et des Pays de l'Adour, Laboratoire de Thermique, Energetique et Procédés (LaTEP), Rue Jules Ferry BP 7511, Pau Cedex 64075, France</footnote>

<footnote><sup>b</sup> Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR, Pau UMR5150, France</footnote>

<footnote><sup>c</sup> Toulouse University, IMT Mines Albi, UMR CNRS 5302, RAPSODEE Center, Jarlard Campus, Albi Cedex 09 F-81013, France</footnote>

<footnote><sup>d</sup> Arkema, High Performance Polymers, Route du Rilsa, Serquigny 27470, France</footnote>

<footnote>* Corresponding author at: Université de Pau et des Pays de l'Adour, Laboratoire de Thermique, Energetique et Procédés (LaTEP), Rue Jules Ferry BP 7511, Pau Cedex 64075, France.</footnote>

<footnote>E-mail address: salaheddine.chabab@univ-pau.fr (S. Chabab).</footnote>
## Nomenclature

a Parameter of the equation of state (energy parameter, [J.m$^3$.mol$^{-2}$])

A Molar Helmholtz free energy, [J.mol$^{-1}$]

B second virial coefficient

b Parameter of the equation of state (molar co volume parameter, [m$^3$.mol$^{-1}$])

c EoS parameter

d EoS parameter

D Wong-Sandler term

Λ Wong-Sandler parameter (denoted “C” in the original paper of WS)

F Objective function

G Molar Gibbs energy, [J. mol$^{-1}$]

k$_{ij}$ Binary interaction parameter

N Number of components

n Number of mole

P Pressure, [bar]

Q Wong-Sandler term

R Perfect gas constant, [J.mol$^{-1}$.K$^{-1}$]

T Temperature, [K]

ν molar volume, [m$^3$/mol]

x Liquid mole fraction (or composition in a phase)

y Vapor mole fraction

Z Compressibility factor

## Greek letters

α$_{ij}$ NRTL model parameter

τ$_{ij}$ NRTL model binary interaction parameter, dimensionless or [J. mol$^{-1}$]

ω Acentric factor

μ Molar chemical potential, [J.mol$^{-1}$]

Φ Fugacity coefficient

## Superscript

E Excess property

R Residual

* Perfect gas

## Subscripts

C Critical property

cal Calculated property

exp Experimental property

i, j Molecular species

L Liquid

Low Low-pressure reference state

m Mixture

r Reduced

RA Rackett equation

shift Used in volume translation parameter

∞ Infinite pressure reference state

# 1. Introduction

The design and optimization of chemical and energy processes require knowledge of the thermodynamic properties of the fluids (pure or mixtures) involved [1]. For this purpose, suitable thermodynamic models (activity coefficient models and/or equations of state) are generally used. Such models are usually optimized on the basis of experimental phase equilibrium data, allowing control of the fluid states (phase diagrams), the distribution of the various constituents in the different phases, and prediction of other volumetric (density, speed of sound) and energetic (enthalpy, heat capacity) properties [2]. For instance, the simulation of thermodynamic cycles relies on the use of equations of state (EoS), whose parameters are generally adjusted on Vapor-Liquid Equilibrium (VLE) properties alone [3], or on both VLE and energetic properties [1,4]. Moreover, accurate representation of the mixture's critical point is essential to ensure correct sizing of trans-critical thermodynamic cycles.

Cubic EoS are the most widely used EoS in chemical and energy process and oil & gas reservoir simulators, due to their simplicity and speed of computation. The most popular cubic equations of state are Soave-Redlich-Kwong (SRK) [5,6] and Peng Robinson (PR) [7]. They are very useful to estimate the pure component vapor pressure and quite accurate to predict the phase diagrams of mixtures. Numerous researches have been done to improve their capability to calculate the pure component vapor pressure by the development of alpha functions (e.g. Coquelet et al. [8], Le Guennec et al. [9]) and for mixtures, the development of mixing rules which involves or not activity coefficient models. Several advanced mixing rules have been developed extending the applicability of cubic equations of state to the prediction of high-temperature high-pressure VLE for polar systems through the use of excess Gibbs free energy (GE) models. In 1979, Huron and Vidal [10] introduced the excess free energy model into cubic equation of state and consequently, several EoS G$^E$ model appear such as MHV1 [11], MHV2 [12], PSRK [13], LCVM [14] and Wong Sandler [15] mixing rules (WS MR). Wong and Sandler have developed mixing rules that combines an EoS with an excess Gibbs energy model but produces the desired EOS

behavior at both low and high densities without being density dependent. These mixing rules are based on several observations including the second virial development and the pressure insensitivity of the excess Helmholtz energy of mixing.

The main problem of the SRK and PR EoS is that they are inaccurate regarding the density prediction, particularly liquid and supercritical phase densities. Privat and Jaubert [16] have shown that for 2-parameter cubic EoS, Peng Robinson was the most accurate cubic EoS in term of liquid density prediction, particularly at saturation. In order to overcome this problem, two main solutions were adopted.

The first one was introduced by Peneloux [17] and improved by Privat et al. [18] and Le Guennec et al. [19]. It consists of the application of a volume translation parameter $c_{shift}$. Le Guennec et al. have demonstrated that the best results were obtained if $c_{shift}$ is calculated at a reduced temperature equal to 0.8. If no experimental data are available, $c_{shift}$ can be calculated using an improved Rackett equation (Eq. 1).

$$c_{shift} = v_L^{Sat.EoS}(T_r = 0.8) - \frac{RT_c Z_{RA}^{[1+0.2^{2/7}]} }{P_c} \quad (1)$$

The calculation of liquid density is strongly improved using SRK and PR EoS except at the critical point as the value of the critical compressibility factor is fixed. Moreover, the critical compressibility factor value is overestimated and do not correspond to the main compressibility factor values of the different chemical species (Ji and Lempe [20]). In the case of application to mixture, Privat et al. recommend to use a linear mixing rule for the $c_{shift}$ parameter. They also recommend to consider temperature independent volume-translation parameter to avoid inconsistency in the calculation of the thermodynamic properties.

The second solution was the development of 3-parameters cubic EoS. It consists of adding a supplementary parameter properly adjusted to correlate simultaneously the saturated densities and the vapor pressure of pure components ([21]). The main advantage is that the critical compressibility factor is not fixed and improved strongly the calculation of saturated liquid and critical densities (Coquelet et al. [22]). The most popular 3-parameter EoS is the Patel-Teja (PT) EoS [23] and its
generalized form (Patel and Teja [24], Valderrama [25]). Like with SRK and PR EoS, classical mixing rules have been considered for the representation of mixtures phase diagrams and densities. For the 3-parameter cubic EoS classical mixing rules and Huron-Vidal type mixing rules (Ji and Lempe [20]) have been developed and used. In 1997, Yang et al. [26] have applied the WS MR to the 3-parameter Patel-Teja (PT) EoS to improve in the near critical region and for highly asymmetric system. The main difference between the WS MR developed for two and three parameters cubic EoS is that the EoS parameter (denoted “C” in the original paper of WS [15]) is constant for SRK and PR EoS (for example, -0.623225 for the PR EoS) but not for Patel Teja EoS as it depends of Zc value. In this paper a methodology is proposed to extend the WS MR to other EoS and particularly in the perspective of the work of Coquelet et al., 2016 [22] and El Abbadi 2016 [27] who have proposed a modified PT EoS special for refrigerants.

Starting from the generalized expression of the Helmholtz energy for the different cubic EoS, a general expression of the WS MR applicable to the 2- and 3-parameter cubic EoS is proposed in Section 2. In Section 3, the mixing rule is applied to the PR EoS, the PT EoS and two modified versions of the PT EoS (Forero and Velásquez 2013 [28] and Coquelet et al., 2016 [22]). The results obtained are then compared with those obtained with classical (vdW) mixing rules in terms of phase diagram representation particularly close to the critical point and prediction of mixture densities.

## 2. Methodology

Huron-Vidal type mixing rules require parameters dependent on the EoS used and their values are constants. This is no longer the case for equations of state which use 3 or 4 parameters. As with Wong and Sandler, the development of the mixing rules begins with the calculation of the Helmholtz free energy.

### 2.1. Expression of the Helmholtz free energy

The selected generalized cubic EoS (Eq. 2) is the one proposed by Zielke and Lempe 1997 [29], which can be transformed from or to the well-known generalized form proposed by Schmidt and Wenzel 1980 [30] (see Appendix A).

$$P = \frac{RT}{v-b} - \frac{a}{(v+c)(v+d)} \quad (2)$$

The excess Helmholtz free energy is calculated using the following equation [31].

$$\underline{A}^E = \underline{A}_m - \sum_i x_i \underline{A}_i - RT \sum_i x_i \ln(x_i) \quad (3)$$

Moreover, $\underline{A}^R = \underline{A} - \underline{A}^*$, so

$$\underline{A}^E = \underline{A}_m^R + \underline{A}_m^* - \sum_i x_i (\underline{A}_i^R + \underline{A}_i^*) - RT \sum_i x_i \ln(x_i) \quad (4)$$

and then

$$\underline{A}^E = \underline{A}_m^R - \sum_i x_i \underline{A}_i^R \quad (5)$$

From basic thermodynamic relationships and Eq. 2, the residual Helmholtz free energy function can be derived

$$\begin{aligned} \underline{A}^R &= \left( - \int_{v=\infty}^{v} P dv \right) - \left( - \int_{v=\infty}^{v=RT/P} \frac{RT}{v} dv \right) \\ &= -RT \ln \left[ \frac{P(v-b)}{RT} \right] + \frac{a}{d-c} \ln \left( \frac{v+c}{v+d} \right) \end{aligned} \quad (6)$$

Under infinite pressure v=b and Eq. 6 became (if d≠c)

$$\underline{A}_{\infty}^{R} = \frac{a}{d-c} \ln\left(\frac{b+c}{b+d}\right) \quad (7)$$

For mixture

$$\underline{A}_{\infty,m}^{R} = \frac{a_m}{d_m - c_m} \ln\left(\frac{b_m + c_m}{b_m + d_m}\right) \quad (8)$$

Eqs. (7) and (8) can be expressed by

$$\underline{A}_{\infty}^{R} = \frac{a}{b} \Lambda \quad (9)$$

where

$$\Lambda = \frac{b}{d-c} \ln\left(\frac{b+c}{b+d}\right) \text{ if } d \neq c \quad (10)$$

$$\Lambda = -\frac{b}{b+c} \quad \text{if } d=c \quad (11)$$

so Eq. 5 is expressed, under infinite pressure by

$$\underline{A}_{\infty}^{E} = \frac{a_m}{b_m} \Lambda - \sum_{i} x_{i} \frac{a_{i}}{b_{i}} \Lambda_{i} \quad (12)$$

The characteristic of the parameter Λ is analyzed for several EoSs in Appendix A.

### 2.2. The second virial coefficients

The van der Waals type EoS can be explained as a virial series.

$$Z = \frac{Pv}{RT} = 1 + \frac{BP}{RT} \quad (13)$$

where

$$B = b - \frac{a}{RT} \quad (14)$$

In a mixture of N components,

$$B = \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j B_{ij}(T) \quad (15)$$

So

$$b_m - \frac{a_m}{RT} = \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j \left( b - \frac{a}{RT} \right)_{ij} \quad (16)$$

Eq. 16 can be related to those of the pure components by

$$\left(b - \frac{a}{RT}\right)_{ij} = \frac{\left(b - \frac{a}{RT}\right)_i + \left(b - \frac{a}{RT}\right)_j}{2} (1 - k_{ij}) \quad (17)$$

where kij is the second virial coefficient binary interaction parameter.

### 2.3. Definition of the generalized mixing rule

To develop the generalized mixing rules, it has been decided to keep the same definitions as in the original publication by Wong and Sandler 1992 [15] for 2-parameter cubic EoS while adopting simplifications similar to those of Zhao 1995 [32] and Yang et al., 1997 [26] for cubic EoS with three or more parameters. The generalized mixing rules presented herein satisfy the same conditions considered by Wong and Sandler 1992 [15], i.e.: 1) infinite pressure Helmholtz energy is equal to low pressure excess gibbs free energy; 2) quadratic dependence of the second virial coefficient. Using Eqs. (12) and (16), the generalized Wong-Sandler Mixing Rule (g-WS MR) is defined as follows:
$$ 
\begin{aligned}
b_m &= \frac{Q}{1 - D} \\
a_m &= \frac{QDRT}{1 - D}
\end{aligned}
\qquad (18) $$

with

$$ 
\begin{aligned}
Q &= \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j \left(b - \frac{a}{RT}\right)_{ij} \\
D &= \frac{\frac{A_{\infty}^{E}}{RT} + \sum_{i=1}^{N} x_i \frac{a_i}{RT b_i} \Lambda_i}{\Lambda}
\end{aligned}
\qquad (19) $$

To simplify the development, the c and d parameters can be expressed as a function of the co-volume b as:

– for pure compounds:

$$ c_i = b_i \times f_i \text{ and } d_i = b_i \times g_i \qquad (20) $$

e.g. the coefficients f and g for Peng-Robinson EoS are (from Table A1): $f_i = (1 - \sqrt{2})$ and $g_i = (1 + \sqrt{2})$

– for mixtures:

$$ c_m = b_m \times f_m \text{ and } d_m = b_m \times g_m \qquad (21) $$

In this equation, the classical quadratic mixing rules are applied to $f_m$ and $g_m$:

$$ f_m = \sum_i \sum_j x_i x_j \text{ and } g_m = \sum_i \sum_j x_i x_j g_{ij} \qquad (22) $$

where

$$ f_{ij} = \frac{f_i + f_j}{2} \text{ and } g_{ij} = \frac{g_i + g_j}{2} \qquad (23) $$

Parameter $\Lambda$ calculated through Eqs. (10) and (11) is then obtained using Eqs. (24–27):

– for pure compounds:

If $c_i \neq d_i$

$$ \Lambda_i = \frac{1}{g_i - f_i} \ln \left( \frac{1 + f_i}{1 + g_i} \right) \qquad (24) $$

If $c_i = d_i$ (or $f_i = g_i$)

$$ \Lambda_i = \frac{-1}{1 + f_i} \qquad (25) $$

for mixtures:

If $c_m \neq d_m$

$$ \Lambda = \frac{1}{g_m - f_m} \ln \left( \frac{1 + f_m}{1 + g_m} \right) \qquad (26) $$

If $c_m = d_m$ (or $f_m = g_m$)

$$ \Lambda = \frac{-1}{1 + f_m} \qquad (27) $$

The approximation given by Wong-Sandler concerning the quasi-equivalence of $G^E$ and $A^E$ (generally the $V^E$ is very small and the $G^E$ models are generally used at low pressure) is maintained also in this work:

$$ \underline{G}^E(T, x, P = \text{low}) = \underline{A}^E(T, x, P = \text{low}) = \underline{A}^E(T, x, P = \infty) \qquad (28) $$

As explain in Le Guennec’s thesis [33], if we consider a linear mixing rule on co-volume b (like in this work), it is desirable to have $\lim_{P \to \infty} A_{res}^{E, EoS} = \underline{A}_{res}^{E, \gamma}$, with subscript res which corresponds to the residual part. In order to respect this constraint, it is recommended of using activity coefficient model with no combinatorial part. Consequently, the NRTL model developed by Renon and Prausnitz [34] is used in this work.

For the calculation of phase equilibria, the fugacity coefficient obtained from the EoS is used together with conventional equilibrium resolution algorithms (e.g. Michelsen’s [35] or Sandler’s [36] book). The fugacity coefficient (Eq. 30) as previously developed by Zielke and Lempe [29] is obtained by integrating the generalized equation of state (Eq. 2) into the thermodynamic definition of the fugacity coefficient (Eq. 29).

$$ \ln \Phi_i = \int_{v}^{\infty} \left[ \frac{1}{RT} \left( \frac{\partial P}{\partial n_i} \right)_{T, v, n_{i+1}} - \frac{1}{v} \right] dv - \ln Z \qquad (29) $$

So,

$$ \ln \Phi_i = -\ln \left[ \frac{P(v-b)}{RT} \right] + Z - 1 + \frac{n}{v-b} \frac{\partial b}{\partial n_i} - \frac{a+n}{RT(c-d)} \ln \left( \frac{v+c}{v+d} \right) + \frac{a}{RT(c-d)^2} \times \left[ n \frac{\partial c}{\partial n_i} \left( \ln \left( \frac{v+c}{v+d} \right) - \frac{c-d}{v+c} \right) + n \frac{\partial d}{\partial n_i} \left( \frac{c-d}{v+d} - \ln \left( \frac{v+c}{v+d} \right) \right) \right] \qquad (30) $$

The derivatives of the various parameters in Eq. 30 are given in Appendix B.

## 2.4. Calculation procedure

Based on the formulations of the g-WS MR, the following procedure is employed to calculate the mixture parameters $a_m$, $b_m$, $c_m$ and $d_m$, as well as the fugacity coefficients of the generalized EoS (Eq. 2):

i. Set initial conditions (temperature, pressure, mixture composition) and binary interaction parameters ($\tau_{ij}$, $\tau_{ji}$ and $k_{ij}$)

ii. Using Tc, Pc and $\omega$, calculate pure compound parameters $a_i$ and $b_i$ (and $c_i^*$ for 3-parameter EoS) according to the chosen EoS

iii. Compute parameters $c_i$ and $d_i$ (from Table A1 or Zielke and Lempe [29]), $f_i = c_i/b_i$ and $g_i = d_i/b_i$ (Eq. 20), and $\Lambda_i$ (Eqs. 24–25)

iv. Compute $f_{ij}$, $g_{ij}$, $f_m$ and $g_m$ (Eqs. 22–23), and $\Lambda_m$ (Eqs. 26–27)

v. Compute excess Helmholtz free energy using NRTL model (Eq. B6, Appendix B)

vi. Compute Q and D (Eqs. 17 and 19)

vii. Compute $a_m$ and $b_m$ (Equation 18), and $c_m$ and $d_m$ (Eq. 21)

viii. Compute the derivatives (Appendix B) of $a_m$, $b_m$, $c_m$ and $d_m$ with respect to $n_i$ (Eqs. B9-B12) using the following sequence of equations: B22, B23, B20 (or B21), B7, B18, B17, B13-B16, and then B9-B12

ix. Solve the cubic equation (e.g. Cardan method) to obtain the molar volumes $v$ (or compressibility factors) of each phase

x. Calculate the fugacity coefficients of each component in each phase (Eq. 30)

# 3. Results

## 3.1. Validation of the generalized Wong-Sandler mixing rules

In order to validate (by comparison) the development and implementation of the generalized equations (denoted by g-WS MR) presented in the previous section, more or less "asymmetric" systems modeled in the literature by two-parameter (Peng-Robinson) or three-parameter (Patel-Teja) cubic equations of state with WS rules were considered.

Regarding the Peng-Robinson EoS with WS mixing rules (PR-WS), the parameters optimized recently by Yang et al., 2019 [37] for the
<figure><img src="image_5.png" /><figcaption>Fig. 1. Validation of the generalized Wong-Sandler mixing rules using Peng-Robinson EoS on the VLE of the CO₂+H₂O system at 323 K. Literature modeling data from: Yang et al. 2019 [37].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2. Validation of the generalized Wong-Sandler mixing rules using Patel-Teja EoS on the VLE of the N₂-CO₂ (left) and CO₂+H₂O (right) systems. Literature modeling data from: Yang 1997 [26].</figcaption></figure>

CO₂-H₂O system were used for comparison. These parameters were used with the original PR-WS model [15] as well as the generalized version proposed in this work. As shown in Fig. 1, the VLE results are perfectly reproduced by the generalized equations, making them valid for 2-parameter cubic EoS.

Similarly, the 3-parameter Patel-Teja (PT) EoS with g-WS MR was evaluated. The parameters proposed by Yang 1997 [26] for the CO₂-H₂O and N₂-CO₂ systems were used with the g-WS MR and compared with the original PT-WS results (Yang 1997 [26]). Fig. 2 illustrates the very similar VLE results obtained with the original PT-WS model and the generalized mixing rules. The very small deviations are due to the fact that the data calculated by Yang 1997 (PT-WS) and shown in Fig. 2 are not "real" data but rather have been extracted from their graphs.

## 3.2. Applications and discussion

In this section, the g-WS MR has been tested on VLE and density data for various classes of mixtures for different types of use (hydrogen, CCUS, refrigeration). We have tested the mixing rules on binary system where both VLE and density data are available. The NRTL parameter $\alpha_{ji}$ was set equal to 0.3 in all cases examined. When available, $\tau_{ij}$, $\tau_{ji}$ and $k_{ij}$ parameters were taken from the literature. Otherwise, the optimal values for NRTL parameters $\tau_{12}$ and $\tau_{21}$ were determined simultaneously with the binary interaction parameter $k_{ij}$ from binary VLE data using the objective function:

$$F = \frac{100}{N_{\text{exp}}} \left[ \sum_{1}^{N_{\text{exp}}} \left( \frac{x_{\text{exp}} - x_{\text{cal}}}{x_{\text{exp}}} \right)^2 + \left( \frac{y_{\text{exp}} - y_{\text{cal}}}{y_{\text{exp}}} \right)^2 \right] \quad (31)$$
<figure><img src="image_7.png" /><figcaption>Fig. 3. VLE (left) and density (right) of the CO₂-R1234yf system: Dotted lines: NEOS with vdW MR; dashed lines: GERG-2008 EoS (REFPROP); solid lines: NEOS with g-WS MR.</figcaption></figure>

where $N_{exp}$ is the number of data points, and $x_{exp}$ (or $x_{cal}$) and $y_{exp}$ (or $y_{cal}$) are the experimental (or calculated) liquid and vapor compositions, respectively.

The first system tested is the refrigerant mixture CO₂-R1234yf, which is of great industrial interest with zero ozone depletion potential and very low global warming potential. To model this binary system, the 3-parameter cubic EoS (NEOS) developed by Coquelet et al., 2016 [22] for refrigerant systems was selected. The parameters of the pure compounds (CO₂ and R1234yf) provided in their work were used. This binary system was modeled using the g-WS MR, and also the classical van der Waals (vdW) mixing rules using the binary interaction parameters ($k_{ij}$) from El Abbadi 2016 [27]. The binary interaction parameters ($\tau_{12}$, $\tau_{21}$ and $k_{12}$) of the g-WS MR were fitted to the VLE data from Juntarachat et al., 2014 [38], and the results are shown in Fig. 3. The VLE data are perfectly reproduced by the g-WS MR (ARD(x)=0.6% and ARD(y)=1.4%; with

$$ARD = \frac{100}{N_{exp}} \sum_{i=1}^{N_{exp}} \left| \frac{x_{(or y)}^{cal} - x_{(or y)}^{exp}}{x_{(or y)}^{exp}} \right| \text{ is the Average Relative Deviation) } $$

and significantly improve the results when compared with the vdW MR (ARD(x)=1.9% and ARD(y)=5.2%). In the same figure, the prediction results of the mixture density as a function of pressure and at different compositions are compared with experimental data from the literature [39] and with the GERG-2008 multiparameter model (available in

REFPROP [40]). The prediction results using the g-WS MR (ARD(ρ)<0.7%) are very comparable to those of GERG-2008 under the different conditions and without any direct adjustment on density.

The second system studied is the H₂-CO₂ binary mixture, which is of particular interest for several industrial applications (production of hydrogen carrier molecules [41] and e-fuels [42], H₂ and CO₂ separations in syngas processes [43], CCUS and underground H₂ storage [44], etc.). To model this binary system, the three-parameter Patel-Teja [24] cubic EoS was selected. The original version of Patel-Teja [24] poorly describes the vapor densities and enthalpy of vaporization of H₂, hence the modified version by Forero and Velásquez [28] (denoted by mPT), which solves this problem by using an exponential alpha function, was used. In addition, the Peng-Robinson EoS with classical mixing rules (denoted by PR-vdW) and with g-WS MR (denoted by PR-WS) are used for comparison with the mPT+g-WS MR model (denoted by mPT-WS). The $k_{12}$ for vdW MR as well as the binary interaction parameters of g-WS MR ($\tau_{12}$, $\tau_{21}$ and $k_{12}$) were optimized on experimental VLE data from Fandino et al. [43] and Tsang and Street [45], and the results are shown in Fig. 4. In contrast to the PR-vdW model (ARD(x)=12.9% and ARD(y)=9.5%), the PR-WS model (ARD(x)=5.3% and ARD(y)=2.2%), and even better the mPT-WS model (ARD(x)=4.5% and ARD(y)=1.5%), reproduce the phase diagram perfectly at different temperatures

<figure><img src="image_8.png" /><figcaption>Fig. 4. VLE (left) and density (right) of the CO₂-H₂ system: Dotted lines: Peng-Robinson with vdW MR; dashed lines: Peng-Robinson with g-WS MR; solid lines: modified Patel-Teja with g-WS MR.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** VLE (left) and density (right) of the CH₄-CO₂ system: Dotted lines: Peng-Robinson with vdW MR; dashed lines: Peng-Robinson with g-WS MR; solid lines: original Patel-Teja with g-WS MR.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Density of the ternary 52.32 %CH₄-31.48 %CO₂-16.20 %N₂ system: Dotted lines: original Patel-Teja with vdW MR; dashed lines: GERG-2008; solid lines: original Patel-Teja with g-WS MR.</figcaption></figure>

and pressures, and across the entire composition range, including near the mixture's critical point. Furthermore, the mPT-WS model accurately (ARD(ρ)<0.3 %) predicts the density of the CO₂-H₂ mixture at different conditions, with an improvement (at high pressure) over the other two models (PR-vdW and PR-WS).

The third system studied is the CH₄-CO₂ binary, which is the main constituent of natural gas and biogas, and has a wide range of uses. The three-parameter cubic equation of state tested for this system is the original Patel-Teja model with g-WS MR (denoted as PT-WS). Again, for comparison, the PR-WS and PR-vdW models are also used. The binary interaction parameters of the PR-vdW and PR-WS models are taken from Nazeri et al. [46] and Theveneau et al. [47], respectively. For the PT-WS model, the binary interaction parameters (τ₁₂, τ₂₁ and k₁₂) were optimized in this work on literature experimental VLE data [48–50], and the results are shown in Fig. 5. Compared with PR-vdW (ARD(x)=4.7 % and ARD(y)=1.9 %) and PR-WS (ARD(x)=9.6 % and ARD(y)=3.9 %), the PT-WS model (ARD(x)=2.5 % and ARD(y)=2.6 %) gives the best representation of the VLE data at different temperatures and pressures, and over almost the entire range of liquid and vapor compositions. Density

predictions are also considerably improved by the PT-WS model (ARD (ρ)<0.6 %) for pressures up to 1000 bar and at different mixture compositions (20 % and 90 % CH₄ in CO₂).

Moreover, Fig. 6 illustrates the results of the density prediction of the CH₄-CO₂-nC₄ ternary mixture with the PT-vdW, PT-WS and GERG-2008 models. The binary interaction parameters of the PT-vdW and PT-WS models were taken from Pan et al. [51] and Yang et al. [26], respectively. At low temperature (294.55 K), the PT-WS model is the only one of the three capable of perfectly (ARD(ρ)<0.7 %) reproducing densities over the entire pressure range, while the PT-vdW and GERG-2008 models reproduce density with an ARD(ρ) of 4.9 % and 7 %, respectively. At higher temperature (333.15 K), the same model correctly captures the mixture density up to 150 bar, then diverges at higher pressure (ARD(ρ)<4.3 %), unlike the other two models (ARD(ρ)<2.8 %). This can be due to a problem with the temperature dependency of the BIP especially at high temperature.

Finally, to get a better view of the performance of the g-WS MR at the critical point, the VLE critical locus of two binary systems were computed. The procedure used for calculating mixtures' critical point lines is that developed by Heidemann and Khalil [52] and Michelsen and Heidemann [53] and generalized by Stockfleth and Dohrn [54]. The latter proposed to calculate the first and second derivatives of the fugacity with respect to the mole number rather numerically (four-point differencing scheme) than analytically, to facilitate access to this method for all types of EoS and mixing rules. Fig. 7 shows the results of calculating the VLE critical locus of the CO₂-R1234yf and CH₄-CO₂ systems using the different EoS+MR combinations previously used to model their VLEs. For the CO₂-R1234yf (nonpolar–polar binary) system, in contrast to the results obtained with the NEOS model combined with classical MR (ARD(T)<1.6 % and ARD(P)<4.1 %), the critical locus is perfectly reproduced by the NEOS model combined with g-WS MR (ARD(T)<0.5 % and ARD(P)<2.3 %) with temperature-independent binary interaction parameters. For the CH₄-CO₂ system, the difference in results is less significant (as these are non-polar binary systems), with a slightly better result from the EoSs (PR and PT) combined with g-WS MR (ARD(T and P)<2 %).

In terms of predicting mixtures' critical curves, to date, the use of cubic EoS but with advanced mixing rules like WS remains the most simple and appropriate, especially if compared with Helmholtz energy-explicit EoS like GERG-2008 which are excellent for pure compounds but can give unphysical critical curves as mentioned recently in 2020 by the works of Deiters and Bell [55,56].
<figure><img src="image_11.png" /><figcaption>Fig. 7. VLE critical locus (P-T diagram) of CO₂-R1234yf (left) and CH₄-CO₂ (right) binary systems. Solid lines: EoS with g-WS MR; Dotted lines: EoS with vdW MR; Dashed lines: PR with g-WS MR.</figcaption></figure>

## 4. Conclusions

Cubic EoSs are still the most widely used models because of their simplicity and speed of calculation. While their use generally leads to reasonable results, the representation of particularly liquid and supercritical phase densities of pure compounds is not precise. Moreover, for asymmetrical mixtures, and in some cases even for symmetrical mixtures, they fail to reproduce accurately the VLE, especially around the critical point. These weaknesses have increased the interest in using 3-parameter cubic EoS with advanced mixing rules combining the advantages of activity coefficient (G⁵) models and EoSs.

In this work, the well-known Wong-Sandler mixing rules have been generalized (g-WS MR), making them usable with all existing 2- and 3-parameter cubic EoSs. The expressions, their derivatives and the sequence of calculations are presented. The g-WS MR developed has been validated by comparing and obtaining exactly the same results as the two existing versions in the literature (PR-WS and PT-WS).

The g-WS MRs were then tested on PR EoS and three 3-parameter cubic EoSs: original Patel-Teja (PT), its modified version (mPT) by Forero and Velásquez 2013, and the EoS recently developed by Coquelet et al., 2016 (NEOS) for refrigerants. VLE and densities of the binary mixtures CO₂-R1234yf, CO₂-H₂, and CH₄-CO₂ and the density of the ternary mixture CH₄-CO₂-nC₄ were studied. The findings clearly demonstrate the advantages of employing such mixing rules, particularly near the critical point where traditional ones often fall short. The results show that the best compromise for the representation of VLE and mixture density at near-critical and supercritical conditions, respectively, is the use of three-parameter cubic EoS with g-WS MR, followed by two-parameter cubic EoS with g-WS MR and finally EoS with classical vdW mixing rules. These types of calculations are essential for oil&gas applications, modelling of trans critical fluid heat pump cycles or

supercritical fluid extraction.

In the perspective of this work, it would be interesting to extensively investigate mixtures for different industrial uses (e.g. refrigerant mixtures by NEOS-WS, or H₂ and CO₂ mixtures by mPT-WS) in order to provide binary interaction parameters, which could encourage the use of three-parameter cubic EoS with advanced mixing rules like WS.

## CRediT authorship contribution statement

Salaheddine Chabab: Writing – original draft, Supervision, Software, Methodology, Conceptualization. Christophe Coquelet: Writing – review & editing, Validation, Supervision, Investigation, Conceptualization. Fabien Rivollet: Writing – review & editing, Visualization, Validation, Investigation.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data Availability

Data will be made available on request.

## Acknowledgments

Financial support from the "HYDR" project partners is gratefully acknowledged. The S.C. "HYDR" Chair is funded by the ISITE E2S program, supported by ProSim SA, ANR PIA, IFPEN and Région Nouvelle-Aquitaine.

## Appendix A. Characteristic of parameter Λ

Eq. 2 can be transformed into Schmidt and Wenzel's [30] well-known generic cubic EoS:

$$P = \frac{RT}{v-b} - \frac{a}{v^2 + ubv + wb^2} \quad (A1)$$

with
$$ 
\begin{aligned}
u &= \frac{c+d}{b} \\
w &= \frac{cd}{b^2}
\end{aligned}
\qquad (A2) $$

The parameter $\Lambda$ (denoted “C” in the original paper of WS), which does not have any physical dimension, can be expressed as a function of critical compressibility factor for different types of EoS. To calculate this parameter using Eqs. 10 or 11, the values of c and d parameters as a function of the EoS are listed in Table A1 for a few EoS selected in this work (for other EoS please refer to the paper of Zielke and Lempe 1997 [29]).

<caption>Table A1</caption>

Values of parameters c and d in function of 5 selected EoS.

<table><thead><tr><th>EoS</th><th>Parameter c</th><th>Parameter d</th></tr></thead><tbody><tr><td>van der Waals</td><td>0</td><td>0</td></tr><tr><td>SRK</td><td>0</td><td>b</td></tr><tr><td>PR</td><td>b(1 - √2)</td><td>b(1 + √2)</td></tr><tr><td>Patel Teja</td><td>b + c* / 2 + √(c* b + (b + c*)² / 2)</td><td>b + c* / 2 - √(c* b + (b + c*)² / 2)</td></tr><tr><td>Coquelet et al.</td><td>b(1 + c* / b - √(1 + c* / b)(5 + c* / b))</td><td>b(1 + c* / b + √(1 + c* / b)(5 + c* / b))</td></tr></tbody></table>

Parameters b and c* in Table A1 are expressed as follows:

$$ 
\begin{aligned}
b &= \Omega_b \frac{RT_c}{P_c} \\
c^* &= \Omega_c \frac{RT_c}{P_c}
\end{aligned}
\qquad (A3) $$

with

$$ 
\begin{aligned}
u\Omega_b &= 1 + \Omega_b - 3Z_c \\
\Omega_b^3 + ((1 - 3Z_c) + (u + w))\Omega_b^2 + 3Z_c^2\Omega_b - Z_c^3 &= 0
\end{aligned}
\qquad (A4) $$

For different EoS and for different values of $Z_c$, some values of parameter $\Lambda$ have been calculated. It is also possible to classify the different EoS as a function of the values of u and w [57] as mentioned in Eq. A1. Table A2 presents the list of some published equations and the corresponding u and w parameters. Different groups of EoS can be distinguished: for example, Patel Teja [24], Harmens and Knapp [58], Heyen [59], Schmidt and Wenzel [30], SRK [5,6] and PR [7] are classified with EoS where u+w=1. Van der Waals [60] is classified with EoS where u+w=0. As the $\Lambda$ parameter is not constant, unlike vdW, PR and SRK EoS, Fig. A1 presents the behavior of the parameter $\Lambda$ as a function of $Z_c$. A remark can be made about Harmens and Knapp's EoS [58], being the only one that is classified with EoS where u+w=1 and whose parameter $\Lambda$ does not respect that as illustrated in Figure A1. The main reason is probably due its third parameter c which does not have any physical sense and is unitless contrary to the one of Patel and Teja EoS.

<caption>Table A2</caption>

Equations of state and their corresponding parameters (Eq. A1)

<table><thead><tr><th>EoS</th><th>Parameters (u and w)</th></tr></thead><tbody><tr><td>van der Waals [60]</td><td>u=0; w=0</td></tr><tr><td>Redlich and Kwong [5,6]</td><td>u=1; w=0</td></tr><tr><td>Peng and Robinson [7]</td><td>u=2; w=-1</td></tr><tr><td>Patel and Teja [24], Harmens and Knapp [58], Heyen [59], Schmidt and Wenzel [30]</td><td>u+w=1</td></tr><tr><td>Coquelet et al. [22]</td><td>u+w=0</td></tr><tr><td>Yu and Lu [61]</td><td>u-w=3</td></tr><tr><td>Twu et al. [62]</td><td>u-w=4</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Values of parameter Λ as a function of critical compressibility factor Zc for the different categories of EoS. Dash-dotted line: u+w=1; dotted line: u+w=0. van der Waals (u+w=0): Λ = -1; PR (u+w=1): Λ = -1/√2 ln(√2 + 1); SRK (u+w=1): Λ = -ln(2).</figcaption></figure>

For 2-parameter cubic EoS (vdW, SRK, PR, etc.) the parameters Λ and Λᵢ are equal and constant. However, for cubic equations of state with more than two parameters, these two parameters are not necessarily equal and constant. Figure A1 shows also that it is possible to develop EoS-specific correlations of Λ=f(Zc). In this study, we have preferred to calculate the Λ parameter for each case study using Eqs. 10 or 11.

## Appendix B. Expressions and derivatives of g-WS MR

The calculation of parameters a and b for the mixture is done with

$$ 
\begin{aligned}
b_m &= \frac{Q}{1-D} \\
a_m &= \frac{QDRT}{1-D}
\end{aligned} 
\qquad (B1) $$

Where Q and D are defined as follows

$$ 
\begin{aligned}
Q &= \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j \left(b - \frac{a}{RT}\right)_{ij} \\
D &= \frac{\frac{A_{\infty}^{E}}{RT} + \sum_{i=1}^{N} x_i \frac{a_i}{RT b_i} \Lambda_i}{\Lambda}
\end{aligned} 
\qquad (B2) $$

The calculation of parameters c and d for the mixture is done with

$$ c_m = b_m \times f_m \quad \text{and} \quad d_m = b_m \times g_m \qquad (B3) $$

In this equation, the classical mixing rules are applied to f_m and g_m:

$$ f_m = \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j f_{ij} \quad \text{and} \quad g_m = \sum_{i=1}^{N} \sum_{j=1}^{N} x_i x_j g_{ij} \qquad (B4) $$

Where

$$ f_{ij} = \frac{f_i + f_j}{2} \quad \text{and} \quad g_{ij} = \frac{g_i + g_j}{2} \qquad (B5) $$

The coefficients fᵢ and gᵢ are specific to the selected cubic EoS (see Table A1, or Zielke and Lempe 1997 [29]); e.g. for Peng-Robinson EoS: fᵢ = (1 - √2) and gᵢ = (1 + √2).

The expressions for the excess Helmholtz free energy and the corresponding activity coefficient of the NRTL model (Renon and Prausnitz [34]) are as follows:
$$ \frac{A_{\infty}^{E}}{RT} = \sum_{i} x_i \left( \frac{\sum_{j} x_j \tau_{ji} g_{ji}}{\sum_{k} x_k g_{ki}} \right) \qquad (B6) $$

$$ \ln \gamma_{\infty,i} = \frac{\sum_{j} x_j \tau_{ji} g_{ji}}{\sum_{k} x_k g_{ki}} + \sum_{j} \frac{x_j g_{ij}}{\sum_{k} x_k g_{kj}} \left( \tau_{ij} - \frac{\sum_{l} x_l \tau_{lj} g_{lj}}{\sum_{k} x_k g_{kj}} \right) \qquad (B7) $$

with

$$ g_{ij} = \exp(-\alpha_{ij} \tau_{ij}) \qquad (B8) $$

The partial derivatives in Eq. 30 were derived analytically and validated using Maple software and numerical symmetric derivatives. Here are the developed analytical expressions:

$$ n \frac{\partial b_m}{\partial n_i} = \frac{\partial n b_m}{\partial n_i} - b_m \qquad (B9) $$

$$ n \frac{\partial a_m}{\partial n_i} = \frac{\partial n a_m}{\partial n_i} - a_m \qquad (B10) $$

$$ n \frac{\partial c_m}{\partial n_i} = \frac{\partial n c_m}{\partial n_i} - c_m \qquad (B11) $$

$$ n \frac{\partial d_m}{\partial n_i} = \frac{\partial n d_m}{\partial n_i} - d_m \qquad (B12) $$

Then,

$$ \frac{\partial n b_m}{\partial n_i} = \frac{1}{1-D} \left( \frac{1}{n} \frac{\partial (n^2 Q)}{\partial n_i} \right) - \frac{Q}{(1-D)^2} \left( 1 - \frac{\partial (nD)}{\partial n_i} \right) \qquad (B13) $$

$$ \frac{\partial n a_m}{\partial n_i} = RTD \times \left( n \frac{\partial b_m}{\partial n_i} - n \frac{\partial \Lambda}{\partial n_i} \times \frac{b_m}{\Lambda} \right) + \frac{b_m}{\Lambda} \times \left( \frac{a_i}{b_i} \Lambda_i + RT \times \ln \gamma_{\infty,i} \right) \qquad (B14) $$

$$ \frac{\partial n c_m}{\partial n_i} = \frac{\partial (n f_m b_m)}{\partial n_i} = f_m b_m + b_m n \frac{\partial f_m}{\partial n_i} + f_m n \frac{\partial b_m}{\partial n_i} \qquad (B15) $$

$$ \frac{\partial n d_m}{\partial n_i} = \frac{\partial (n g_m b_m)}{\partial n_i} = g_m b_m + b_m n \frac{\partial g_m}{\partial n_i} + g_m n \frac{\partial b_m}{\partial n_i} \qquad (B16) $$

with

$$ \left( \frac{1}{n} \frac{\partial (n^2 Q)}{\partial n_i} \right) = 2 \sum_{j}^{N} x_j \left( b - \frac{a}{RT} \right)_{ij} \qquad (B17) $$

$$ \frac{\partial (nD)}{\partial n_i} = \frac{1}{\Lambda RT} \left( \frac{a_i}{b_i} \Lambda_i + RT \times \ln \gamma_{\infty,i} - n \frac{\partial \Lambda}{\partial n_i} \times \frac{\left( \sum_{l}^{N} \left( x_{l} \frac{a_{l}}{b_{l}} \Lambda_{l} \right) + A_{\infty}^{E} \right)}{\Lambda} \right) \qquad (B18) $$

$$ \frac{\partial n A_{\infty}^{E}}{\partial n_i} = RT \times \ln \gamma_{\infty,i} \qquad (B19) $$

and from Eqs. 26 and 27:

if $d_m \neq c_m$

$$ \left( n \frac{\partial \Lambda}{\partial n_i} \right) = - \frac{n \times \left( \frac{\partial g_m}{\partial n_i} - \frac{\partial f_m}{\partial n_i} \right)}{(g_m - f_m)^2} \ln \left( \frac{1 + f_m}{1 + g_m} \right) + \frac{1}{f_m - g_m} \times n \times \left( \frac{\frac{\partial g_m}{\partial n_i}}{1 + g_m} - \frac{\frac{\partial f_m}{\partial n_i}}{1 + f_m} \right) \qquad (B20) $$

If $d_m = c_m$ ($\text{orf}_m = g_m$)

$$ \left( n \frac{\partial \Lambda}{\partial n_i} \right) = \frac{n \frac{\partial f_m}{\partial n_i}}{(1 + f_m)^2} \qquad (B21) $$

where
$$n \frac{\partial f_m}{\partial n_i} = 2 \times \left( \sum_j (x_j \times f_{ij}) - f_m \right) \qquad (B22)$$

$$n \frac{\partial g_m}{\partial n_i} = 2 \times \left( \sum_j (x_j \times g_{ij}) - g_m \right) \qquad (B23)$$

## References

[1] C. Coquelet, A. Valtz, P. Théveneau, Experimental determination of thermophysical properties of working fluids for ORC applications, Org. Rank. Cycles Waste Heat. Recovery-Anal. Appl. (2019).

[2] C. Feng, C. Guo, J. Chen, S. Tan, Y. Jiang, Thermodynamic analysis of a dual-pressure evaporation high-temperature heat pump with low GWP zeotropic mixtures for steam generation, Energy 294 (2024) 130964.

[3] V. Koulocheris, M. Koliva, V. Loui, E. Voutsas, Modelling of pure refrigerant thermodynamic properties and vapor-liquid equilibrium of refrigerant mixtures with the UMR-PRU model, Fluid Phase Equilib. (2024) 114052.

[4] S. Lasala, A.-P. Martinez, J.-N. Jaubert, A predictive equation of state to perform an extending screening of working fluids for power and refrigeration cycles, Org. Rank. Cycles Waste Heat. Recovery-Anal. Appl. (2020).

[5] O. Redlich, J.N. Kwong, On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions, Chem. Rev. 44 (1949) 233-244.

[6] G. Soave, Equilibrium constants from a modified Redlich-Kwong equation of state, Chem. Eng. Sci. 27 (1972) 1197-1203.

[7] D.-Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundam. 15 (1976) 59-64.

[8] C. Coquelet, A. Chapoy, D. Richon, Development of a new alpha function for the Peng-Robinson equation of state: comparative study of alpha function models for pure gases (natural gas components) and water-gas systems, Int. J. Thermophys. 25 (2004) 133-158.

[9] Y. Le Guennec, S. Lasala, R. Privat, J.-N. Jaubert, A consistency test for α-functions of cubic equations of state, Fluid Phase Equilib. 427 (2016) 513-538.

[10] M.-J. Huron, J. Vidal, New mixing rules in simple equations of state for representing vapour-liquid equilibria of strongly non-ideal mixtures, Fluid Phase Equilib. 3 (1979) 255-271.

[11] M.L. Michelsen, A modified Huron-Vidal mixing rule for cubic equations of state, Fluid Phase Equilib. 60 (1990) 213-219.

[12] S. Dahl, M.L. Michelsen, High-pressure vapor-liquid equilibrium with a UNIFAC-based equation of state, AIChE J. 36 (1990) 1829-1836.

[13] T. Holderbaum, J. Gmehling, PSRK: a group contribution equation of state based on UNIFAC, Fluid Phase Equilib. 70 (1991) 251-265.

[14] C. Boukouvalas, N. Spiliotis, P. Coutsikos, N. Tzouvaras, D. Tassios, Prediction of vapor-liquid equilibrium with the LCVM model: a linear combination of the Vidal and Michelsen mixing rules coupled with the original UNIF, Fluid Phase Equilib. 92 (1994) 75-106.

[15] D.S.H. Wong, S.I. Sandler, A theoretically correct mixing rule for cubic equations of state, AIChE J. 38 (1992) 671-680.

[16] R. Privat, J.-N. Jaubert, The state of the art of cubic equations of state with temperature-dependent binary interaction coefficients: from correlation to prediction, Fluid Phase Equilib. 567 (2023) 113697.

[17] A. Péneloux, E. Rauzy, F. Fréze, A consistent correction for Redlich-Kwong-Soave volumes, Fluid Phase Equilib. 8 (1982) 7-23.

[18] R. Privat, J.-N. Jaubert, Y. Le Guennec, Incorporation of a volume translation in an equation of state for fluid mixtures: which combining rule? Which effect on properties of mixing? Fluid Phase Equilib. 427 (2016) 414-420.

[19] Y. Le Guennec, R. Privat, J.-N. Jaubert, Development of the translated-consistent tc-PR and tc-RK cubic equations of state for a safe and accurate prediction of volumetric, energetic and saturation properties of pure compounds in the sub-and super-critical domains, Fluid Phase Equilib. 429 (2016) 301-312.

[20] W.-R. Ji, D. Lempe, A systematic study of cubic three-parameter equations of state for deriving a structurally optimized PVT relation, Fluid Phase Equilib. 147 (1998) 85-103.

[21] H. Segura, D. Seiltgens, A. Mejía, F. Lovell, L.F. Vega, An accurate direct technique for parameterizing cubic equations of state: Part II. Specializing models for predicting vapor pressures and phase densities, Fluid Phase Equilib. 265 (2008) 155-172.

[22] C. Coquelet, J. El Abbadi, C. Houriez, Prediction of thermodynamic properties of refrigerant fluids with a new three-parameter cubic equation of state, Int. J. Refrig. 69 (2016) 418-436.

[23] P.N. Ghoderao, M. Narayan, V.H. Dalvi, H.S. Byun, Patel-Teja cubic equation of state-a review of modifications and applications till 2022, Fluid Phase Equilib. 567 (2023) 113707.

[24] N.C. Patel, A.S. Teja, A new cubic equation of state for fluids and fluid mixtures, Chem. Eng. Sci. 37 (1982) 463-473.

[25] J.O. Valderrama, A generalized Patel-Teja equation of state for polar and nonpolar fluids and their mixtures, J. Chem. Eng. Jpn. 23 (1990) 87-91.

[26] T. Yang, G.-J. Chen, W. Yan, T.-M. Guo, Extension of the Wong-Sandler mixing rule to the three-parameter Patel-Teja equation of state: application up to the near-critical region, Chem. Eng. J. 67 (1997) 27-36.

[27] J. El Abbadi, Etude des propriétés thermodynamiques des nouveaux fluides frigorigènes, in, PSL Research University, 2016.

[28] L.A. Forero G, J.A. Velasquez J, A modified Patel-Teja cubic equation of state: part I-generalized model for gases and hydrocarbons, Fluid Phase Equilib. 342 (2013) 8-22.

[29] F. Zielke, D. Lempe, Generalized calculation of phase equilibria by using cubic equations of state, Fluid Phase Equilib. 141 (1997) 63-85.

[30] G. Schmidt, H. Wenzel, A modified van der Waals type equation of state, Chem. Eng. Sci. 35 (1980) 1503-1512.

[31] J. Smith, H. Van Ness, M. Abbott, M. Swihart. Introduction to Chemical Engineering Thermodynamics, 8th Ed, McGraw Hill, 2018.

[32] E. Zhao, Equations of state for nonelectrolyte and electrolyte solutions, University of Ottawa (Canada), 1995.

[33] Y.Le Guennec, Développement d'équations d'état cubiques adaptées à la représentation de mélanges contenant des molécules polaires (eau, alcools, amines...) et des hydrocarbures, in, Université de Lorraine, 2018.

[34] H. Renon, J.M. Prausnitz, Local compositions in thermodynamic excess functions for liquid mixtures, AIChE J. 14 (1968) 135-144.

[35] M.L. Michelsen, J. Mollerup, Thermodynamic Modelling: Fundamentals and Computational Aspects, Tie-Line Publications, 2004.

[36] S.I. Sandler, Chemical, Biochemical, and Engineering Thermodynamics, John Wiley & Sons, 2017.

[37] Z.-L. Yang, H.-Y. Yu, Z.-W. Chen, S.-Q. Cheng, J.-Z. Su, A compositional model for CO2 flooding including CO2 equilibria between water and oil using the Peng-Robinson equation of state with the Wong-Sandler mixing rule, Pet. Sci. 16 (2019) 874-889.

[38] N. Juntarachat, A. Valtz, C. Coquelet, R. Privat, J.-N. Jaubert, Experimental measurements and correlation of vapor-liquid equilibrium and critical data for the CO2+ R1234yf and CO2+ R1234ze (E) binary mixtures, Int. J. Refrig. 47 (2014) 141-152.

[39] S. Peng, E. Wang, K. Qing, Z. Yang, Y. Duan, Gaseous pvTx properties for binary and ternary mixtures of R1234yf, R32, and CO2, J. Chem. Eng. Data 68 (2023) 1624-1637.

[40] E. Lemmon, I.H. Bell, M. Huber, M. McLinden, NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0, National Institute of Standards and Technology, (2018).

[41] M. Niermann, A. Beckendorff, M. Kaltschmitt, K. Bonhoff, Liquid organic hydrogen carrier (LOHC)-assessment based on chemical and economic properties, Int. J. Hydrog. Energy 44 (2019) 6631-6654.

[42] F. Ueckerdt, C. Bauer, A. Dirnaichner, J. Everall, R. Sacchi, G. Luderer, Potential and risks of hydrogen-based e-fuels in climate change mitigation, Nat. Clim. Change 11 (2021) 384-393.

[43] O. Fandino, J.M. Trusler, D. Vega-Maza, Phase behavior of (CO2+ H2) and (CO2+ N2) at temperatures between (218.15 and 303.15) K at pressures up to 15 MPa, Int. J. Greenh. Gas. Control 36 (2015) 78-92.

[44] M. Kanaani, B. Sedaee, M. Asadian-Pakfar, Role of cushion gas on underground hydrogen storage in depleted oil reservoirs, J. Energy Storage 45 (2022) 103783.

[45] C. Tsang, W. Street, Phase equilibria in the H2/CO2 system at temperatures from 220 to 290 K and pressures to 172 MPa, Chem. Eng. Sci. 36 (1981) 993-1000.

[46] M. Nazeri, A. Chapoy, R. Burgass, B. Tohidi, Measured densities and derived thermodynamic properties of CO2-rich mixtures in gas, liquid and supercritical phases from 273 K to 423 K and pressures up to 126 MPa, J. Chem. Thermodyn. 111 (2017) 157-172.

[47] P. Theveneau, X. Xu, O. Baudouin, J.-N. Jaubert, P. Ceragioli, C. Coquelet, Vapor-liquid equilibria of the CH4+ CO2+ H2S ternary system with two different global compositions: experiments and modeling, J. Chem. Eng. Data 65 (2020) 1802-1813.

[48] T. Al-Sahhaf, A. Kidnay, E.D. Sloan, Liquid+ vapor equilibria in the N2+ CO2+ CH4 system, Ind. Eng. Chem. Fundam. 22 (1983) 372-380.

[49] N. Xu, J. Dong, Y. Wang, J. Shi, High pressure vapor liquid equilibria at 293 K for systems containing nitrogen, methane and carbon dioxide, Fluid Phase Equilib. 81 (1992) 175-186.

[50] M.S.-W. Wei, T.S. Brown, A.J. Kidnay, E.D. Sloan, Vapor+ liquid equilibria for the ternary system methane+ ethane+ carbon dioxide at 230 K and its constituent binaries at temperatures from 207 to 270 K, J. Chem. Eng. Data 40 (1995) 726-731.

[51] H.-Q. Pan, T. Yang, T.-M. Guo, Measurement and prediction of the bubble/dew point locus in the near-critical region and of the compressed fluid densities of the methane[ ] carbon dioxide-n-butane ternary system, Fluid Phase Equilib. 105 (1995) 259-271.

[52] R.A. Heidemann, A.M. Khalil, The calculation of critical points, AIChE J. 26 (1980) 769-779.
[53] M.L. Michelsen, R.A. Heidemann, Calculation of critical points from cubic two-constant equations of state, AI Ch. EJ;(U. S. 27 (1981).

[54] R. Stockfleth, R. Dohrn, An algorithm for calculating critical points in multicomponent mixtures which can easily be implemented in existing programs to calculate phase equilibria, Fluid Phase Equilib. 145 (1998) 43–52.

[55] U.K. Deiters, I.H. Bell, Unphysical critical curves of binary mixtures predicted with GERG models, Int. J. Thermophys. 41 (2020) 169.

[56] U.K. Deiters, I.H. Bell, Calculation of critical curves of fluid mixtures through solution of differential equations, Ind. Eng. Chem. Res. 59 (2020) 19062–19076.

[57] A. Pina-Martinez, R. Privat, S. Lasala, G. Soave, J.-N. Jaubert, Search for the optimal expression of the volumetric dependence of the attractive contribution in cubic equations of state, Fluid Phase Equilib. 522 (2020) 112750.

[58] A. Harmens, H. Knapp, Three-parameter cubic equation of state for normal substances, Ind. Eng. Chem. Fundam. 19 (1980) 291–294.

[59] G. Heyen, A cubic equation of state with extended range of application, Chemical Engineering Thermodynamics, SA Newman, 175 (1980).

[60] J. Van der Waals, Over de Continuïteit van den Gas-en Vloestoftoestand.(Über die Kontinuittät des Gas-und Flüssigkeitszustands), 1873, in, Dissertation, Universität Leiden, Niederlande, deutsche Übersetzung, Leipzig, 1899.

[61] J.-M. Yu, B.C.-Y. Lu, A three-parameter cubic equation of state for asymmetric mixture density calculations, Fluid Phase Equilib. 34 (1987) 1–19.

[62] C.H. Twu, J.E. Coon, J.R. Cunningham, A new cubic equation of state, Fluid Phase Equilib. 75 (1992) 65–79.