# Let us rethink advanced mixing rules for cubic equations of state

Romain Privat, Jean-Noël Jaubert, Georgios M Kontogeorgis

► To cite this version:

Romain Privat, Jean-Noël Jaubert, Georgios M Kontogeorgis. Let us rethink advanced mixing rules for cubic equations of state. Fluid Phase Equilibria, 2025, 596, pp.114455. 10.1016/j.fluid.2025.114455. hal-05097975

HAL Id: hal-05097975

https://hal.univ-lorraine.fr/hal-05097975v1

Submitted on 4 Jun 2025

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>

## Perspective Article

# Let us rethink advanced mixing rules for cubic equations of state

Romain Privat<sup>a, *</sup>, Jean-Noël Jaubert<sup>a, *</sup> ID, Georgios M. Kontogeorgis<sup>b, *</sup>

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

## ARTICLE INFO

### Keywords:

- Cubic equations of state
- Advanced mixing rules
- Excess Gibbs energy
- Huron-Vidal
- MHV1

## ABSTRACT

One of the most challenging aspects of using a Van der Waals type equation of state for mixtures is determining the appropriate expressions for the coefficients *a* (attractive parameter) and *b* (covolume) involved in this equation. It has been 45 years since Huron and Vidal first proposed the “EoS/gE” advanced mixing rules. By equating under infinite reference pressure, the mathematical expression of the excess Gibbs energy gE derived from an equation of state, to the same quantity issued from an explicit activity coefficient model, they deduced an expression for the ratio *a*/*b*. A decade after Huron and Vidal’s initial proposal, building upon the original proposal, Michelsen subsequently derived the “zero reference pressure” (ZRP) approach and proposed the approximate ZRP mixing rules MHV1 and MHV2. Throughout the 1990’s and 2000’s, the Huron-Vidal and ZRP approaches were subject, often empirically, to multiple revisions in order to remedy some of their well identified shortcomings. It would appear that the debates surrounding advanced mixing rules are now over, with the latest conclusions proposed in the 2000s enjoying a degree of consensus.

The objective of this article is to reopen the debate in light of the scientific insights gained from our recent research on advanced mixing rules for cubic equations of state. The concept of deriving mixing rules by equating the excess Gibbs energy expressed from an equation of state to the same quantity expressed from an activity coefficient model (this equality is called a “matching equation”) was undoubtedly an appealing one. However, experience has shown that such a matching equation is not without its limitations, particularly due to the lack of sufficient constraints. We have reached the conclusion that the only way to derive advanced mixing rules that are free from shortcomings is to ensure that not only are the complete expressions of gE from an equation of state and from an activity coefficient model equal but also that the three separate contributions that make it up (i.e., combinatorial, residual and the product of the pressure by the excess volume P·vE) are equal. As discussed in this paper, achieving this objective is challenging and we conclude that the best and unique solution for developing safe mixing rules is to modify the matching equation proposed by Huron-Vidal and Michelsen and to only equate the residual contributions.

Based on this observation, we demonstrate how the demonstrations of the ZRP and HV mixing rules can be reworked to arrive at a unique and universal (independent of the reference pressure) mixing rule, called UHVMI (Unified Huron Vidal Michelsen) mixing rule. We hope that this result will encourage new thinking on mixing rules for cubic equations.

## 1. Introduction

In 1873, Van der Waals introduced his famous equation and tested its capacity to model the behaviour of pure substances. The question of its applicability to mixtures was subsequently addressed. In the speech that he gave on the occasion of being awarded the Nobel Prize in 1910, he said: "When I first conceived the idea of utilizing my equation to study the properties of binary mixtures I can no longer say. But even 20 years ago, at

the insistence of my friend Kamerlingh Onnes, I was able to publish a complete theory for binary mixtures" [1]. As reported by Rowlinson [2], Lorentz demonstrated as early as 1881 that the covolume *b* of a binary mixture can be expressed as a quadratic function of the amounts of each molecule. Inspired by this result, Van der Waals opened his "Molecular theory of a substance made up of two different constituents" with *a* and *b*, the two characteristic parameters of the Van der Waals equation, expressed as quadratic functions of *x*<sub>i</sub>, the molar fractions of components

<footnote><sup>a</sup> Université de Lorraine, CNRS, LRGP, F-54000 Nancy, France</footnote>

<footnote><sup>b</sup> Center for Energy Resources Engineering (CERE), Department of Chemical and Biochemical Engineering, Building 229, Technical University of Denmark, Denmark</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: romain.privat@univ-lorraine.fr (R. Privat), jean-noel.jaubert@univ-lorraine.fr (J.-N. Jaubert), gk@kt.dtu.dk (G.M. Kontogeorgis).</footnote>
i present in the mixture [2]. Still in his Nobel Prize acceptance speech, Van der Waals commented: "Nevertheless, the consideration that even with b constant my theory had not been unimportant in the case of a single substance gave me the hope that many phenomena would be explained qualitatively provided that a suitable value for a and b were introduced for mixtures". These proposals were, in fact, the first attempts to answer the big question facing equation-of-state developers in the 20th century: what are the optimal mixing rules for the parameters a and b of the cubic equations of state? Further details regarding the intriguing history of Van der Waals and his followers, particularly with regard to the utilisation of the Van der Waals model to elucidate and comprehend the phase behaviour of mixtures, can be found in the excellent book How Fluids Unmix [3].

Until the late 1970's, model developers considered exclusively the original "Van der Waals mixing rules" which express parameters a and b either as quadratic functions of the composition or a quadratic function for a and a linear function for b, in accordance with Van der Waals' practice [2]. To bring more flexibility to the model, a binary interaction parameter, namely $k_{ij}$, was added in the expression of a. The numerous methodologies for expressing and predicting $k_{ij}$ are extensively discussed in reference [4] and chapter 3 of [5].

In 1979, Huron and Vidal [6] observed that the traditional mixing rules were unable to accurately represent highly polar systems. To address this deficiency, they proposed advanced mixing rules aimed at rectifying this issue. Their seminal idea was to deduce the mixing rules from activity-coefficient (AC) models, a category of well-established thermodynamic models that are recognized for their ability to represent liquid mixtures of polar and associating molecules, albeit with a limitation in their range of applicability to low-pressure conditions. In other words, the objective was to develop a novel model that integrates the strengths of both cubic equations of state (CEoS) and AC models. In order to derive the new mixing rule (MR), Huron and Vidal equated the expressions of the excess Gibbs energies ($g^E$) from the CEoS (designated as $g^{E,EoS}$) and from the AC model (designated as $g^{E,γ}$). They then solved this equation for the ratio $a/b$.

Given that CEoS do not assume pressure independence of liquid properties, while AC models do, the $g^{E,EoS}$ expression depends on three variables (the temperature $T$, pressure $P$ and composition characterized by the vector of molar fractions $\mathbf{x}$) and one unknown parameter ($a/b$), whereas the $g^{E,γ}$ expression depends solely on temperature and composition. Consequently, the property $g^{E,EoS}$ is estimated at a specific reference pressure $P_{ref}$ that has to be specified.

$$g^{E,EoS}(T, P_{ref}, \mathbf{x}, \frac{a}{b}) = g^{E,γ}(T, \mathbf{x}) \quad (1)$$

This equation will henceforth be referred to as the "matching equation".

Huron and Vidal proposed the use of $P_{ref} = +\infty$ [6] which requires the selection of a linear mixing rule for $b$ (see below for further details). The mixing rule for parameter $a$ can therefore be readily derived from the expression of $a/b$ and the mixing rule for $b$. Subsequently, Michelsen [7] developed a parallel theory based on the selection of $P_{ref} = 0$ which led to the exact "Zero Reference Pressure" (ZRP) mixing rule for the ratio $a/b$. Such a mixing rule was however implicit which means that an iterative procedure is needed for calculating $a/b$. In order to obtain an explicit mixing rule and to address some intrinsic limitations of the exact ZRP mixing rule at high temperatures, Michelsen developed a number of approximate ZRP mixing rules known as MHV1 [8] and MHV2 [9,10] (where MHV is the abbreviation for "Modified Huron-Vidal"). These fascinating research works have resulted in the formation of two principal schools of thought regarding mixing rules for CEoS, namely, Huron-Vidal (HV) and ZRP, respectively. The detailed expressions of these mixing rules are introduced in the next section.

A significant number of mixing rules for CEoS were proposed throughout the 1990s and 2000s. Some of these rely on the infinite reference pressure, while the majority of them refers to the ZRP. The

most notable examples are the Wong-Sandler approach [11], the combinations of Soave-Redlich-Kwong (SRK) or Peng-Robinson (PR) EoS with the MHV1 or MHV2 versions of the ZRP mixing rules [9,10], the predictive SRK model [12–16], as well as the LCVM EoS which employs a mixing rule for $a/b$ that is a simple linear combination of the HV and MHV1 mixing rules [17]. In the early 2000's, advanced mixing rules were reworked – often empirically – in the light of some shortcomings identified along the 1990's, giving birth to new models, namely, the VTPR EoS (proposed by the developers of the PSRK EoS) [18,19], the UMR-PRU EoS (proposed by the developers of the LCVM EoS) [20,21], the PPR78 model [22], the NRTL-PR EoS [23] etc.

Despite the substantial progress that has been made in the field over the past 45 years, a number of questions remain:

- What is the optimal matching equation? In other words, which property estimated from the EoS should be equated to the same property expressed from an activity coefficient (AC) model: the complete excess Gibbs energy $g^E$, or only some of its three separate contributions (combinatorial part, residual part or $P \cdot v^E$, i.e., the product of the pressure by the excess molar volume)?

- What is the most appropriate reference pressure ($P_{ref} = 0$ or $P_{ref} = +\infty$) for the development of advanced mixing rules?

- Can existing model parameters for $g^{E,γ}$ (determined at low pressure) be used in EoS whose advanced mixing rules incorporate a $g^{E,γ}$?

- What insights can be gained from the multitude of studies that have been conducted on the relative merits and drawbacks of HV and ZRP mixing rules for the development of advanced mixing rules that are devoid of shortcomings?

The aforementioned questions can also be considered in light of the “double objective” sometimes associated with the $EoS/g^E$ mixing rules: i) to extend cubic equations of state to polar systems by replacing the classical Van der Waals one fluid mixing rules with the local composition concept (implied in many explicit $g^{E,γ}$ models like NRTL and UNIQUAC) and ii) to develop (in addition to (i)) predictive successful cubic EoS for polar systems. The latter sometimes implies that existing parameters of explicit $g^{E,γ}$ models can be incorporated into the final model and no new parameter estimation should be carried out. As a “conditio sine qua non” for the proper modelling of polar systems, target (i) is considered more important, although in the early days of the development of $EoS/g^E$ models, target (ii) has been considered very important as well.

The objective of this paper is twofold. Firstly, it addresses the aforementioned questions that have preoccupied model developers from the 1990s to the present day. Secondly, it proposes a novel approach for advanced mixing rules that reconciles the HV and ZRP approaches.

By reworking the original demonstrations proposed by Huron & Vidal (for the HV approach), and Michelsen (for the ZRP approach) and, in particular, by rethinking the terms to be incorporated in the matching equation (with special focus on the combinatorial and residual terms of $g^E$), it will be demonstrated that the HV and MHV mixing rules can be unified through a unique mixing rule for $a/b$, called UHVM (Unified Huron Vidal Michelsen) mixing rule.

In the light of this result, it will be possible to propose an a posteriori analysis of some of the most commonly used $EoS/g^E$ mixing rules, with a view to identifying their respective strengths and weaknesses.

## 2. General presentation of the $EoS/g^E$ mixing rules: definition and key lessons learnt over the past 30 years

### 2.1. Prologue

As stated in the introduction, the approach called “$EoS/g^E$ mixing rules” is a way to express the mixture property $a/b$ by equating the excess Gibbs energy $g^{E,γ}$ coming from an AC model with the same quantity $g^{E,EoS}$ coming from an EoS at a given reference pressure (see the
<table><caption>**Table 1**
Universal constants $r_i$ of the PR and SRK CEoS.</caption><thead><tr><th></th><th>SRK</th><th>PR</th></tr></thead><tbody><tr><td>r<sub>1</sub></td><td>-1</td><td>-1 - &radic;2</td></tr><tr><td>r<sub>2</sub></td><td>0</td><td>-1 + &radic;2</td></tr></tbody></table>

matching Eq. (1)).

The two main CEoS used in practice are the Soave-Redlich-Kwong (SRK) [24] and Peng-Robinson (PR) [25,26] EoS that can be written in the following general form [27]:

$$P = \frac{RT}{v-b} - \frac{a}{(v-r_1 b)(v-r_2 b)} \quad (2)$$

where the $r_i$ are two universal constants characterizing the SRK and PR CEoS (see Table 1).

In order to gain a comprehensive understanding of the fundamental principles underlying the derivation of EoS/g$^E$ mixing rules, the excess Gibbs energy is divided into three contributions: one reflecting combinatorial, or entropic, effects (the molar combinatorial excess Helmholtz energy denoted $a_{comb}^E$), another reflecting residual, or enthalpic, effects (the molar residual excess Helmholtz energy denoted $a_{res}^E$) and the last being a $Pv^E$ term resulting from the relationship between the Gibbs energy and the Helmholtz energy, i.e., $g = a + Pv$. Eq. (3) is then obtained.

$$g^E = a_{comb}^E + a_{res}^E + Pv^E \quad (3)$$

At this stage we should mention that “entropic” and “enthalpic” effects, as mentioned above, follow the terminology often associated with activity coefficient models. These effects and corresponding terms of the models (activity coefficients and EoS) are equivalent but not identical to the excess entropy and excess enthalpy terms of these models. More details on these differences are presented in [28] and this is not discussed further in this manuscript.

The molar excess Gibbs energy from a CEoS for a *n*-component mixture is given in Eq. (4) [28].

$$ \left\{ \begin{aligned} \frac{g^{E,EoS}}{RT} &= \frac{a_{comb}^{E,EoS}}{RT} + \frac{a_{res}^{E,EoS}}{RT} + \frac{Pv^{E,EoS}}{RT} \\ \frac{a_{comb}^{E,EoS}}{RT} &= \underbrace{\sum_{i=1}^{n} x_i \ln \left( \frac{v_{pure \ i} - b_i}{v - b} \right)}_{a_{Flory free volume}^{E,EoS} / RT} = \underbrace{\sum_{i=1}^{n} x_i \ln \left( \frac{b_i}{b} \right)}_{a_{Flory}^{E,EoS} / RT} + \underbrace{\sum_{i=1}^{n} x_i \ln \left( \frac{u_{pure \ i} - 1}{u - 1} \right)}_{a_{NonFlory}^{E,EoS} / RT} \\ \frac{a_{res}^{E,EoS}}{RT} &= \frac{1}{r_1 - r_2} \sum_{i=1}^{n} x_i \alpha_{pure \ i} \ln \left( \frac{u_{pure \ i} - r_2}{u_{pure \ i} - r_1} \right) - \frac{\alpha}{r_1 - r_2} \ln \left( \frac{u - r_2}{u - r_1} \right) \\ Pv^{E,EoS} &= v - \sum_{i=1}^{n} x_i v_{pure \ i} \\ \text{with: } \alpha &= \frac{a}{bRT}, \alpha_{pure \ i} = \frac{a_{pure \ i}}{b_i RT}, u = \frac{v}{b} \text{ and } u_{pure \ i} = \frac{v_{pure \ i}}{b_i} \end{aligned} \right. \quad (4) $$

The reader is invited to take note of the form of the combinatorial term from the EoS, which is designated “Flory free volume” and which can be expressed as the sum of a classical “Flory” term and a complement named “Non Flory”. By Flory term, it is referred to the excess Helmholtz energy ($a_{Flory}^E$) of an athermal solution composed of molecules of different sizes and shapes but showing the same molecular interactions. According to Flory’s theory, $a_{Flory}^E/(RT)$ is simply the sum over all the components of terms $x_i \ln(\phi_i/x_i)$ where $\phi_i = \frac{x_i v_{pure \ i}}{v}$ denotes the volume fraction of component *i* and $x_i$ is its molar fraction. In the context of CEoS, the ratio $\phi_i/x_i$ can be approximated by the ratio $b_i/b$ of the covolume $b_i$ of pure *i* and the covolume *b* of the mixture (rigorously speaking, this equality is satisfied at infinite pressure where the volume is equal to the co-volume).

At this stage, if the two terms $a_{Non Flory}^{E,EoS}$ and $a_{res}^{E,EoS}$ are grouped together, the resulting expression can be written in terms of a Q function, as highlighted in Eq. (5).

$$ \left\{ \begin{aligned} & \frac{a_{Non Flory}^{E,EoS}}{RT} + \frac{a_{res}^{E,EoS}}{RT} = Q(\alpha, \beta) - \sum_{i=1}^{n} x_i Q(\alpha_i, \beta_i) \\ & \text{with: } \beta = \frac{P_{ref}b}{RT} \text{ and } \beta_i = \frac{P_{ref}b_i}{RT} \end{aligned} \right. \quad (5) $$

The formal definition of the Q function is given in Eq. (6).

$$Q(z_1, z_2) = -1 - \ln[u(z_1, z_2) - 1] - \frac{z_1}{r_1 - r_2} \ln \left[ \frac{u(z_1, z_2) - r_2}{u(z_1, z_2) - r_1} \right] \quad (6)$$

Where *u* represents the ratio of the molar volume over the covolume, which is the solution to the third-degree Eq. (7), derived from Eq. (2):

$$ \begin{aligned} z_2 u^3 &+ [z_2(r_1 + r_2 + 1) + 1] u^2 + [z_2(r_1 r_2 + r_1 + r_2) + r_1 + r_2 + z_1] u \\ &- r_1 r_2 (z_2 + 1) - z_1 \\ &= 0 \end{aligned} \quad (7) $$

Eventually, the general expression of $g^{E,EoS}$ can be written in the form of Eq. (8) which is derived from the combination of Eqs. (4) and (5):

$$\frac{g^{E,EoS}}{RT} = Q(\alpha, \beta) - \sum_{i=1}^{n} x_i Q(\alpha_i, \beta_i) + \sum_{i=1}^{n} x_i \ln \left( \frac{b_i}{b} \right) + \frac{Pv^{E,EoS}}{RT} \quad (8)$$

## 2.2. Introduction of the Huron-Vidal (HV) approach

As mentioned above, the HV approach is derived by setting the reference pressure $P_{ref}$ at $+\infty$ in the matching Eq. (1) and by expressing the mixture covolume *b* as a linear function of the molar fractions (this requirement for the covolume expression is justified in Section 2.4).

Without going into all the details, let us analyse how the expression
of the quantity $g^{E,EoS}$ (see Eq. (4)) is modified when the pressure becomes infinite.

$$ \left\{ \begin{array}{l} \lim_{P \to +\infty} \frac{a^{E,EoS}_{\text{comb}}}{RT} = 0 \\[1.5em] \lim_{P \to +\infty} \left( \frac{P \cdot v^{E,EoS}}{RT} \right) = \lim_{P \to +\infty} \left[ \frac{P}{RT} \left( v - \sum_{i=1}^{n} x_i v_{\text{pure } i} \right) \right] = \lim_{P \to +\infty} \left[ \frac{P}{RT} \left( b - \sum_{i=1}^{n} x_i b_i \right) \right] = 0 \quad \text{if } b = \sum_{i=1}^{n} x_i b_i \\[1.5em] \lim_{P \to +\infty} \frac{a^{E,EoS}_{\text{res}}}{RT} = C_{\text{EoS}}^{\text{HV}} \left( \underbrace{\frac{a}{bRT}}_{\alpha} - \sum_{i=1}^{n} \underbrace{x_i a_{\text{pure } i}}_{a_{\text{pure } i}} \right) \end{array} \right. \qquad (9) $$

The definition of the universal constant $C_{EoS}^{HV}$ is:

$$ C_{EoS}^{HV} = \frac{-1}{r_1 - r_2} \ln \left( \frac{1 - r_2}{1 - r_1} \right) = \begin{cases} -1 & \text{for the Van der Waals CEoS} \\ -\frac{\sqrt{2}}{2} \ln(1 + \sqrt{2}) & \approx -0.623 \text{ for the PR EoS} \\ -\ln(2) & \approx -0.693 \text{ for the SRK CEoS} \end{cases} \qquad (10) $$

By combining Eqn 1, 3 and 9, the general expression of the HV mixing rule is obtained:

$$ \begin{aligned} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} &= \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{g^{E,\gamma}(T, \mathbf{x})}{C_{EoS}^{HV}} \\ b(\mathbf{x}) &= \sum_{i=1}^{n} x_i b_i \end{aligned} \qquad (11) $$

Despite many good results, the HV mixing rule was criticized for various reasons:

* As it is believed that the excess Gibbs energies at infinite and low pressures are not the same, existing model parameters for $g^{E,γ}$ may not be suitable for the HV mixing rule using $g^{E,γ}$.

* The second virial coefficient is not a quadratic function of the molar fractions, whereas the quadratic form is imposed by statistical thermodynamics.

* It is very restrictive to always use the linear mixing rule for the covolume.

* Confusion arose as to which $g^{E,γ}$ should be used in Eq. (11). Could we use any external activity coefficient model? To answer this question, we refer to Jean Vidal himself, who wrote in the book he published in 2003 [29]: "the $g^{E,γ}$ excess Gibbs energy expression can be chosen as a function of the problem, and the nature of the mixture, as we do when we apply a heterogeneous ($\gamma - \varphi$) method such as the Wilson equation, the NRTL models, UNIQUAC or UNIFAC." It therefore seems clear that, according to Jean Vidal, any expression for $g^{E,γ}$ could be combined with the HV mixing rules. The same has been stated by many others, e.g., Michelsen and Mollerup [30] in their book used the HV mixing rule with the Wilson equation. And they attributed the different results obtained by HV/Wilson and Wilson to the inability of the HV model in being used with existing parameters of activity coefficient models. However, proof will be given in the next section that claiming that any $g^{E,γ}$ model can be used in the context of HV mixing rule is controversial.

## 2.3. Introduction of the zero reference pressure (ZRP) approach

In the ZRP approach, the reference pressure $P_{ref}$ is set to zero in the matching Eq. (1). Unlike the HV approach, this approach does not constrain the expression of the covolume of the mixture. There are three different versions of the ZRP approach: (i) the exact mixing rule, applicable at low temperature only, (ii) the MHV1 version and (iii) the MHV2 version.

### (i) Exact ZRP mixing rule

As previously stated, the quantity $g^{E,EoS}$ is given by Eq. (8). By combining Eqs. (1), (8) and setting $P_{ref} = 0$, the nonlinear Eq. (12) is obtained [30]. The exact ZRP mixing rule can be extracted by solving Eq. (12) whose unknown is $\alpha = a/(bRT)$, at specified values of T and x.

$$ \underbrace{Q(\alpha, 0)}_{Q_0(\alpha)} - \sum_i x_i \underbrace{Q(\alpha_i, 0)}_{Q_0(\alpha_i)} = \frac{g^{E,\gamma}(T, \mathbf{x})}{RT} + \sum_i x_i \ln \left[ \frac{b(\mathbf{x})}{b_i} \right] \qquad (12) $$

For the sake of brevity, functions $Q(\alpha, 0)$ and $Q(\alpha_i, 0)$ will be simply denoted $Q_0(\alpha)$ and $Q_0(\alpha_i)$ in what follows. Note that the expression of the $Q_0$ function was first introduced by Michelsen [8].

In the expression of $Q_0$, the reduced density $u_0$ at $P_{ref} = 0$ is calculated by solving Eq. (7) and by setting $z_2 = 0$. In this case, Eq. (7) is no longer a cubic equation but a quadratic equation. When a real liquid root exists, its expression is given by Eq. (13).

$$ \underbrace{u(z_1, 0)}_{u_0(z_1)} = \frac{2(r_1 r_2 + z_1)}{r_1 + r_2 + z_1 + \sqrt{(r_1 + r_2 + z_1)^2 - 4(r_1 r_2 + z_1)}} \qquad (13) $$

The condition of existence of the real liquid root is:

$$ z_1 \ge -r_1 - r_2 + 2 + 2\sqrt{1 - r_1 - r_2 + r_1 r_2} \text{ with } z_1 \in \{\alpha, \alpha_i\} \qquad (14) $$

In practice, the consequence of this limitation is the impossibility to apply the exact ZRP mixing rule at high temperature.

### (ii) MHV1 mixing rule

To circumvent this issue, the MHV1 approach uses a linear approximation of the $Q_0$ vs. $\alpha = \frac{a}{bRT}$ function (see Eq. (15)) thereby eliminating the temperature limitation from the application range of the exact ZRP approach.

$$ Q_0(\alpha) \approx k + C_{EoS}^{\text{MHV1}} \cdot \alpha \qquad (15) $$

By combining Eqs. (12) and (15), the MHV1 mixing rule (16) is obtained.
$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{g^{E,\gamma}(T, \mathbf{x}) - RT \sum_{i} x_i \ln \left[ \frac{b_i}{b(\mathbf{x})} \right]}{C_{EoS}^{MHV1}} \quad (16) $$
$b(\mathbf{x}): \text{free choice}$

The universal constant $C_{EoS}^{MHV1}$ involved in Eq. (16) was estimated by Michelsen for the SRK and PR EoS [8].

$$ C_{EoS}^{MHV1} = \begin{cases} -0.530 & \text{for the PR CEoS} \\ -0.593 & \text{for the SRK CEoS} \end{cases} \quad (17) $$

The values assigned to $C_{EoS}^{MHV1}$, along with the methodology employed to obtain them, are discussed in the following paragraphs.

It is noticeable that Eq. (11) (the HV mixing rule) and Eq. (16) (the MHV1 mixing rule) fit into the common formalism of Eq. (18) for the expression of $a/b$:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{E(T, \mathbf{x})}{C_{EoS}} \quad (18) $$

$E(T, \mathbf{x})$ is an excess function involving an AC model $g^{E,\gamma}(T, \mathbf{x})$. The relationship between $E$ and $g^{E,\gamma}$ depends on the used reference pressure (infinite or zero) as discussed above. $C_{EoS}$ is a universal constant that depends only on the nature of the EoS and the type of mixing rule (HV or MHV1).

### (i) MHV2 mixing rule

The MHV2 approach consists in approximating the $Q_0$ function as a quadratic function of $\alpha = \frac{a}{bRT}$. For the present study, using the MHV2 approach would not alter the general conclusions we will draw from the MHV1 approach. For this reason, only the MHV1 version of the ZRP mixing rule will be considered thereafter.

Note that it is beyond the scope of this paper to discuss in details the consequences of using an approximate (MHV1, MHV2) or exact version of the ZRP mixing rule. We can however state that, contrary to the exact ZRP mixing rule, the approximate ZRP rules do not fully match the underlying explicit $g^{E,\gamma}$ model at zero pressure, which was the original target by Michelsen.

## 2.4. On the mixing rule for the covolume in the HV approach

Having considered the mixing rules for the parameters $a/b$, we now turn our attention to the selection of a mixing rule for the covolume $b(\mathbf{x})$. While the ZRP approach affords a free selection of the mixing rule for $b$, the HV approach imposes to use a linear mixing rule for $b$ in order to ensure that $g^{E,EoS}$ remains finite at infinite pressure and can thus be identified with the finite quantity $g^{E,\gamma}$ [6]. Indeed, the $g^{E,EoS}$ function involves a “$P \cdot v^E$” contribution (see Eq. (4)) that diverges at infinite

$$ \begin{aligned} \lim_{P \to +\infty} v^{E,EoS}(T, P, \mathbf{x}) &= \lim_{P \to +\infty} \left( v(T, P, \mathbf{x}) - \sum_{i} x_i v_{\text{pure}}(T, P) \right) \\ &= b(\mathbf{x}) - \sum_{i} x_i b_i = b^E(\mathbf{x}) \end{aligned} \quad (19) $$

Consequently, selecting a linear mixing rule for $b$ ensures that $v^{E,EoS}$ is zero at infinite pressure and thus, that the “$Pv^{E,EoS}$” term vanishes in such conditions (this is what was reported in Eq. (9)).

## 2.5. The Wong-Sandler mixing rule: an alternative to the HV mixing rule

Before closing this section, let us mention also the mixing rule proposed by Wong and Sandler [11], based on the HV approach. Wong and Sandler considered AC models as excess Helmholtz energy ($a^{E,\gamma}$) models instead of excess Gibbs energy models ($g^{E,\gamma}$). Indeed, as pointed out in Eq. (20), classical AC models are all pressure independent and therefore, they predict that excess molar volumes are all zero. Consequently, $a^{E,\gamma}$ and $g^{E,\gamma}$ are strictly equal.

$$ \left\{ \begin{aligned} v^{E,\gamma} &= \left( \frac{\partial g^{E,\gamma}}{\partial P} \right)_{T,x} = 0 \\ a^{E,\gamma} &= g^{E,\gamma} - P \cdot v^{E,\gamma} = g^{E,\gamma} \end{aligned} \right. \quad (20) $$

For this reason, Wong and Sandler based their mixing rule on the matching equation, shown in Eq. (21):

$$ a^{E,\gamma}(T, \mathbf{x}) = a^{E,EoS}(T, P \to +\infty, \mathbf{x}, \alpha) \quad (21) $$

Because $\lim_{P \to +\infty} a^{E,EoS}$ remains finite [contrary to $\lim_{P \to +\infty} g^{E,EoS} = \lim_{P \to +\infty} (a^{E,EoS} + Pv^{E,EoS})$], they did not need to impose the selection of a specific mixing rule for $b$ that satisfies $\lim_{P \to +\infty} v^{E,EoS} = 0$ and makes the diverging term $Pv^{E,EoS}$ vanishing. Therefore, they could select freely the mixing rule for $b$. This new degree of freedom in the model formulation allowed them to choose a specific form of the mixing rule for $b$ such that the second virial coefficient ($B = b - \frac{a}{RT}$) is a quadratic function of the composition, as predicted by statistical mechanics. However, by doing so, the covolume becomes temperature-dependent which is known to lead to unphysical behaviour [31–34].

In short, the mixing rule for $a/b$ is identical to the one proposed by Huron and Vidal. The Wang-Sandler and HV models however differ in the mixing rule for $b$. The equations that describe the Wong-Sandler model are summarized in Eq. (22). Unlike the other EoS/g$^E$ models, it contains one more binary interaction parameter ($l_{ij}$) in the cross second virial coefficient ($B_{ij}$). The common practice [35–38] is to determine the binary parameter $l_{ij}$ in order the equation of state reproduces as closely as possible the experimental $g^{E,\gamma}$ curve obtained from the activity coefficient model.

$$ \frac{a(T, \mathbf{x})}{b(T, \mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{g^{E,\gamma}(T, \mathbf{x})}{C_{EoS}^{HV}} \quad \text{with: } C_{EoS}^{HV} = \begin{cases} -1 & \text{for the VdW EoS} \\ -0.623 & \text{for the PR EoS} \\ -0.693 & \text{for the SRK EoS} \end{cases} \\ b(T, \mathbf{x}) = \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} x_i x_j B_{ij}}{1 - \sum_{i} x_i RT b_i - \frac{g^{E,\gamma}(T, \mathbf{x})}{RTC_{EoS}^{HV}}} \quad \text{with: } B_{ij} = \frac{1}{2} \left[ \left(b_i - \frac{a_i}{RT}\right) + \left(b_j - \frac{a_j}{RT}\right) \right] (1 - l_{ij}) \quad (22) $$

pressure. Since molar volumes tend to covolumes at infinite pressure, the excess volume term at infinite pressure can be expressed as an excess covolume term:

## 2.6. Successes and limitations of the HV and MHV1 mixing rules

We here review the major successes and limitations of HV and MHV1
mixing rules which are well documented in several review articles [5,30, 39–41]. For each of the mixing rules, a theoretical interpretation of the observed weaknesses is given, which allows a new derivation of advanced mixing rules in the next section.

### (i) HV mixing rule

The HV mixing rule is first addressed. In the paper introducing their mixing rule, Huron and Vidal decided to identify the complete excess Gibbs energies $g^{E,EoS}$ and $g^{E,γ}$ at infinite reference pressure (see the matching Eq. (1) with $P_{ref} = +\infty$). In doing so, they obtained the mixing rules reported in Eq. (11) which were employed with a variety of AC models, including UNIQUAC, Van Laar, and NRTL. The primary finding was that accurate correlation of the data was achieved when purely residual AC models, such as NRTL or Van Laar, were utilized [5, 30, 39–41]. In their original paper, Huron and Vidal indeed used NRTL without justifying their choice. To understand the reasons behind this important conclusion, an analysis of the contributions to $g^E$ incorporated in the matching equation, is required.

As mentioned earlier, at infinite pressure, the expression of the excess Helmholtz energy coming from the CEoS contains only a residual part since, as recalled in Eq. (9), the combinatorial part vanishes [5,28, 42]. Therefore, the matching Eq. (1) can be rewritten in the form of Eq. (23) (it is recalled that a linear mixing rule is used for b).

$$ \left\{ \begin{array}{l} \text{Matching equation: } (g^{E,EoS})_{p \to +\infty} = g^{E,\gamma} \\ (g^{E,EoS})_{p \to +\infty} = a_{res}^{E,EoS} + \underbrace{a_{comb}^{E,EoS}}_{=0} + \underbrace{Pv^{E,EoS}}_{=0} \\ g^{E,\gamma} = a_{res}^{E,\gamma} + \underbrace{a_{comb}^{E,\gamma}}_{=0} + \underbrace{Pv^{E,\gamma}}_{=0} \end{array} \right. \Rightarrow \underbrace{a_{res}^{E,EoS} = a_{res}^{E,\gamma} + a_{comb}^{E,\gamma}}_{=0} \quad (23) $$

Eq. (23) demonstrates an unbalance in the matching between the $g^E$ of the AC model and CEoS, as the combinatorial contribution is present solely on the right-hand side of the equality. This unbalanced equation is the underlying cause of the observed reduction in modelling accuracy when a complete AC model is incorporated in Eq. (11). The unique possibility to remove the unbalance is to set $a_{comb}^{E,\gamma} = 0$ in Eq. (23). It is thus immediately possible to conclude, in agreement with Kontogeorgis and Folas [5] that: “Although any activity coefficient model could, in principle, be used, it is clear from the derivation that a ‘solely energetic’ model should be preferred, i.e., a model which does not have an explicit combinatorial term”. In [6], Huron and Vidal did not report the problem pointed out by Eq. (23) and circumvented it by incorporating the purely residual NRTL AC model in their mixing rule. From the above, the HV mixing rules should be written in the revised form proposed in Eq. (24):

$$ \left\{ \begin{array}{l} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{a_{res}^{E,\gamma}(T, \mathbf{x})}{C_{EoS}^{HV}} \\ b(\mathbf{x}) = \sum_{i} x_i b_i \end{array} \right. \quad (24) $$

We can also conclude that, although inconsistent, Eq. (23), should in principle be used to accurately correlate the properties of regular solutions characterized by $s^E = 0$, i.e. by $a_{comb}^{E,\gamma} = 0$. Proof of this important conclusion is given in Fig. 1 and Fig. 2 where we decided to correlate the VLE data of the strongly non-ideal binary system acetone (C$_3$H$_6$O) + chloroform (CHCl$_3$) [43] under different conditions. The two molecules are approximately the same size so that the excess entropy, $s^E$ (or equivalently $a_{comb}^{E,\gamma}$), although not zero, is small. The complete Wilson AC model which contains both a combinatorial and a residual part [44,45] was selected. The numerical values of the two binary interaction parameters ($A_{12} = -190$ K and $A_{21} = -80$ K) were determined in order to minimize the deviations between experimental [46,47] and calculated bubble pressures. Fig. 1 shows that, the model accurately reproduces the

<figure><img src="image_4.png" /><figcaption>Fig. 1. Correlation at P/bar = 0.94 and along the bubble curve of the excess Gibbs energy for the acetone (1) + chloroform (2) system with the Peng–Robinson EoS and the HV MR. The complete Wilson ACM was used in the MR. The symbols are the experimental data points and the full curve shows the result of the correlation.</figcaption></figure>

experimental values of $g^E$ and that, as expected, the combinatorial contribution to $g^E$ is small in comparison to the residual one. In addition and as expected, Fig. 2 clearly shows that the original HV MRs (Eq. (11)), coupled with a complete $g^E$ model (having both an entropic and an enthalpic term) are able to correlate nearly regular solutions.

Symmetrically, the inconsistencies linked to the use of an unbalanced matching equation are further accentuated when examining the modelling of size-asymmetric mixtures (nearly athermal mixtures like, e.g., a mixture of ethane and n-eicosane) for which $a_{comb}^{E,\gamma}$ takes large values. For such systems, we can classically assume that there is no enthalpic effects, so that $a_{res}^{E,\gamma} = 0$. Therefore, the original HV mixing rules (Eq. (11)) combined with a complete $g^{E,\gamma}$ model take the form of Eq. (25):

$$ \left\{ \begin{array}{l} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{a_{comb}^{E,\gamma}(T, \mathbf{x})}{C_{EoS}^{HV}} \\ b(\mathbf{x}) = \sum_{i} x_i b_i \end{array} \right. \quad (25) $$

<figure><img src="image_5.png" /><figcaption>Fig. 2. Correlation of the acetone (1) + chloroform (2) system with the Peng–Robinson EoS and the HV MR. The complete Wilson ACM was used in the MR. The symbols are the experimental data points and the full curve shows the result of the correlation.</figcaption></figure>
Which is known to produce markedly inaccurate predictions of thermodynamic properties in size-asymmetric systems [41,48]. To convince the reader of this weakness, Fig. 3 shows the poor correlation of the vapor-liquid equilibrium data for the binary system ethane (1) + eicosane (2) at three different temperatures [49,50] when the original mixing rule of HV is coupled with the Wilson AC model. For such an athermal system, we have imposed zero values on the binary interaction parameters ($A_{12} = A_{21} = 0$).

The inconsistency of Eq. (25) is even more obvious if one writes the corresponding matching equation (obtained by setting $a_{res}^{E,y} = 0$ in Eq. (23)):

$$a_{res}^{E,EoS} = a_{comb}^{E,y} \qquad (26)$$

Eq. (26) highlights that $a_{res}^{E,EoS}$, which should be equal to zero for this type of system is in fact nonzero in order to compensate for the presence of the combinatorial $a_{comb}^{E,y}$ term from the AC model. This is definitely inconsistent.

If the revised HV mixing rules (Eq. (24)) are applied to size-asymmetric systems, Eq. (27) is obtained.

$$\left\{ \begin{array}{l} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} \\ b(\mathbf{x}) = \sum_{i} x_{i} b_{i} \end{array} \right. \qquad (27)$$

As shown in several papers [51,52], Eq. (27) leads to very satisfactory predictions for size-asymmetric systems, the modelling of which only requires a combinatorial contribution. We can thus conclude that the combinatorial free-volume term of cubic EoS is sound.

**Remark:** under infinite pressure, the mixing rule $a/b = \sum_i x_i a_i / b_i$ can be straightforwardly derived from $(a_{res}^{E,EoS})_{P \to +\infty} = 0$.

As a first conclusion regarding the HV approach, it is imperative to impose that only identical contributions of $g^E$ ($a_{comb}^E$, $a_{res}^E$ and $P \cdot v^E$) from the EoS and the AC model are equated in order to obtain a correct and consistent advanced mixing rule. The use of an unbalanced matching equation should be avoided at all costs.

<figure><img src="image_6.png" /><figcaption>Fig. 3. Correlation of the ethane (1) + eicosane (2) system with the Peng-Robinson EoS and the HV MR. The Wilson ACM with null binary interaction parameters was used in the MR. The symbols are the experimental data points at T/K = 300, 440, 572.85 and the full curve shows the result of the correlation.</figcaption></figure>

### (ii) MHV1 mixing rule

The ZRP mixing rule is now considered and we will use the approximate MHV1 version as an illustration of this family of EoS/g^E models. The MHV1 mixing rule expression given in Eq. (16) can be reformulated by simply expressing the AC model $g^E$ as the sum of a combinatorial and a residual term.

$$\frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \frac{a_{res}^{E,y}(T, \mathbf{x}) + \left[ a_{comb}^{E,y}(T, \mathbf{x}) - RT \sum_{i} x_{i} \ln \left[ \frac{b_{i}}{b(\mathbf{x})} \right] \right]}{C_{EoS}^{\text{MHV1}}} \qquad (28)$$

It is immediately evident that such an equation introduces a so-called “double combinatorial term” (DCT), which reflects the difference between two combinatorial terms: one coming from the AC model and the other one coming from the EoS [the Flory term in Eq. (4)]. Eq. (28) can be rewritten in the concise form of Eq. (29).

$$\frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \frac{a_{res}^{E,y}(T, \mathbf{x}) + \text{DCT}}{C_{EoS}^{\text{MHV1}}} \qquad (29)$$

It is imperative to emphasise that the DCT does not encompass the entire combinatorial contribution from the EoS (which is, as shown in Eq. (4), a Flory free-volume term), but only incorporates the classical EoS Flory term (defined in Eq. (4)).

A number of articles have highlighted that the presence of this DCT is the root cause of the issues observed when utilising MHV1 [5,30,39–41, 53]. It has been shown empirically that removing the DCT from the MHV1 mixing rule enables experimental phase-equilibrium data to be better correlated by the CEoS [5,17,30,39–41,51]. However, we would like to qualify this assertion and provide further insight into the reasons why the modelling of some systems may be inaccurate when the complete MHV1 mixing rule is used.

To begin, it is important to recall that such a DCT arises naturally in the derivation of the MHV1 mixing rule, which is based on the following matching equation:

$$\underbrace{a_{\text{comb}}^{\text{E,EoS}}(T, P=0, \mathbf{x}) + a_{\text{res}}^{\text{E,EoS}}(T, P=0, \mathbf{x}, \alpha)}_{g^{\text{E,EoS}}(T,P=0,\mathbf{x})} = \underbrace{a_{\text{comb}}^{\text{E,y}}(T, \mathbf{x}) + a_{\text{res}}^{\text{E,y}}(T, \mathbf{x})}_{g^{\text{E,y}}(T, \mathbf{x})} \qquad (30)$$

In other words, imposing the equality of the excess Gibbs energies $g^E$ and $g^{E,EoS}$ as a prerequisite for deriving the mixing rule for $a/b$ irretrievably generates the presence of the DCT in this mixing rule.

In the early 2000s, the modelling of size-asymmetric systems by EoS based on the MHV1 mixing rule (and more generally, by EoS relying on the ZRP approach) [53] revealed a major shortcoming in the MHV1 mixing rule. The term DCT was identified as the culprit. To support this statement, let us cite the Section 6.6 in the book by Kontogeorgis and Folas [5]: “an inspection of the MHV1-type mixing rules[...] shows that in reality, when they are combined with the UNIQUAC/UNIFAC models, as typically done, the resulting models contain two different combinatorial terms”. “The two combinatorial terms [...] must be of similar magnitude [...] but this is generally not the case. Moreover, the difference increases with asymmetry [...]. Problems are, as might be anticipated, observed for size-asymmetric systems”.

Although correct, such sentences could give – at first sight – the impression that there are issues with the derivation of the MHV1 mixing rule and that the DCT should be removed from Eq. (29). However, as previously stated this term is necessary as long as it is desired to ensure that the matching equation $g^{E,EoS} = g^{E,y}$ is satisfied. This apparent contradiction invites us to clarify the situation and reopen the debate on the DCT.

Firstly, it is important to recall that incorporating a complete AC model (e.g., UNIQUAC) into the MHV1 mixing rule has the following implications (proof is given below), which are responsible for MHV1’s weaknesses:

- The DCT is never zero at the reference pressure $P_{ref} = 0$ so that the combinatorial contribution of the AC model ($a_{comb}^{E,y}$) is always
different from the combinatorial contribution of the EoS ($a_{comb}^{E,EoS}$). As discussed below, this will cause a problem for the correlation of non-regular systems for which the combinatorial contribution is significant ($s^E \neq 0$, i.e. $a_{comb}^{E,\gamma} \neq 0$) and more particularly for size-asymmetric systems, characterized by large values of $s^E$.

- However, since the matching equation $g^{E,EoS} = g^{E,\gamma}$ is satisfied at $P_{ref} = 0$, the residual contribution of the AC model ($a_{res}^{E,\gamma}$) is also always different from the residual contribution of the EoS ($a_{res}^{E,EoS}$) at $P_{ref} = 0$. Near-regular solutions, should not be strongly affected by this drawback since for such systems, the combinatorial contributions are extremely weak ($a_{comb}^{E,\gamma} \approx a_{comb}^{E,EoS} \approx 0$).

From the above, it is noticeable that the HV and the MHV1 MRs suffer from the same issue. Indeed, for both, it is impossible to satisfy $a_{comb}^{E,\gamma} = a_{comb}^{E,EoS}$. As a direct consequence, size-asymmetric systems are poorly correlated while regular solutions, whose modelling does not require a combinatorial contribution, are not strongly affected by this drawback and should be properly correlated.

To illustrate this last point, Fig. 4 shows that the VLE data of the strongly non-ideal binary system acetone (1) + chloroform (2) are accurately correlated with the MHV1 MR coupled with the complete

Wilson model. These MRs were also used to correlate the experimental $g^E$ values shown in Fig. 1. The correlation results are very similar to those reported in Fig. 1 using HV MRs and are therefore not shown here.

Let us now return to the issue of size-asymmetric (i.e., near athermal) systems in order to gain a deeper understanding of the reasons behind the poor modelling of such systems. As previously highlighted with the HV mixing rule, such systems are characterized by $a_{res}^{E,\gamma} = 0$ so that the MHV1 matching equation, Eq. (30), can be written in this context:

$$a_{comb}^{E,EoS} + a_{res}^{E,EoS} = a_{comb}^{E,\gamma} \qquad (31)$$

with:

$$a_{comb}^{E,EoS} = RT \sum_{i=1}^{n} x_i \ln \left( \frac{v_{pure \ i} - b_i}{v - b} \right) \qquad (32)$$

As previously discussed, the combinatorial part of $g^{E,EoS}$ (see Eq. (32)) is expressed as a free-volume Flory term (i.e., a classical Flory contribution where the volumetric fraction is expressed with free volumes ($v-b$) instead of volumes $v$). The mathematical expressions of the combinatorial terms stemming from well-known AC models (Wilson, UNIQUAC, NRTL or Van Laar) are provided in Eq. (33) for illustration.

$$ \left\{ \begin{array}{l} a_{comb}^{E, Wilson} = RT \sum_{i=1}^{n} x_i \ln \left( \frac{v_{pure \ i}}{v_{liq}} \right) \quad \text{with: } v_{liq} = \sum_{i=1}^{n} x_i v_{pure \ liq \ i} \\ \\ a_{comb}^{E, UNIQUAC} = RT \sum_{i=1}^{n} x_i \ln \left( \frac{r_{pure \ i}}{r_{liq}} \right) + \text{Staverman-Guggenheim term: } r_{liq} = \sum_{i=1}^{n} x_i r_{pure \ liq \ i} \\ \\ r_{pure \ i} \quad \text{is the molecular volume (Van der Waals volume) of component i.} \\ \\ a_{comb}^{E, NRTL} = a_{comb}^{E, Van \ Laar} = 0 \end{array} \right. \qquad (33) $$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Correlation of the acetone (1) + chloroform (2) system with the Peng-Robinson EoS and the MHV1 MR. The complete Wilson ACM was used in the MR. The symbols are the experimental data points and the full curve shows the result of the correlation.</figcaption></figure>

While CEoS produce a combinatorial Helmholtz energy $a_{comb}^{E,EoS}$ that contains a Flory term expressed in terms of free volumes (see Eq. (32)), it is worth noting that the combinatorial term of AC models is either zero (as, for example, in the NRTL and Van Laar models) or a Flory term involving fractions of volumes instead of free volumes (as, for example, in the Wilson, UNIQUAC and UNIFAC models). In the UNIQUAC/UNIFAC cases, the complete combinatorial term consists of the Flory contribution augmented by a Staverman-Guggenheim (SG) term.

It is evident that regardless of the selected AC model incorporated in the matching equation, the combinatorial term derived from the EoS (Eq. (32)) will never align with the combinatorial term derived from an AC model (Eq. (33)). In order to satisfy the matching equation (Eq. (31)), despite the fact that the combinatorial contributions are different, it is imperative to have a non-zero $a_{res}^{E,EoS}$. In other words, a non-zero $a_{res}^{E,EoS}$ has to be introduced to satisfy the matching equation: $a_{comb}^{E,\gamma} = a^{E,EoS}$ at $P_{ref} = 0$ meaning that the non-zero $a_{res}^{E,EoS}$ value corrects the unsatisfactory $a_{comb}^{E,EoS}$ value in order to align $a_{comb}^{E,EoS} + a_{res}^{E,EoS}$ with $a_{comb}^{E,\gamma}$. This appears to be unfortunate.

To summarize, when modelling size-asymmetric systems, the only way to satisfy the matching equation (Eq. (31)) is to introduce a non-zero $a_{res}^{E,EoS}$. However, this is not a physically sound solution for athermal systems which are not supposed to have a residual contribution and it leads to a poor correlation of phase behaviour data. This drawback is a simple consequence of the impossible alignment of $a_{comb}^{E,EoS}$ with $a_{comb}^{E,\gamma}$. In other words, rather than attributing the unsatisfactory correlation of athermal systems to the presence of a DCT which is the strict
consequence of the matching equation expressed in terms of $g^E$ (see Eq. (30)), we believe the root cause lies in the need to add of a non-zero EoS residual term to satisfy at all cost the matching equation $a_{comb}^{E,γ} = a^{E,EoS}$. The situation with the MHV1 MR, where an unbalanced matching equation exists between the AC and EoS models ($a_{comb}^{E,γ} = a_{comb}^{E,EoS} + a_{res}^{E,EoS}$), is analogous to what was observed with the HV mixing rules. The very poor correlation of size-asymmetric systems with the MHV1 MR is illustrated in Fig. 5, where the binary system ethane (1) + eicosane (2) has been chosen as a representative example. As was done previously, the MHV1 MR is coupled with the Wilson AC model in which we have set the binary interaction parameters to zero ($A_{12} = A_{21} = 0$) because the solution is athermal.

As a second conclusion regarding the ZRP approach, we would like to draw the reader's attention to the fact that if the DCT is removed from the mixing rule to improve the correlation of size-asymmetric systems, it is necessary to bear in mind that the MHV1 matching equation (Eq. (31)) is no more satisfied so that the initial objective of aligning the excess Gibbs energies from the EoS and the AC model is not reached. And indeed, this is the case with LCVM and similar models. Matching the underlying $g^{E,γ}$ is sacrificed, but improved results are obtained, also for size-asymmetric systems.

It can therefore be stated that the fundamental weakness of the MHV1 approach is, in fact, the absence of constraints on the matching equation on which it is based. While MHV1 requires that the complete $g^E$ from the AC model and the EoS are the same at the reference pressure $P_{ref} = 0$, we are unable to accurately model athermal solutions that only require, strictly speaking, a combinatorial contribution. This is due to the fact that the residual and combinatorial contributions from the AC model and the EoS will never align in pairs. In conclusion, the MHV1 weaknesses come from the impossibility that the distinct contributions of $g^E$ ($a_{comb}^{E}$ and $a_{res}^{E}$) from the EoS and the AC model are equal. If we impose $g^{E,EoS} = g^{E,γ}$, we will always have in return:

$$ \begin{cases} a_{comb}^{E,EoS}(T, P=0, \mathbf{x}) \neq a_{comb}^{E,\gamma}(T, \mathbf{x}) \\ a_{res}^{E,EoS}(T, P=0, \mathbf{x}, \alpha) \neq a_{res}^{E,\gamma}(T, \mathbf{x}) \end{cases} \quad (34) $$

<figure><img src="image_8.png" /><figcaption>Fig. 5. Correlation of the ethane (1) + eicosane (2) system with the Peng- Robinson EoS and the MHV1 MR. The Wilson ACM with null binary interac- tion parameters was used in the MR. The symbols are the experimental data points at T/K = 300, 440, 572.85 and the full curve shows the result of the correlation.</figcaption></figure>

The following section will present the necessary modifications to the ZRP approach (and, in particular, to the MHV1 mixing rule) in order to ensure a proper modelling of athermal and other systems.

### 3. Formulation of a universal mixing rule unifying the HV and ZRP mixing rules

In this section, in order to derive advanced mixing rules without shortcomings, we propose to reformulate the HV and ZRP mixing rules.

To begin, let us recapitulate the key lessons from the previous section:

- If the only constraint imposed to define an EoS/g$^E$ mixing rule is the equality of the complete $g^E$ (see the matching Eq. (1)), modelling issues will arise, especially in systems containing size-asymmetric compounds. This observation reveals the existence of a fundamental problem with this way of deriving mixing rules.

- To avoid the observed shortcomings, we evidenced that the three separate contributions (combinatorial, residual and $Pv^E$) from the EoS should be equal to the corresponding contributions from the AC model. Ideally, at a specified reference pressure, the objectives are:

$$ g^{E,EoS} = g^{E,\gamma} \text{ with } \begin{cases} a_{comb}^{E,EoS} = a_{comb}^{E,\gamma} \\ a_{res}^{E,EoS} = a_{res}^{E,\gamma} \\ P\gamma^{E,EoS} = P\gamma^{E,\gamma} \end{cases} \quad (35) $$

Any matching equation that reflects an impossible objective (i.e. where it is not possible to ensure strict equality between the contribution of the EoS and the corresponding contribution of the AC model) should be considered *unbalanced* and will result in modelling issues, as stated above.

Two examples of unbalanced matching equations were provided in the previous section:

- In Eq. (23), a residual contribution from the EoS is equated to both a combinatorial and a residual contribution from the AC model.

- A second example is Eq. (31) where the opposite occurs, as both a combinatorial and a residual contribution from the EoS are equated with only a combinatorial term from the AC model.

To guarantee safe correlation of regular ($s^E = 0$), athermal ($h^E = 0$) or mixed ($s^E \neq 0$ and $h^E \neq 0$) solutions, it is absolutely essential to respect the division of $g^E$ in three separate contributions (see Eq. (35)). However, as discussed in the previous section, it is never possible to get $a_{comb}^{E,γ} = a_{comb}^{E,EoS}$, neither in the HV approach (since $a_{comb}^{E,EoS} = 0$ at $P_{ref} = +\infty$ while $a_{comb}^{E,γ} \neq 0$ in the general case), nor in the MHV1 approach (a free volume term will never align with the combinatorial term of an AC model).

In addition, except in very special cases (e.g., if the pressure $P$ is zero), it is never possible to satisfy $P\gamma^{E,EoS} = P\gamma^{E,\gamma}$ since $\gamma^{E,\gamma}$ is always zero.

There is thus no solution enabling to satisfy Eq. (35) and we have to conclude that the unique remaining possibility to develop – drawback free – mixing rules is thus to use as matching equation:

$$ a_{res}^{E,EoS}(T, P_{ref}, \mathbf{x}, \alpha) = a_{res}^{E,\gamma}(T, \mathbf{x}) \quad (36) $$

This is the basis of the approach we develop in the next section. In some sense, by accepting this, we abandon one of the original targets of the EoS/g$^E$ models (to use existing parameters of the underlying activity coefficient models) and we focus on developing a sound cubic EoS model with mixing rules permitting extension of the model to polar systems. However, there are cases where, even under the conditions implied in
Eq. (36), existing parameters of external activity coefficient models can be used. We return to this point later.

It is worth noting that although the matching equation contains a residual contribution only, the expression of the excess Gibbs energy from the EoS still contains combinatorial and $Pv^E$ terms, the expression of which are:

This is in agreement with the common modelling practices for such systems [5,41,51,52].

## 3.2. The revisited MHV1 (RMHV1) mixing rule

The matching Eq. (36) is now applied with $P_{ref} = 0$.

To understand how we revisited the original MHV1 mixing rule, let us follow the same steps as in the demonstration originally proposed by Michelsen [7,8], now applied to the new matching equation.

The mathematical expression of $a_{res}^{E,EoS}(T, P = 0, \mathbf{x}, \alpha)$ is first recalled:

$$ \left\{ \begin{array}{l} \frac{a_{res}^{E,EoS}(T, P = 0, \mathbf{x}, \alpha)}{RT} = \sum_{i} x_{i} \frac{\alpha_{pure \ i}}{r_{1} - r_{2}} \ln \left( \frac{u_{0}(\alpha_{pure \ i}) - r_{2}}{u_{0}(\alpha_{pure \ i}) - r_{1}} \right) - \frac{\alpha}{r_{1} - r_{2}} \ln \left( \frac{u_{0}(\alpha) - r_{2}}{u_{0}(\alpha) - r_{1}} \right) \\ \\ \text{With : } u_{0}(z) = \frac{2(r_{1}r_{2} + z)}{r_{1} + r_{2} + z + \sqrt{(r_{1} + r_{2} + z)^{2} - 4(r_{1}r_{2} + z)}} \quad \text{and } \alpha_{pure \ i} = \frac{a_{pure \ i}}{b_{i}RT} \end{array} \right. \qquad (42) $$

$$ a_{comb}^{E,EoS} = RT \sum_{i=1}^{n} x_i \ln \left( \frac{v_{pure \ i} - b_i}{v - b} \right) \qquad (37) $$

$$ P v^{E,EoS} = P \left( \nu - \sum_{i=1}^{n} x_i v_{pure \ i} \right) \qquad (38) $$

As previously stated, this combinatorial free-volume term from the EoS is known to be accurate [42] and sound [5,39–41,51].

## 3.1. The revisited HV (RHV) mixing rule

In this approach, we specify $P_{ref} = +\infty$. The mixing rule for $a/b$ is derived from the new matching Eq. (36) that is combined with the expression of the EoS residual term under infinite pressure (Eq. (9)). One obtains:

$$ \left\{ \begin{array}{l} a_{res}^{E,Y}(T, \mathbf{x}) = C_{EoS}^{RHV} \left[ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} - \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} \right] \\ \\ \text{With : } C_{EoS}^{RHV} = \frac{-1}{r_{1} - r_{2}} \ln \left( \frac{1 - r_{2}}{1 - r_{1}} \right) \end{array} \right. \qquad (39) $$

The expression of the revisited HV (RHV) mixing rule is then:

$$ \left\{ \begin{array}{l} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \frac{a_{res}^{E,Y}(T, \mathbf{x})}{C_{EoS}^{RHV}} \quad \text{with:} \quad \begin{cases} C_{PR}^{RHV} = -0.623 \\ C_{SRK}^{RHV} = -0.693 \end{cases} \\ b(\mathbf{x}): \text{free choice} \end{array} \right. \qquad (40) $$

Compared to the original version of the HV mixing rule (see Eq. (11)), we immediately notice that:

* The complete $g^{E,Y}$ is replaced by its residual contribution,

* The derivation of the mixing rule for $a/b$ does not require $b$ to be a linear function of $\mathbf{x}$. The mixing rule for $b$ can be freely chosen.

* For size-asymmetric systems, characterized by $a_{res}^{E,Y} = 0$, the mixing rule is simply:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} \qquad (41) $$

At this step, in order to have a more compact mathematical expression, it is desirable – following Michelsen – to introduce a function that is similar to the $Q_0$ function that was defined in Eqs. (6) and (12). The new function is denoted $Q_{0,RMHV1}$ and is defined as:

$$ Q_{0,RMHV1}(z) = \frac{-z}{r_1 - r_2} \ln \left[ \frac{u_0(z) - r_2}{u_0(z) - r_1} \right], \quad z \in \{\alpha, \alpha_i\} \qquad (43) $$

By combining Eqs. (42) and (43), one obtains:

$$ \frac{a_{res}^{E,EoS}(T, P = 0, \mathbf{x}, \alpha)}{RT} = Q_{0,RMHV1}(\alpha) - \sum_{i} x_{i}Q_{0,RMHV1}(a_{i}) \qquad (44) $$

The function $Q_{0,RMHV1}$ is plotted in Fig. 6(a) for the PR EoS. Still following the way Michelsen derived the MHV1 mixing rule, it is observed that for $\alpha > 9$, the plot of the function can be approximated by a straight line and consequently, $Q_{0,RMHV1}(\alpha)$ can be approximated by an affine function (see Eq. (45)).

$$ Q_{0,RMHV1}(\alpha) \approx k_{RMHV1} + \alpha \cdot C_{EoS}^{RMHV1} \qquad (45) $$

$C_{EoS}^{RMHV1}$ is thus simply the slope of the blue straight line plotted in Fig. 6 (a).

The combination of Eqs. (44) and (45) provides the final expression of $a_{res}^{E,EoS}$:

$$ \frac{a_{res}^{E,EoS}}{RT} \approx C_{EoS}^{RMHV1} \left( \alpha - \sum_{i} x_{i} \alpha_{i} \right) \qquad (46) $$

By injecting the expression for $a_{res}^{E,EoS}$ in the matching Eq. (36), the final expression of the RMHV1 mixing rule is obtained:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \frac{a_{res}^{E,Y}(T, \mathbf{x})}{C_{EoS}^{RMHV1}} \qquad (47) $$

In order to fully define this mixing rule, it is necessary to define the value of the universal constant $C_{EoS}^{RMHV1}$. When he developed the MHV1 mixing rule, Michelsen set $C_{EoS}^{RMHV1}$ to -0.53 for the PR EoS. This value was obtained as the first derivative of the $Q_0$ function linearized over a range of low $\alpha$ values, specifically between $\alpha = 9$ and $\alpha = 14$ [30]. We can also state that $C_{EoS}^{RMHV1}$ is the slope of the straight line obtained by plotting $Q_0$ as a function of $\alpha$. However, the value of $C_{EoS}^{RMHV1}$ [the slope of the green curve in Fig. 6(a)] is very sensitive to the range of $\alpha$ values considered for the linearization of $Q_0$. To illustrate this statement, $C_{EoS}^{RMHV1}$ was estimated as a function of the upper bound $\alpha_{upper bound}$ of the $\alpha$ range
<figure><img src="image_9.png" /><figcaption>Fig. 6. (a) Plots of the Q₀ and Q₀,rmhv1 functions on the α range [7.5, 30] for the Peng-Robinson EoS. (b) Slope of the linear regression line (i.e., C<sup>MHV1</sup><sub>PR</sub> and C<sup>RMHV1</sup><sub>PR</sub> values) versus the upper bound of the α range used to linearize the Q₀ and Q₀,rmhv1 functions.</figcaption></figure>

[9; α<sub>upper bound</sub>] used to linearize the Q₀ function. This upper bound was varied between α<sub>upper bound</sub> = 10 and α<sub>upper bound</sub> = 10<sup>4</sup>. The plot of C<sup>MHV1</sup><sub>EoS</sub> versus α<sub>upper bound</sub> is reported in Fig. 6(b) for the PR EoS (refer to the thin green line). It can be observed that C<sup>MHV1</sup><sub>PR</sub> exhibits a range between -0.48 (for α<sub>upper bound</sub> = 10) and C<sup>HV</sup><sub>PR</sub> = -0.623 (for α<sub>upper bound</sub> → +∞). Furthermore, it can be noted that the specific value C<sup>MHV1</sup><sub>PR</sub> = -0.53, selected by Michelsen, is associated with α<sub>upper bound</sub> ≈ 14 [see the black circle in Fig. 6(b)], as expected. Eq. (48) below gives the range of variation of C<sup>MHV1</sup><sub>EoS</sub> for both the PR and the SRK CEoS. The values selected, somewhat arbitrarily, by Michelsen are also recalled.

$$
\left\{
\begin{array}{l}
-0.480 \le C_{PR}^{\text{MHV1}} \le (C_{PR}^{\text{HV}} = -0.623) \\
\text{Michelsen selected } C_{PR}^{\text{MHV1}} = -0.530 \text{ without strong justification} \\
\\
-0.562 \le C_{SRK}^{\text{MHV1}} \le (C_{SRK}^{\text{HV}} = -0.693) \\
\text{Michelsen selected } C_{SRK}^{\text{MHV1}} = -0.593 \text{ without strong justification}
\end{array}
\right.
\qquad (48)
$$

It may be surprising to readers that the universal constant of the HV mixing rule is the limit value of C<sup>MHV1</sup><sub>EoS</sub> when the Q<sub>0</sub> function is linearized over the range [9: + ∞]; a mathematical justification of this result is given just below.

The same methodology applied to the constant C<sup>RMHV1</sup><sub>EoS</sub> leads to the
conclusion that this parameter is also a sensitive function of α<sub>upper bound</sub>,
as illustrated in Fig. 6(b). The C<sup>RMHV1</sup><sub>PR</sub> parameter (slope of the blue line in

Fig. 6(a)) ranges from -0.694 to C<sup>HV</sup><sub>PR</sub> = -0.623. In a manner analogous to the behaviour observed for C<sup>MHV1</sup><sub>EoS</sub>, C<sup>RMHV1</sup><sub>EoS</sub> tends to C<sup>HV</sup><sub>EoS</sub> when the function Q<sub>0,RMHV1</sub> is linearized over the α range [9; + ∞].

To justify why C<sup>MHV1</sup><sub>EoS</sub>→C<sup>HV</sup><sub>EoS</sub> and C<sup>RMHV1</sup><sub>EoS</sub>→C<sup>HV</sup><sub>EoS</sub> when the Q<sub>0</sub> and Q<sub>0,RMHV1</sub> functions are linearized over the α range [9; + ∞], it is enough to demonstrate that these functions are asymptotically equivalent to the function (C<sup>HV</sup><sub>EoS</sub>·α) when α tends to + ∞. For example, let us consider the function Q<sub>0,RMHV1</sub> that is defined in Eqs. (43) (and (42) for the u<sub>0</sub> function involved in the expression of Q<sub>0,RMHV1</sub>(α)). By passing these expressions to the limit when α→ + ∞, the result is immediately obtained:

$$
\left\{
\begin{aligned}
\lim_{\alpha \to +\infty} u_0(\alpha) &= \lim_{\alpha \to +\infty} \frac{2(r_1 r_2 + \alpha)}{r_1 + r_2 + \alpha + \sqrt{(r_1 + r_2 + \alpha)^2 - 4(r_1 r_2 + \alpha)}} = 1 \\
Q_{0,\text{RMHV1}}(\alpha) &= \frac{-\alpha}{r_1 - r_2} \ln \left[ \frac{u_0(\alpha) - r_2}{u_0(\alpha) - r_1} \right] \underset{\alpha \to +\infty}{\sim} \frac{-\alpha}{r_1 - r_2} \ln \left( \frac{1-r_2}{1-r_1} \right) = \alpha \cdot C_{\text{EoS}}^{\text{HV}}
\end{aligned}
\right.
\qquad (49)
$$

The same demonstration applies to the function Q₀,rmhv1 which is
related to Q₀,rmhv1 by Eq. (50).

$$
Q_{0,\text{RMHV1}}(\alpha) = -1 - \ln[u_0(\alpha) - 1] + Q_{0,\text{RMHV1}}(\alpha) \quad (50)
$$

When $\alpha \to +\infty$; the quantity $-1 - \ln[u_0(\alpha) - 1]$ is dominated by $Q_{0,\text{RMHV1}}$ so that $Q_{0,\text{RMHV1}}(\alpha) \approx Q_{0,\text{RMHV1}}(\alpha)$, thus proving that both functions exhibit analogous limiting behaviours.

Eventually, the expression of the RMHV1 mixing rule is:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} \frac{a_i(T)}{b_i} + \frac{a_{res}^{E,\gamma}(T, \mathbf{x})}{C_{EoS}^{\text{RMHV1}}} \quad \text{with:} \\
\displaystyle \qquad \left\{
\begin{array}{l}
-0.694 \le C_{PR}^{\text{RMHV1}} \le (C_{PR}^{\text{HV}} = -0.623) \\
-0.734 \le C_{SRK}^{\text{RMHV1}} \le (C_{SRK}^{\text{HV}} = -0.693)
\end{array}
\right.
\end{array}
\right.
\quad \text{b}(\mathbf{x}): \text{free choice}
$$

(51)
Compared to the original version of the MHV1 mixing rule (see Eq. (16)), we immediately notice that:

* The complete $g^{E,Y}$ is replaced by its residual contribution,

* For size-asymmetric systems, characterized by $a_{res}^{E,Y} = 0$, the mixing rule for $a/b$ is simply:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} \qquad (52) $$

### 3.3. Towards a universal and unified mixing rule

To conclude Section 3, for both the zero-reference pressure (ZRP) and the infinite reference pressure (HV), the unique matching equation that we propose to consider is:

$$ a_{res}^{E,EoS}(T, P_{ref}, \mathbf{x}, \alpha) = a_{res}^{E,Y}(T, \mathbf{x}) \qquad (53) $$

This matching equation has many key advantages:

1. It is physically consistent, as the quantity from the EoS and the AC model that are equated have the same physical meaning.

2. It is always possible to find an $\alpha$ value that satisfies this equation.

3. It prevents from poor correlation in size-asymmetric systems.

4. Although Eq. (53) contains a residual contribution only, the expression of the excess Gibbs energy from the EoS still contains accurate combinatorial and $Pv^E$ terms, the expression of which are given in Eqs. (37) and (38).

5. Perhaps most surprisingly of all, the mixing rule can be expressed in a universal form, i.e., a form that is independent of the reference pressure considered thus unifying the HV and MHV mixing rules. It is indeed clear that Eqs. (40) and (51) are identical. They differ only in that the constant $C_{EoS}$ takes on a fixed value at infinite pressure, whereas, at zero pressure, it is allowed to vary within a certain range (encompassing the value at infinite pressure). We can thus state that the RMHV1 mixing rule encompasses the RHV mixing rule. One way of introducing flexibility into the proposed mixing rule could therefore be to consider $C_{EoS}$ as an adjustable parameter.

6. There is no restriction regarding the selection of the mixing rule for the covolume $b$: a linear or a quadratic mixing rule can be selected.

In conclusion, Eq. (54) here below [a compilation of Eqs. (40) and (51)], summarizes what we decided to call the UHVM (Unified Huron Vidal Michelsen) mixing rule. For a given mixture, the constant parameter $C_{EoS}^{UHVM}$ may be arbitrarily chosen in the ranges recalled in Eq. (54) or fitted to experimental data in order to determine the value that best correlate the data.

**UHVM mixing rule:**

$$ \left\{ \begin{array}{l} \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \frac{a_{res}^{E,Y}(T, \mathbf{x})}{C_{EoS}^{UHVM}} \quad \text{with:} \\ b(\mathbf{x}): \text{free choice} \end{array} \right. \qquad \left\{ \begin{array}{l} -0.694 \le C_{PR}^{UHVM} \le \underbrace{-0.623}_{\text{maximum value}} \\ -0.734 \le C_{SRK}^{UHVM} \le \underbrace{-0.693}_{\text{maximum value}} \end{array} \right. \qquad (54) $$

The maximum value of $C_{EoS}^{UHVM}$ checks: $C_{EoS}^{UHVM,MAX} = \frac{-1}{r_1 - r_2} \ln(\frac{1 - r_2}{1 - r_1})$

$C_{EoS}^{UHVM}$ can be arbitrarily fixed or adjusted on experimental data.

The PPR78 model [22,54–59] and its enhanced version E-PPR78 [60, 61] rely on such mixing rules in which a linear mixing rule was used for $b$ and a Van Laar type $a_{res}^{E,Y}(T, \mathbf{x})$ model was selected. More recently [42], the proposed mixing rule [Eq. (40)] was combined with the NRTL, residual-UNIQUAC and residual-Wilson AC models and the following quadratic mixing rule was used for the $b$ parameter:

$$ \left\{ \begin{array}{l} b_m = \sum_{i} \sum_{j} z_i z_j b_{ij} \\ b_{ij} = \left( \frac{b_i^{1/s} + b_j^{1/s}}{2} \right)^s \quad \text{with } s = \frac{3}{2} \end{array} \right. \qquad (55) $$

By employing these mixing rules, the PR EoS was able to obtain an excellent correlation of the phase behaviour and energetic properties of 200 binary systems [62]. To give proof of the accuracy of such advanced MRs, they were used to correlate the VLE data of the binary system ethane (1) + eicosane (2), as was done previously with the original HV (see Fig. 3) and MHV1 (see Fig. 5) MRs. The results are shown in Fig. 7

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Correlation of the ethane (1) + eicosane (2) system with the Peng-Robinson EoS and the UHVM MR. The Wilson ACM with null binary interaction parameters was used in the MR. The symbols are the experimental data points at T/K = 300, 440, 572.85 and the full curve shows the result of the correlation.</figcaption></figure>
and speak for themselves: the proposed MRs are indeed capable of perfectly correlating the experimental data. Compared with what was obtained in Fig. 3 and Fig. 5, the improvement is really significant.

We can now discuss a hot topic related to the use in Eq. (54) of existing g<sup>E,γ</sup> parameters, typically obtained from low pressure experimental data and available in the DECHEMA Chemistry Data Series. In case, the external g<sup>E,γ</sup> model only contains a residual contribution like NRTL or Van Laar, and was used to correlate the experimental data of a nearly regular solution characterized by a<sup>E,γ</sup><sub>comb</sub>(T, x) ≈ 0, it should be possible to use existing parameters of g<sup>E,γ</sup> in the proposed UHVM mixing rule. On the other hand, if the residual term of Wilson or UNIQUAC/ UNIFAC is used in Eq. (54), it is not clear if the existing parameters of g<sup>E,γ</sup> could be used in the UHVM mixing rule because they have been fitted together with a combinatorial term which was not a free-volume term. Once again, in the case of a nearly regular solution for which the combinatorial term is small, the use of existing parameters should be possible.

At this stage, we feel it is important to draw the reader's attention to the fact that the binary interaction parameters available in the DECHEMA tables have, in general, been determined independently for each series of experimental data points. For a binary system, if experimental data are available at several temperatures (isothermal data) and/or pressures (isobaric data), the DECHEMA tables will contain as many binary interaction parameters values as there are different series, i.e. different temperatures and pressures. Therefore, when discussing the accuracy of the new mixing rules, it would be unfair to use the binary interaction parameters determined at one temperature or pressure and to predict experimental data at other conditions.

## 4. Is the constant packing fraction assumption compatible with the MHV1 and RMHV1 mixing rules?

In the 1990's, some authors introduced an alternative way to demonstrate the MHV1 mixing rule based on the assumption that the packing fraction (u = v/b) remains constant during the mixing process at a fixed temperature and pressure [37,53,63–66]. Briefly, their demonstration can be summarized as follows:

- According to Eq. (4), the exact expression of g<sup>E,EoS</sup>/RT under zero pressure is:

$$ \frac{g^{E,EoS}(P=0)}{RT} = \frac{a_{\text{comb}}^{E,EoS}}{RT} + \frac{a_{\text{res}}^{E,EoS}}{RT} \\ = \sum_{i=1}^{n} x_i \ln\left(\frac{b_i}{b}\right) + \sum_{i=1}^{n} x_i \ln\left(\frac{u_{\text{pure } i} - 1}{u - 1}\right) + \frac{1}{r_1 - r_2} \sum_{i=1}^{n} x_i \alpha_{\text{pure } i} \ln\left(\frac{u_{\text{pure } i} - r_2}{u_{\text{pure } i} - r_1}\right) - \frac{\alpha}{r_1 - r_2} \ln\left(\frac{u - r_2}{u - r_1}\right) \quad (56) $$

- Assuming the constant packing fraction approximation which implies $u_{\text{pure } i} = u$, Eq. (56) in which the constant packing fraction under zero pressure is denoted $u_0$ yields:

$$ \frac{g^{E,EoS}(P=0, u_{\text{constant}})}{RT} = \sum_{i=1}^{n} x_i \ln\left(\frac{b_i}{b}\right) \\ - \frac{1}{r_1 - r_2} \ln\left(\frac{u_0 - r_2}{u_0 - r_1}\right) \left(\alpha - \sum_{i=1}^{n} x_i \alpha_{\text{pure } i}\right) \quad (57) $$

By introducing the MHV1 matching equation: g<sup>E,EoS</sup>(P = 0) = g<sup>E,γ</sup>, we immediately get:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_i \frac{a_i(T)}{b_i} + \frac{g^{E,\gamma}(T, \mathbf{x}) - RT \sum_{i} x_i \ln\left[\frac{b_i}{b(\mathbf{x})}\right]}{-\frac{1}{r_1 - r_2} \ln\left(\frac{u_0 - r_2}{u_0 - r_1}\right)} \quad (58) $$

By denoting:

$$ C_{EoS}^{\text{MHV1}} = -\frac{1}{r_1 - r_2} \ln \left( \frac{u_0 - r_2}{u_0 - r_1} \right) \quad (59) $$

The expression of the MHV1 mixing rule appears (compare Eqs. (58) and (59) with Eq. (16)). Note that Eq. (59) connects the C<sup>MHV1</sup><sub>EoS</sub> and the u<sub>0</sub> constant. The expression of the constant packing fraction deduced from Eq. (59) is then:

$$ u_0 = \frac{r_1 - r_2 e^{-C_{\text{EoS}}^{\text{MHV1}}(r_1 - r_2)}}{1 - e^{-C_{\text{EoS}}^{\text{MHV1}}(r_1 - r_2)}} \quad (60) $$

For the PR EoS, C<sup>MHV1</sup><sub>PR</sub> = -0.53 (see Eq. (17)) and the value of the constant packing fraction is ≈ 1.228.

Remarkably, this demonstration makes it possible to derive the MHV1 mixing rule although the assumptions on which it is based are not the same as those used by Michelsen (see Michelsen's demonstration in Section 2.3(ii)).

For clarity, the differences between the two approaches are here summarized:

* Michelsen never assumed a constant packing fraction and the equations he derived led to values of u<sub>0</sub>(α) (refer to Eq. (13)) that actually vary between 1 and 2.4 when the Peng-Robinson EoS is considered and when α varies between α<sub>min</sub> ≈ 6.8 and an infinite value. α<sub>min</sub> = 4 + 2√2 ≈ 6.8 is the lowest value of α at which the PR EoS has a liquid root under zero pressure (value calculated from Eq. (14)). As expected, Fig. 8, in which u<sub>0</sub>, obtained from the Peng-Robinson EoS, is plotted as a function of α, clearly demonstrates that the constant packing fraction assumption is not valid.

* Furthermore, Michelsen showed that the function Q<sub>0</sub>(α), recalled in Eq. (61), can be well approximated by a straight line.

$$ \begin{aligned} Q_0(\alpha) &= -1 - \ln[u_0(\alpha) - 1] - \frac{\alpha}{r_1 - r_2} \ln\left[\frac{u_0(\alpha) - r_2}{u_0(\alpha) - r_1}\right] \\ &\approx k + C_{EoS}^{\text{MHV1}} \cdot \alpha \end{aligned} \quad (61) $$

* On the other hand, when the constant packing fraction assumption is used, u<sub>0</sub> is obviously constant (it becomes α-independent) so that the Q<sub>0</sub>(α) function recalled in Eq. (61) is rigorously linear:

$$ Q_0(\alpha) = \underbrace{-1 - \ln[u_0 - 1]}_{\text{constant}} + \underbrace{-\frac{1}{r_1 - r_2} \ln\left[\frac{u_0 - r_2}{u_0 - r_1}\right]}_{(C_{\text{EoS}}^{\text{MHV1}}) \text{ constant packing fraction assumption}} \cdot \alpha \quad (62) $$
<figure><img src="image_11.png" /><figcaption>Fig. 8. Graphical representation of the function $u_0(\alpha)$, defined in Eq. (13). The Peng-Robinson EoS is here considered.</figcaption></figure>

Let us now investigate further why so large differences between the two approaches lead to the same mixing rule. As indicated above, if the packing fraction is not assumed constant, function $Q_0$ defined by Michelsen can be considered linear to a first approximation. Its graphical representation is a straight line the slope of which can be approximated by the derivative of $Q_0(\alpha)$ with respect to $\alpha$ (denoted $Q'_0(\alpha)$).

This derivative is expressed in Eq. (63):

$$ 
\begin{aligned}
Q'_{0}(\alpha) &= \frac{dQ_{0}}{d\alpha} \\
&= -u'_{0}(\alpha) \left( \frac{1}{u_{0}(\alpha) - 1} - \frac{\alpha}{(u_{0}(\alpha) - r_{1})(u_{0}(\alpha) - r_{2})} \right) \\
&\quad - \frac{1}{r_{1} - r_{2}} \ln \left( \frac{u_{0}(\alpha) - r_{2}}{u_{0}(\alpha) - r_{1}} \right)
\end{aligned}
\qquad (63) $$

By noticing that the quantity $\left(\frac{1}{u_0(\alpha)-1} - \frac{\alpha}{(u_0(\alpha)-r_1)(u_0(\alpha)-r_2)}\right)$ is actually the ratio $\frac{bP}{RT}$ and remembering that the pressure is assumed to be zero, the final expression of $Q'_0(\alpha)$ is obtained:

$$ Q'_{0}(\alpha) = \frac{dQ_{0}}{d\alpha} = -\frac{1}{r_{1} - r_{2}} \ln \left( \frac{u_{0}(\alpha) - r_{2}}{u_{0}(\alpha) - r_{1}} \right) \qquad (64) $$

Therefore, $Q'_0(\alpha)$ has the same expression as the slope of the rigorously linear $Q_0(\alpha)$ function obtained under the constant packing fraction assumption (compare Eqs. (59) and (64)). At this step, we want to insist on the fact that the particular expression of Eq. (64) was never obtained by assuming that the packing fraction $u_0(\alpha)$ under zero pressure is constant, or, in other words, by postulating that $u'_0(\alpha) = 0$; in Eq. (63), the term $u'(\alpha) \cdot \left(\frac{1}{u_0(\alpha)-1} - \frac{\alpha}{(u_0(\alpha)-r_1)(u_0(\alpha)-r_2)}\right)$ is equal to zero, not because the first factor $u'(\alpha)$ is zero but because the second term (proportional to the pressure) is equal to zero. In other words, the reason why the proof of MHV1 based on the idea that the packing fraction is constant and Michelsen's original proof of MHV1 lead to the same expression can be considered as a coincidence and cannot be ascribed to the fact that the packing fraction is constant to a first approximation.

Let us now answer the question: could a demonstration based on the constant packing fraction be used to derive the RMHV1 mixing rule?

As explained below, the RMHV1 mixing rule uses a matching equation based on the equality of $a_{res}^E$ contributions. The expression of this mixing rule is recalled in Eq. (65).

$$ \frac{a_{\text{res}}^{E/\gamma}(T, \mathbf{x})}{RT} = \sum_{i} x_{i} \frac{\alpha_{\text{pure } i}}{r_{1} - r_{2}} \ln \left( \frac{u_{0}(\alpha_{\text{pure } i}) - r_{2}}{u_{0}(\alpha_{\text{pure } i}) - r_{1}} \right) + \underbrace{\frac{-Q_{0,\text{RMHV1}}(a_{\text{pure } i})}{Q_{0,\text{RMHV1}}(\alpha)}}_{\text{constant packing fraction assumption}} \qquad (65) $$

Assuming a constant packing fraction which implies $u_0(\alpha_{pure i}) = u_0(\alpha)$, the $Q_{0,\text{RMHV1}}(\alpha)$ function is rigorously linear and Eq. (65) yields:

$$ \frac{a_{\text{res}}^{E/\gamma}(T, \mathbf{x})}{RT} = \frac{-1}{r_1 - r_2} \ln \left( \frac{u_0 - r_2}{u_0 - r_1} \right) \cdot \left( \alpha - \sum_{i} x_{i} \alpha_{i} \right) \qquad (66) $$

So that:

$$ \frac{a(T, \mathbf{x})}{b(\mathbf{x})} = \sum_{i} x_{i} \frac{a_{i}(T)}{b_{i}} + \underbrace{\frac{a_{\text{res}}^{E/\gamma}(T, \mathbf{x})}{\frac{1}{r_{1} - r_{2}} \ln \left( \frac{u_{0} - r_{2}}{u_{0} - r_{1}} \right)}}_{\text{constant packing fraction assumption}} \qquad (67) $$

Eq. (67) is the expression of the RMHV1 mixing rule (similar to Eq. (47)) and the $C_{EoS}^{\text{RMHV1}}$ would then have the following expression:

$$ (C_{EoS}^{\text{RMHV1}})_{\text{constant packing fraction assumption}} = \frac{-1}{r_1 - r_2} \ln \left( \frac{u_0 - r_2}{u_0 - r_1} \right) \qquad (68) $$

This expression is strictly identical to the expression of $C_{EoS}^{\text{RMHV1}}$ (see Eq. (59)) found by assuming a constant packing fraction to demonstrate the MHV1 mixing rule.

Let us now forget the constant packing fraction assumption. In this case, $Q_{0,\text{RMHV1}}(\alpha)$ can be considered linear to a first approximation and its graphical representation is a straight line the slope of which can be approximated by the derivative of $Q_{0,\text{RMHV1}}(\alpha)$ with respect to $\alpha$ (denoted $Q'_{0,\text{RMHV1}}(\alpha)$). This derivative is expressed in Eq. (69).

$$ 
\begin{aligned}
Q'_{0,\text{RMHV1}}(\alpha) &= \frac{dQ_{0,\text{RMHV1}}}{d\alpha} = -\frac{1}{r_1 - r_2} \ln \left( \frac{u_0(\alpha) - r_2}{u_0(\alpha) - r_1} \right) + u'_{0}(\alpha) \frac{b}{v - b} \\
&= \underbrace{Q'_{0}(\alpha)}_{(C_{EoS}^{\text{RMHV1}})_{\text{constant packing fraction assumption}}} + u'_{0}(\alpha) \frac{b}{v - b}
\end{aligned}
\qquad (69) $$

It appears that the expression of $Q'_{0,\text{RMHV1}}(\alpha)$ is far different from the expression of $C_{EoS}^{\text{RMHV1}}$ obtained from the constant packing fraction assumption. To support this statement, we provide some estimations of $u_0(\alpha)$, $u'_0(\alpha)$, $Q'_0(\alpha)$ and $Q'_{0,\text{RMHV1}}(\alpha)$ obtained for $\alpha = 14$ (value arbitrarily selected) and the Peng-Robinson EoS, in Eq. (70).

$$ 
\begin{cases}
\alpha = 14 \\
u_0(\alpha) = 1.2042 \\
u'_{0}(\alpha) = -0.0219 \\
Q'_{0}(\alpha) = (C_{EoS}^{\text{RMHV1}})_{\text{constant packing fraction assumption}} = -0.5380 \\
Q'_{0,\text{RMHV1}}(\alpha) = Q'_{0}(\alpha) + u'_{0}(\alpha) \frac{b}{v-b} = -0.6423
\end{cases}
\qquad (70) $$

Two conclusions arise:

- It is not possible to obtain the expression of the RMHV1 mixing rule by assuming constant packing fraction,

- The constant packing fraction assumption is not rigorously supported by the mathematical analysis proposed in this section and therefore, should not be used to derive the MHV1 or RMHV1 mixing rules.
## 5. Study of the literature models based on the HV and ZRP mixing rules in the light of our analysis grid

This section reviews some of the main models that rely on the HV and ZRP mixing rules and that were proposed in the open literature. Most of them are obtained from empirical modifications of the MHV1 mixing rule. However, as previously mentioned, the derivation of the MHV1 mixing rule is rigorous (the resulting solution is unique) and flawless. Consequently, when researchers empirically modified the mixing rule to overcome some drawbacks, they inevitably lost the matching equation and the associated reference pressure. Consequently, most of the existing EoS/g^E mixing rules available on the open literature do not rely on a clearly-established matching equation and reference pressure.

### 5.1. The original HV-NRTL model

This model, already discussed in Section 2, is the model used by Huron and Vidal when they decided to incorporate the purely residual NRTL model in their mixing rule. The starting matching equation is Eq. (1) with an infinite reference pressure. However as previously stated, under infinite pressure, a^E,EoS^comb = 0 and simultaneously, the NRTL model does not contain a combinatorial contribution. In addition, P^E,EoS = 0, since a linear mixing rule is selected for b. In addition, P^E,NRTL = 0. We thus immediately notice that the adoption of a purely residual AC model and a linear mixing rule for b in the HV mixing rules satisfies all the objectives defined in Eq. (35). This is the reason why such a model is able to accurately correlate many types of VLE data and must be considered as physically consistent. However, Vidal himself never carried out such an analysis, and he never gave a clear explanation as to why he chose NRTL. Without wishing to minimize the impact of the research conducted by Huron & Vidal, it can be assumed that they proceeded by trial and error.

Since a^E,EoS^comb = a^E,NRTL^comb = P^E,EoS = P^E,NRTL = 0, the considered matching equation can thus be simply re-written as follows:

$$ \lim_{P \to +\infty} \frac{a_{res}^{E,EoS}(T, P, \mathbf{x}, \alpha)}{RT} = \frac{a_{res}^{E,NRTL}(T, \mathbf{x})}{RT} \qquad (71) $$

And the final mixing rules are:

$$ \left\{ \begin{array}{l} \alpha^{\text{HV-NRTL}}(T, \mathbf{x}) = \sum_{i} x_{i} \alpha_{i}(T) + \frac{a_{\text{res}}^{\text{E,NRTL}}(T, \mathbf{x})}{C_{\text{EoS}}^{\text{HV}}} \\ \\ b = \sum_{i} x_{i} b_{i} \end{array} \right. \qquad (72) $$

Note that this model imposes the adoption of a linear mixing rule for b, and a fixed value for C^HV_EoS. It is believed (see Eq. (54)) that these 2 constraints are not useful and limit the flexibility and applicability of the model.

### 5.2. The Wong-Sandler-UNIFAC approach

Unlike most developers of EoS mixing rules, Wong and Sandler considered as matching equation the equality of excess Helmholtz energies under infinite pressure [see Eq. (21)]. Eventually, they obtained the same mixing rule for α as Huron and Vidal but, as previously stated, could select an expression for b leading to a quadratic composition dependence for the second virial coefficient [refer to Eq. (22)].

The Wong-Sandler approach was combined, by Orbey et al. [36] (among other authors) with the complete UNIFAC model (including both combinatorial and residual contributions). Following our analysis, and remembering that the EoS combinatorial contribution vanishes at infinite pressure, using an activity coefficient model that embeds a combinatorial term is not consistent with the derivation of the model at infinite pressure. Problems are bound to arise when modeling size-asymmetric systems. Following the discussion with the HV mixing

rule, the Wong-Sandler mixing rule should be used with an energetic activity coefficient model like NRTL or Van Laar.

### 5.3. The MHV1 (or MHV2)-UNIFAC model

The MHV1-UNIFAC or MHV2-UNIFAC [9,10] models are implemented in the majority of commercial simulation software and are appreciated for their predictive nature. The initial matching equation is:

$$ g^{E,EoS}(T, P=0, \mathbf{x}, \frac{\alpha}{b}) = g^{E,\gamma}(T, \mathbf{x}) \qquad (73) $$

Under a zero pressure, P^E,EoS = P^E,γ = 0, so that Eq. (73) becomes:

$$ a_{\text{Flory}}^{\text{E,EoS}} + a_{\text{NonFlory}}^{\text{E,EoS}} + a_{\text{res}}^{\text{E,EoS}} = a_{\text{Flory}}^{\text{E,UNIFAC}}(T, \mathbf{x}) + a_{\text{SG}}^{\text{E,UNIFAC}}(T, \mathbf{x}) + a_{\text{res}}^{\text{E,UNIFAC}}(T, \mathbf{x}) \qquad (74) $$

Or equivalently:

$$ \frac{a_{\text{res}}^{\text{E,EoS}}(T, P=0, \mathbf{x}, \alpha)}{RT} = \underbrace{\left( \frac{a_{\text{Flory}}^{\text{E,UNIFAC}}(T, \mathbf{x})}{RT} - \frac{a_{\text{Flory}}^{\text{E,EoS}}(T, P=0, \mathbf{x})}{RT} \right)}_{\text{DCT}} + \underbrace{\left( \frac{a_{\text{SG}}^{\text{E,UNIFAC}}(T, \mathbf{x})}{RT} - \frac{a_{\text{Non Flory}}^{\text{E,EoS}}(T, P=0, \mathbf{x})}{RT} \right)}_{\text{DCT}_2} + \frac{a_{\text{res}}^{\text{E,UNIFAC}}(T, \mathbf{x})}{RT} \qquad (75) $$

It is noticeable that the matching equation introduces two differences, noted DCT and DCT2, between combinatorial contributions. The first one is the classical DCT, the second one is the difference between the "non-Flory" term originating from the EoS and a Staverman-Guggenheim contribution (denoted DCT2 below).

Eq. (75) is consistent but it suffers from the MHV1 shortcoming discussed earlier. Indeed, the combinatorial contributions from the EoS and from the UNIFAC model do not match and consequently, as shown in literature [5,41,51,53], athermal mixtures and more specifically size-asymmetric systems will not be properly correlated.

The final mixing rules can be written:

$$ \left\{ \begin{array}{l} \alpha^{\text{MHV1-UNIFAC}}(T, \mathbf{x}) = \sum_{i} x_{i} \alpha_{i}(T) + \frac{\text{DCT} + \text{DCT}_2 + \frac{a_{\text{res}}^{\text{E,UNIFAC}}(T, \mathbf{x})}{RT}}{C_{\text{EoS}}^{\text{MHV1}}} \\ \\ b = \text{free choice} \end{array} \right. \qquad (76) $$

We conclude that this type of model is problematic, in particular if it is used to correlate athermal size-asymmetric systems, e.g., mixtures of alkanes with different asymmetries.

### 5.4. The PSRK model

The famous PSRK model was introduced in 1991 by Holderbaum and Gmehling [12]. The main difference between this model and the MHV1-UNIFAC model is the choice of the mixing rule constant:

$$ \left\{ \begin{array}{l} a^{\text{PSRK}}(T, \mathbf{x}) = \sum_{i} x_{i} \alpha_{i}(T) + \frac{\text{DCT} + \text{DCT}_2 + \frac{a_{\text{res}}^{\text{E,UNIFAC}}(T, \mathbf{x})}{RT}}{C_{\text{PSRK}}} \\ \\ C_{\text{PSRK}} = -0.64663 \\ \\ b = \sum_{i} x_{i} b_{i} \end{array} \right. \qquad (77) $$

Therefore, the PSRK model involves undesirable DCT and DCT2 terms that are detrimental to the correlation of size-asymmetric systems. The authors provide an empirical justification for the change in the constant
value $C_{PSRK}$ compared to the value of "-0.593" proposed by Michelsen for the SRK EoS. They wrote: "This value is changed to C = -0.64663 in the PSRK equation, which yields better results at higher pressure" [12]. Interestingly, the new value of $C_{PSRK}$ they introduced is in the acceptable range of $C_{SRK}^{MHV1}$ values defined in Eq. (48); it is closer to the value of $C_{SRK}^{HV} \approx -0.693$ and thus corresponds to the first derivative of the $Q_0$ function linearized over a large range of $\alpha$ values.

Some years later, Fischer and Gmehling proposed an a posteriori justification of their mixing rule [67] based on the main assumption that the packing fraction is constant upon mixing: $\frac{b}{v} = \frac{b_i}{v_{pure i}}$. Inherited from Guggenheim's quasi lattice theory, this assumption is discussed in detail in a paper by Péneloux et al. published in 1989 [68] which summarizes it as follows: "the lattice on which the pure components are placed undergoes no distortion during the mixing process". Fischer and Gmehling considered the average value $u = 1.1$ for deriving the PSRK mixing rule (the "-0.593 value, proposed by Michelsen, would correspond to a packing fraction $u = 1.234$). In 1996, Novenario et al. [69], using diverse data for a wide range of systems, proposed for the Peng-Robinson EoS a value of $u = 1.15$. As discussed in Section 4, we considered that the constant packing fraction assumption is not supported by a rigorous mathematical analysis and should be disregarded.

To conclude this section devoted to the PSRK model, let us recall that in 2002, Chen et al. [14] published a new version in which they empirically skipped the DCT and DCT2 terms. They also introduced a quadratic mixing rule for the covolume with the generalized formulation of the combining rule proposed by Lorentz for $b_{ij}$ [see Eq. (78)].

$$ \left\{ \begin{aligned} b_m &= \sum_i \sum_j z_i z_j b_{ij} \\ b_{ij} &= \left( \frac{b_i^{1/s} + b_j^{1/s}}{2} \right)^s \quad \text{with } s = \frac{4}{3} \end{aligned} \right. \qquad (78) $$

This new version of the PSRK (derived empirically) model is thus totally compatible with the mixing rules we mathematically derived in this paper [refer to Eq. (51)].

### 5.5. The VTPR model

VTPR was introduced in the early 2000's by the group that developed the PSRK model [18,19,70]. Their authors decided to remove the DCT term present in the original MHV1 mixing rule in order to improve the correlation of athermal systems. They also used a quadratic function of the composition $x$ for the parameter $b$. Eq. (79) was obtained.

$$ \left\{ \begin{aligned} \alpha^{\text{VTPR}}(T, \mathbf{x}) &= \sum_i x_i \alpha_i(T) + \frac{\left(\frac{a_{\text{res}}^{\text{UNIFAC}}(T, \mathbf{x})}{RT}\right)}{C_{\text{PR}}^{\text{MHV1}}}, & C_{\text{PR}}^{\text{MHV1}} &= -0.530 \\ b &= \sum_i \sum_j x_i x_j \left(\frac{b_i^{3/4} + b_j^{3/4}}{2}\right)^{4/3} \end{aligned} \right. \qquad (79) $$

At this step, it seems difficult to define the matching equation on which such a model relies. Indeed, by equating the complete $g^E$, as Michelsen did, the DCT would not disappear. By equating the residual contributions only, as we did in this paper, a value lower than or equal to $C_{EoS}^{HV} = -0.623$ should be preferred for the universal constant of the VTPR mixing rule.

Following the demonstration presented in [70], it clearly appears that these authors had in mind to equate the residual term of the AC model with the EoS residual term. However, as discussed above (see the text located below Eq. (29)), the DCT is not the difference between the combinatorial term from the AC model and the one from the EoS. The DCT only includes the Flory part of the free-volume combinatorial contribution stemming from the EoS. By removing the DCT, the Flory terms are removed from the matching equation but not the "Non Flory" combinatorial terms.

In conclusion, the matching equation of the VTPR model can be written:

$$ a_{\text{Non Flory}}^{\text{EoS}} + a_{\text{res}}^{\text{EoS}} = a_{\text{res}}^{\text{UNIFAC}}(T, \mathbf{x}) \qquad (80) $$

The magnitude of the non-Flory part of the combinatorial part is generally small, but strictly speaking, this equation is unbalanced and therefore not entirely consistent.

However, as also shown from the authors, this problem is corrected via the inevitable parameter estimation in the VTPR model. The parameters of the UNIFAC residual term must be re-estimated (and indeed this is what is done in literature for VTPR). We cannot use the original parameters of UNIFAC, even though they only appear in the residual term, as they have been fitted for the complete UNIFAC model containing both combinatorial and residual terms. This parameter estimation can accommodate for the "small inconsistency" implied by Eq. (80).

### 5.6. The NRTL-PR model

Neau et al. have proposed a model that is similar to the VTPR approach: it combines the PR EoS and the MHV1 mixing rule in which the DCT was set to zero [23,71–73]. The main difference lies in the choice of the AC model: VTPR's authors used the residual term of the UNIFAC model while Neau et al. used the NRTL model which is a purely residual AC model. Like the authors of VTPR, Neau et al. intended to identify the residual EoS contribution with the residual NRTL model but actually used a matching equation similar to Eq. (80) that is unbalanced because it involves a "Non Flory" contribution on the left-hand side, and nothing equivalent on the right-hand side.

### 5.7. The UMR-PRU model

The UMR-PRU model is a predictive EoS based on a volume-translated version of the PR EoS, combined with the UNIFAC AC model through an empirically modified version of the MHV1 mixing rule [20,21,74]. This model represents a hybrid between PSRK and VTPR. Their authors decided to delete the DCT contribution from the PSRK model but to keep the $DCT_2$ (see the definition of this term in Eq. (75)) and simultaneously to keep the universal constant ($-0.53$) proposed by Michelsen.

The resulting mixing rule is:

$$ \left\{ \begin{aligned} \alpha^{\text{UMR-PRU}}(T, \mathbf{x}) &= \sum_i x_i \alpha_i(T) + \frac{\mathbf{DCT}_2 + \left(\frac{a_{\text{res}}^{\text{UNIFAC}}(T, \mathbf{x})}{RT}\right)}{-0.53} \\ \text{with : } \mathbf{DCT}_2 &= \frac{a_{\text{SG}}^{\text{UNIFAC}}(T, \mathbf{x})}{RT} - \frac{a_{\text{Non Flory}}^{\text{EoS}}(T, P = 0, \mathbf{x})}{RT} \\ b &= \sum_i \sum_j x_i x_j \left(\frac{b_i^{1/2} + b_j^{1/2}}{2}\right)^2 \end{aligned} \right. \qquad (81) $$

Consequently, this mixing rule suffers from the $DCT_2$ issue. The magnitude of the $DCT_2$ term is however small, compared to DCT. It should be mentioned that the UMR-PRU model has been used with the existing UNIFAC parameters and no new parameterization has taken place. This further confirms the importance of using only the residual part of the external activity coefficient models in the EoS/g$^E$ mixing rules (or EoS/$\alpha^E$ models, as they should more accurately be termed). The same conclusion is obtained with the results using the LCVM model which, as shown in the literature, also, in a rather empirical way, eliminates the DCT and, like UMR-PRU, uses existing UNIFAC parameters without refitting.
<table><caption>**Table 2**
Synthesis of the characteristics of the many advanced mixing rules published since 1980.</caption><thead><tr><th>Model</th><th>Type of mixing rule:<br>a/b = &Sigma;x<sub>i</sub>a<sub>i</sub>/b<sub>i</sub> + E</th><th>Matching equations (note that for MHV1-type mixing rule while P<sub>ref</sub> = +&infin; for HV type mixing rule)</th><th>Presence of a double comb. term DCT<sub>1</sub> in the matching equations<sup>(1)</sup></th><th>Presence of a double comb. term DCT<sub>2</sub> in the matching equations<sup>(2)</sup></th><th>Type of mixing rule for b</th><th>Capacity to model size-asymmetric mixtures</th></tr><tr><th></th><th>HV type:<br>E = g<sup>E,Y</sup>/RT / C</th><th>MHV1 type:<br>E = g<sup>E,Y</sup>/RT - &Sigma;x<sub>i</sub>ln(b<sub>i</sub>/b)</th><th>a<sup>E,Y</sup> = g<sup>E,EoS</sup> at (T, P<sub>ref</sub>, x)</th><th>a<sup>E,Y</sup><sub>comb</sub> = a<sup>E,EoS</sup><sub>comb</sub> at (T, P<sub>ref</sub>, x)</th><th>a<sup>E,Y</sup><sub>res</sub> = a<sup>E,EoS</sup><sub>res</sub> at (T, P<sub>ref</sub>, x)</th><th></th><th></th></tr></thead><tbody><tr><td>PR + HV + NRTL</td><td>{ PR EoS<br>C = -0.62</td><td></td><td>&#10003; (P<sub>ref</sub> = +&infin;)</td><td>&#10003; (comb. terms are zero)</td><td>No</td><td>No</td><td>Linear</td><td>Rather good</td></tr><tr><td>PR + Wong Sandler + UNIFAC</td><td>{ PR EoS<br>C = -0.62</td><td></td><td>&#10003; (P<sub>ref</sub> = +&infin;)</td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003;</td><td>Complex (neither Linear, nor quadratic)</td><td>Low</td></tr><tr><td>PR + MHV1 or MHV2 + UNIFAC</td><td></td><td>{ PR EoS<br>C = -0.53</td><td>&#10003; (P<sub>ref</sub> &approx; 0)</td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003;</td><td>Linear or quadratic</td><td>Low</td></tr><tr><td>PSRK</td><td></td><td>{ SRK EoS<br>C = -0.64</td><td>&#10003; (P<sub>ref</sub> &approx; 0)<sup>(3)</sup></td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003;</td><td>Linear</td><td>Low</td></tr><tr><td>VTPR</td><td>{ PR EoS<br>C = -0.53</td><td></td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003; (P<sub>ref</sub> &approx; 0)*</td><td>No</td><td>&#10003;</td><td>Quadratic</td><td>Rather good</td></tr><tr><td>NRTL-PR</td><td>{ PR EoS<br>C = -0.53</td><td></td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003; (P<sub>ref</sub> &approx; 0)*</td><td>No</td><td>&#10003;</td><td>Linear</td><td>Rather good</td></tr><tr><td>UMR-PRU</td><td>{ PR EoS<br>C = -0.53</td><td></td><td>Not satisfied</td><td>Not satisfied</td><td>&#10003; (P<sub>ref</sub> &approx; 0)*</td><td>No</td><td>&#10003;</td><td>Quadratic</td><td>Rather good</td></tr><tr><td>E-PPR78 (PR + HV + Van Laar)</td><td>{ PR EoS<br>C = -0.62</td><td></td><td>&#10003; (P<sub>ref</sub> = +&infin;)</td><td>&#10003; (comb. terms are zero)</td><td>&#10003; (P<sub>ref</sub> = +&infin;)</td><td>No</td><td>No</td><td>Linear</td><td>Rather good</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table 2 (continued)**</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Type of mixing rule:</th><th colspan="3">Matching equations (note that for MHV1-type mixing rule while P<sub>ref</sub> = +∞ for HV type mixing rule)</th><th rowspan="2">Presence of a double comb. term DCT<sub>1</sub> in the matching equations<sup>(1)</sup></th><th rowspan="2">Presence of a double comb. term DCT<sub>2</sub> in the matching equations<sup>(2)</sup></th><th rowspan="2">Type of mixing rule for b</th><th rowspan="2">Capacity to model size-asymmetric mixtures</th></tr><tr><th>a/b = &sum;<sub>i</sub>x<sub>i</sub>a<sub>i</sub>/b<sub>i</sub> + E</th><th>MHV1 type:</th><th>a<sup>E<sub>Y</sub></sup> = g<sup>E,EoS</sup> at (T, P<sub>ref</sub>, x)</th><th>a<sup>E<sub>Y</sub></sup><sub>comb</sub> = a<sup>E,EoS</sup> at (T, P<sub>ref</sub>, x)</th><th>a<sup>E<sub>Y</sub></sup><sub>res</sub> = a<sup>E,EoS</sup> at (T, P<sub>ref</sub>, x)</th></tr></thead><tbody><tr><td></td><td>HV type:<br>E = g<sup>E,Y</sup>/RT<br>C</td><td>MHV1 type:<br>E =<br>g<sup>E,Y</sup>/RT - &sum;<sub>i</sub>x<sub>i</sub>ln(b<sub>i</sub>/b)<br>C</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>tc-PR (PR + HV + residual part of g<sup>E,Y</sup>)</td><td>{ PR EoS<br>C = -0.62</td><td></td><td>✓ (P<sub>ref</sub> = +∞)</td><td>✓ (comb. terms are zero)</td><td>✓ (P<sub>ref</sub> = +∞)</td><td>No</td><td>No</td><td>Quadratic</td><td>Rather good</td></tr></tbody></table>

(1) DCT₁ denotes a nonzero difference in the matching equation between the Flory term coming from the EoS and the Flory term coming from the AC model.

(2) (1) DCT₂ denotes a nonzero difference in the matching equation between the non-Flory term coming from the EoS and the non-Flory term (e.g., Sta-verman–Guggenheim) coming from the AC model.

(3) As discussed earlier, the PSRK model involves a C constant in the mixing rule which is compatible with the MHV1 approach and a linearization of the Q₀ function over a large range of α value.

## 5.8. The tc-PR model

The tc-PR EoS for mixtures associates the translated-consistent PR EoS and the universal mixing rule proposed in Section 3.3 [42,75] [see Eq. (40)]. The covolume is expressed as a quadratic function of the composition (Eq. (82)).

$$b = \sum_{i} \sum_{j} x_{i}x_{j} \left( \frac{b_{i}^{2/3} + b_{j}^{2/3}}{2} \right)^{3/2} \quad (82)$$

Different AC models were combined to this mixing rule in order to identify the one leading to optimal performances. In ref. [42], the performances of the three tc-PR versions (tc-PR EoS + unified mixing rule defined in Eq. (40) + residual part of an AC model), that are the tc-PR NRTL, tc-PR residual-UNIQUAC and tc-PR residual-Wilson models, were compared against a high-quality reference database [62] containing VLE, LLE, enthalpy of mixing, heat capacity of mixing, critical data for 200 binary systems. The tc-PR-residual-Wilson EoS was identified as the most promising of the three tc-PR EoS versions making it one of the most accurate EoS available today.

## 6. Conclusion

This work constitutes a comprehensive acknowledgment of the significance of the contributions made by Huron-Vidal and Michelsen regarding the formulation of mixing rules for cubic EoS. They introduced the fundamental concept of the “matching equation”, which enables the derivation of a mixing rule for the a/b (or α = a/(bRT)) parameter. Table 2 attempts to summarize the characteristics of the many advanced mixing rules published since 1980 and the main question we addressed in this article is: are the different versions of the Huron-Vidal and MHV mixing rules fully consistent from a physical viewpoint? Our main conclusions are as follows:

- The Huron and Vidal (HV) proposition is physically sound provided that only the residual part of an activity coefficient model is incorporated inside this mixing rule. This is precisely what they did when they combined the SRK EoS with the purely residual NRTL model. As this paper points out, incorporating a complete activity coefficient model that includes both a combinatorial and a residual part in the

HV mixing rules leads to an unbalanced matching equation that will cause many problems when correlating/predicting phase behaviour (is it always problematic to equate quantities – from the EoS and the AC model – that do not have the same physical meaning). In addition, the derivation proposed by Huron and Vidal imposes the adoption of a linear mixing rule for b, which we believe is not a needed constraint, which moreover limits the flexibility of the overall model.

- The ZRP mixing rule, as introduced by Michelsen, is fully consistent from a physical point of view. He used a matching equation to align the complete g^E from the EoS and from the AC model, which are quantities with the same physical significance. This paper however highlights that it is not sufficient to satisfy the unique equation (g^E,EoS = g^E,Y) at a reference pressure in order to satisfactorily predict or correlate the phase behaviour and thermodynamic properties of all types of binary systems (athermal, regular or mixed-effect solutions). We demonstrated that a suitable definition of an EoS/g^E mixing rule should be based on the requirement that the three individual contributions (combinatorial, residual and PV^E) from the EoS are equal to the corresponding contributions of the AC model. Failure to comply with this rule will result in inconsistencies when correlating some systems, particularly those that are size-asymmetric for which the concern is known as the double combinatorial term issue in the literature.

After reviewing the pioneering works of Huron, Vidal and Michelsen, our ambition was to develop "drawback free" advanced mixing rules. We soon realised, however that, although desirable, it will never be possible to simultaneously equate the three separate contributions of g^E,EoS (combinatorial, residual and excess volume) to those of an AC model. This is simply because none of the typically-used AC models has a combinatorial part which is a Flory free-volume expression and because the excess volume from an AC model is always zero. The unique remaining possibility to develop "drawback free" mixing rules is thus to only equate the residual contributions of the EoS and AC model, and this is what we proposed here.

As a remarkable result, the mixing rule for a/b we derived after selecting the reference pressure P<sub>ref</sub> = +∞ has the same expression as the mixing rule associated with the reference pressure P<sub>ref</sub> = 0. The resulting reference-pressure-independent mixing rule which unifies the HV and MHV mixing rules was called the UHVM (Unified Huron Vidal
Michelsen) mixing rule.

CRediT authorship contribution statement

**Romain Privat:** Writing – review & editing, Writing – original draft, Formal analysis, Conceptualization.
**Jean-Noël Jaubert:** Writing – review & editing, Writing – original draft, Formal analysis, Conceptualization.
**Georgios M. Kontogeorgis:** Writing – review & editing, Formal analysis, Conceptualization.

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

No data was used for the research described in the article.

References

[1] J.D. Van der Waals, The equation of state for gases and liquids (December 12, 1910). Nobel Lectures Including Presentation Speeches and Laureates' Biographies, Elsevier, Amsterdam, 1967, pp. 254–265. Physics, 1901-1921, https://www.nobelprize.org/uploads/2018/06/waals-lecture.pdf.

[2] J.D. Van der Waals, J.S. Rowlinson, On the Continuity of the Gaseous and Liquid States, Dover Publications, Mineola, N.Y., 2004.

[3] J.M.H. Levelt Sengers, How Fluids Unmix: Discoveries by the School of Van der Waals and Kamerlingh Onnes, Koninklijke Nederland Akademie van Wetenschappen, Amsterdam, 2002.

[4] 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, https://doi.org/10.1016/j.fluid.2022.113697.

[5] G.M. Kontogeorgis, G.K. Folas, Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories, 1st ed., Wiley, 2010 https://doi.org/10.1002/9780470747537.

[6] 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, https://doi.org/10.1016/0378-3812(79)80001-1.

[7] M.L. Michelsen, A method for incorporating excess gibbs energy models in equations of state, Fluid. Phase Equilib. 60 (1990) 47–58, https://doi.org/10.1016/0378-3812(90)85042-9.

[8] M.L. Michelsen, A modified Huron-Vidal mixing rule for cubic equations of state, Fluid. Phase Equilib. 60 (1990) 213–219, https://doi.org/10.1016/0378-3812(90)85053-D.

[9] S. Dahl, M.L. Michelsen, High-pressure vapor-liquid equilibrium with a UNIFAC-based equation of state, AIChE J. 36 (1990) 1829–1836, https://doi.org/10.1002/aic.690361207.

[10] S. Dahl, A. Fredenslund, P. Rasmussen, The MHV2 model: a UNIFAC-based equation of state model for prediction of gas solubility and vapor-liquid equilibria at low and high pressures, Ind. Eng. Chem. Res. 30 (1991) 1936–1945, https://doi.org/10.1021/ie00056a041.

[11] D.S.H. Wong, S.I. Sandler, A theoretically correct mixing rule for cubic equations of state, AIChE J. 38 (1992) 671–680, https://doi.org/10.1002/aic.690380505.

[12] T. Holderbaum, J. Gmehling, PSRK: a group contribution equation of state based on UNIFAC, Fluid. Phase Equilib. 70 (1991) 251–265, https://doi.org/10.1016/0378-3812(91)85038-V.

[13] S. Horstmann, K. Fischer, J. Gmehling, PSRK group contribution equation of state: revision and extension III, Fluid. Phase Equilib. 167 (2000) 173–186, https://doi.org/10.1016/S0378-3812(99)00333-7.

[14] J. Chen, K. Fischer, J. Gmehling, Modification of PSRK mixing rules and results for vapor-liquid equilibria, enthalpy of mixing and activity coefficients at infinite dilution, Fluid. Phase Equilib. 200 (2002) 411–429, https://doi.org/10.1016/S0378-3812(02)00048-1.

[15] J. Gmehling, Potential of thermodynamic tools (group contribution methods, factual data banks) for the development of chemical processes, Fluid. Phase Equilib. 210 (2003) 161–173, https://doi.org/10.1016/S0378-3812(03)00177-8.

[16] S. Horstmann, A. Jabloniec, J. Krafcyk, K. Fischer, J. Gmehling, PSRK group contribution equation of state: comprehensive revision and extension IV, including critical constants and α-function parameters for 1000 components, Fluid. Phase Equilib. 227 (2005) 157–164, https://doi.org/10.1016/j.fluid.2004.11.002.

[17] 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 UNIFAC and the t-mPR equation of state, Fluid. Phase Equilib. 92 (1994) 75–106, https://doi.org/10.1016/0378-3812(94)80043-X.

[18] J. Ahlers, J. Gmehling, Development of an universal group contribution equation of state, Fluid. Phase Equilib. 191 (2001) 177–188, https://doi.org/10.1016/S0378-3812(01)00626-4.

[19] B. Schmid, J. Gmehling, Present status of the group contribution equation of state VTPR and typical applications for process development, Fluid. Phase Equilib. 425 (2016) 443–450, https://doi.org/10.1016/j.fluid.2016.06.042.

[20] E. Voutsas, K. Magoulas, D. Tassios, Universal mixing rule for cubic equations of state applicable to symmetric and asymmetric systems: results with the Peng–Robinson equation of state, Ind. Eng. Chem. Res. 43 (2004) 6238–6246, https://doi.org/10.1021/ie049580p.

[21] E. Voutsas, V. Loui, C. Boukouvalas, K. Magoulas, D. Tassios, Thermodynamic property calculations with the universal mixing rule for EoS/GE models: results with the Peng–Robinson EoS and a UNIFAC model, Fluid. Phase Equilib. 241 (2006) 216–228, https://doi.org/10.1016/j.fluid.2005.12.028.

[22] J.-N. Jaubert, F. Mutelet, VLE predictions with the Peng-Robinson equation of state and temperature-dependent $k_{ij}$ calculated through a group contribution method, Fluid. Phase Equilib. 224 (2004) 285–304, https://doi.org/10.1016/j.fluid.2004.06.059.

[23] E. Neau, C. Nicolas, J.-N. Jaubert, F. Mutelet, The generalized NRTL model associated with the Peng-Robinson equation of state to predict liquid-liquid equilibria between hydrocarbons, water and ethylene glycol, Pol. J. Chem. 80 (2006) 27–35.

[24] G. Soave, Equilibrium constants from a modified Redlich-Kwong equation of state, Chem. Eng. Sci. 27 (1972) 1197–1203, https://doi.org/10.1016/0009-2509(72)80096-4.

[25] D.-Y. Peng, D.B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fund. 15 (1976) 59–64, https://doi.org/10.1021/i160057a011.

[26] D.B. Robinson, H.-J. Ng, A.D. Leu, Behavior of $CH_4-CO_2-H_2S$ mixtures at sub-ambient temperatures (RR-47), Res. Rep. GPA (1981) 1–38.

[27] A. Pina-Martinez, R. Privat, J.-N. Jaubert, D.-Y. Peng, Updated versions of the generalized Soave α-function suitable for the Redlich-Kwong and Peng-Robinson equations of state, Fluid. Phase Equilib. 485 (2019) 264–269, https://doi.org/10.1016/j.fluid.2018.12.007.

[28] G.M. Kontogeorgis, R. Privat, J.-N. Jaubert, Taking another look at the Van der Waals equation of state–almost 150 years later, J. Chem. Eng. Data 64 (2019) 4619–4637, https://doi.org/10.1021/acs.jcbed.9b00264.

[29] J. Vidal, Thermodynamics: Applications in Chemical Engineering and the Petroleum Industry, Editions Technip, Paris, 2003.

[30] M.L. Michelsen, J.M. Mollerup, Thermodynamic Models: Fundamentals & Computational Aspects, 2nd ed., Tie-Line Publications, Holte, 2007.

[31] O. Pfohl, Evaluation of an improved volume translation for the prediction of hydrocarbon volumetric properties, Fluid. Phase Equilib. 163 (1999) 157–159, https://doi.org/10.1016/S0378-3812(99)00199-5.

[32] P. Ungerer, H.B. De Sant'Ana, Reply to the letter to the editor by O. Pfohl about the paper "evaluation of an improved volume translation for the prediction of hydrocarbon volumetric properties" [FPE 154, 193–204 (1999)], Fluid. Phase Equilib. 163 (1999) 161–162, https://doi.org/10.1016/S0378-3812(99)00198-3.

[33] V. Kalikhman, D. Kost, I. Polishuk, About the physical validity of attaching the repulsive terms of analytical EOS models by temperature dependencies, Fluid. Phase Equilib. 293 (2010) 164–167, https://doi.org/10.1016/j.fluid.2010.03.003.

[34] J.-N. Jaubert, R. Privat, Y. Le Guennec, L. Coniglio, Note on the properties altered by application of a Péneloux-type volume translation to an equation of state, Fluid. Phase Equilib. 419 (2016) 88–95, https://doi.org/10.1016/j.fluid.2016.03.012.

[35] D.S.H. Wong, H. Orbey, S.I. Sandler, Equation of state mixing rule for nonideal mixtures using available activity coefficient model parameters and that allows extrapolation over large ranges of temperature and pressure, Ind. Eng. Chem. Res. 31 (1992) 2033–2039, https://doi.org/10.1021/ie00008a027.

[36] H. Orbey, S.I. Sandler, D.S.H. Wong, Accurate equation of state predictions at high temperatures and pressures using the existing UNIFAC model, Fluid. Phase Equilib. 85 (1993) 41–54, https://doi.org/10.1016/0378-3812(93)80003-6.

[37] H. Orbey, S.I. Sandler, On the combination of equation of state and excess free energy models, Fluid. Phase Equilib. 111 (1995) 53–70, https://doi.org/10.1016/0378-3812(95)02764-6.

[38] H. Orbey, S.I. Sandler, A comparison of various cubic equation of state mixing rules for the simultaneous description of excess enthalpies and vapor-liquid equilibria, Fluid. Phase Equilib. 121 (1996) 67–83, https://doi.org/10.1016/0378-3812(96)03030-0.

[39] P. Coutsikos, N.S. Kalospiros, D.P. Tassios, Capabilities and limitations of the Wong-Sandler mixing rules, Fluid. Phase Equilib. 108 (1995) 59–78, https://doi.org/10.1016/0378-3812(94)02675-Q.

[40] G.M. Kontogeorgis, I.G. Economou, Equations of state: from the ideas of Van der Waals to association theories, J. Supercrit. Fluids. 55 (2010) 421–437, https://doi.org/10.1016/j.supflu.2010.10.023.

[41] G.M. Kontogeorgis, P. Coutsikos, Thirty years with EoS/gE models—What have we learned? Ind. Eng. Chem. Res. 51 (2012) 4119–4142, https://doi.org/10.1021/ie2015119.

[42] A. Pina-Martinez, R. Privat, I.K. Nikolaidis, I.G. Economou, J.-N. Jaubert, What is the optimal activity coefficient model to be combined with the translated-consistent Peng–Robinson equation of state through advanced mixing rules? Cross-comparison and grading of the Wilson, UNIQUAC, and NRTL aE models against a benchmark database involving 200 binary systems, Ind. Eng. Chem. Res. 60 (2021) 17228–17247, https://doi.org/10.1021/acs.iecr.1c30003.

[43] J.-N. Jaubert, R. Privat, N. Juntarachat, General reflection on critical negative azeotropy and upgrade of the Bancroft's rule with application to the acetone+ chloroform binary system, J. Supercrit. Fluids 94 (2014) 17–29, https://doi.org/10.1016/j.supflu.2014.06.014.
[44] R. Privat, J.-N. Jaubert, G.M. Kontogeorgis, The secret of the Wilson equation, Fluid. Phase Equilib. 579 (2024) 114018, https://doi.org/10.1016/j.fluid.2023.114018.

[45] R. Privat, J.-N. Jaubert, G.M. Kontogeorgis, Can liquid-liquid equilibria be predicted by the combination of a cubic equation of state and a gE model not suitable for liquid-liquid equilibria? Fluid. Phase Equilib. (2024) 114249 https://doi.org/10.1016/j.fluid.2024.114249.

[46] H. Segura, A. Mejía, R. Reich, J. Wisniak, S. Loras, Isobaric vapor-liquid equilibria and densities for the binary systems oxolane + ethyl 1,1-dimethylethyl ether, oxolane + 2-propanol and propan-2-one + trichloromethane, Phys. Chem. Liquids 41 (2003) 283–301, https://doi.org/10.1080/0031910021000044456.

[47] I.B. Rabinovich, P.N. Nikolaev, Isotopic effect in the thermodynamic properties of liquid mixtures. I. The isotopic shift in the vapor pressure-composition diagram of mixtures of deuterochloroform with acetone, Zhur. Fiz. Khim. 34 (1960). https://www.osti.gov/biblio/4104884 (Accessed 22 April 2025).

[48] G. Soave, Application of equations of state and the theory of group solutions to phase equilibrium prediction, Fluid. Phase Equilib. 87 (1993) 23–35, https://doi.org/10.1016/0378-3812(93)85016-F.

[49] S.H. Huang, H.M. Lin, K.C. Chao, Solubility of carbon dioxide, methane, and ethane in n-eicosane, J. Chem. Eng. Data 33 (1988) 145–147, https://doi.org/10.1021/je00052a026.

[50] C.J. Peters, J.L. De Roo, R.N. Lichtenthaler, Measurements and calculations of phase equilibria of binary mixtures of ethane + eicosane. Part I: vapour + liquid equilibria, Fluid. Phase Equilib. 34 (1987) 287–308, https://doi.org/10.1016/0378-3812(87)80037-7.

[51] G.M. Kontogeorgis, P. Coutsikos, Comments on “predictions of activity coefficients of nearly athermal binary mixtures using cubic equations of state, Ind. Eng. Chem. Res. 44 (2005) 3374–3375, https://doi.org/10.1021/ie0500063.

[52] J.-N. Jaubert, R. Privat, SAFT and cubic EoS: type of deviation from ideality naturally predicted in the absence of BIPs. Application to the modelling of athermal mixtures, Fluid. Phase Equilib. 533 (2021) 112924, https://doi.org/10.1016/j.fluid.2020.112924.

[53] G.M. Kontogeorgis, P.M. Vlamos, An interpretation of the behavior of EoS/S<sup>E</sup> models for asymmetric systems, Chem. Eng. Sci. 55 (2000) 2351–2358, https://doi.org/10.1016/S0009-2509(99)00472-8.

[54] J.-N. Jaubert, L. Coniglio, The group contribution concept: a useful tool to correlate binary systems and to predict the phase behavior of multicomponent systems involving supercritical CO<sub>2</sub> and fatty acids, Ind. Eng. Chem. Res. 38 (1999) 5011–5018, https://doi.org/10.1021/ie990544d.

[55] S. Vitu, R. Privat, J.-N. Jaubert, F. Mutelet, Predicting the phase equilibria of CO<sub>2</sub> + hydrocarbon systems with the PPR78 model (PR EoS and kij calculated through a group contribution method), J. Supercrit. Fluids 45 (2008) 1–26, https://doi.org/10.1016/j.supflu.2007.11.015.

[56] J.-W. Qian, R. Privat, J.-N. Jaubert, Predicting the phase equilibria, critical phenomena, and mixing enthalpies of binary aqueous systems containing alkanes, cycloalkanes, aromatics, alkenes, and gases (N<sub>2</sub>, CO<sub>2</sub>, H<sub>2</sub>S, H<sub>2</sub>) with the PPR78 equation of state, Ind. Eng. Chem. Res. 52 (2013) 16457–16490, https://doi.org/10.1021/ie402541h.

[57] J.-N. Jaubert, R. Privat, F. Mutelet, Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach, AIChE J. 56 (2010) 3225–3235, https://doi.org/10.1002/aic.12232.

[58] J.-W. Qian, R. Privat, J.-N. Jaubert, P. Duchet-Suchaux, Enthalpy and heat capacity changes on mixing: fundamental aspects and prediction by means of the PPR78 cubic equation of state, Energy Fuels 27 (2013) 7150–7178, https://doi.org/10.1021/ef401605c.

[59] J.-N. Jaubert, R. Privat, Relationship between the binary interaction parameters (kij) of the Peng-Robinson and those of the Soave-Redlich-Kwong equations of state: application to the definition of the PR2SRK model, Fluid. Phase Equilib. 295 (2010) 26–37, https://doi.org/10.1016/j.fluid.2010.03.037.

[60] X. Xu, S. Lasala, R. Privat, J.-N. Jaubert, E-PPR78: a proper cubic EoS for modelling fluids involved in the design and operation of carbon dioxide capture and storage (CCS) processes, Int. J. Greenh. Gas Control 56 (2017) 126–154, https://doi.org/10.1016/j.ijggc.2016.11.015.

[61] J.-N. Jaubert, J.-W. Qian, S. Lasala, R. Privat, The impressive impact of including enthalpy and heat capacity of mixing data when parameterising equations of state. Application to the development of the E-PPR78 (Enhanced-Predictive-Peng-Robinson-78) model, Fluid. Phase Equilib. 560 (2022) 113456, https://doi.org/10.1016/j.fluid.2022.113456.

[62] J.-N. Jaubert, Y. Le Guennec, A. Pina-Martinez, N. Ramirez-Velez, S. Lasala, B. Schmid, I.K. Nikolaidis, I.G. Economou, R. Privat, Benchmark database containing binary-system-high-quality-certified data for cross-comparing thermodynamic models and assessing their accuracy, Ind. Eng. Chem. Res. 59 (2020) 14981–15027, https://doi.org/10.1021/acs.iecr.0c01734.

[63] K. Tochigi, K. Kojima, T. Sako, Prediction of vapor-liquid equilibria in polymer solutions using EOS-group contribution model consistent with the second virial coefficient condition, Fluid. Phase Equilib. 117 (1996) 55–60, https://doi.org/10.1016/0378-3812(95)02936-2.

[64] H. Orbey, S.I. Sandler, Analysis of excess free energy based equations of state models, AIChE J. 42 (1996) 2327–2334, https://doi.org/10.1002/aic.690420822.

[65] J. Mollerup, A note on the derivation of mixing rules from excess gibbs energy models, Fluid. Phase Equilib. 25 (1986) 323–327, https://doi.org/10.1016/0378-3812(86)80007-3.

[66] J. Mollerup, A note on excess Gibbs energy models, equations of state and the local composition concept, Fluid. Phase Equilib. 7 (1981) 121–138, https://doi.org/10.1016/0378-3812(81)85017-0.

[67] K. Fischer, J. Gnehling, Further development, status and results of the PSRK method for the prediction of vapor-liquid equilibria and gas solubilities, Fluid. Phase Equilib. 112 (1995) 1–22, https://doi.org/10.1016/0378-3812(95)02792-D.

[68] A. Peneloux, W. Abdoul, E. Rauzy, Excess functions and equations of state, Fluid. Phase Equilib. 47 (1989) 115–132, https://doi.org/10.1016/0378-3812(89)80172-4.

[69] C.R. Novenario, J.M. Caruthers, K.-C. Chao, A mixing rule to incorporate solution model into equation of state, Ind. Eng. Chem. Res. 35 (1996) 269–277, https://doi.org/10.1021/ie950203+.

[70] J. Ahlers, J. Gnehling, Development of a universal group contribution equation of state. 2. Prediction of vapor–liquid equilibria for asymmetric systems, Ind. Eng. Chem. Res. 41 (2002) 3489–3498, https://doi.org/10.1021/ie020047o.

[71] J. Escandell, E. Neau, C. Nicolas, A new formulation of the predictive NRTL-PR model in terms of kij mixing rules. Extension of the group contributions for the modeling of hydrocarbons in the presence of associating compounds, Fluid. Phase Equilib. 301 (2011) 80–97, https://doi.org/10.1016/j.fluid.2010.11.009.

[72] E. Neau, J. Escandell, C. Nicolas, Modeling of highly nonideal systems: 2. Prediction of high pressure phase equilibria with the group contribution NRTL-PR EoS, Ind. Eng. Chem. Res. 49 (2010) 7589–7596, https://doi.org/10.1021/ie101266x.

[73] E. Neau, J. Escandell, C. Nicolas, Modeling of highly nonideal systems: 1. A generalized version of the NRTL equation for the description of low-pressure equilibria, Ind. Eng. Chem. Res. 49 (2010) 7580–7588, https://doi.org/10.1021/ie100121c.

[74] E.C. Voutsas, G.D. Pappa, K. Magoulas, D.P. Tassios, Vapor liquid equilibrium modeling of alkane systems with equations of state: “simplicity versus complexity, Fluid. Phase Equilib. 240 (2006) 127–139, https://doi.org/10.1016/j.fluid.2005.12.010.

[75] 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, https://doi.org/10.1016/j.fluid.2016.09.003.