A General Model for Thermodynamic Properties of Fluid Mixtures Based on Helmholtz Energy Formulations for the Components. Virial Expansion and Reduction to van der Waals Mixing Rules

Abstract

Over the recent decades, Helmholtz energy formulations became available for a broad range of fluids. These multiparameter equations of state (R. Span, Springer 2000) allow computation of thermodynamic properties essentially within the experimental errorbars. Corresponding states-based model by Lemmon and Tillner-Roth (Fluid Phase Equilib 165:1, 1999) enabled construction of Helmholtz energy formulations for mixtures. However, we show that this model generates a non-physical dependence of virial coefficients on composition, which can be strong when the components are dissimilar. We propose a new mixture model that overcomes this deficiency. It has two main ingredients: (i) Quadratic mixing of “Helmholtz volumities”. This quantity with units of molar volume is introduced as a ratio of the molar residual Helmholtz energy to a product of gas constant, thermodynamic temperature, and molar density. It reduces to the second virial coefficient in the zero-density limit. Helmholtz volumities are considered for components and “cross-components”, hypothetical fluids representing the binary interactions. (ii) Replacing the variables—reduced reciprocal temperatures and reduced densities—with temperature and density scaling functions. Different scaling functions can be used for different components and cross-components, thus providing a highly flexible framework for representing the properties of mixtures. The scaling functions must be expandable into Taylor series in terms of molar concentrations in the zero-density limit. For the proposed mixture model, we develop formulas for computing virial coefficients up to the fourth order. Furthermore, we show that when the proposed mixture model is applied to a cubic equation of state, the conventional van der Waals mixing rules can be retrieved. These findings allow to consider the new model as a viable alternative to the corresponding states method of modeling thermodynamic properties of fluid mixtures.

1 Introduction

Multiparameter equations of state enable computation of thermodynamic properties with uncertainties comparable with the uncertainties of the underlying experimental data. They are implemented in thermophysical property databases such as REFPROP [1], TREND [2], and CoolProp [3]. Various aspects of multiparameter equations of state were analyzed by Span [4]. More recent developments were reviewed by Thol and Bell [5]. Tillner-Roth [6] and Lemmon [7] independently developed a corresponding states method for modeling mixtures of fluids which, as pure components, are described by multiparameter equations of state in the form of Helmholtz energy [8]. This concept was successfully used to model mixtures of air gases [9], natural gas mixtures [10, 11], and mixtures of refrigerants [12].

The present work is motivated by computations performed to support and analyze measurements of homogeneous nucleation of water droplets in various gases and gaseous mixtures [13,14,15,16]. The analysis required computing isentropic expansions and chemical potentials in supersaturated humid gases with a goal to determine the nucleation temperature and the composition of critical clusters [17]. A mixture model for the relevant systems was elaborated by Gernert and Span [18]. Although this model reproduces favorably experimental data including phase equilibria, the application to supersaturated humid gas lead to nonphysical predictions. An appropriate model was obtained by combining a virial equation of state with the ideal-gas parts of the multiparameter equations of state. Application of quantum mechanically computed cross second virial coefficients [19,20,21,22] enabled an accurate analysis of the homogeneous nucleation experiments which were restricted to pressures up to 1 MPa. The reason of the failure of mixture model [7] in this range of conditions is the fact that it does not properly reproduce the rigorous dependence of the virial coefficients on the composition of the mixture.

Second virial coefficient expresses the first correction to the ideal gas equation of state at medium densities, where forces between the molecules forming the gas cannot be neglected [23]. Virial coefficient $B_{11}$ represents interactions of two molecules of kind 1, cross second virial coefficient $B_{12}$ represents interactions of one molecule of kind 1 and one molecule of kind 2. These coefficients can be expressed as sums (in quantum mechanical computations) or integrals (in the semi-classical approximation) over the configuration space of the two molecules with no reference to the concentrations of the components. Consequently, coefficients $B_{ij}$ only depend on temperature. Second virial coefficient for a mixture of gases is rigorously given by quadratic form 14 in terms of molar fractions of the components, in which $B_{ij}$ form the matrix of coefficients. This fact can be traced to statistical mechanics; however, it can also be seen as a consequence of a phenomenologically justified assumption that, at given temperature, pressure is an analytical function of concentrations of the components. Rigorous mixing rules can also be found for the third and higher mixture virial coefficients, as we explain in Sect. 2.1.

The problem of composition dependence of virial coefficients as computed from the corresponding states-based mixture model [8] is known to the community active in the development of multiparameter equations of state. Jäger et al. [24] provided an extensive study of the composition dependence of the cross second virial coefficient $B_{12}$ computed from models of binary mixtures. They expressed $B_{12}$ from rigorous quadratic mixing rule 18 based on virial coefficients of the components $B_{11}$ and $B_{22}$, and on the virial coefficient B of the mixture, as given by the corresponding states model:

$$\begin{aligned} B_{12}=\frac{B(T,x_1)-x_1^2 B_{11}(T)-x_2^2 B_{22}(T)}{2x_1 x_2}. \end{aligned}$$

(1)

Here, $x_1$ and $x_2$ are molar fractions of components 1 and 2. The result of formula 1 should depend on temperature T only. This actually happens, e.g., when a cubic equation of state with standard van der Waals mixing rules is used. However, as shown by Jäger et al. [24], the corresponding states model [7] produces $B_{12}$ depending on composition.

We address this problem from a different perspective. (Cross) second virial coefficients $B_{ij}$, (cross) third virial coefficients $C_{ijk}$, etc., for a system of I components can be identified with zero-density limits of derivatives of the residual pressure with respect to molar concentrations $\varvec{\rho }=\{\rho _1,\dots ,\rho _I\}$, as given in Eq. 10, or, equivalently, to the zero-density limits of derivatives of the residual Helmholtz energy density, Eq. 27. A condition of existence of these limits is that they are invariant to the direction they are approached, i.e., invariant with respect to the mixture composition. It turns out that this condition is not satisfied for the corresponding states mixture model [7] not only for the cross terms $B_{ij}$, but also for diagonal (“pure component”) terms $B_{ii}$. (Later demonstrated in Figs. 12, 3 and 4). This testing approach can be extended to multi-component cases and to higher virial coefficients.

The origin of the problem is in the composition dependent reducing volume and reducing temperature, Eqs. 45 and 46. It concerns the fundamentals of the corresponding states approach. Jäger et al. [24], expressed their idea how the problem could be fixed: “...to have a noncorresponding states-based mixing rule in the limit of zero density and keeping the corresponding states character of the model at higher densities.” They also suggested that in the low density limit, the model should reduce to quadratic mixing of reduced Helmholtz energies $\alpha _{ij}$ of the pure components ($i=j$) and interaction terms ($i\ne j$)

$$\begin{aligned} \alpha ^\textrm{r}(T,\rho ,{\textbf{x}})=\sum _{i,j=1}^I x_i x_j \alpha ^\textrm{r}_{ij}(T,\rho ). \end{aligned}$$

(2)

Equation 2 is a quadratic form in molar fractions. The currently used expression 42 can be re-arranged into form 2 by defining $\alpha ^\textrm{r}_{ij}$ as given in Eq. 47. Using the quadratic form 2 directly appears advantageous, because it more closely corresponds to the nature of intermolecular interactions. In addition, the interaction function $\alpha ^\textrm{r}_{ij}$ might be estimated as an average of the residual reduced Helmholtz energies of components i and j, as suggested by Jäger et al. [24]. More importantly, reduced parameters $\tau$ and $\delta$ are abandoned in Eq. 2. This leads immediately to the correct mixing rule 14 for the second virial coefficient. Composition dependence of the third and higher virial coefficient is also correct, although it is restricted: e.g., coefficients $C_{112}$ and $C_{221}$ cannot be adjusted independently. The quadratic mixing can be considered as an acceptable approximation for higher virial coefficients in some cases [25]. However, form 2 needs to be modified in some (sophisticated) manner in order to capture the behavior of mixtures at high densities.

Cubic equations of state with van der Waals mixing rules are a classical approach to mixture modeling. They allow to model complex phase behavior with a few adjustable parameters. We show that the cubic equations provide the correct composition dependence of the mixture second, third, and fourth, virial coefficients. Cubic equations of course cannot approach the level of accuracy of the multiparameter equations of state. Yet it appears desirable that the physically sound van der Waals mixing rules are included as a special case in a general mixture model applicable to multiparameter equations of state.

In this article, we propose a route to mixture modeling which, by design, provides a correct composition dependence of the virial coefficients. We provide formulas for second, third, and fourth, virial coefficients. When applied to a cubic equation of state and the model parameters are chosen appropriately, the new model reduces to the classical van der Waals mixing.

2 Preliminaries

In this section we introduce three standard equations of state applicable to mixtures: the virial equation of state, cubic equations of state (represented by Peng–Robinson equation), and the corresponding states model applied to multiparameter equations of state. We review the facts important for the development of the new model.

Virial and cubic equations of state yield pressure as function of temperature T, density $\rho$, and composition represented by a vector of molar fractions ${\textbf{x}}=\{x_1,\dots ,x_I\}$. Reduced residual Helmholtz energy $\alpha ^\textrm{r}$ can be computed using thermodynamic relation

$$\begin{aligned} \alpha ^\textrm{r}=\frac{a^\textrm{r}}{RT} =\frac{1}{RT}\int _0^\rho \frac{p^\textrm{r}(T,\rho ,{\textbf{x}})}{\rho ^2} \,{\mathrm d}\rho , \end{aligned}$$

(3)

where we introduced the residual pressure $p^\textrm{r}$

$$\begin{aligned} p^\textrm{r}(T,\rho ,{\textbf{x}})=p(T,\rho ,{\textbf{x}})- RT\rho . \end{aligned}$$

(4)

Integration 3 is performed at constant temperature and composition.

2.1 Virial Equation of State

The virial equation of state for a pure fluid at temperature T and molar density $\rho$ gives residual pressure $p^\textrm{r}$ as

$$\begin{aligned} p^\textrm{r}(T,\rho )=RT\left[ B(T)\rho ^2+C(T)\rho ^3+D(T)\rho ^4+\cdots \right] . \end{aligned}$$

(5)

Here, B, C, and D, are, respectively, the second, third, and fourth virial coefficients. The virial coefficients strongly depend on temperature. Relation 5 with infinite number of terms is usually called the virial expansion. Mathematically, it can be understood as a Taylor expansion of function $p^\textrm{r}(T,\rho )/RT$ for a constant T at point $\rho =0$. For a mixture of I components with molar concentrations $\varvec{\rho }=\{\rho _1,\dots ,\rho _I\}$, the virial equation is generalized as

$$\begin{aligned} p^\textrm{r}(T,\varvec{\rho }){} & {} =RT\biggl [\,\sum _{i,j=1}^I B_{ij}(T)\rho _i\rho _j +\sum _{i,j,k=1}^I C_{ijk}(T)\rho _i \rho _j \rho _k \nonumber \\{} & {} \quad +\sum _{i,j,k,\ell =1}^I D_{ijk\ell }(T)\rho _i \rho _j\rho _k\rho _\ell +\dots \biggr ]. \end{aligned}$$

(6)

In Eq. 6, $B_{ii}$, $C_{iii}$, and $D_{iiii}$ are virial coefficients of pure fluid i. $B_{ij}$ with $i\ne j$ is the cross-second virial coefficient. Similarly, $C_{ijk}$, and $D_{ijk\ell }$ are, respectively, the third, and fourth, cross virial coefficients. $B_{ij}$ represents interactions of two molecules of kind i and j, $C_{ijk}$ represents interactions of three molecules of kinds i, j, and k, etc. The (cross) virial coefficients are invariant with respect to any permutation (changing the order) of the indices and they depend on temperature only.

Equation 6 is a polynomial expression with a structure corresponding to a multivariate Taylor expansion of function $p^\textrm{r}(T,\varvec{\rho })$ at point $\varvec{\rho }={\textbf{0}}$ (all concentrations equal to zero) and at given T, with no constant and linear terms:

$$\begin{aligned} p^\textrm{r}(T,\varvec{\rho }){} & {} =\frac{1}{2}\sum _{i,j=1}^I p^\textrm{r}_{ij}(T,{\textbf{0}})\,\rho _i\rho _j +\frac{1}{6}\sum _{i,j,k=1}^I p^\textrm{r}_{ijk}(T,{\textbf{0}})\,\rho _i \rho _j \rho _k \nonumber \\{} & {} \quad +\frac{1}{24}\sum _{i,j,k,\ell =1}^I p^\textrm{r}_{ijk\ell }(T,{\textbf{0}})\,\rho _i \rho _j\rho _k\rho _\ell +\dots . \end{aligned}$$

(7)

Here we denoted partial derivatives of residual pressure with respect to concentrations as

$$\begin{aligned} p^\textrm{r}_{ij}(T,\varvec{\rho })=\dfrac{\partial ^2{p^\textrm{r}(T,\varvec{\rho })}}{\partial {\rho _i}\partial {\rho _j}},\quad p^\textrm{r}_{ijk}(T,\varvec{\rho })=\frac{\partial ^3 p^\textrm{r}(T,\varvec{\rho })}{\partial \rho _i \partial \rho _j \partial \rho _k}, \quad \dots , \end{aligned}$$

(8)

and their zero-density limits we denoted as

$$\begin{aligned} p^\textrm{r}_{ij}(T,{\textbf{0}})=\lim _{\varvec{\rho }\rightarrow {\textbf{0}}} p^\textrm{r}_{ij}(T,\varvec{\rho }),\quad p^\textrm{r}_{ijk}(T,{\textbf{0}})=\lim _{\varvec{\rho }\rightarrow {\textbf{0}}} p^\textrm{r}_{ijk}(T,\varvec{\rho }),\quad \dots . \end{aligned}$$

(9)

By comparing Eqs. 6 and 7, we can identify the virial coefficients with limits of derivatives 9,

$$\begin{aligned} B_{ij}(T)=\frac{p^\textrm{r}_{ij}(T,{\textbf{0}})}{2RT},\quad C_{ijk}(T)=\frac{p^\textrm{r}_{ijk}(T,{\textbf{0}})}{6RT},\quad D_{ijk\ell }(T)=\frac{p^\textrm{r}_{ijk\ell }(T,{\textbf{0}})}{24RT}. \end{aligned}$$

(10)

We note that a condition of the existence of the Taylor expansion 7 and, correspondingly, of the virial expansion 6, is the existence of limits 9. This condition is satisfied in physical systems, but it is not the case in the corresponding states mixture model as discussed later.

The concentrations are expressed in terms of molar fractions ${\textbf{x}}$ and molar density $\rho$ as

$$\begin{aligned} \rho _i=x_i \rho ,\quad \rho =\sum _{i=1}^I \rho _i. \end{aligned}$$

(11)

When expression 11 is substituted in Eq. 6, we obtain a virial equation in a form analogous to the pure fluid Eq. 5:

$$\begin{aligned} p^\textrm{r}(T,{\textbf{x}},\rho )=RT\left[ B(T,{\textbf{x}})\rho ^2+C(T,{\textbf{x}})\rho ^3 +D(T,{\textbf{x}})\rho ^4+\cdots \right] , \end{aligned}$$

(12)

At constant composition, virial Eq. 12 is analogous to a virial equation of a pure fluid. The mixture virial coefficients can be identified with zero-density limits of derivatives with respect to molar density,

$$\begin{aligned} B(T,{\textbf{x}})= & {} \frac{1}{2RT}\lim _{\rho \rightarrow 0} \dfrac{\partial ^2 {p^\textrm{r}(\rho ,T,{\textbf{x}})}}{\partial {\rho }^2}, \end{aligned}$$

(13a)

$$\begin{aligned} C(T,{\textbf{x}})= & {} \frac{1}{6RT}\lim _{\rho \rightarrow 0} \frac{\partial ^3 p^\textrm{r}(\rho ,T,{\textbf{x}})}{\partial \rho ^3}. \end{aligned}$$

(13b)

The mixture virial coefficients depend on temperature and composition. By substituting in Eq. 6 for concentrations from Eq. 11 and comparing with the virial expansion in form of Eq. 12, we obtain mixing rules for the virial coefficients,

$$\begin{aligned} B(T,{\textbf{x}})= & {} \sum _{i,j=1}^I B_{ij}(T)\; x_i x_j, \end{aligned}$$

(14)

$$\begin{aligned} C(T,{\textbf{x}})= & {} \sum _{i,j,k=1}^I C_{ijk}(T)\; x_i x_j x_k, \end{aligned}$$

(15)

$$\begin{aligned} D(T,{\textbf{x}})= & {} \sum _{i,j,k,\ell =1}^I D_{ijk\ell }(T)\; x_i x_j x_k x_\ell . \end{aligned}$$

(16)

These relations can be understood as a direct consequence of the existence of the multivariate Taylor expansion 6, or, in other words, a consequence of the analyticity of function $p^\textrm{r}(T,\varvec{\rho })$ at point $\varvec{\rho }={\textbf{0}}$.

As already mentioned, cross virial coefficients differing only in the order of subscripts are equal. For example, the (cross) second virial coefficients $B_{ij}$ form a symmetric matrix with pure-component second virial coefficients $B_{ii}$ forming its diagonal. Therefore, sums 14–16 contain many equal terms. We will find the number of occurrences of the equal virial terms. The order (second, third, fourth,...) of the virial coefficient we denote as N. If all the indices are different (the interacting molecules are all different kinds), the number of occurrences is equal to N! (factorial). If the cross-virial coefficient represents interactions of $N_1$ molecules of kind 1, $N_2$ molecules of kind 2, and so on up to $N_I$ molecules of kind I, so that $N=\sum _{i=1}^I N_i$, the number of occurrences of the respective N-th order cross-virial coefficient in the virial expansion 6 is given by the multinomial coefficient

$$\begin{aligned} \left( {\begin{array}{c}N\\ N_1,N_2,\dots ,N_I\end{array}}\right) =\frac{N!}{N_1!\times N_2!\times \ldots \times N_I!}. \end{aligned}$$

(17)

Let us consider a binary mixture ($I=2$). Mixing rule 14 for the second virial coefficient becomes with help of 17

$$\begin{aligned} B(T,{\textbf{x}})=x_1^2 B_{11}+2x_1 x_2 B_{12}+x_2^2 B_{22}. \end{aligned}$$

(18)

Here, $B_{11}$ and $B_{22}$ are virial coefficients of pure components 1 and 2, and $B_{12}$ is the cross second virial coefficient. Representative cross virial coefficients are selected such that their indices form a non-decreasing sequence. The temperature variable of $B_{ij}$ is not shown for brevity. The third virial coefficient 15 reduces for a binary mixture to

$$\begin{aligned} C(T,{\textbf{x}})=x_1^2 C_{111}+3x_1^2 x_2 C_{112}+3x_1 x_2^2 C_{122}+x_2^3 C_{222}, \end{aligned}$$

(19)

and for the fourth virial coefficient 16 of a binary mixture we obtain

$$\begin{aligned} D(T,{\textbf{x}})=D_{1111}x_1^4+4 D_{1112}x_1^3x_2+6 D_{1122}x_1^2 x_2^2+4D_{1222}x_1 x_2^3+D_{2222}x_2^4. \end{aligned}$$

(20)

Reduced residual Helmholtz energy $\alpha ^\textrm{r}$ for the virial equation of state can be obtained by integration 3 with the residual pressure expressed by Eq. 12:

$$\begin{aligned} \alpha ^\textrm{r}=B(T,{\textbf{x}})\rho +\frac{1}{2}C(T,{\textbf{x}})\rho ^2 +\frac{1}{3}D(T,{\textbf{x}})\rho ^3+\dots . \end{aligned}$$

(21)

This expression can be considered as a Taylor expansion with respect to density at point $\rho =0$ at constant temperature and composition. Correspondingly, the mixture virial coefficients can be identified with limits

$$\begin{aligned} B(T,{\textbf{x}})= & {} \lim _{\rho \rightarrow 0} \dfrac{\partial {\alpha ^\textrm{r}(\rho ,T,{\textbf{x}})}}{\partial {\rho }}, \end{aligned}$$

(22a)

$$\begin{aligned} C(T,{\textbf{x}})= & {} \lim _{\rho \rightarrow 0} \dfrac{\partial ^2 {\alpha ^\textrm{r}(\rho ,T,{\textbf{x}})}}{\partial {\rho }^2}. \end{aligned}$$

(22b)

For the present analysis it is advantageous to introduce reduced Helmholtz energy density

$$\begin{aligned} \alpha ^\textrm{d}=\frac{A}{VRT}=\frac{a\rho }{RT}=\alpha \rho . \end{aligned}$$

(23)

Ratio of the Helmholtz energy of the system A (in joules) to the system volume V (in cubic meters) is the density of Helmholtz energy (J·m⁻³). The unit of $\alpha ^\textrm{d}$ is m⁻³. We use function $\alpha ^\textrm{d}(\theta ,\varvec{\rho })$ as a thermodynamic potential. In Supplementary Information (SI), we review relations needed to compute thermodynamic quantities based on this function. Because the reciprocal temperature

$$\begin{aligned} \theta =\frac{1}{T}, \end{aligned}$$

(24)

acts as natural variable there, we will use $\theta$ rather than T in subsequent analysis. Reduced Helmholtz energy density $\alpha ^\textrm{d}$ can be separated into ideal and residual parts. The residual reduced Helmholtz energy density for the virial equation of state 21 is

$$\begin{aligned} \alpha ^\textrm{dr}=\alpha ^\textrm{r}\rho =B(\theta ,{\textbf{x}})\rho ^2+\frac{1}{2}C(\theta ,{\textbf{x}})\rho ^3 +\frac{1}{3}D(\theta ,{\textbf{x}})\rho ^4+\dots . \end{aligned}$$

(25)

In the superscript, “d” stands for “density”, and “r” stands for “residual”. When the virial coefficients are expressed with Eqs. 14–16 and the molar fractions are expressed in terms of concentrations, we obtain a virial expansion of the residual reduced Helmholtz energy density

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho }){} & {} =\sum _{i,j=1}^I B_{ij}(\theta )\rho _i\rho _j +\frac{1}{2}\sum _{i,j,k=1}^I C_{ijk}(\theta )\rho _i \rho _j \rho _k \nonumber \\{} & {} \quad +\frac{1}{3}\sum _{i,j,k,\ell =1}^I D_{ijk\ell }(\theta )\rho _i \rho _j\rho _k\rho _\ell +\dots . \end{aligned}$$

(26)

By comparing virial expansion 26 with the Taylor expansion of function $\alpha ^\textrm{dr}(\theta ,\varvec{\rho })$ at $\varvec{\rho }={\textbf{0}}$ and constant temperature, the (cross) virial coefficients can be identified with limits of derivatives as

$$\begin{aligned} B_{ij}(\theta )=\tfrac{1}{2}\alpha ^\textrm{dr}_{ij}(\theta ,{\textbf{0}}),\; C_{ijk}(\theta )=\tfrac{1}{3}\alpha ^\textrm{dr}_{ijk}(\theta ,{\textbf{0}}),\; D_{ijk\ell }(\theta )=\tfrac{1}{8}\alpha ^\textrm{dr}_{ijk\ell }(\theta ,{\textbf{0}}). \end{aligned}$$

(27)

The notation for derivatives of the residual reduced Helmholtz energy density with respect to molar concentrations is analogous to that used for derivatives of the residual pressure, Eqs. 8 and 9.

2.2 Cubic Equations of State

Cubic equations of state such as (Soave)–Redlich–Kwong [26, 27] provide a consistent (albeit less accurate) model of mixtures in gas and liquid phases. Bell and Jäger [28] provided thermodynamic equations for a class of cubic equations. For the sake of concreteness we chose the Peng–Robinson equation [29], giving for a pure component i pressure p as function of temperature T and molar volume $v=1/\rho$ as

$$\begin{aligned} p=\frac{RT}{v-{\bar{b}}} -\frac{{\bar{a}}}{v^2+2{\bar{b}}v -{\bar{b}}^2}. \end{aligned}$$

(28)

The notation with a “bar” over ${\bar{a}}$ and ${\bar{b}}$ is used to distinguish from other quantities introduced later. For a mixture, the two parameters ${\bar{a}}$ and ${\bar{b}}$ are given, respectively, by linear form 29 and quadratic form 30 in terms of molar fractions $x_i$ of the components as

$$\begin{aligned} {\bar{b}}({\textbf{x}})= & {} \sum _{i=1}^I {\bar{b}}_i x_i, \end{aligned}$$

(29)

$$\begin{aligned} {\bar{a}}(T,{\textbf{x}})= & {} \sum _{i,j=1}^I {\bar{a}}_{ij}(T)\, x_i x_j, \end{aligned}$$

(30)

where ${\bar{b}}_i$ and ${\bar{a}}_{ii}$ are parameters of the pure components and ${\bar{a}}_{ij}$, $i\ne j$, capture interactions of unlike molecules. Relations 29 and 30 are generally referred to as the van der Waals mixing rules. Parameter ${\bar{b}}_i$ is related to the size of molecule i and it is taken as a substance-specific constant. Parameter ${\bar{a}}_{ii}$ in the attraction term depends on temperature and on the substance-specific acentric factor by a function designed by original authors [29] or later modifications. Interaction parameters ${\bar{a}}_{ij}$ are estimated as a geometric average of the pure component parameters ${\bar{a}}_{ii}$ and ${\bar{a}}_{jj}$, possibly modified by an interaction coefficient. These details are not relevant to the present analysis; it is sufficient to consider ${\bar{a}}_{ij}$ for like and dislike interactions as general functions of temperature as indicated in Eq. 30. To shorten formulas and to make the subsequent derivations easier, it is suitable to introduce a modified attraction parameter

$$\begin{aligned} {\mathcal {A}}(T,{\textbf{x}})=\sum _{ij}{\mathcal {A}}_{ij}(T)\,x_i x_j,\quad {\mathcal {A}}_{ij}(T)=\frac{{\bar{a}}_{ij}}{RT}. \end{aligned}$$

(31)

We express the residual pressure defined by Eq. 4 for the Peng–Robinson equation 28 as function of density,

$$\begin{aligned} p^\textrm{r}=RT\left( \frac{{\bar{b}}\rho ^2}{1-{\bar{b}}\rho } -\frac{{\mathcal {A}}\rho ^2}{1+2{\bar{b}}\rho -{\bar{b}}^2\rho ^2}\right) . \end{aligned}$$

(32)

The residual reduced Helmholtz energy is obtained by integration 3,

$$\begin{aligned} \alpha ^\textrm{r}=-\ln (1-{\bar{b}}\rho ) +\frac{{\mathcal {A}}}{2\sqrt{2}\,{\bar{b}}} \ln \frac{1+\left( 1-\sqrt{2}\right) {\bar{b}}\rho }{1+\left( 1+\sqrt{2}\right) {\bar{b}}\rho }. \end{aligned}$$

(33)

The residual pressure $p^\textrm{r}$ or the residual reduced Helmholtz energy density $\alpha ^\textrm{dr}$ can be Taylor-expanded in terms of density $\rho$ at point $\rho =0$. By comparing the expansion of $p^\textrm{r}$ with Eq. 12, or by comparing expansion of $\alpha ^\textrm{dr}$ with Eq. 21, we obtain expressions for the second, third, and fourth virial coefficient

$$\begin{aligned} B={\bar{b}}-{\mathcal {A}},\quad C={\bar{b}}^2+2{\bar{b}}{\mathcal {A}},\quad D={\bar{b}}^3-5{\bar{b}}^2{\mathcal {A}}. \end{aligned}$$

(34)

To obtain the composition dependence of the mixture virial coefficients, we use Eqs. 29 and 31 in expressions 34. The resulting multiple sums need to be re-arranged as described in SI. As a result, we recover mixing rules 14, 15, and 16, with (cross) virial coefficients given as

$$\begin{aligned} B_{ij}= & {} \tfrac{1}{2}({\bar{b}}_i+{\bar{b}}_j)-{\mathcal {A}}_{ij}, \end{aligned}$$

(35)

$$\begin{aligned} C_{ijk}= & {} \tfrac{1}{3}\left( {\bar{b}}_i{\bar{b}}_j +{\bar{b}}_i {\bar{b}}_k+{\bar{b}}_j{\bar{b}}_k\right) +\tfrac{2}{3}\bigl ({\bar{b}}_i {\mathcal {A}}_{jk} +{\bar{b}}_j {\mathcal {A}}_{ik}+{\bar{b}}_k {\mathcal {A}}_{ij}\bigl ), \end{aligned}$$

(36)

$$\begin{aligned} D_{ijk}= & {} \tfrac{1}{4}\bigl ({\bar{b}}_i {\bar{b}}_j {\bar{b}}_k +{\bar{b}}_i {\bar{b}}_j {\bar{b}}_\ell +{\bar{b}}_i {\bar{b}}_k {\bar{b}}_\ell +{\bar{b}}_j {\bar{b}}_k {\bar{b}}_\ell \bigr )\nonumber \\{} & {} -\tfrac{5}{6} \bigl ({\bar{b}}_i{\bar{b}}_j {\mathcal {A}}_{k\ell } +{\bar{b}}_i{\bar{b}}_k {\mathcal {A}}_{j\ell } +{\bar{b}}_i{\bar{b}}_\ell {\mathcal {A}}_{jk} +{\bar{b}}_j{\bar{b}}_k {\mathcal {A}}_{i\ell } +{\bar{b}}_j{\bar{b}}_\ell {\mathcal {A}}_{ik} +{\bar{b}}_k{\bar{b}}_\ell {\mathcal {A}}_{ij}\bigr ). \end{aligned}$$

(37)

The obtained (cross) second virial coefficients $B_{ij}$, (cross) third virial coefficients $C_{ijk}$, and (cross) fourth virial coefficients $D_{ijk\ell }$ depend only on temperature, not on composition. In this way, we have verified that the Peng–Robinson equation, as a typical representative of the family of cubic equations, exhibits a correct limiting behavior.

2.3 Multiparameter Equations of State and the Corresponding States Mixture Model

Modern multiparameter equations of state of pure fluids are generally developed in form of reduced (or dimensionless) Helmholtz energy $\alpha$

$$\begin{aligned} \alpha (\tau ,\delta )=\frac{a}{RT} =\alpha ^\circ (\tau ,\delta )+\alpha ^\textrm{r}(\tau ,\delta ), \end{aligned}$$

(38)

where a is molar Helmholtz energy, R is the universal gas constant, and T is thermodynamic temperature. Variables of the reduced Helmholtz energy are the reduced reciprocal temperature

$$\begin{aligned} \tau =\frac{T_\textrm{r}}{T}=T_\textrm{r}\theta , \end{aligned}$$

(39)

and reduced density

$$\begin{aligned} \delta =v_\textrm{r}\rho , \end{aligned}$$

(40)

where $T_\textrm{r}$ is the reducing temperature (chosen as the critical temperature of the fluid), and $v_\textrm{r}$ is the reducing volume (chosen as the critical volume of the fluid). As indicated by Eq. 38, the reduced Helmholtz energy is divided in the ideal gas part $\alpha ^\circ$, and the residual part $\alpha ^\textrm{r}$. Modern multiparameter equations of state [30,31,32,33,34,35,36,37,38,39,40], express the residual reduced Helmholtz energy as a sum of ${K_\textrm{pol}}$ polynomial terms, ${K_\textrm{exp}}$ exponential terms, and ${K_\textrm{gbs}}$ Gaussian bell-shaped terms:

$$\begin{aligned}{} & {} \alpha ^\textrm{r}(\tau ,\delta ) =\sum _{k=1}^{K_\textrm{pol}}n_k \delta ^{d_k}\tau ^{t_k} +\sum _{k={K_\textrm{pol}}+1}^{{K_\textrm{pol}}+{K_\textrm{exp}}}n_k\delta ^{d_k}\tau ^{t_k} \exp \left( -l_k\delta ^{p_k}\right) \nonumber \\{} & {} \quad +\sum _{k={K_\textrm{pol}}+{K_\textrm{exp}}+1}^{{K_\textrm{pol}}+{K_\textrm{exp}}+{K_\textrm{gbs}}} n_k\delta ^{d_k}\tau ^{t_k} \exp \left[ -\eta _k\;(\delta -\varepsilon _k)^2 -\beta _k\;(\tau -\gamma _k)^2\right] . \end{aligned}$$

(41)

While the density exponents $d_k$ are generally natural numbers, temperature exponents $t_k$ are typically rational or optimized irrational numbers. The polynomial terms correspond to the virial equation of state. The exponential terms, originally introduced in the Benedict–Webb–Rubin equation of state [41], enable to mathematically decouple the vapor and liquid regions. The Gaussian terms [30] are used for the critical region. Some earlier formulations [42, 43] contain additional “non-analytical” terms attempting to mimic the behavior of fluids in close vicinity of the critical point.

Reduced Helmholtz energy of a mixture can be separated into an ideal gas part and a residual part. The ideal gas part is given by rigorous thermodynamics. Difficult is the residual part. Residual reduced Helmholtz energy $\alpha ^\textrm{r}$ of a mixture of I components is modeled [8, 11] as

$$\begin{aligned} \alpha ^\textrm{r}(T,\rho ,{\textbf{x}})=\sum _{i=1}^I x_i \alpha _i^\textrm{r}(\tau ,\delta ) +\sum _{i=1}^{I-1} \sum _{j=i+1}^I x_i x_j F_{ij}\Delta \alpha ^\textrm{r}_{ij}(\tau ,\delta ), \end{aligned}$$

(42)

where dimensionless density is defined as

$$\begin{aligned} \delta =\rho v_\textrm{r}({\textbf{x}}), \end{aligned}$$

(43)

and dimensionless reciprocal temperature is defined as

$$\begin{aligned} \tau =\frac{T_\textrm{r}({\textbf{x}})}{T}. \end{aligned}$$

(44)

Reducing volume $v_\textrm{r}({\textbf{x}})$ and reducing temperature $T_\textrm{r}({\textbf{x}})$ depend on the mixture composition. These functions are constructed such that in the limit $x_i\rightarrow 1$ they approach, respectively, the reducing molar volume $v_{\textrm{r}i}$ and the reducing temperature $T_{\textrm{r}i}$ of component i. Their contemporary functional forms [10, 11] are Lorentz-Berthelot formulas modified to capture asymmetric mixture behavior,

$$\begin{aligned} v_\textrm{r}({\textbf{x}})=\sum _{i,j=1}^I x_i x_j \beta _{\textrm{v}ij} \gamma _{\textrm{v}ij}\frac{x_i+x_j}{\beta _{\textrm{v}ij}^2 x_i+x_j} \cdot \frac{1}{8} \left( v_{\textrm{r}i}^{1/3}+v_{\textrm{r}j}^{1/3}\right) ^3, \end{aligned}$$

(45)

$$\begin{aligned} T_\textrm{r}({\textbf{x}})=\sum _{i,j=1}^I x_i x_j \beta _{\textrm{T}ij} \gamma _{\textrm{T}ij}\frac{x_i+x_j}{\beta _{\textrm{T}ij}^2 x_i+x_j} \cdot \left( T_{\textrm{r}i}\cdot T_{\textrm{r}j}\right) ^{1/2}. \end{aligned}$$

(46)

The adjustable parameters in Eqs. 45 and 46 satisfy symmetries $\beta _{\textrm{T}ji}=1/\beta _{\textrm{T}ij}$, $\beta _{\textrm{v}ji}=1/\beta _{\textrm{v}ij}$, $\gamma _{\textrm{T}ji}=\gamma _{\textrm{T}ij}$, and $\gamma _{\textrm{v}ji}=\gamma _{\textrm{v}ij}$. Further in Eq. 42, $\Delta \alpha ^\textrm{r}_{ij}(\tau ,\delta )$ is an interaction function whose structure is analogous to the functional form of the reduced Helmholtz energy of a pure component, Eq. 41, and $F_{ij}$ is an adjustable parameter. Function $\Delta \alpha ^\textrm{r}_{ij}$ is denoted here with the added “$\Delta$” in order to distinguish from $\alpha ^\textrm{r}_{ij}$ introduced in Eq. 2. Both quantities can be related. Expression 42 can be re-arranged into quadratic form 2 by defining

$$\begin{aligned} \alpha ^\textrm{r}_{ij}=\tfrac{1}{2}\bigl (\alpha ^\textrm{r}_i+\alpha ^\textrm{r}_j +F_{ij}\Delta \alpha ^\textrm{r}_{ij}\bigr )\, \end{aligned}$$

(47)

with

$$\begin{aligned} F_{ji}\Delta \alpha ^\textrm{r}_{ji}=F_{ij}\Delta \alpha ^\textrm{r}_{ij},\quad F_{ii}\Delta \alpha ^\textrm{r}_{ii}=0. \end{aligned}$$

(48)

2.3.1 Virial Expansion Applied to the Corresponding States Mixture Model

Relations 22 can be written in terms dimensionless variables 43 and 44 as

$$\begin{aligned} B(T,{\textbf{x}})=v_\textrm{r}({\textbf{x}})\lim _{\delta \rightarrow 0} \alpha ^\textrm{r}_{\delta }(\tau ,\delta ,{\textbf{x}}),\quad C(T,{\textbf{x}})=v_\textrm{r}^2({\textbf{x}})\lim _{\delta \rightarrow 0} \alpha ^\textrm{r}_{\delta \delta }(\tau ,\delta ,{\textbf{x}}). \end{aligned}$$

(49)

Here, subscript $\delta$ indicates derivative with respect to $\delta$ at constant $\tau$ and ${\textbf{x}}$. By application of relations 49 to the corresponding states formulation Eq. 42 we find

$$\begin{aligned} B(T,{\textbf{x}}){} & {} = v_\textrm{r}({\textbf{x}})\Biggl \{\sum _{i=1}^I \frac{x_i}{v_{\textrm{r}i}} \,B_i\!\left[ \frac{T\cdot T_{\textrm{r}i}}{T_\textrm{r}({\textbf{x}})}\right] \nonumber \\{} & {} \quad +\sum _{i=1}^{I-1} \sum _{j=i+1}^I x_i x_j F_{ij}\,\Delta \alpha ^\textrm{r}_{ij,\delta }\! \left[ \frac{T_\textrm{r}({\textbf{x}})}{T},0\right] \Biggr \}, \end{aligned}$$

(50)

where we used for the second virial coefficient of component i

$$\begin{aligned} B_i=v_{\textrm{r}i} \alpha ^\textrm{r}_{i,\delta }(\tau ,0) . \end{aligned}$$

(51)

In Eq. 50, second virial coefficients of the components $B_i$ are evaluated at temperatures given by the expression in square brackets, differing from the actual temperature T. As it can be seen from expression 50, the corresponding states contribution to the second virial coefficients generally does not obey the quadratic mixing rule 14. Only in some special cases the quadratic mixing is recovered. One such case is that (i) the reducing temperature $T_\textrm{r}$ is constant (composition independent), (ii) zero correction term $\Delta \alpha _{ij,\delta }^\textrm{r}(\tau ,0)=0$, and (iii) a linear form for the mixture reducing volume is assumed,

$$\begin{aligned} v_\textrm{r}({\textbf{x}})=\sum _{j=1}^I x_j v_{\textrm{r}j}. \end{aligned}$$

(52)

For this case, we obtain a special case of quadratic scaling

$$\begin{aligned} B(T,{\textbf{x}}) =\sum _{i,j=1}^I x_i x_j \frac{1}{2}\left[ \frac{v_{\textrm{r}j}}{v_{\textrm{r}i}}B_i(T) +\frac{v_{\textrm{r}i}}{v_{\textrm{r}j}} B_j(T)\right] . \end{aligned}$$

(53)

Another special case is when (i) the reducing temperature $T_\textrm{r}$ is constant (composition independent), and, consequently, $\Delta \alpha _{ij,\delta }^\textrm{r}(\tau ,0)$ does not depend on composition, (ii) the reducing volume $v_\textrm{r}$ is constant, and, consequently, $v_\textrm{r}=v_{\textrm{r}i}$. In this case, Eq. 50 reduces to

$$\begin{aligned} B(T,{\textbf{x}}) =\sum _{i,j=1}^I x_i x_j \tfrac{1}{2}\left[ B_i(T)+B_j(T)+F_{ij}\Delta \alpha ^\textrm{r}_{ij,\delta }(T)\right] , \end{aligned}$$

(54)

where we assumed conditions 48. These conditions are not met by any physical system, but for a mixture of similar gases the deviations from the rigorous quadratic mixing rule might be insignificant.

In the Introduction we described the method of testing the composition dependence of second virial coefficient of a binary mixture introduced by Jäger et al. [24], which is based on Eq. 1. In this work we introduce a test based on Eq. 10, which identifies the (cross) second virial coefficients $B_{ij}$ (and higher coefficients as well) to the zero-density limits of the derivatives $p^\textrm{r}_{ij}$ of the residual pressure with respect to molar concentrations $\rho _i$ and $\rho _j$ at constant temperature. We note that relations of the (cross) virial coefficients to derivatives of the residual reduced Helmholtz energy density, Eq. 27, could be used with the same results. Using Eqs. 9 and 10, the second (cross) virial coefficient should be expressible as a limit

$$\begin{aligned} B_{ij}(T)=\frac{1}{2RT}\lim _{\varvec{\rho }\rightarrow {\textbf{0}}} p^\textrm{r}_{ij}(T,\varvec{\rho }). \end{aligned}$$

(55)

Point $\varvec{\rho }={\textbf{0}}$ can be approached from various directions in the I-dimensional space. A condition of the existence of limit 55 is that the same value is reached for any direction of approach. The direction of approach can be labeled with the vector of molar fractions ${\textbf{x}}$, and the distance from point ${\textbf{0}}$ can be measured as molar density $\rho$. Thus, the studied limit can be written as

$$\begin{aligned} B_{ij}(T,{\textbf{x}})=\frac{1}{2RT}\lim _{\rho \rightarrow 0} p^\textrm{r}_{ij}(T,{\textbf{x}}\rho )=\frac{1}{2RT}\lim _{\rho \rightarrow 0} p^\textrm{r}_{ij}(T,\rho ,x_1), \end{aligned}$$

(56)

where the right-most expression is meant for the case of a binary mixture. The dependence of $B_{ij}(T,{\textbf{x}})$ on composition is nonphysical. Yet we need to consider it here to evaluate how strong is this nonphysical effect in the corresponding states mixture model.

In the present test, we computed the (cross) second virial coefficients $B_{ij}(T,{\textbf{x}})$ in two ways, as explained below and detailed in SI. Both methods gave identical results, within the numerical accuracy.

In the first method, we computed the residual pressure using REFPROP10 software [1] which implements multiparameter equations of state and mixture models for the tested binary systems. Second-order partial derivatives $p^\textrm{r}_{ij}$ were computed numerically using standard differentiation schemes as explained in SI. We computed these derivatives along lines of constant composition with molar density as a parameter. The computed values $p^\textrm{r}_{ij}$ were extrapolated towards zero density as exemplified in Fig. 1.

Fig. 1

Composition dependent second virial coefficient of nitrogen $B_{11}$ in the mixture with water vapor. (An artifact of the corresponding states mixture model [8, 18]). Left: lines $p^\textrm{r}_{11}(T,x_1,\rho )/2RT$ [18], $T=600$ K, circle: $B_{11}$ [33], cross: $B_{11}$ [44]. Right: isotherms $B_{11}(T,{\textbf{x}})$ [18] (solid lines) and $B_{11}(T)$ [33] (dashed lines)

In the second method, we computed $B_{ij}(T,{\textbf{x}})$ based on derivatives of the mixture virial coefficient $B(T,{\textbf{x}})$, as given by the corresponding states mixture model, with respect to molar fraction $x_1$; $B_x$ denotes the first derivative, and $B_{xx}$ denotes the second derivative. In SI, we derived formulas

$$\begin{aligned} B_{11}(T,x_1)= & {} B(T,x_1) +x_2 B_x(T,x_1) +\frac{x_2^2}{2}\,B_{xx}(T,x_1), \end{aligned}$$

(57)

$$\begin{aligned} B_{12}(T,x_1)= & {} B(T,x_1)+\frac{x_2-x_1}{2} B_x(T,x_1) -\frac{ x_1 x_2}{2}\,B_{xx}(T,x_1), \end{aligned}$$

(58)

$$\begin{aligned} B_{22}(T,x_1)= & {} B(T,x_1)-x_1 B_x(T,x_1) +\frac{x_1^2}{2} B_{xx}(T,x_1). \end{aligned}$$

(59)

Two binary systems were studied as certain extreme cases: nitrogen-oxygen (rather similar gases) and nitrogen-water (gas and vapor dissimilar in critical parameters and other aspects).

Corresponding states model for mixtures of nitrogen [33] (component 1) and water [43] (component 2) was developed by Gernert and Span [18]. In Fig. 1 (left), expression $p^\textrm{r}_{11}/(2RT)$ is shown as function of density for various compositions of the gaseous mixture at 600 K. The composition dependent (artifact) coefficient $B_{11}(T,{\textbf{x}})$ is obtained as the zero-density limit (intersection of the lines with the left axis). Also shown is the second virial coefficient of nitrogen as given by the used multiparameter equation of state [33] and as given by recent quantum mechanical computations [44]. It can be seen that lines for high nitrogen content approach the second virial coefficient of nitrogen, whereas other lines do not (although they should). The composition dependent (artifact) second virial coefficient of nitrogen $B_{11}(T,{\textbf{x}})$ is shown in Fig. 1 (right). The lines are flat in the region of high nitrogen fraction (low humidity), but they widely diverge in mixtures containing more water. Similarly, in Fig. 2 (right) we see that lines representing the composition dependent (artifact) second virial coefficient of water $B_{22}(T,{\textbf{x}})$ approach the correct (composition independent) values $B_{22}(T)$ when water vapor prevails in the mixture. Figure 2 (left) shows that the composition dependent (artifact) cross-second virial coefficient $B_{12}(T,{\textbf{x}})$ approaches the correct value $B_{12}(T)$ in the region of low humidity. This corresponds to the fact that the largest part the vapor-liquid equilibrium data is in the region of low humidity and this region is correctly represented by the mixture model.

Fig. 2

Solid lines: composition dependent cross second virial coefficient of nitrogen and water $B_{12}(T,{\textbf{x}})$ (left) and composition dependent second virial coefficient of water in the mixture with nitrogen $B_{22}(T,{\textbf{x}})$ (right). (Artifacts of mixture model [8, 18]). Dashed lines: cross second virial coefficient of nitrogen and water $B_{12}(T)$ [21] (left) and second virial coefficient for water $B_{22}(T)$ [43]

Nitrogen [33] (component 1) and oxygen [45] (component 2) have similar second virial coefficients, and also the cross second virial coefficient [46] is of similar magnitude. As shown in Fig. 3, the nonphysical composition dependence of second virial coefficients for nitrogen $B_{11}(T,{\textbf{x}})$ and oxygen $B_{22}(T,{\textbf{x}})$ as obtained from the mixture model by Lemmon et al. [9] is almost negligible for this system. The same can be said for the cross virial coefficient $B_{12}(T,{\textbf{x}})$ shown in Fig. 4.

Fig. 3

Solid lines: composition dependent second virial coefficient of nitrogen $B_{11}(T,{\textbf{x}})$ (left) and oxygen $B_{22}(T,{\textbf{x}})$ (right). (Artifacts of mixture model [8, 9]). Dashed lines: second virial coefficient of nitrogen $B_{11}(T)$ [33] (left) and second virial coefficient for oxygen $B_{22}(T)$ [45]

Fig. 4

Solid lines: composition dependent cross second virial coefficient of nitrogen and oxygen $B_{12}(T,{\textbf{x}})$ (Artifact of mixture model [8, 9]). Dashed lines: cross second virial coefficient of nitrogen and oxygen $B_{12}(T)$ [46]

Finally, we note that when the (by artifact) composition dependent (cross) second virial coefficients (56) are used in the quadratic mixing rule, a correct mixture second virial coefficient is obtained:

$$\begin{aligned} B(T,{\textbf{x}})=\sum _{ij} B_{ij}(T,{\textbf{x}})\,x_i x_j . \end{aligned}$$

(60)

For a binary system, this fact can be proved by substituting Eqs. 57–59 into Eq. 60. It can be seen as a “generalization” of the rigorous mixing rule 14. We conclude that for the corresponding-states based mixture models, the quadratic mixing rule is valid only locally (in the vicinity of given ${\textbf{x}}$), whereas in physical systems it captures the complete composition dependence of the second virial coefficient.

3 Suggested Mixture Model

The purpose of the suggested formulation is to provide a complete scheme for modeling thermodynamic properties of mixtures starting form Helmholtz energy formulations for its components. Such a scheme is provided by the corresponding states mixture model. The model needs to be sufficiently general to embrace various non-ideal mixtures. In constructing the model, we aim at accurately reproducing the rigorous mixing rules for virial coefficients. This is a weak point of the corresponding states model, as discussed in the previous section.

3.1 Quadratic Mixing of “Helmholtz Volumities”

We introduce a quantity $\alpha ^\textrm{vr}$ as

$$\begin{aligned} \alpha ^\textrm{vr}=\frac{a^\textrm{r}}{\rho R T} =\frac{\alpha ^\textrm{r}}{\rho } =\frac{\alpha ^\textrm{dr}}{\rho ^2}. \end{aligned}$$

(61)

Quantity $\alpha ^\textrm{vr}$ has units of molar volume. Consequently, it appears appropriate to call it Helmholtz volumity. In the superscript, “v” stands for “volumity” and “r” stands for “residual”. This quantity will be used exclusively for the residual part of the Helmholz energy. Therefore, we skip the adjective “residual” in its name.

Let us first consider a pure fluid. With help of a reducing volume $v_\textrm{r}$ it is possible to define dimensionless reduced Helmholz volumity

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}=\frac{\alpha ^\textrm{vr}}{v_\textrm{r}}. \end{aligned}$$

(62)

With respect to Eqs. 40 and 61, the reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}$ of a pure component can be related to the residual reduced Helmholtz energy $\alpha ^\textrm{r}$ as

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}(\tau ,\delta )=\frac{\alpha ^\textrm{r}(\tau ,\delta )}{\delta }. \end{aligned}$$

(63)

We note that function ${\hat{\alpha }}^\textrm{vr}(\tau ,\delta )$ can be obtained analytically by simply lowering by one the powers of $\delta$ in individual terms of the reduced Helmholtz energy $\alpha$ as given by thermodynamic property formulations; these powers are generally greater or equal to one [31,32,33,34,35,36,37,38,39,40, 43], There are exceptions—in [30] some of the “modified two-dimensional Gaussian terms” appear with power $\delta ^0$, so that special measures need to be taken when evaluating ${\hat{\alpha }}^\textrm{vr}$ in the zero-density limit.

From Eq. 25 we see that the Helmholtz volumity has virial expansion

$$\begin{aligned} \alpha ^\textrm{vr}= \frac{\alpha ^\textrm{dr}}{\rho ^2} =B+\frac{1}{2}C\rho +\frac{1}{3}D\rho ^2+\dots . \end{aligned}$$

(64)

We recall that the virial coefficients depend on temperature and, for mixtures, also on composition as given by Eqs. 14–16. In the low-density limit, Helmoltz volumity $\alpha ^\textrm{vr}$ approaches the second virial coefficient. Inspired by this finding, we construct a function $\alpha ^\textrm{vr}$ for a mixture based on a quadratic mixing rule analogous to mixing rule 14 for the second virial coefficient:

$$\begin{aligned} \alpha ^\textrm{vr}(\theta ,{\textbf{x}},\rho )\approx \sum _{i,j=1}^I x_i x_j \alpha ^\textrm{vr}_{ij}(\theta ,\rho ). \end{aligned}$$

(65)

By design, Eq. 25 is accurate in the zero-density limit. We can hope that it will allow, at least to some extent, to represent the mixture properties also at higher densities. This actually happens, although a substantial refinement of the model, introduced in Sect. 3.2, is needed to make the suggested formulation sufficiently general. Formulation 65 treats interactions between like and unlike components in a symmetric form. Functions $\alpha ^\textrm{vr}_{ii}$ are Helmholtz volumities for individual components. Functions $\alpha ^\textrm{vr}_{ij}$ for $i\ne j$ are to some extent analogous to the functions for components. It appears appropriate in this context to introduce a term cross-component to denote a hypothetical fluid formed by particles interacting by i-j interaction forces. At the end of this subsection we provide a short discussion of the cross-components. Similarly as the (cross) second virial coefficients, functions $\alpha ^\textrm{vr}_{ij}$ form a symmetric matrix:

$$\begin{aligned} \alpha ^\textrm{vr}_{ji}(\theta ,\rho )=\alpha ^\textrm{vr}_{ij}(\theta ,\rho ). \end{aligned}$$

(66)

Multiplication of Eq. 23 with $\rho ^2$ yields an expression in terms of the residual reduced Helmholtz energy density (23),

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho })\approx \sum _{i,j=1}^I \rho _i \rho _j \alpha ^\textrm{vr}_{ij}(\theta ,\rho ). \end{aligned}$$

(67)

The reduced Helmholtz energies $\alpha ^\textrm{r}_{ii}$ for pure components as given by thermodynamic property formulations are generally functions of reduced density

$$\begin{aligned} \delta _{ii}=v_{ii}\rho , \end{aligned}$$

(68)

and reduced reciprocal temperature

$$\begin{aligned} \tau _{ii}=\frac{T_{ii}}{T}=T_{ii} \theta , \end{aligned}$$

(69)

where the reducing volume $v_\textrm{r}=v_{ii}$ is chosen as the critical volume of pure fluid i, and the reducing temperature $T_\textrm{r}=T_{ii}$ is chosen as the critical temperature of pure fluid i. Analogously we define reduced reciprocal temperature and reduced density for the cross-components,

$$\begin{aligned} \tau _{ij}= & {} \frac{T_{ij}}{T}=T_{ij} \theta , \end{aligned}$$

(70)

$$\begin{aligned} \delta _{ij}= & {} v_{ij}\rho . \end{aligned}$$

(71)

The reducing parameters form symmetric matrices,

$$\begin{aligned} T_{ji}=T_{ij},\quad v_{ji}=v_{ij}. \end{aligned}$$

(72)

Reduced Helmholz volumities as defined by Eqs. 62 and 63 can be expressed for all components and cross-components:

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ij}(\tau _{ij},\delta _{ij}) =\frac{\alpha ^\textrm{vr}_{ij}}{v_{ij}} =\frac{\alpha ^\textrm{r}(\tau _{ij},\delta _{ij})}{\delta _{ij}}. \end{aligned}$$

(73)

The quadratic expression 67 then becomes

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho }) \approx \sum _{i,j=1}^I \rho _i \rho _j v_{ij}{\hat{\alpha }}^\textrm{vr}_{ij}[\tau _{ij}(\theta ),\delta _{ij}(\rho _{ij})]. \end{aligned}$$

(74)

As we demonstrate in Sect. 3.5, the Helmholtz volumity for a cross-component (i, j) can be assumed in a form corresponding to a hypothetical pure fluid whose molecules interact through a force field physically existing between molecules of types i and j. This allows its construction based on equations of state for pure components i and j (or other pure fluids) and justifies the term cross-component. However, at conditions where more complex interactions and critical phenomena play a role, the shape of function ${\hat{\alpha }}^\textrm{vr}_{ij}(\tau _{ij},\delta _{ij})$ will differ from a pure component case. One reason is that thermodynamic properties in the critical region of mixtures are strongly influenced by large composition fluctuations, which are absent in pure fluids. To some extent there is an analogy with functions $\Delta \alpha ^\textrm{r}_{ij}(\tau ,\delta )$ of the corresponding states model 42, for which a similar functional form is used as for the residual reduced Helmholtz energy of a pure fluid, Eq. 41, however including special exponential-Gaussian terms [11, 18].

3.2 General Density and Temperature Scaling

Function 74 provides a correct virial expansion and a considerable freedom in designing functions ${\hat{\alpha }}^\textrm{vr}_{ij}(\tau _{ij},\delta _{ij})$ for the cross-components. However, it is still not sufficiently general. In the gas phase, it does not allow to independently adjust the third and higher order cross virial coefficients. In the liquid phase, it is not capable of modeling mixtures of molecules with significantly different sizes (different partial volumes) and non-ideal mixtures. In fact, expression 74 corresponds to the formulation 42 for the case that the corresponding states procedure is not adopted.

Instead of the corresponding states approach, we refine the formulation in the following way. At place of the reduced reciprocal temperature $\tau _{ij}$ defined by Eqs. 69 and 70 we substitute a temperature scaling function ${\mathcal {T}}_{ij}(\theta ,\varvec{\rho })$, and at place of the reduced density $\delta _{ij}$ defined by Eqs. 68 and 71 we substitute a density scaling function ${\mathcal {D}}_{ij}(\theta ,\varvec{\rho })$. Both classes of scaling functions depend on reciprocal temperature $\theta$ and on the vector of concentrations $\varvec{\rho }$. In this way, the residual reduced Helmholtz energy density is expressed in a general form

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho })=\sum _{i,j=1}^I \rho _i \rho _j v_{ij}{\hat{\alpha }}^\textrm{vr}_{ij}[ {\mathcal {T}}_{ij}(\theta ,\varvec{\rho }), {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })]. \end{aligned}$$

(75)

Equation 75 is the central formulation of the proposed model. We note that the temperature scaling functions ${\mathcal {T}}_{ij}(\theta ,\varvec{\rho })$ and the density scaling functions ${\mathcal {D}}_{ij}(\theta ,\varvec{\rho })$ are dimensionless functions of dimensional arguments $\theta$ and $\varvec{\rho }$. These arguments can be reduced, e.g., using Eqs. 70 and 71. We do not proceed in this direction here in order to avoid unnecessary complexity of the resulting relations.

By introducing the generalized scaling functions, we obtain an extreme level of flexibility. Note that different scaling functions can be designed for different components and cross-components. However, the temperature and density scaling functions must satisfy certain conditions in order that the model is physically sound. We will investigate these conditions.

Because the scaling functions adjust the representation of interactions of two components which are symmetric, we expect also the scaling functions to be symmetric,

$$\begin{aligned} {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })= {\mathcal {T}}_{ji}(\theta ,\varvec{\rho }),\quad {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })= {\mathcal {D}}_{ji}(\theta ,\varvec{\rho }). \end{aligned}$$

(76)

When the composition reduces to a single component i, the temperature scaling function ${\mathcal {T}}_{ii}$ must be equal to the reduced inverse temperature,

$$\begin{aligned} {\mathcal {T}}_{ii}(\theta ,\{0,\dots ,\rho _i,\dots ,0\}) =\tau _{ii}=T_{ii}\theta , \end{aligned}$$

(77)

and also the density scaling function ${\mathcal {D}}_{ii}$ must become equal to the reduced density,

$$\begin{aligned} {\mathcal {D}}_{ii}(\theta ,\{0,\dots ,\rho _i,\dots ,0\}) =\delta _{ii}=v_{ii}\rho _i. \end{aligned}$$

(78)

Relations 77 and 78 can conveiently be called “boundary conditions” of the scaling functions. For scaling functions of the cross-components ($i\ne j$), direct analogy of these boundary conditions is not available.

We turn attention to the derivatives of the scaling functions. We denote them as follows:

$$\begin{aligned} {\mathcal {T}}_{ij,k}(\theta ,\varvec{\rho })= & {} \dfrac{\partial { {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })}}{\partial {\rho _k}} ,\quad {\mathcal {T}}_{ij,k\ell }(\theta ,\varvec{\rho }) =\dfrac{\partial ^2{ {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })}}{\partial {\rho _k}\partial {\rho _\ell }} , \end{aligned}$$

(79)

$$\begin{aligned} {\mathcal {D}}_{ij,k}(\theta ,\varvec{\rho })= & {} \dfrac{\partial { {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })}}{\partial {\rho _k}},\quad {\mathcal {D}}_{ij,k\ell }(\theta ,\varvec{\rho }) =\dfrac{\partial ^2{ {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })}}{\partial {\rho _k}\partial {\rho _\ell }}. \end{aligned}$$

(80)

As given by Eq. 76, scaling functions ${\mathcal {T}}_{ij}$ and ${\mathcal {D}}_{ij}$ are symmetric with respect to an interchange of indices. Therefore, also their derivatives are symmetric,

$$\begin{aligned} {\mathcal {T}}_{ji,k}(\theta ,\varvec{\rho })= & {} {\mathcal {T}}_{ij,k}(\theta ,\varvec{\rho }) ,\quad {\mathcal {T}}_{ji,k\ell }(\theta ,\varvec{\rho })= {\mathcal {T}}_{ij,k\ell }(\theta ,\varvec{\rho }) , \end{aligned}$$

(81)

$$\begin{aligned} {\mathcal {D}}_{ji,k}(\theta ,\varvec{\rho })= & {} {\mathcal {D}}_{ij,k}(\theta ,\varvec{\rho }) ,\quad {\mathcal {D}}_{ji,k\ell }(\theta ,\varvec{\rho })= {\mathcal {D}}_{ij,k\ell }(\theta ,\varvec{\rho }) . \end{aligned}$$

(82)

Furthermore, because the sequence of taking derivatives with respect to various variables is irrelevant, the second-order derivatives are symmetric with respect to indices denoting derivatives:

$$\begin{aligned} {\mathcal {T}}_{ij,\ell k}(\theta ,\varvec{\rho }) = {\mathcal {T}}_{ij,k\ell }(\theta ,\varvec{\rho }),\quad {\mathcal {D}}_{ij,\ell k}(\theta ,\varvec{\rho }) = {\mathcal {D}}_{ij,k\ell }(\theta ,\varvec{\rho }). \end{aligned}$$

(83)

Scaling functions need to be designed such that they enable capturing detailed features of the mixture thermodynamics at higher densities, but they must not affect the proper virial expansion. This can be achieved if the scaling functions can be Taylor-expanded in terms of molar concentrations as

$$\begin{aligned} {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })= & {} T_{ij} \theta +\sum _{k=1}^I {\mathcal {T}}_{ij,k}^<(\theta )\; \rho _k +\frac{1}{2}\sum _{k,\ell =1}^I {\mathcal {T}}_{ij,k\ell }^<(\theta )\; \rho _k \rho _\ell \dots , \end{aligned}$$

(84)

$$\begin{aligned} {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })= & {} \sum _{k=1}^I {\mathcal {D}}_{ij,k}^<(\theta )\;\rho _k +\frac{1}{2}\sum _{k,\ell =1}^I {\mathcal {D}}_{ij,k\ell }^<(\theta ) \; \rho _k \rho _\ell . \end{aligned}$$

(85)

In the zero-density limit, the temperature scaling function ${\mathcal {T}}_{ij}$ becomes equal to the reduced reciprocal temperature $\tau _{ij}=T_{ij}\theta$, while the density scaling function approaches zero as a linear function of $\varvec{\rho }$.

Coefficients of Taylor series 84 and 85 are zero-density limits of derivatives 79 and 80. We denoted these limits as follows:

$$\begin{aligned} {\mathcal {T}}_{ij,k}^<(\theta )= & {} \lim _{\varvec{\rho }\rightarrow {\textbf{0}}} {\mathcal {T}}_{ij,k}(\theta ,\varvec{\rho }) , \quad {\mathcal {T}}_{ij,k\ell }^<(\theta ) =\lim _{\varvec{\rho }\rightarrow {\textbf{0}}} {\mathcal {T}}_{ij,k\ell }(\theta ,\varvec{\rho }) , \end{aligned}$$

(86)

$$\begin{aligned} {\mathcal {D}}_{ij,k}^<(\theta )= & {} \lim _{\varvec{\rho }\rightarrow {\textbf{0}}} {\mathcal {D}}_{ij,k}(\theta ,\varvec{\rho }), \quad {\mathcal {D}}_{ij,k\ell }^<(\theta ) =\lim _{\varvec{\rho }\rightarrow {\textbf{0}}} {\mathcal {D}}_{ij,k\ell }(\theta ,\varvec{\rho }) . \end{aligned}$$

(87)

Note that the zero-density limits of the derivatives of the scaling functions may depend on (reciprocal) temperature only. They must not depend on composition. Zero density limits of second and higher order derivatives are defined analogously to the limits of the first derivatives.

Equations 81 and 82 are valid also in the zero-density limit:

$$\begin{aligned} {\mathcal {T}}_{ij,k}^<= {\mathcal {T}}_{ji,k}^<,\quad {\mathcal {T}}_{ij,k\ell }^<= {\mathcal {T}}_{ji,k\ell }^<,\quad {\mathcal {D}}_{ij,k}^<= {\mathcal {D}}_{ji,k}^<,\quad {\mathcal {D}}_{ij,k\ell }^<= {\mathcal {D}}_{ji,k\ell }^<. \end{aligned}$$

(88)

In addition, Eqs. 83 hold in the zero-density limit,

$$\begin{aligned} {\mathcal {T}}_{ij,k\ell }^<= {\mathcal {T}}_{ij,\ell k}^<,\quad {\mathcal {D}}_{ij,k\ell }^<= {\mathcal {D}}_{ij,\ell k}^<. \end{aligned}$$

(89)

Boundary condition 77 can be expressed in terms of series 84 as

$$\begin{aligned} T_{ii} \theta + {\mathcal {T}}_{ii,i}^<\; \rho _i +\tfrac{1}{2} {\mathcal {T}}_{ii,ii}^<\; \rho _i^2+\ldots =T_{ii}\theta . \end{aligned}$$

(90)

In order that the Eq. 90 is satisfied for any $\rho _i$, the coefficients of the linear, quadratic, and higher terms must vanish:

$$\begin{aligned} {\mathcal {T}}_{ii,i}^<=0,\; {\mathcal {T}}_{ii,ii}^<=0,\; \;\text {etc.} \end{aligned}$$

(91)

Similarly, boundary condition78 can be expressed in terms of series 85 as

$$\begin{aligned} {\mathcal {D}}_{ii,i}^<\;\rho _i +\tfrac{1}{2}{\mathcal {D}}_{ii,ii}^< \; \rho _i^2+\ldots =v_{ii}\rho _i. \end{aligned}$$

(92)

In order that this condition holds for any $\rho _i$, the coefficients must be

$$\begin{aligned}{} & {} {\mathcal {D}}_{ii,i}^<=v_{ii},\; \end{aligned}$$

(93)

$$\begin{aligned}{} & {} {\mathcal {D}}_{ii,ii}^<=0,\; \;\text {etc.} \end{aligned}$$

(94)

It can be convenient to set a condition analogous to Eq. 93 also for the cross-components. We suggest a straightforward generalization of Eq. 93:

$$\begin{aligned} \tfrac{1}{2} {\mathcal {D}}_{ij,i}^<+\tfrac{1}{2} {\mathcal {D}}_{ij,j}^<=v_{ij}. \end{aligned}$$

(95)

Also, Eq. 94 can be generalized for a cross-component scaling function:

$$\begin{aligned} {\mathcal {D}}_{ij,ij}^<=0 . \end{aligned}$$

(96)

Conditions 95 and 96 appear be suitable for reducing extra degrees of freedom of the scaling functions. However, their application is not necessary.

The temperature scaling is assumed to substantially enhance the flexibility of the model, in particular for nonideal mixtures showing local composition effects. For simpler systems it is possible that an isothermal model with density scaling only is sufficient. In the isothermal model, Eq. 75 simplifies to

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho })=\sum _{i,j=1}^I \rho _i \rho _j v_{ij}{\hat{\alpha }}^\textrm{vr}_{ij}[ \tau _{ij}(\theta ), {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })], \end{aligned}$$

(97)

where $\tau _{ij}(\theta )$ is just a transformation of the reciprocal temperature $\theta$ to the dimensionless reciprocal temperatures according to Eq. 70.

3.2.1 Examples of the Density and Temperature Scaling Functions

The basic cases of the density scaling function are linear and quadratic.

Linear density scaling function is a linear function of concentrations $\rho _k$,

$$\begin{aligned} {\mathcal {D}}_{ij}(\theta ,\varvec{\rho }) =\sum _{k=1}^I b_{ijk}(\theta )\; \rho _k . \end{aligned}$$

(98)

Coefficients $b_{ijk}$ are generally functions of temperature. In the low density limit, we can easily identify

$$\begin{aligned} {\mathcal {D}}_{ij,k}^<=b_{ijk} . \end{aligned}$$

(99)

Therefore, all relations developed for the derivatives of the density scaling function in the low-density limit apply to coefficients $b_{ijk}$.

Quadratic scaling includes linear and quadratic terms,

$$\begin{aligned} {\mathcal {D}}_{ij}(\theta ,\varvec{\rho }) =\sum _{k=1}^I b_{ijk}(\theta )\; \rho _k +\frac{1}{2}\sum _{k,\ell =1}^I b_{ijk\ell }(\theta )\; \rho _k \rho _\ell . \end{aligned}$$

(100)

By comparison with Taylor expansion 85 we see that the linear term corresponds to the linear part of the series as given by Eq. 99, and also the quadratic terms correspond,

$$\begin{aligned} {\mathcal {D}}_{ij,k\ell }^<=b_{ijk\ell }. \end{aligned}$$

(101)

Linear and quadratic temperature scaling functions can be developed. The quadratic temperature scaling function can be written as

$$\begin{aligned} {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })=T_{ij} \theta +\sum _{k=1}^I c_{ijk}(\theta )\; \rho _k +\frac{1}{2}\sum _{k,\ell =1}^I c_{ijk\ell }(\theta )\; \rho _k \rho _\ell . \end{aligned}$$

(102)

By comparison with Eq. 84 we see that coefficients $c_{ijk}$ and $c_{ijk\ell }$ can be identified with the zero-density limits of the derivatives of the temperature scaling function.

3.3 Virial Expansion for the Proposed formulation

We first develop an expansion of the reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}(\tau ,\delta )$ of a pure fluid. Using Eqs. 63 and 64 we obtain a virial expansion in dimensionless form

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}(\tau ,\delta )={\bar{B}}(\tau )+\frac{1}{2}{\bar{C}}(\tau ) \delta +\frac{1}{3}{\bar{D}}(\tau ) \delta ^2+\dots . \end{aligned}$$

(103)

Here we defined reduced second, third, and fourth, virial coefficients

$$\begin{aligned} {\bar{B}}=\frac{B}{v_\textrm{r}},\quad {\bar{C}}=\frac{C}{v_\textrm{r}^2},\quad {\bar{D}}=\frac{D}{v_\textrm{r}^3}. \end{aligned}$$

(104)

By comparing expression 103 with a Taylor expansion of function ${\hat{\alpha }}^\textrm{vr}(\tau ,\delta )$ at $\delta =0$, $\tau$ being fixed, we find for the reduced virial coefficients

$$\begin{aligned} {\bar{B}}(\tau )={\hat{\alpha }}^\textrm{vr}(\tau ,0),\; {\bar{C}}(\tau )=2{\hat{\alpha }}^\textrm{vr}_\delta (\tau ,0),\; {\bar{D}}(\tau )=\frac{3}{2}{\hat{\alpha }}^\textrm{vr}_{\delta \delta }(\tau ,0). \end{aligned}$$

(105)

We also consider derivatives of the reduced virial coefficients ${\bar{B}}$ and ${\bar{C}}$ with respect to $\tau$. We denote these derivatives with primes and we relate them to the derivatives of the reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}$ in the zero-density limit:

$$\begin{aligned} {\bar{B}}'(\tau ){} & {} =\dfrac{\,{\mathrm d}{{\bar{B}}}}{\,{\mathrm d}{\tau }} ={\hat{\alpha }}^\textrm{vr}_\tau (\tau ,0),\; {\bar{B}}''(\tau )=\dfrac{\,{\mathrm d}^2 {{\bar{B}}}}{\,{\mathrm d}{\tau }^2} ={\hat{\alpha }}^\textrm{vr}_{\tau \tau }(\tau ,0),\;\nonumber \\{} & {} \quad {\bar{C}}'(\tau )=\dfrac{\,{\mathrm d}{{\bar{C}}}}{\,{\mathrm d}{\tau }} =2{\hat{\alpha }}^\textrm{vr}_{\delta \tau }(\tau ,0). \end{aligned}$$

(106)

We can now develop a two-dimensional second-order Taylor expansion of function ${\hat{\alpha }}^\textrm{vr}$ at point $(\tau ,0)$:

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}(\tau +\Delta \tau ,\delta ){} & {} ={\hat{\alpha }}^\textrm{vr}(\tau ,0)+{\hat{\alpha }}^\textrm{vr}_\delta (\tau ,0)\,\delta +{\hat{\alpha }}^\textrm{vr}_\tau (\tau ,0)\,\Delta \tau \nonumber \\{} & {} \quad +\frac{1}{2}{\hat{\alpha }}^\textrm{vr}_{\delta \delta }(\tau ,0)\,\delta ^2 +{\hat{\alpha }}^\textrm{vr}_{\delta \tau }(\tau ,0)\,\delta \,\Delta \tau +\frac{1}{2}{\hat{\alpha }}^\textrm{vr}_{\tau \tau }(\tau ,0)\,(\Delta \tau )^2. \end{aligned}$$

(107)

Using Eqs. 105 and 106 we replace the derivatives of function ${\hat{\alpha }}^\textrm{vr}(\tau ,\delta )$ with reduced virial coefficients and their temperature derivatives,

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}(\tau +\Delta \tau ,\delta ){} & {} ={\bar{B}}(\tau )+\frac{1}{2}{\bar{C}}(\tau )\,\delta +{\bar{B}}'(\tau )\,\Delta \tau \nonumber \\{} & {} \quad +\frac{1}{3}{\bar{D}}(\tau )\,\delta ^2 +\frac{1}{2}{\bar{C}}'(\tau )\,\delta \,\Delta \tau +\frac{1}{2}{\bar{B}}''(\tau )\,(\Delta \tau )^2. \end{aligned}$$

(108)

We apply expansion 108 to all components and cross-components,

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ij}{} & {} ={\bar{B}}_{ij}+\frac{1}{2}{\bar{C}}_{ij}\,\delta _{ij} +{\bar{B}}_{ij}'\,\Delta \tau _{ij}\nonumber \\{} & {} \quad +\frac{1}{3}{\bar{D}}_{ij}\,\delta _{ij}^2 +\frac{1}{2}{\bar{C}}_{ij}'\,\delta _{ij}\,\Delta \tau _{ij} +\frac{1}{2}{\bar{B}}_{ij}''\,(\Delta \tau _{ij})^2. \end{aligned}$$

(109)

Because of the symmetry of the dimensionless Helmholtz energy with respect to the inversion of indices, also the virial coefficients for the cross-components are symmetric,

$$\begin{aligned} {\bar{B}}_{ji}={\bar{B}}_{ij},\quad {\bar{C}}_{ji}={\bar{C}}_{ij},\quad {\bar{D}}_{ji}={\bar{D}}_{ij}. \end{aligned}$$

(110)

The symmetry of course holds also for temperature derivatives 106,

$$\begin{aligned} {\bar{B}}_{ji}'={\bar{B}}_{ij}',\quad {\bar{B}}_{ji}''={\bar{B}}_{ij}'',\quad {\bar{C}}_{ji}'={\bar{C}}_{ij}'. \end{aligned}$$

(111)

Now we perform the scaling step—at place of $\delta _{ij}$ we substitute in 109 the density scaling function ${\mathcal {D}}_{ij}(\theta ,\varvec{\rho })$ in form of Taylor expansion 85. For $\tau _{ij}$ we substitute the temperature scaling function ${\mathcal {T}}_{ij}(\theta ,\varvec{\rho })$ in form of Taylor expansion 84. For $\Delta \tau _{ij}$ we substitute

$$\begin{aligned} {\mathcal {T}}_{ij}(\theta ,\varvec{\rho })-\tau _{ij} ={\mathcal {T}}_{ij}(\theta ,\varvec{\rho })-T_{ij} \theta =\sum _{k=1}^I {\mathcal {T}}_{ij,k}^<(\theta )\; \rho _k +\frac{1}{2}\sum _{k,\ell =1}^I {\mathcal {T}}_{ij,k\ell }^<(\theta )\; \rho _k \rho _\ell \dots . \end{aligned}$$

(112)

Performing the substitutions and neglecting higher than quadratic terms we obtain expansion for reduced Helmholtz volumity of (cross)component ij including the temperature and density scaling,

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ij}{} & {} ={\bar{B}}_{ij}+\sum _{k=1}^I \left( \frac{1}{2}{\bar{C}}_{ij}{\mathcal {D}}_{ij,k}^< +{\bar{B}}_{ij}' {\mathcal {T}}_{ij,k}^<\right) \, \rho _k +\sum _{k,\ell =1}^I \biggl (\frac{1}{4}{\bar{C}}_{ij} {\mathcal {D}}_{ij,k\ell }^< +\frac{1}{2}{\bar{B}}_{ij}' {\mathcal {T}}_{ij,k\ell }^<\nonumber \\{} & {} \quad +\frac{1}{3}{\bar{D}}_{ij}{\mathcal {D}}_{ij,k}^<{\mathcal {D}}_{ij,\ell }^< +\frac{1}{2}{\bar{C}}_{ij}' {\mathcal {D}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^< +\frac{1}{2}{\bar{B}}_{ij}'' {\mathcal {T}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^<\biggr )\, \rho _k\rho _\ell \dots . \end{aligned}$$

(113)

In order to simplify consequent derivations, we symmetrize the terms in the double summation. It is possible to write

$$\begin{aligned} \sum _{k,\ell =1}^I {\mathcal {D}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^< \,\rho _k\rho _\ell{} & {} =\sum _{k,\ell =1}^I {\mathcal {D}}_{ij,\ell }^<{\mathcal {T}}_{ij,k}^< \,\rho _k\rho _\ell \nonumber \\{} & {} \quad =\sum _{k,\ell =1}^I \left( \frac{1}{2}{\mathcal {D}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^< +\frac{1}{2}{\mathcal {D}}_{ij,\ell }^<{\mathcal {T}}_{ij,k}^<\right) \,\rho _k\rho _\ell . \end{aligned}$$

(114)

Using this result, we re-arrange expression 113 as

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ij}{} & {} ={\bar{B}}_{ij}+\sum _{k=1}^I \left( \frac{1}{2}{\bar{C}}_{ij}{\mathcal {D}}_{ij,k}^< +{\bar{B}}_{ij}' {\mathcal {T}}_{ij,k}^<\right) \, \rho _k\nonumber \\{} & {} \quad +\sum _{k,\ell =1}^I \biggl [\frac{1}{4}{\bar{C}}_{ij} {\mathcal {D}}_{ij,k\ell }^< +\frac{1}{2}{\bar{B}}_{ij}' {\mathcal {T}}_{ij,k\ell }^< +\frac{1}{3}{\bar{D}}_{ij}{\mathcal {D}}_{ij,k}^<{\mathcal {D}}_{ij,\ell }^<\nonumber \\{} & {} \quad +\frac{1}{4}{\bar{C}}_{ij}' \left( {\mathcal {D}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^<+{\mathcal {D}}_{ij,\ell }^<{\mathcal {T}}_{ij,k}^<\right) +\frac{1}{2}{\bar{B}}_{ij}'' {\mathcal {T}}_{ij,k}^<{\mathcal {T}}_{ij,\ell }^<\biggr ]\, \rho _k\rho _\ell \dots . \end{aligned}$$

(115)

We substitute Eq. 115 into expression 97 for the residual reduced Helmholtz energy density. After a re-arrangement, we obtain expansion

$$\begin{aligned} \alpha ^\textrm{dr}(\theta ,\varvec{\rho }){} & {} =\sum _{i,j=1}^I {\hat{B}}_{ij}\rho _i \rho _j +\frac{1}{2}\sum _{i,j,k=1}^I {\hat{C}}_{ijk} \;\rho _i\rho _j\rho _k\nonumber \\{} & {} \quad +\frac{1}{3}\sum _{i,j,k,\ell =1}^I {\hat{D}}_{ijk\ell }\;\rho _i\rho _j\rho _k\rho _\ell +\dots , \end{aligned}$$

(116)

where we introduced coefficients

$$\begin{aligned} {\hat{B}}_{mn}= & {} v_{mn} {\bar{B}}_{mn}, \end{aligned}$$

(117)

$$\begin{aligned} {\hat{C}}_{mnp}= & {} v_{mn}\left( {\bar{C}}_{mn} {\mathcal {D}}_{mn,p}^< +2{\bar{B}}'_{mn} {\mathcal {T}}_{mn,p}^<\right) , \end{aligned}$$

(118)

$$\begin{aligned} {\hat{D}}_{mnpq}= & {} v_{mn}\biggl [{\bar{D}}_{mn} {\mathcal {D}}_{mn,p}^< {\mathcal {D}}_{mn,q}^< +\frac{3}{4}{\bar{C}}_{mn} {\mathcal {D}}_{mn,pq}^< +\frac{3}{4}{\bar{C}}_{mn}'\Bigl ({\mathcal {D}}_{mn,p}^< {\mathcal {T}}_{mn,q}^<\nonumber \\{} & {} +{\mathcal {T}}_{mn,p}^<{\mathcal {D}}_{mn,q}^<\Bigr ) +\frac{3}{2}{\bar{B}}_{mn}'' {\mathcal {T}}_{mn,p}^< {\mathcal {T}}_{mn,q}^< +\frac{3}{2}{\bar{B}}_{mn}' {\mathcal {T}}_{mn,pq}^<\biggr ]. \end{aligned}$$

(119)

Coefficients ${\hat{B}}_{mn}$, ${\hat{C}}_{mn}$, and ${\hat{D}}_{mn}$ depend only on temperature. Due to the symmetry of the scaling volume 72, symmetries of derivatives of the scaling functions 88 and 89, and symmetries of the virial coefficients of the cross-components 110, the coefficients ${\hat{B}}_{mn}$, ${\hat{C}}_{mn}$, and ${\hat{D}}_{mn}$, also show some symmetries:

$$\begin{aligned} {\hat{B}}_{ij}= & {} {\hat{B}}_{ji}, \end{aligned}$$

(120)

$$\begin{aligned} {\hat{C}}_{ijk}= & {} {\hat{C}}_{jik}, \end{aligned}$$

(121)

$$\begin{aligned} {\hat{D}}_{ijk\ell }= & {} {\hat{D}}_{ij\ell k} ={\hat{D}}_{jik\ell }={\hat{D}}_{ji\ell k}. \end{aligned}$$

(122)

However, coefficients ${\hat{C}}_{ijk}$ and ${\hat{D}}_{ijk\ell }$ are not completely symmetric as it is the case of the virial coefficients $C_{ijk}$ an $D_{ijk\ell }$.

We collect the terms in 116 into groups with the same combination of indices (and the same combination of exponents of concentrations). We further use symmetries 120, 121, and 122. This procedure results in a rather long expression (SI.2) given in the SI. By comparison of terms in Eq. SI.41 with virial expansion (SI.2), we can express virial and cross-virial coefficients in terms of the present model as follows.

Pure-component virial coefficients are straightforward,

$$\begin{aligned} B_{ii}= & {} {\hat{B}}_{ii}, \end{aligned}$$

(123a)

$$\begin{aligned} C_{iii}= & {} {\hat{C}}_{iii}, \end{aligned}$$

(123b)

$$\begin{aligned} D_{iiii}= & {} {\hat{D}}_{iiii}. \end{aligned}$$

(123c)

Next we consider coefficients related to two different molecular kinds i and j. We will not show variants of the formulas obtainable by re-naming the indices. Second cross-virial coefficient is simple:

$$\begin{aligned} B_{ij}={\hat{B}}_{ij}=v_{ij} {\bar{B}}_{ij}. \end{aligned}$$

(124)

We used Eq. 117 for the right-most expression. Third cross-virial coefficients for two kinds of molecules are

$$\begin{aligned} C_{iij}=\frac{1}{3}\left( {\hat{C}}_{iij}+2{\hat{C}}_{iji}\right) . \end{aligned}$$

(125)

Fourth order cross-virial coefficient describing interactions of three, two, or one, molecules of kind i and one, two, or three molecules of kind j can be expressed, respectively, as

$$\begin{aligned} D_{iiij} =\frac{1}{2}\left( {\hat{D}}_{iiij}+{\hat{D}}_{ijii}\right) ,\quad D_{iijj}=\frac{1}{6} \left( {\hat{D}}_{iijj}+4{\hat{D}}_{ijij}+{\hat{D}}_{jjii}\right) . \end{aligned}$$

(126)

We consider interactions involving three molecular kinds i, j, and k. By comparing terms with powers $\rho _i \rho _j \rho _k$, in Eqs. SI.2 and SI.41, we obtain expression

$$\begin{aligned} C_{ijk} =\frac{1}{3}\left( {\hat{C}}_{ijk}+{\hat{C}}_{ikj}+{\hat{C}}_{jki}\right) . \end{aligned}$$

(127)

Equation 125 for two different kinds, as well as Eq. 123b for three-body interactions of a single molecular kind, can be obtained as special cases of Eq. 127. By comparing terms with powers $\rho _i^2 \rho _j \rho _k$ in Eqs. SI.2 and SI.41, we obtain expression

$$\begin{aligned} D_{iijk}=\frac{1}{6} \left( {\hat{D}}_{iijk}+2{\hat{D}}_{ijik}+2{\hat{D}}_{ikij} +{\hat{D}}_{jkii}\right) . \end{aligned}$$

(128)

Finally, we consider interactions of four different molecules. By comparing terms with product $\rho _i \rho _j \rho _k \rho _\ell$ in Eqs. SI.2 and SI.41, we obtain expression

$$\begin{aligned} D_{ijk\ell } =\frac{1}{6}\biggl ( {\hat{D}}_{ijk\ell } +{\hat{D}}_{ikj\ell } +{\hat{D}}_{i\ell jk} +{\hat{D}}_{jki\ell } +{\hat{D}}_{j\ell ik} +{\hat{D}}_{k\ell ij} \biggr ). \end{aligned}$$

(129)

Relations for cross-virial coefficients involving four-body interactions of three kinds of molecules 128, two kinds of molecules 126), and a single molecular kind 123c, can be obtained as special cases of Eq. 129.

We expressed the second, third, and fourth, (cross) virial coefficients of a mixture $B_{ij}$, $C_{ijk}$, and $D_{ijk\ell }$, in terms of the parameters of the proposed model: the virial coefficients of the components and cross-components ${\bar{B}}_{ij}$, ${\bar{C}}_{ij}$, ${\bar{D}}_{ij}$, and the zero-density limits of first and second derivatives of the density scaling functions ${\mathcal {D}}_{ij,k}^<$ and ${\mathcal {D}}_{ij,k\ell }^<$, and the temperature scaling functions ${\mathcal {T}}_{ij,k}^<$ and ${\mathcal {T}}_{ij,k\ell }^<$. This was accomplished through Eqs. 117, 118, and 119, giving, respectively, coefficients ${\hat{B}}_{mn}$, ${\hat{C}}_{mnp}$, and ${\hat{D}}_{mnpq}$, in terms of ${\bar{B}}_{ij}$, ${\bar{C}}_{ij}$, ${\bar{D}}_{ij}$, ${\mathcal {D}}_{ij,k}^<$, ${\mathcal {D}}_{ij,k\ell }^<$, ${\mathcal {T}}_{ij,k}^<$, and ${\mathcal {T}}_{ij,k\ell }^<$, and through Eqs. 124, 127, and 129, giving, respectively, the (cross) virial coefficients $B_{ij}$, $C_{ijk}$, and $D_{ijk\ell }$, in terms of coefficients ${\hat{B}}_{mn}$, ${\hat{C}}_{mnp}$, and ${\hat{D}}_{mnpq}$.

3.4 Why NOT to Use a Quadratic form in Terms of Residual Reduced Helmholtz Energy

We suggested formulation 75 as a quadratic form in terms of reduced Helmholtz volumities ${\hat{\alpha }}^\textrm{vr}_{ij}$, which are related to the residual reduced Helmholtz energies $\alpha ^\textrm{r}_{ij}$ by Eq. 63 as ${\hat{\alpha }}^\textrm{vr}_{ij}(\tau _{ij},\delta _{ij}) =\alpha ^\textrm{r}_{ij}(\tau _{ij},\delta _{ij})/\delta _{ij}$. Here we would like to explain, why it is not wise to use a quadratic form in terms of residual reduced Helmholtz energies, as suggested by Jäger et al. [24] in Eq. 2. We note that this form can be considered as convenient rearrangement of the contemporarily used formulation 42, as we have shown by Eqs. 47 and 48. The important difference comes when we attempt to introduce the temperature and density scaling functions into Eq. 2:

$$\begin{aligned} \alpha ^\textrm{r}=\sum _{i,j=1}^I x_i x_j \alpha ^\textrm{r}_{ij}[ {\mathcal {T}}_{ij}(\theta ,\varvec{\rho }), {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })]. \end{aligned}$$

(130)

We multiply Eq. 130 by $\rho ^2$ and find that the left-hand side is equal to the product of molar density $\rho$ and the reduced density of residual Hemlholtz energy $\alpha ^\textrm{dr}$,

$$\begin{aligned} \rho ^2\alpha ^\textrm{r}=\rho \alpha ^\textrm{dr}=\sum _{i,j=1}^I \rho _i \rho _j \alpha ^\textrm{r}_{ij}[ {\mathcal {T}}_{ij}(\theta ,\varvec{\rho }), {\mathcal {D}}_{ij}(\theta ,\varvec{\rho })]. \end{aligned}$$

(131)

We expand the reduced residual Helmholtz energy analogously to Eq. 109, however keeping fewer terms:

$$\begin{aligned} \alpha _{ij}^\textrm{r}={\bar{B}}_{ij}\,\delta _{ij} +\frac{1}{2}{\bar{C}}_{ij}\,\delta _{ij}^2 +{\bar{B}}_{ij}'\,\delta _{ij}\Delta \tau _{ij}+\dots . \end{aligned}$$

(132)

We perform the scaling step analogously to Eq. 113, but keeping only the first term, which is sufficient for the argument:

$$\begin{aligned} \alpha _{ij}^\textrm{r}={\bar{B}}_{ij}\,\left[ \sum _{k=1}^I {\mathcal {D}}_{ij,k}^<(\theta )\;\rho _k \right] +\dots . \end{aligned}$$

(133)

We substitute expansion 133 into the right-hand side of Eq. 131, and in the left-hand side we apply virial expansion 26:

$$\begin{aligned} \sum _{k=1}^I\rho _k \cdot \sum _{i,j=1}^I B_{ij}\rho _i \rho _j+\dots =\sum _{i,j=1}^I \rho _i \rho _j {\bar{B}}_{ij}\,\left[ \sum _{k=1}^I {\mathcal {D}}_{ij,k}^<\;\rho _k\right] +\dots . \end{aligned}$$

(134)

The dots on both sides stay for fourth and higher powers of molar concentrations. Since equality 134 must hold for all values of concentrations $\rho _i,\rho _j,\rho _k$, and all the coefficients are independent of concentrations (only depending on temperature), we find that this can only happen if ${\mathcal {D}}_{ij,k}^<$ is the same number for given i and j (independent of k). The choice compatible with boundary conditions 93 and 95 is

$$\begin{aligned} {\mathcal {D}}_{ij,k}^<=v_{ij}, \end{aligned}$$

(135)

where $v_{ij}$ is the reducing volume. In this case, Eq. 134 simplifies to

$$\begin{aligned} \sum _{i,j=1}^I B_{ij}\rho _i \rho _j+\dots =\sum _{i,j=1}^I \rho _i \rho _j {\bar{B}}_{ij}v_{ij}+\dots . \end{aligned}$$

(136)

With respect to the definition of reduced virial coefficients Eq. 104 we find that Eq. 136 is correct. However, condition 135 makes it impossible to choose various values ${\mathcal {D}}_{ij,k}^<$ for various k, that is, various degrees of influence of different components k on the scaling of (cross) component (i, j). For limits of higher derivatives, it can be deduced in a similar manner that formulation 130 with boundary condition 94 enforce ${\mathcal {D}}_{ij,k\ell }=0$. The same result will be found for higher derivatives. Therefore, we conclude that formulation 130 does not allow introducing the scaling functions.

3.5 Application of the Suggested Formulation to a Cubic Equation of State

We will show that conventional van der Waals mixing rules 29 and 30 for parameters of cubic equations of state represent a special case of the suggested method for modeling mixtures based on Helmholtz energy equations of state.

Reduced Helmholz energy for a pure component i can be obtained as a special case of Eq. 33

$$\begin{aligned} \alpha ^\textrm{r}_{ii}=-\ln (1-{\bar{b}}_i \rho ) +\frac{{\mathcal {A}}_{ii}(T)}{2\sqrt{2}\,{\bar{b}}_i} \ln \frac{1+\left( 1-\sqrt{2}\right) {\bar{b}}_i\rho }{1+\left( 1+\sqrt{2}\right) {\bar{b}}_i\rho }. \end{aligned}$$

(137)

We convert density into a dimensionless variable

$$\begin{aligned} \delta _{ii}=v_{ii}\rho , \end{aligned}$$

(138)

where $v_{ii}$ is the critical volume of component i. The considered scaling will be isothermal; parameter ${\mathcal {A}}_{ii}(T)$ will be evaluated at physical temperature T and introduction of a reduced temperature is needless in this context. Reduced Helmholtz energy 137 can be expressed in terms of the reduced density 138 as

$$\begin{aligned} \alpha ^\textrm{r}_{ii} =-\ln (1-\beta \delta _{ii}) +\frac{{\mathcal {A}}_{ii}}{2\sqrt{2}\,v_{ii}\beta } \ln \frac{1+\left( 1-\sqrt{2}\right) \beta \delta _{ii}}{1+\left( 1+\sqrt{2}\right) \beta \delta _{ii}}. \end{aligned}$$

(139)

We note hat the critical volume $v_{ii}$ and critical temperature $T_{ii}$ are related to parameters ${\bar{a}}_{ii}$ (or ${\mathcal {A}}_{ii}$) and ${\bar{b}}_i$. In Eq. 139, a constant $\beta$ was introduced

$$\begin{aligned} \beta =\frac{{\bar{b}}_i}{v_{ii}} =0.253\dots . \end{aligned}$$

(140)

This numerical value of constant $\beta$ is specific for the Peng–Robinson equation and it is different for other cubic equation of state.

Equation 139 is an equation of state of a pure component i in the form of reduced Helmholtz energy depending on reduced density and reduced reciprocal temperature, as it is usual for the multiparameter equations of state. We can now start to model the mixture based on the present method. First, we express reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}_{ii}$ according to Eq. 63

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ii}=\frac{\alpha ^\textrm{r}_{ii}}{ \delta _{ii}} =-\frac{1}{ \delta _{ii}}\ln (1-\beta \delta _{ii}) +\frac{{\mathcal {A}}_{ii} }{2\sqrt{2} \,v_{ii}\beta \delta _{ii}} \ln \frac{1+\left( 1-\sqrt{2}\right) \beta \delta _{ii}}{1+\left( 1+\sqrt{2}\right) \beta \delta _{ii}}. \end{aligned}$$

(141)

As the second step we suggest a suitable equation for the reduced Helmholtz volumity of the cross-components ${\hat{\alpha }}^\textrm{vr}_{ij}$, $i\ne j$. We chose a form analogous to the equation for a pure component 141:

$$\begin{aligned} {\hat{\alpha }}^\textrm{vr}_{ij}= -\frac{1}{ \delta _{ij}}\ln (1-\beta \delta _{ij}) +\frac{{\mathcal {A}}_{ij} }{2\sqrt{2} \,v_{ij}\beta \delta _{ij}} \ln \frac{1+\left( 1-\sqrt{2}\right) \beta \delta _{ij}}{1+\left( 1+\sqrt{2}\right) \beta \delta _{ij}}. \end{aligned}$$

(142)

Here, ${\mathcal {A}}_{ij}$ is the attractive term of the original Peng–Robinson equation modified according to Eq. 31, and $v_{ij}$ is the critical volume of the hypothetical cross-component fluid; Eq. 140 holds also for ratio ${\bar{b}}_{ij}/v_{ij}$. Reduced density for the cross-component is defined as

$$\begin{aligned} \delta _{ij}=v_{ij}\rho . \end{aligned}$$

(143)

As the third step, we replace reduced densities $\delta _{ij}$ in Eqs. 141 and 142 with density scaling functions ${\mathcal {D}}_{ij}(\theta ,\varvec{\rho })$. The scaling functions are so far unspecified. We construct the modified residual Helmholtz energy density $\alpha ^\textrm{dr}(\theta , \varvec{\rho })$ as given by Eq. 97:

$$\begin{aligned}{} & {} \alpha ^\textrm{dr}(\theta ,\varvec{\rho }) =-\sum _{i,j=1}^I \rho _i \rho _j \frac{v_{ij}}{ {\mathcal {D}}_{ij}}\ln (1-\beta {\mathcal {D}}_{ij})\nonumber \\{} & {} \quad +\sum _{i,j=1}^I \rho _i \rho _j \frac{{\mathcal {A}}_{ij}}{2\sqrt{2}\, \beta {\mathcal {D}}_{ij}} \ln \frac{1+\left( 1-\sqrt{2}\right) \beta {\mathcal {D}}_{ij}}{1+\left( 1+\sqrt{2}\right) \beta {\mathcal {D}}_{ij}}. \end{aligned}$$

(144)

As the final step, we have to choose the density scaling functions. We chose linear scaling 98 ,

$$\begin{aligned} {\mathcal {D}}_{ij} =\sum _{k=1}^I b_{ijk}\,\rho _k . \end{aligned}$$

(145)

From condition 93 it follows

$$\begin{aligned} {\mathcal {D}}_{ii,i}^<=b_{iii}=v_{ii}, \end{aligned}$$

(146)

and from condition 95 we have for $i\ne j$

$$\begin{aligned} \tfrac{1}{2} {\mathcal {D}}_{ij,i}^<+\tfrac{1}{2} {\mathcal {D}}_{ij,j}^< =\tfrac{1}{2}b_{iji}+\tfrac{1}{2}b_{ijj}=v_{ij}. \end{aligned}$$

(147)

These conditions still leave much freedom in shaping the scaling functions. In order to end up with the standard mixing rules for the cubic equations of state, we choose

$$\begin{aligned} b_{ijk}=v_{kk}. \end{aligned}$$

(148)

With this choice, ${\mathcal {D}}_{ij}$ as given by Eq. 145 is the same for all components and cross-components:

$$\begin{aligned} {\mathcal {D}}_{ij}=\sum _{k=1}^I v_{kk} \rho _k={\mathcal {D}}. \end{aligned}$$

(149)

We removed the indices i and j to indicate that the value ${\mathcal {D}}$ is, for a given $\varvec{\rho }$, the same for all components and cross-components. From Eq. 147 it follows that choice 148 requires a special value of the scaling volume of a cross-component as the arithmetic average of critical volumes of the respective pure components,

$$\begin{aligned} v_{ij}=\tfrac{1}{2}v_{ii}+\tfrac{1}{2}v_{jj}. \end{aligned}$$

(150)

For a pure fluid i at critical density equal to $1/v_{ii}$, the density scaling function ${\mathcal {D}}$ given by Eq. 149 equals to unity. For an equimolar binary mixture $x_1=x_2=1/2$ at density equal to the hypothetical critical density $1/v_{12}$ of the cross-component, the density scaling function 149 is also equal to unity.

Using Eqs. 149 and 150 in expression 144 we obtain

$$\begin{aligned} \alpha ^\textrm{dr}{} & {} =-\sum _{i,j=1}^I \rho _i \rho _j \frac{v_{ij}}{{\mathcal {D}}} \ln (1-\beta {\mathcal {D}}) +\sum _{i,j=1}^I \rho _i \rho _j \frac{{\mathcal {A}}_{ij}}{2\sqrt{2}\,\beta \delta } \ln \frac{1+\left( 1-\sqrt{2}\right) \beta {\mathcal {D}}}{1+\left( 1+\sqrt{2}\right) \beta {\mathcal {D}}}\nonumber \\{} & {} \quad = -\rho \ln (1-\beta {\mathcal {D}}) + \frac{1}{2\sqrt{2}\,\beta {\mathcal {D}}} \ln \frac{1+\left( 1-\sqrt{2}\right) \beta {\mathcal {D}}}{1+\left( 1+\sqrt{2}\right) \beta {\mathcal {D}}} \sum _{i,j=1}^I \rho _i \rho _j {\mathcal {A}}_{ij}. \end{aligned}$$

(151)

The last result can be considered as a final expression for the residual reduced Helmholtz energy density, which can be used to compute all thermodynamic properties of the mixture. We will now transform it to the form of reduced Helmholtz energy in order to compare it with expression 33 obtained using the conventional mixing rules. Using Eq. 61 we obtain reduced residual Helmholtz energy $\alpha ^\textrm{r}$ of the mixture

$$\begin{aligned} \alpha ^\textrm{r}{} & {} =\frac{\alpha ^\textrm{dr}}{\rho }\nonumber \\{} & {} \quad = -\ln (1-\beta {\mathcal {D}}) + \frac{1}{2\sqrt{2}\,\beta {\mathcal {D}}\rho } \ln \frac{1+\left( 1-\sqrt{2}\right) \beta {\mathcal {D}}}{1+\left( 1+\sqrt{2}\right) \beta {\mathcal {D}}} \sum _{i,j=1}^I \rho _i \rho _j {\mathcal {A}}_{ij}. \end{aligned}$$

(152)

With help of Eqs. 140, 149, and 11 we find

$$\begin{aligned} \beta {\mathcal {D}}=\sum _{k=1}^I \beta v_{kk} \rho _k =\sum _{k=1}^I {\bar{b}}_k \rho _k =\left( \sum _{k=1}^I {\bar{b}}_k x_k\right) \;\rho ={\bar{b}}\rho . \end{aligned}$$

(153)

In the right-most equality we recover the linear van der Waals mixing rule for parameter ${\bar{b}}$, Eq. 29. Further, using relation 11 we see that the modified attraction parameter ${\mathcal {A}}$ for the mixture is given by quadratic mixing rule 31, which is a straightforward modification of the van der Waals mixing rule 30:

$$\begin{aligned} \sum _{i,j=1}^I \rho _i \rho _j {\mathcal {A}}_{ij} =\rho ^2\sum _{i,j=1}^I {\mathcal {A}}_{ij} x_i x_j =\rho ^2{\mathcal {A}}. \end{aligned}$$

(154)

Using results 153 and 154 in Eq. 152 we recover expression 33 for the reduced residual Helmholtz energy of the mixture described by the Peng–Robinson equation of state.

We conclude that the present method of mixture modeling becomes identical with application of the conventional mixing rules 29 and 30 for the parameters of cubic equations of state when the following choices are adopted. (i) The cross-components are represented by the same cubic equation of state as the components. (ii) The reducing volume (hypothetical critical volume) of cross-component ij is taken as an arithmetic average of critical volumes of the pure components i and j, Eq. 150. The hypothetical critical temperature $T_{ij}$ of the cross-component is unconstrained; it corresponds to the choice of parameter ${\bar{a}}_{ij}$. (iii) The dimensionless volumities ${\hat{\alpha }}^\textrm{vr}_{ij}$ of the pure components and cross-components are evaluated for the same value of reduced density ${\mathcal {D}}$ given by linear scaling function 149 and the same temperature (isothermal scaling). However, the present method is much more general—a plentitude of other, physically sound mixture models can be developed when the model is not constrained as described in points (i) to (iii).

4 Conclusions

The contemporary corresponding states approach to modeling thermodynamic properties of mixtures [8] based on the equations of state of the components in form of Helmholtz energy is generally very successful, but it fails to reproduce the rigorous concentration dependence rules for the mixture virial coefficients [24]. To analyze this aspect deeper, we computed the second virial coefficients for pairs of like molecules and the cross virial coefficients for pairs of unlike molecules based on their definition as second order partial derivatives of the residual pressure. We developed formulas for computing these (cross) second virial coefficients based on the mixture virial coefficient and its derivatives with respect to molar fraction. The (cross) second virial coefficients should be independent of composition of the mixture. This rigorous requirement cannot be satisfied by the corresponding states mixture model. Yet for mixtures of similar gases (given example of nitrogen-oxygen), the deviations from the correct constant value can be negligible. For mixtures of fluids widely differing in critical parameters and other properties, as exemplified by nitrogen-water system, the (cross) virial coefficients may be sufficiently accurate in a region satisfactorily covered by experimental data, but strongly deviate in other regions.

We suggest a mixture model which, by design, provides correct composition dependences of second and higher virial coefficients. The concept of the model can be summarized in the following points.

• (1) Residual thermodynamic properties of the mixture are modeled based on Helmholtz energy models of the components with arbitrary mathematical structure, with reduced density and reduced reciprocal temperature as independent variables. (The same framework as for the corresponding states model.)

• (2) Quadratic mixing is used as a basis. This leads to the introduction of “cross-components” which reflect the interactions of unlike molecules. The cross-component is represented by a Helmholtz energy model analogous to a model of a component. (See the last paragraph of Sect. 3.1 for a brief discussion of this analogy and its limitations.)

• (3) The quadratically mixed property is not the residual reduced Helmholtz energy $\alpha ^\textrm{r}$, but rather a ratio $\alpha ^\textrm{r}/\rho$. This ratio has a dimension of volume. Therefore, we call it “Helmholtz volumity”. We also introduce reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}=\alpha ^\textrm{r}/\delta$, which is dimensionless.

• (4) Functions for the reduced Helmholtz volumity ${\hat{\alpha }}^\textrm{vr}_{ij}(\tau _{ij},\delta _{ij})$ for the components $i=j$ and cross-components $i\ne j$ are considered; arguments of these functions are the reduced reciprocal temperature $\tau _{ij}=T_{ij}/T$ and the reduced density $\delta _{ij}=v_{ij}\rho$, where $T_{ij}$ and $v_{ij}$ are, respectively, the reducing temperature and reducing volume of the (cross) component.

• (5) Reduced reciprocal temperatures $\tau _{ij}$ are substituted for with temperature scaling functions ${\mathcal {T}}_{ij}(\theta ,\varvec{\rho })$. Reduced densities $\delta _{i,j}$ are substituted for with density scaling functions ${\mathcal {D}}_{ij}(\theta ,\varvec{\rho })$. Variables of the scaling functions are reciprocal temperature $\theta =1/T$ and vector of molar concentrations $\varvec{\rho }=\{\rho _1,\dots ,\rho _I\}$. The scaling functions must satisfy certain boundary conditions in order that pure component properties are reproduced by the mixture model. In addition, multivariate Taylor expansions with respect to concentrations $\varvec{\rho }$ at point $\varvec{\rho }={\textbf{0}}$ must exist, such that the lowest term of any temperature scaling function expansion is a constant (only depends on $\theta$), and the lowest term of any density scaling function expansion is a linear form of $\varvec{\rho }$ with coefficients depending on $\theta$ only.

• (6) A variety of temperature and density scaling functions can be proposed. The basic ones are linear in $\varvec{\rho }$ and quadratic in $\varvec{\rho }$.

Possibly the linear model already provides enough flexibility for fitting the thermodynamic surfaces of simpler mixtures. We note that the coefficients $b_{ijk}$ and others are assumed as temperature functions. The model can be further simplified by keeping the temperature variable unscaled (${\mathcal {T}}_{ij}=\tau _{ij}$). Such a model can be called “isothermal”. We have shown that the van der Waals mixing rules for cubic equations of state can be considered as a special case of an isothermal model. On the other hand, for non-ideal mixtures which show local composition and local structuring effects, having the two variables (scaled temperature and scaled density) per binary interaction might be important. The two “tweaked” variables correspond to the two main parameters of the interaction - energy and entropy. Local ordering leads to lowering the residual entropy and this could be modeled by scaling down the temperature for a component and/or a cross-component.

Mathematically, the present approach is to some extent similar to the contemporary corresponding state method. It is also a “hybrid” model combining a quadratic form in concentrations or molar fractions with scaling of the variables of the functions forming the matrix of the quadratic form. The main differences can be summarized as follows:

• (a) The functions forming the matrix of coefficients are not the residual reduced Helmholtz energies but rather the Helmholtz volumities introduced in this article. The reason for this choice was explained in Sect. 3.4.

• (b) Density and temperature scalings are used for more subtle purposes in the present approach. The scalings fade away in the zero density limit. In the corresponding states approach, the scaling of temperature and density do the “brutal” job of transforming the properties of the mixture components from completely different physical conditions where there were measured for the pure fluids to the actual conditions of the mixture.

• (c) While for the corresponding states method the properties of all components and also the interaction terms (corresponding to the cross-components in language of this article) are evaluated at the same reduced temperature and density, in the present approach each component and cross-component has individual temperature and density scaling functions.

In this article we elaborated some basic properties of the proposed model. In particular we focused on the virial expansion, which is the weakness of the corresponding states method. We have shown how the virial coefficients can be obtained up to the fourth order. Further, we have shown that the van der Waals mixing rules can be considered as special case of the present model when applied to a cubic equation of state. These findings allow to consider the present model as a viable alternative to the corresponding states method of modeling thermodynamic properties of fluid mixtures. The model can be developed at several levels of detail: from a simple variant where the cross-component is modeled as a pure fluid and a few scaling parameters are adjusted to binary mixture data, over a complete fit of the cross-component functions to binary mixture data, up to highly detailed models enabling adjustment to ternary or quaternary mixture data.