Fluid Phase Equilibria 427 (2016) 513e538

Contents lists available at ScienceDirect

Fluid Phase Equilibria

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / fl u i d

A consistency test for a-functions of cubic equations of state

€l Jaubert*
Yohann Le Guennec, Silvia Lasala, Romain Privat, Jean-Noe
Universit 
e de Lorraine, Ecole Nationale Sup
erieure des Industries Chimiques, Laboratoire R
eactions et G
enie des Proc
ed
es (UMR CNRS 7274), 1rue
Grandville, 54000 Nancy, France

a r t i c l e i n f o a b s t r a c t

Article history: This study highlights that a-functions used in cubic equations of state must be of class C2 i.e. that their
Received 8 June 2016 first ðda=dTÞ and second ðd2 a=dT 2 Þ derivatives must exist and must be continuous, positive ða > 0Þ,
Received in revised form monotonically decreasing ðda=dT < 0Þ, convex ðd2 a=dT 2 > 0Þ and also verify d3 a=dT 3 < 0, for any value of
2 July 2016
the temperature T. Our proposed “consistency test for a-functions” gathers all these conditions. The non-
Accepted 27 July 2016
Available online 4 August 2016
respect of one of them can entail low-accuracy prediction of binary phase diagrams involving at least one
supercritical compound (this statement is illustrated through the case-studies of the CO2-argon and CO2-
decane systems) as well as improper variations of pure-component supercritical properties (h and cP )
Keywords:
Alpha-function
with respect to the temperature. Finally, an extensive study of the mostly used a-functions described in
Equation of state the open literature is performed and shows that all of them fail this test. Some component-dependent a-
Consistency test functions may however pass this test but only if mathematical constraints are added to their parameters.
Soave © 2016 Elsevier B.V. All rights reserved.
Twu

1. Introduction Since the Van der Waals proposal, even the most successful two-
parameter CEoS still expresses their attractive and repulsive terms
Numerous modifications to the Van der Waals cubic equation of by introducing a parameter a which is a measure of the attractive
state (CEoS) have been presented over the years [1], with the aim of forces (“energy parameter”) between molecules, and a parameter b
representing with increasingly accuracy the properties of complex which is a measure of the size (“intrinsic volume” or “co-volume”)
systems. As a matter of facts, most of these improvements do not of the molecules. It results the general formulation of a two-
incorporate a substantial difference in the functional combination parameter CEoS:
of attractive and repulsive forces, with respect to the original
formulation introduced by Van der Waals. As a rare exception, let us RT aðTÞ
mention the Cubic-Plus-Association model (CPA) [2], which com- PðT; vÞ ¼ þ (2)
v  b ðv  r1 bÞðv  r2 bÞ
bines the Soave-Redlich-Kwong attractive and repulsive terms [3]
with an association term introduced by Wertheim [4]. Indeed, where r1 and r2 are two universal constants which depend on the
analogously to the Van der Waals expression of the pressure selected CEoS.
equation [5], a generic CEoS may be written as the sum of a The most successful theoretically-based modification to the Van
repulsive and an attractive term: der Waals CEoS consisted in the recognition by Clausius in 1879 [7]
of the temperature-dependency of the energy parameter. In 1949,
P ¼ Prep þ Patt (1) Redlich and Kwong [8] proposed the first temperature-dependent
Although the analytical expressions of the repulsive and formulation which resulted to be sufficiently accurate to gain
attractive terms proposed in literature do not correctly quantify the popularity in the industry. Starting from this outcome, many re-
actual repulsive and attractive contributions to pressure, their sum searches have been focused on the improvement of the functional
results in a quantitative representation of fluid properties being form of this term, which soon proved to strongly affect the
sufficiently accurate to make their combination the “cornerstone of modelling capability of a CEoS in fluid equilibrium calculations. The
the generalized Van der Waals theory” [6]. attractive term of any CEoS can be written as:

aðTÞ ¼ ac ,aðTÞ (3)

* Corresponding author.
E-mail address: [email protected] (J.-N. Jaubert). i.e. as the product of the value of the attractive term at the critical

http://dx.doi.org/10.1016/j.fluid.2016.07.026
0378-3812/© 2016 Elsevier B.V. All rights reserved.
514 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

temperature ðac Þ multiplied by a so-called a-function. Such a-

functions have a strong impact on the accuracy of CEoS and must be aðTc Þ ¼ 1 (7)
selected with caution. Not only the properties of pure compounds Therefore, a-function should be considered as a deviation factor
in the supercritical region are highly affected by the mathematical of the energy parameter of a CEoS to its value at the critical
expression of the a-function but also the vapor-liquid equilibrium temperature.
calculations of multicomponent systems. From our experience,
when applied to mixtures, the accuracy of CEoS is equally affected
by the choice of the mixing rules and by the expression of the a- 2.2. Historical background and general classification of the a-
function. functions
The first aim of this paper is thus to present the requirements
that an a-function should absolutely fulfil in order to guarantee safe One of the most influential contributions in the prediction of
predictions of the vapor-liquid equilibrium and of derived ther- vapor pressures of non-polar and slightly polar pure compounds
modynamic properties at all temperatures. Emphasis will be made has been the introduction in 1972 of the Soave a-function [3]. Over
on the incapacity of a-functions regressed on subcritical properties the years, further modifications have been introduced to this
to correctly extrapolate in the supercritical temperature domain functional form, mainly aimed at improving the correlation of va-
without addition of supplementary constraints during the regres- por pressures. Higher accuracy was achieved by the addition of
sion procedure. Such constraints ensure, among others, the proper more parameters [10] [11], but as stated by Poling et al. [12], most of
curvature of the a-function. Similarly to the consistency tests the a-functions presented in literature have been developed
developed to certify the quality of experimental vapor-liquid without evaluating their predictive capability in the supercritical
equilibrium data, it is proposed to develop a consistency test for domain. Indeed, the Soave a-function diverges at very high tem-
a-functions, i.e., to derive a list of consistent constraints applicable perature, leading to unrealistic calculations of thermodynamic
to any a-function. Such a test aims at identifying which a-functions properties. Therefore, a-functions inspired by the Soave formula-
are (or not) thermodynamically consistent and should (or not) be tion, called polynomial a-functions as they are polynomials of the
used. In this paper, this consistency test will be applied to twelve a- square-root of the reduced temperature, should not be used at high
functions issued from the literature. reduced temperature. To overcome this limitation, some authors
proposed exponential a-functions [13e17] which are positive and
decreasing functions on the whole temperature range.
2. Definition, historical background and general classification Aware of the benefits of using a polynomial a-function in the
of the a-functions subcritical temperature domain and using an exponential a-function
in the supercritical temperature domain, some authors proposed
2.1. Definition of the a-function piecewise a-functions with different mathematical expressions for
the a-function depending on whether the temperature is above or
As for the original Van der Waals equation of state, the critical below the critical temperature, leading to possible discontinuities
attractive parameter (ac ), the covolume ðbÞ and the critical molar of the a-function derivatives at the critical temperature. Although
volume ðvc Þ of any CEoS can be determined by applying the so- this approach seems to be a good compromise between the
called critical constraints: complexity of the a-function expression and the accuracy of results,
it was objected by Coquelet et al. [18] that continuity of the first and
8  exp EoS 
>
> P Tc ; vc ¼ Pexp
c second derivatives of the a-function at the critical temperature
>
>   must be enforced. If not, calculated residual enthalpies and heat
>
> 
>
> vP 
< ¼0 capacities show discontinuities at the critical temperature. This
vv T T¼Texp ;v¼vEoS (4) abnormal behavior of the piecewise a-functions was also high-
c
>
c
>
> ! 
>
> 2  lighted by Boston et al. [19], Twu et al. [16,17] when they respec-
> v P 
> ¼0 tively proposed their own a-functions and deeply analyzed by Neau
>
: vv2 
exp
T T¼Tc ;v¼vc EoS et al. [20,21].
As reported by Valderrama [1], numerous a-functions can be
where Tcexp and Pcexp are the experimental critical temperature and found in the literature and are compared in some recent papers
pressure and where the critical molar volume ðvEoS
c Þ is calculated by [22,23]. Enumerating all the existing a-functions is not of practical
solving the CEoS under the conditions T ¼ Tcexp and P ¼ Pcexp . The interest as only a dozen of them are used in the commercial process
resulting critical energy parameter ðac Þ and covolume ðbÞ are usu- simulation software. Therefore, a general classification of the a-
ally expressed in terms of the CEoS-dependent coefficients Ua and functions based on their mathematical form (polynomial or expo-
Ub : nential) and on the unicity of their formulation (T-overall formu-
8 lation or piecewise) is presented in Table 1. Moreover, we found
 2
>
>
> R2 Tcexp desirable to distinguish two kinds of a-functions: component-
> ac ¼ Ua
< exp dependent a-functions and generalized ones. Parameters of
P c
(5) component-dependent a-functions have to be regressed, compo-
>
>
>
> RTcexp nent by component, on experimental data while coefficients of
: b ¼ Ub exp
Pc generalized a-functions are determined with correlations needing
the acentric factor ðuÞ as input parameter. Therefore, generalized a-
While a covolume is necessarily temperature-independent [9],
functions can be applied to pure compounds for which no experi-
the energy parameter aðTÞ is written as the product of its value at
mental vapor-liquid equilibrium data, except u, are known. This
the critical point, ac , and of a temperature-dependent a-function:
advantage is counterbalanced by a less accurate prediction of the
thermodynamic properties.
aðTÞ ¼ ac ,aðTÞ (6)
In the rest of the article, emphasis will be made on two a-
Doing so, it follows that aðTÞ is a non-dimensional factor which functions having a T-overall formulation: the model proposed by
becomes unity at the critical temperature Soave (1972) (from now on, called Soave) and the one introduced
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 515

Table 1
Classification of a-functions commonly found in process simulation software.

T-overall formulation (unique temperature formulation) Piecewise

Component-dependent a-functions
Polynomial Soave (1984) [24] Stryjek-Vera [11]
  pffiffiffiffiffi
a ¼ 1 þ ð1  Tr Þ a þ Tbr a ¼ ½1 þ kð1  pTrffiffiffiffiÞffi 2
k ¼ k0 þ k1 ð1 þ Tr Þð0:7  Tr Þ
k0 ¼ 0:378893 þ 1:4897153u  0:17131848u2 þ 0:0196554u3
k1 ðT  0:7Tc Þs0
k1 ðT > 0:7Tc Þ ¼ 0
Gibbons-Laughton [25] Mathias-Copeman [10]
pffiffiffiffiffi "
a ¼ 1 þ XðTr  1Þ þ Yð Tr  1Þ pffiffiffiffiffi pffiffiffiffiffi #2
 Tr Þ þ c2 ð1  Tr Þ2
aSUB ¼ aðTr  1Þ ¼ 1 þ c1 ð1 p ffiffiffiffiffi 3
þc3 ð1  Tr Þ
pffiffiffiffiffi
aSUP ¼ aðTr  1Þ ¼ ½1 þ mð1  Tr Þ2
Exponential Twu (1988) [14]
a ¼ Tr2ðM1Þ exp½Lð1  Tr2M Þ
Twu (1991) [15]
e Not any e
a ¼ Trd exp½Lð1  Trg Þ
Trebble-Bishnoi [26]
a ¼ exp½Lð1  Tr Þ
Polynomial & Exponential Melhem-Saini-Goodwin [27] Boston-Mathias [19]
pffiffiffiffiffi 8 pffiffiffiffiffi 2
a ¼ exp½mð1  Tr Þ þ nð1  Tr Þ2  <a ¼ aðT  1Þ ¼ a ¼ ½1 þ mð1  Tr Þ
SUB r SOAVE
:a ¼ aðTr > 1Þ ¼ exp½Lð1  Tr Þ g
SUP
8
>
> 2m
>L ¼
< 1þm
with
>
>
>
: m
g¼1þ
2
Mathias [28]
8 " pffiffiffiffiffi #2
>
> 1 þ mðuÞð1  Tr Þ
< aSUB ¼ aðTr  1Þ ¼
pð1  Tr Þð0:7  Tr Þ
>
>
:
aSUP ¼ aðTr > 1Þ ¼ ½exp½Lð1  Trg Þ2
8
> m
>
> L ¼ 1 þ þ 0:3p
< 2
with
>
> L  1
>
:g ¼
L
Generalized a-functions
Polynomial Soave (1972) [3]
8 pffiffiffiffiffi
< a ¼ ½1 þ mðuÞð1  Tr Þ2 e Not any e
m ðuÞ ¼ 0:480 þ 1:574u  0:176u2
: RK
mPR ðuÞ ¼ 0:37464 þ 1:54226u  0:26992u2
Exponential Trebble-Bishnoi [29] Twu (1995) [16,17]
8
a ¼ exp½Lð1  Tr Þ > ð0Þ ð1Þ
 að0Þ 
For alkanes; water and carbon dioxide : < a ¼ a þð0Þu½a
>
ð0Þ N ½M ð0Þ 1 ð0Þ ð0Þ
if u  0:4: L ¼ 0:418 þ 1:58u  0:580u2 a ¼ Tr exp½Lð0Þ ð1  TrM N Þ (see appendix C for more details)
>
>
if u  0:4: L ¼ 0:212 þ 2:20u  0:831u2 : að1Þ ¼ T Nð1Þ ½Mð1Þ 1 exp½Lð1Þ ð1  T Mð1Þ Nð1Þ Þ
r r

by Twu (1991) (from now on, called Twu91). Their expressions are
reported below.

 The Soave a-function [3]

8
< aðTr Þ ¼ h1 þ mðuÞ1  pffiffiffiffi
Tr
ffi i2

(8)
: mðuÞ ¼ 0:480 þ 1:574u  0:176u
2
for the Redlich  Kwong CEoS
mðuÞ ¼ 0:37464 þ 1:54226u  0:26992u2 for the Peng  Robinson CEoS

h  i
This a-function is polynomial, generalized and it has a T-overall aðTr Þ ¼ TrNðM1Þ exp L 1  TrMN (9)
formulation. Having regard to its simplicity and to the accuracy of
results obtained by using it, this a-function should be regarded as a This a-function is exponential, component-specific and it has a
reference before deriving any new a-function. These reasons justify, T-overall formulation. Parameters L, M and N have to be optimized,
in fact, its widespread popularity. component by component, over experimental data so that better
accuracy is expected in comparison to the Soave a-function.
 The Twu91 a-function [15] Component-dependent a-functions are usually parameterized
516 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

over vapor-liquid equilibrium data, which are also the most re- the residual molar enthalpy, the first derivative of a appears in the
ported data in the literature. Therefore, most of the a-functions will mathematical expression of the residual molar enthalpy and en-
have a similar behavior in the subcritical temperature domain. tropy and the second derivative of a appears in the mathematical
However, the literature lacks of specific, theoretically-based expression of the residual heat capacity at constant volume.
guidelines to derive an a-function that could guarantee consistent

8  
> 1 ac aðTÞ vb
>
> ln f ðT; vÞ ¼ v   1  ln
>
> v  b RTðv  r1 bÞðv  r2 bÞ RT
>
>
>
>  
>
> ac aðTÞ v  r1 b
>
> þ ln
>
> RTbðr1  r2 Þ v  r2 b
>
>
>
>    
< vb ac daðTÞ v  r1 b
sres ðT; vÞ ¼ R ln  ln (10)
>
> v bðr1  r2 Þ dT v  r2 b
>
>  
>
> ac aðTÞv daðTÞ v  r1 b
>
> hres ðT; vÞ ¼ RTb  ac
>
> þ aðTÞ  T ln
>
> v  b ðv  r1 bÞðv  r2 bÞ bðr1  r2 Þ dT v  r2 b
>
>
>
>  
>
> cres ðT; vÞ ¼ Tac d aðTÞ ln v  r1 b
2
>
: v
bðr1  r2 Þ dT2 v  r2 b

predictions, also in the supercritical temperature domain, of ther- We can summarize our observations by:
modynamic properties related to the a-function itself and to its 8
derivatives. With that respect, the purpose of the next section is to >
> aðTÞ  0 and continuous
present requirements that a proper a-function should fulfil. >
>
>
< da
 0 and continuous
For all T : dT (11)
3. Requirements for a consistent a-function
>
>
>
> d2 a
>
: 2 continuous
dT
Firstly, to avoid non-physical positive values of the attractive
term, a-functions should be positive and thus lead to a decrease of These conditions are often cited in the literature as sufficient
the pressure of the system when the attraction between molecules conditions for consistent predictions of thermodynamic properties
increases. As stated by Deiters [30], a negative a-function “only but: are they really sufficient? In the next section, we are going to
results from compensating an insufficient repulsion term” and answer this key question through a case study: the calculation of
should be regarded as a consequence of a maladjusted model but the isothermal phase diagram of the system CO2-argon at 253.1 K,
not as a theoretical valid behavior. with the Peng-Robinson CEoS, associated with two a-functions that
Secondly, a-functions should be constant at the infinite tem- satisfy criteria in (11). Before presenting calculation results over
perature limit. This constant can be positive or zero. Colina et al. this system, it is worth pointing out that it is a complex mixture to
[31] give proof that the limiting value of a at the infinite temper- be modelled since it exhibits a type III phase behavior in the clas-
ature limit would be equal to zero, by assuming that (i) “at infinite sification scheme of Van Konynenburg and Scott [34].
temperature, the cubic EoS reduces to a hard-body EoS, with re-
sidual internal energy identically zero at all densities” and (ii) the
4. Case-study: prediction of the CO2-argon binary phase
Van der Waals repulsive term RT=ðv  bÞ correctly represents hard-
diagram
body repulsions. However, Abbott and Prausnitz [32] and Sandler
[33] leave some uncertainty on whether the a-function at infinite
The CO2-argon binary system is considered as a key system for
temperature should assume a zero rather than a positive finite
CO2-Capture and Storage (CCS) applications [35]. It is desired to
value. Since no unquestionable theoretical evidences have been
predict the phase diagram of this system at 253.1 K with the Peng-
found on this issue, we do not conclude about the constant value of
Robinson CEoS and highly efficient mixing rules which combine the
the a-function at infinite temperature.
expression by Huron and Vidal [36] with the residual part of the
Thirdly, as temperature decreases, the average molecular kinetic
Wilson activity-coefficient model [37]. The two parameters of the
energy is reduced, resulting in a system where molecules are
Wilson model were regressed in order to minimize the following
particularly prone to interacting, upon collision. This explains why
objective function which takes into account the deviations between
the CEoS energy parameter should become increasingly important
calculated and experimental molar fractions of CO2 in the liquid (x)
when temperature decreases. In other words, a has to be a
and gas (y) phases:
decreasing function of temperature, its first derivative has not only
to be negative but has also to cancel out at the infinite temperature !2 !2
X
Nexp exp
xi  xcalc X
Nexp exp
yi  ycalc
limit to ensure a finite a-value. It should be noticed that, if an a- Fobj;1 ¼ i
þ i
(12)
function is both positive and decreasing, it will automatically admit i¼1
sexp
xi i¼1
sexp
yi
a finite value at the infinite temperature limit.
exp
The last point is that not only the a-function but also its first and where sG is the experimental standard deviation associated to the
i
second derivatives with respect to the temperature have to be ith data point of property G.
continuous functions of temperature in order to avoid discontinu- In order to test the influence of the selected a-function on the
ities in the calculated state functions. As highlighted in Eq. (10), a resulting calculations, two different a-functions, that are the Soave
indeed appears in the mathematical expression of the fugacity and and the Twu91 a-functions, were successively considered for both
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 517

pure carbon dioxide and argon. The m-parameter of the Soave a- in Fig. 2) that, at this temperature, the two a-functions return the
function was calculated by applying Eq. (8), while the three pa- same value. On the other hand, at 253.1 K, argon is supercritical
rameters L, M and N of the Twu91 a-function were determined by (Tr ¼ 1.7), that is at a temperature completely outside the range
minimizing the objective function expressed by Eq. (13), which where the parameters were fitted. At this temperature, the two a-
takes into account the deviations between calculated and pseudo- functions take totally different values (aSoave z 0.8 and aTwu z 0.2),
experimental vapor pressures (P sat ), enthalpies of vaporization as shown in Fig. 3. At this step it is thus possible to conclude that the
ðDvap HÞ and saturated liquid heat capacities (csat
P;L ): mismatch in a-values for argon is the unique responsible for the

 
X50  sat  50  
PEoS ðP; Ti Þ  PDIPPR ðTi Þ XDvap HEoS ðP; Ti Þ  Dvap HDIPPR ðTi Þ
sat
Fobj;2 ¼   þ  
 sat
PDIPPR ðTi Þ  Dvap HDIPPR ðTi Þ
i¼1 i¼1
  (13)
50 csat sat 
P  P;L;EoS ðP; Ti Þ  cP;L;DIPPR ðTi Þ
þ  
i¼1  sat
cP;L;DIPPR ðTi Þ 

The pseudo-experimental P sat , Dvap H and csat P;L data were

extracted from the correlations available in the DIPPR database: for
160
each property, 50 equidistant data points were generated in the
Soave
temperature validity range of the correlations. Following this
regression procedure, optimal parameters have been obtained for 140 Twu 91
the Twu91 a-function. Those are reported in Table 2. Exp
A significant difference between the two calculated phase dia- 120
grams of the CO2-argon system is observable in Fig. 1. To better
analyze the cause of such a difference, it is worth recalling that the 100
P / bar

two phase diagrams were modeled using the same CEoS and the
same mixing rules; they thus only differ in the choice of the a-
80
function.
When the Soave a-function is used for both compounds, a nearly
perfect prediction of the phase diagram is achieved (the continuous 60
curve is extremely close to the experimental data points). On the
other hand, very poor results are obtained with the Twu91 a- 40
function. By having a look at the dashed curve in Fig. 1, it could be
erroneously concluded that the Peng-Robinson EoS is totally unable
20
to correlate the CO2-argon system and that it over predicts the bi-
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
nary critical pressure and the solubility of argon in CO2. It is un-
deniably surprising to observe that the generalized Soave a- xAr, yAr
function, which does not contain any adjustable parameter, per-
Fig. 1. Binary phase diagram of the CO2-Ar system at 253.1 K calculated with the Peng-
forms better than the component-dependent Twu91 a-function
Robinson CEoS and two different a-functions parameterized on subcritical vapor-liquid
that involves three adjustable parameters per component. In order equilibrium data.
to better understand the origin of the significant difference be-
tween the two calculated phase diagrams, it is necessary to analyze
more closely each of the two a-functions (Soave and Twu91) which
were designed to accurately reproduce pure component vapor-
liquid equilibrium data. To this aim, a-curves of CO2 and argon
are reported in Figs. 2 and 3. The two figures highlight that, in the
temperature range where the parameters were fitted (let us say
between 0.4Tc and Tc), the Soave and Twu 91 a-functions are
superimposed for both pure components.
At the working temperature of 253.1 K, pure CO2 is subcritical
(Tc ¼ 304.21 K) while pure argon is supercritical (Tc ¼ 150.86 K).
Since for CO2 the working temperature is in the range where the
parameters were fitted, it is not surprising to observe (see x symbol

Table 2
Twu91 a-function parameters of CO2 and argon obtained by minimizing
Eq. (13).

CO2 Argon

L 0.040 0.027
M 0.943 0.968
N 8.538 8.213

Fig. 2. a-curves of CO2 parameterized on subcritical vapor-liquid equilibrium data.

518 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

Fig. 4. Illustration of the unphysical isobars crossing when the second derivative of an
a-function cancels e case of CO2.

uniquely dependent on temperature. It is thus possible to write

Fig. 3. a-curves of argon parameterized on subcritical vapor-liquid equilibrium data.
that, at T * :
   
poor correlation of the CO2-argon vapor-liquid equilibrium data cv T * ; v ¼ cv T * (17)
with the Twu91 a-function at 253.1 K. Such a difference between a-
function values in the supercritical domain, for both CO2 and argon, and whatever the molar volume is (i.e. irrespectively of the pres-
appears to be caused by the presence of two inflection points along sure), cv is the same (that of the perfect gas). Since, at the tem-
the Twu91 a-curve. As highlighted by Figs. 2 and 3, Twu91 a- perature T * , cv results to be pressure independent, all the isobaric
functions of argon and CO2 go from convex to concave to convex curves in a cv -T plane are going to intersect at such a temperature.
with a steep and quick decrease in the concave part of the a-curve. For illustration purpose, the change of cv with respect to T, calcu-
On the other hand, the Soave a-function is strictly convex on the lated for pure CO2 with the Twu91 a-function (see Table 2 for values
whole temperature range with a significantly slower decrease in of L, M and N), is shown in Fig. 4 within the pressure range
50e150 bar. As expected, at the two reduced temperatures (Tr;1 * and
the supercritical domain. Figs. 2 and 3 thus give rise to two ques-
* ) where the a-curve exhibits an inflection point (i.e. where the
Tr;2
tions which are addressed in the next section: are the two inflection
points observed on the Twu91 a-curve the actual responsible for second derivative vanishes) all the isobars (P/bar ¼ 50, 100, 150) in
the poor correlation of the binary vapor-liquid equilibrium data the (cv ,T) plane intersect. Experimental evidences about this
and, above all, are such inflection points thermodynamically peculiar behavior have never been observed and therefore we are
consistent? convinced of its unphysical foundation. Consequently, the presence
of inflection points on a-curves should be regarded as thermody-
namically inconsistent and should, thus, be avoided. In other
5. Are inflection points on a-curves thermodynamically words, the second derivative of a-functions should never cancel
consistent? out.
To avoid inconsistent inflection points on the a-curve we must
To begin, let us recall that at an inflection point, the curvature of impose the convexity of the a-function, for any temperature value:
a function changes sign and hence its second derivative vanishes:
 d2 a
 0 for all T
d2 a
(18)
 ¼0 (14) dT 2
dTr2  That being said, parameters of the Twu91 a-function of CO2 and
Tr;inflection point

argon have thus been re-regressed in this section by imposing the

Among all the thermodynamic properties that can be derived
constraint d2 a=dT 2 > 0 on the second derivative (the mathematical
from a CEoS, the residual heat capacity at constant volume is
aspects are derived in Appendix A). Corresponding values for the L,
directly proportional to the second derivative of the a-function:
M and N parameters are reported in Table 3.
  a-curves obtained with and without the imposition of the
Tac d2 a v  r2 b
cres
v ðT; vÞ ¼ ln (15) convexity of the corresponding a-function (Eq. (18)), are reported
bðr1  r2 Þ dT 2 v  r1 b in Figs. 5 and 6. The second derivative of the convex a-function is
The cv coefficient is given by:

cv ðT; vÞ ¼ cv ðTÞ þ cres

v ðT; vÞ (16) Table 3
Parameters regressed for the Twu91 a-function when convexity of the a-
curve is enforced.
with cv the heat capacity at constant volume of the perfect gas. As a
CO2 Argon
direct consequence, if an a-curve exhibits an inflection points at
temperature T * , the second derivative of a with respect to the L 0.091 0.072
temperature vanishes (see Eq. (14)) at T * , cres
v cancels out and, thus,
M 0.890 0.919
N 3.805 2.667
cv takes the value of that of the perfect gas, cv (see Eq. (16)),
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 519

Fig. 7. Binary phase diagram of the CO2-Ar system at 253.1 K when the convexity
Fig. 5. a-function of pure argon and its second derivative when the convexity
constraint is enforced (d2 a=dT 2 > 0).
constraint is enforced.

objective function value, i.e. leading to the same a-function values

in the subcritical domain, are obtained when Eq. (13) is minimized.
Such sets of parameters however lead to totally different behaviors
in the supercritical region.
Furthermore, Fig. 7 reports the comparison between experi-
mental VLE data and phase diagrams calculated by means of the PR-
EoS either combined with the Soave a-function (as previously
shown in Fig. 1) or with the re-parameterized Twu91 model (L, M, N
parameters collected in Table 3). Comparing Figs. 1 and 7, that only
differ in the VLE calculations performed with the Twu91 model, it is
possible to attest that the imposition of the condition d2 a=dT 2 > 0
to Twu91 parameters (Fig. 7) results in a spectacular improvement
of the accuracy of such a model. Fig. 7 also highlights that a similar
accuracy of prediction is achieved using either the Twu91 a-func-
tion with the convexity constraint or the Soave a-function.
A satisfactory prediction of the binary phase diagram being now
achieved by the imposition of d2 a=dT 2 > 0, our following priority
concern is the investigation of the necessity of adding new con-
straints to properly predict state functions, such as enthalpies and
heat capacities, also in the supercritical domain.

Fig. 6. a-function of pure CO2 and its second derivative when the convexity constraint
is enforced. 6. Relevance of the parameters obtained by imposing d2a/
dT2 > 0 to predict state functions in the supercritical domain

also reported, to enable the reader to attest that, as imposed, it is Pseudo-experimental enthalpy (h) and heat capacity (cP ) were
always positive. generated for pure CO2 and pure argon in the supercritical domain,
As highlighted by Fig. 5 at the working temperature of 253.1 K, in particular for reduced temperatures ranging from 1.2 to 5. The
the value of the non-convex a-function of argon (being supercritical reference Span-Wagner equation [38] was used in this work to
in this condition) is 0.19. Differently, Twu91 L, M and N parameters estimate properties of CO2 at 10, 59.7, 109, 159 and 209 bar, while
optimized under the imposition of the convexity constraint lead to the Tegeler-Span-Wagner reference equation [39] was used for
an a-function value of 0.81, at the same temperature of 253.1 K. On argon at 10, 42, 74, 106 and 138 bar. Fig. 8 shows the comparison
the other hand, at 253.1 K a-value of CO2, (being subcritical at that between properties calculated with these reference models and
temperature), is only slightly affected by the addition of the con- their values determined with the Peng-Robinson EoS combined
vexity constraint (see Fig. 6); it raises from 1.10 to 1.13. with the Twu91 a-function (see Table 3 for values of L, M, and N). As
Moreover, it should be noticed, from both Figs. 5 and 6, that it can be seen from Fig. 8c and d, molar enthalpies (h) are well
major changes on a-curves, deriving from the imposition of represented by the newly regressed a-functions for both CO2 and
constraint in Eq. (18), lie in the supercritical domain, leaving almost Ar: not only the numerical values are accurately predicted but also
unchanged the subcritical part of it. This is a straightforward the change of h with respect to the temperature is properly
consequence of the fact that, as mentioned above, only subcritical reproduced. Totally different conclusions apply for the cP . The
data were used in the optimization procedure. It is also possible to Peng-Robinson EoS with the currently evaluated a-function pre-
conclude that different sets of parameters leading to the same dicts a “wave shape” while the application of reference equations
520 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

Fig. 8. Prediction of supercritical cP of argon (a) and CO2 (b) and prediction of supercritical h of argon (c) and CO2 (d) when convexity constraint (d2 a=dT 2 > 0) only is applied.

Table 4 equations.
Parameters regressed for the Twu91 a-function when both the convexity of To check this hypothesis, the parameters of the Twu91 a-func-
the a-curve ðd2 a=dT 2 > 0Þ and the negativity of the third derivative of the tion of CO2 and argon have been thus re-optimized under the
a-function ðd3 a=dT 3 < 0Þ are enforced.
simultaneous imposition of d2 a=dT 2 > 0 and d3 a=dT 3 < 0. Optimal
CO2 Argon values of parameters are reported in Table 4.
L 0.178 0.123 a-functions that result from the use of these parameters
M 0.859 0.905 (Table 4) are represented in Figs. 9 and 10. These figures also report
N 2.411 1.854 the second derivatives, calculated with respect to temperature,
relative to these a-functions. A major change in the shape of the a-
curve which results from the use of this new set of constraints
results in a weak change of the cP with respect to temperature (a concerns the rate at which a-functions decrease to zero. It can be
minimum is only observed). Such unexpected behavior is particu- seen that the use of parameters reported in Table 3 (optimized over
larly visible in Fig. 8b, relative to CO2. By comparing Figs. 6 and 8b, it the unique constraint d2 a=dT 2 > 0) results in a-functions that equal
is noticeable that both d2 a=dT 2 and cP curves exhibit the same zero at a reduced temperature Tr ~ 4, while the application of pa-
“wave shape” in a similar temperature range. It is thus expected rameters optimized under the simultaneous imposition of
that the cP curve would exhibit a more realistic trend by removing d2 a=dT 2 > 0 and d3 a=dT 3 < 0 generates a-functions that equal zero
the “wave shape” on the d2 a=dT 2 curve. To do that, it is enough to at higher reduced temperatures, Tr > 8. Therefore, for such func-
eliminate (see Fig. 6) the increasing part of such a function. In other tions, the attractive term of the CEoS used for the calculations still
words, it is expected that the imposition of condition d3 a=dT 3 < 0 plays a relevant role at high temperature.
would result in cP variations with temperature, calculated with the Supercritical h and cP of CO2 and argon are represented in Fig. 11.
PR EoS, being similar to those obtained with the reference For both compounds, the quality of prediction of h is similar to the
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 521

Moreover, it is worth noting that P sat of CO2 is better predicted

when the negativity constraint on the third derivative of the a-
function is set and that Dvap H is calculated with the same accuracy
in both cases.
Finally, the CO2-argon binary phase diagram is represented in
Fig. 12. Impressive results are obtained by the application of the full
set of constraints (d2 a=dT 2 > 0 and d3 a=dT 3 < 0), even slightly
better than when only the convexity constraint is enforced. As
shown by Fig. 12, the resulting optimal Twu91 formulation is
characterized by an accuracy similar to that of the Soave a-function.
In conclusion, the addition of a new negativity constraint on the
third derivative of the a-function (i) leads to a thermodynamically
consistent a-function on the whole temperature range, avoiding
non-physical cP changes with respect to temperature, (ii) allows for
a good prediction of subcritical properties used in the regression
procedure and (iii) correctly extrapolates to the supercritical
domain. We can thus claim that a proper a-function must satisfy all
the following constraints:

8
> a  0 and aðT Þ continuous
Fig. 9. a-function of CO2 and it second derivative when constraints d2 a=dT 2 > 0 and >
>
>
>
d3 a=dT 3 < 0 are simultaneously applied. >
>
>
> da da
>
>  0 and continuous
>
> dT dT
>
>
r r
<
For all T : d2 a d2 a (19)
>
>  0 and continuous
>
> 2
>
> dTr dTr2
>
>
>
>
>
>
>
> d3 a
>
: 0
dTr3

The set of conditions reported in Eq. (19) establishes what we

decided to call a “consistency test for an a-function”. In other words,
if one of these constraints is not satisfied, the EoS will possibly
return e depending on the temperature and pressure domain e
inaccurate or unexpected values.
At this step, it could be argued that our derivations rely on the
unique analysis of the Ar þ CO2 binary system which contains two
molecules having the same size. In order to convince the reader
that our conclusions are general, the size-asymmetric CO2 þ n-
decane system is studied in the next section.

7. Study of a size-asymmetric binary mixture with classical

mixing-rules: case of the CO2 þ n-C10 system
Fig. 10. a-function of argon and it second derivative when constraints d2 a=dT 2 > 0 and
d3 a=dT 3 < 0 are simultaneously applied.
In this section, our objective is to test the influence of the
selected a-function on the correlation of the properties of the
one previously obtained when the third derivative of the a-func- CO2 þ n-decane size-asymmetric binary system with the PR EoS. As
tions could be positive. On the other hand, this new set of param- was previously made, three a-functions will be compared: the
eters not only guarantees more accurate cP values but also correct soave a-function, an inconsistent Twu91 a-function (which fails the
trends of cP -variations with temperature, with respect to results proposed consistency test) and a consistent-Twu91 a-function
obtained from reference models. (which passes the proposed consistency test). To treat this system,
It is legitimate now to wonder what was the prize to pay to get classical mixing rules are preferred over complex mixing rules for
such a large improvement of the state functions in the supercritical two reasons: (i) excellent accuracy of prediction is achieved for this
domain. Were subcritical properties damaged? To answer this system by only regressing a kij over experimental phase equilib-
question, deviations on P sat , Dvap H and csat rium data and (ii) discussion will be facilitated when it will come to
P;L for both compounds
have been calculated and reported in Table 5. The content of this compare different optimal kij values. Their expressions are reported
table enables the comparison between results obtained with the below:
addition of the constraint on the third derivative of the a-function 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
> N P
P N 
to the results obtained when only the convexity constraint was >
< aðT; xÞ ¼
> xi xj ai ðTÞaj ðTÞ 1  kij
enforced. As highlighted by Table 5, this huge improvement on the i¼1 j¼1
(20)
prediction of the supercritical properties is marginally counter- >
> P
N
>
: bðxÞ ¼ x i bi
balanced by a small accuracy decrease of csat P;L of CO2 and argon. i¼1
522 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

Fig. 11. Prediction of supercritical cP of argon (a) and CO2 (b) and prediction of supercritical h of argon (c) and CO2 (d) when constraints d2 a=dT 2 > 0 and d3 a=dT 3 < 0 are both
applied.

Table 5 specific constraints and (ii) by simultaneously imposing

Comparison of the errors calculated on subcritical properties for CO2 and argon d2 a=dT 2 > 0 and d3 a=dT 3 < 0 in order to get a consistent a-function.
when only a constraint on the second derivative of the a-function is applied and
At 520 K, CO2 is supercritical ðTr;CO2 ¼ 1:71Þ and as highlighted by
when the constraint on the third derivative is added to the previous one.
Figs. 2 and 9, at such a temperature, the inconsistent and consistent
Constraints applied Error Psat Error DvapH Error csat
p;L Twu91 a-functions have very different values. On the other hand, n-
Argon 2
d a 0 0.27% 2.16% 3.95% C10 is subcritical ðTr;nC10 ¼ 0:84Þ so that the two Twu91 a-func-
dTr2
d2 a d3 a 0.28% 2.30% 4.74%
tions, the parameters of which were fitted in this temperature
dTr2
 0 and dTr3
0
range, have very similar values. Calculations of the phase diagrams
CO2 d2 a 0 0.56% 0.92% 5.52%
dTr2 at 520 K are reported in Fig. 13. In this figure, no difference can be
d2 a  0 and d3 a 0 0.29% 1.01% 6.99% seen between the red full curve (Soave a-function) and the blue
dTr2 dTr3
dashed curve (consistent Twu91 a-function): they are super-
imposed even in the vicinity of the critical pressure. Fitted kij values
In a first step, the isothermal (P,x,y) phase diagram was calcu- to be used with the Soave and the consistent Twu91 a-functions are
lated at 520 K with the 3 different a-functions. For CO2, the pa- also very close with respective values of 0.14 and 0.11 which
rameters are those determined in the previous sections (see correspond to expected values of this parameter. Fig. 13 also
Tables 2 and 4) and the plots of the inconsistent and consistent highlights that the green dashed curve corresponding to the
Twu91 a-functions are shown in Figs. 2 and 9 respectively. For n- inconsistent Twu91 a-function does not correlate badly the
C10, the three parameters of the Twu91 a-function are those which experimental data points although it over predicts the critical
minimize Eq. (13). They were determined twice: (i) without any pressure. However, the major drawback with the use of this
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 523

function does not pass the proposed consistency test. Therefore, a

140 Soave too large kij value is not only an unordinary value but also a warning
Twu 91 against poorly regressed a-function parameters that might lead to
Exp inaccurate predictions of supercritical properties of pure
120 compounds.
In the next few lines, in order to convince the reader of the
100 absolute necessity of using a consistent a-function in their calcu-
lations, the mixing enthalpies (hM) of the CO2 þ n-C10 system are
P / bar

predicted with the three considered a-functions (Soave, consistent

80 and inconsistent Twu91). Calculations are performed at 125 bar and
at 470.11 K and 573.11 K and the results are reported in Fig. 14 and
60 Fig. 15. At 470.11 K and 125 bar (Fig. 14), the inconsistent Twu91 a-
function (that corresponds to kij ¼ 2:43) over predicts by a factor
100 the experimental mixing enthalpies and wrongly predicts a
40 one-phase system whatever the composition is. On the other hand,
coupling a kij of 0.11 with the consistent Twu91 a-function or a kij of
20 0.14 with the Soave a-function leads to an accurate prediction of
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 hM. It is believed that the consistent Twu91 a-function performs
better than the Soave a-function (see Fig. 14b) because its param-
xAr, yAr
eters were not only regressed on vapor pressures but also on
vaporization enthalpies and heat capacities (Soave solely parame-
Fig. 12. Binary phase diagram of the CO2-argon system at 253.1 K when constraints
d2 a=dT 2 > 0 and d3 a=dT 3 < 0 are simultaneously applied. terized his a-function with the objective of reproducing the
experimental vapor pressures). At 573.11 K and 125 bar (Fig. 15) a
noneexisting VLE is predicted and poor results are obtained for the
mixing enthalpies when a kij of 2.34 is coupled with the incon-
sistent Twu91 a-function while the two others a-functions allow
for an excellent adequacy between the experimental points and the
predicted curves.
In conclusion, we can state that the constraints defined by Eq.
(19) apply whatever the binary system. This section also pointed
out that the classical Van der Waals mixing rules, involving a kij , do
not permit to overcome the limitations of a poorly regressed a-
function. Indeed, when a very large kij is calculated to balance the
deficiencies of a non-consistent a-function, not so bad VLE corre-
lation can be achieved but extremely poor predictions of the
derived properties (here illustrated with the mixing enthalpies) is
performed.
Although it is claimed that the consistency test developed in this
paper apply whatever the component, a special attention should be
paid to the so-called “quantic fluids” like He and H2. This is the aim
of the next section.

8. The special case of quantic fluids (H2 and He)

Fig. 13. Prediction of the CO2en-C10 phase diagram at 520 K with different a-functions
He and H2, the critical temperatures of which are 5.20 K and
and classical mixing rules. Red full curve is calculated with the Soave a-function while 33.15 K respectively, are qualified of quantic fluids as their low
blue and green dashed curves are calculated with the Twu91 a-function. Parameters of temperature physical behaviors can only be represented if quantum
the a-function used to calculate the green curve do not satisfy the consistency test effects are taken into account in the modelling procedure. Conse-
while parameters of the blue curve do. (For interpretation of the references to colour in
quently, they must be treated apart from other classical fluids. As an
this figure legend, the reader is referred to the web version of this article.)
example, Gasem et al. [40] and Mahmoodi et al. [22] noticed that
when a-functions of H2 and He are both regressed on vapor pres-
inconsistent a-function is the value of the optimal kij which was sures, they exhibit an increasing ðda=dT > 0Þ and concave portion
found to be equal to 2.43, far beyond an acceptable value. This ðd2 a=dT 2 < 0Þ followed, at least for H2, by a decreasing ðda=dT < 0Þ
abnormal kij value can be explained by the role this parameter plays and convex portion ðd2 a=dT 2 > 0Þ in the supercritical domain.
when the inconsistent Twu91 a-function is used. Indeed, in this Working with quantic fluids and in order to determine the
case, the kij endorses the role of a corrective term which restores proper shape (convex or concave) of a-functions in the subcritical
the quality of prediction of the CEoS by counterbalancing the temperature range, let us recall that cres
v and the second derivative
wrongly estimated a-value of CO2 at 520 K. At this step, one could of the a-function have always the same sign. We indeed can write:
believe that it is not necessary to regress the a-function parameters  
by imposing the constraints listed in Eq. (19) as long as the kij , or Tac v  r2 b d2 a
cres
v ðT; vÞ ¼ ln  (21)
any other parameter from any mixing rule, can compensate the bðr1  r2 Þ v  r1 b dT 2
|ffl{zffl}
|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
defects of the a-function. This reasoning is however not admissible >0 phase
independant
because, as previously demonstrated, heat capacities and en-
thalpies of CO2 at 520 K will be extremely poorly predicted if the a- Along the vaporization curve, i.e. when the pure component is in
524 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

Fig. 14. Calculation of the mixing enthalpies of the CO2en-C10 system at 470.11 K and 125 bar.

Fig. 15. Calculation of the mixing enthalpies of the CO2en-C10 system at 573.11 K and Fig. 16. Residual heat capacities of vapor and liquid phases of saturated H2.
125 bar.

Although the quantum effects cannot be described by a 2-

VLE, the 2 phases have the same temperature i.e. the same d2 a=dT 2 parameter CEoS, it would be valuable to accurately represent the
value and both the liquid and vapor phases cres v calculated from a behavior of He and H2 in some specific temperature ranges due to
CEoS have the same sign (the sign of the second derivative of the a- their use in some chemical process plants. First, it must be noticed
function). Extracted from the NIST database, the quantity cres
v of He that in most applications, the operating temperature is far above
and H2 for both the liquid and vapor phases are shown in Fig. 16 and the respective critical temperatures of these two compounds.
Fig. 17. Therefore, focusing on the subcritical temperature range is not
One can see that for H2 at a temperature lower than 28.5 K, relevant. Secondly, it is expected that at reasonably high temper-
experimental cres
v of the liquid phase is negative while it is positive atures the a-function of each of these pure compounds behaves like
for the vapor phase. Therefore, at these low temperatures, where a classical fluid, i.e. it satisfies, all the criteria defined by Eq. (19).
fluids are quantic, the second derivative of the a-function should be In order to verify whether or not a-functions of He and H2 at
both positive and negative. This contradiction leads us to affirm sufficiently high temperatures behave conventionally, a regression
that the same a-function cannot represent simultaneously the of the Twu91 parameters was performed on both subcritical
liquid and vapor phases. In other words, these fluids do not follow [P sat ðTÞ and Dvap HðTÞ] and supercritical [cres res
v ðT; PÞ and h ðT; PÞ]
the 3-parameter law of corresponding states in the low tempera- data. The same weight was given to all the properties. Subcritical
ture range. It is also possible to conclude that cresv data in the cres
v were excluded from the regression procedure since values for
subcritical temperature range are not relevant data to fit the pa- the liquid and vapor phases are conflicting. The supercritical
rameters of a-functions. experimental data were generated using the NIST database at
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 525

Fig. 17. Residual heat capacities of vapor and liquid phases of saturated He. Fig. 19. Optimized a-functions of H2. Small dashed curve is regressed on subcritical
data, large dashed curve is regressed on supercritical data, full curve is regressed on
both subcritical and supercritical data.

 When the L, M and N parameters are only regressed on

subcritical data, excellent representation of the subcritical
properties is achieved although an increasing a-function is
calculated on this temperature range, demonstrating the un-
usual behavior of fluids in “quantic state”. Meanwhile, super-
critical properties are poorly represented (see Table 6). As seen
in Fig. 18, in the case of He, the a-function goes to zero at a
temperature of 40 K (red dashed curve) and, thus, is unable to
accurately extrapolate at higher temperatures making this a-
function unsuitable for most of the actual processes involving
this molecule.
 When the a-function parameters are regressed on both
subcritical and supercritical properties, a slight accuracy
decrease of prediction is observed for the subcritical properties
while supercritical ones are largely improved. This enhance-
ment is caused by the shape of the a-function at high temper-
ature, which is adjusted to represent the residual heat capacities
and enthalpies. This improvement is obvious when attention is
paid to the change of the He a-function with respect to the case
Fig. 18. Optimized a-functions of He. Small dashed curve is regressed on subcritical
when only subcritical data were regressed. Now the a-function
data, large dashed curve is regressed on supercritical data, full curve is regressed on
both subcritical and supercritical data. is very different from zero at temperatures above 150 K, dividing
by two the error over heat capacities while leaving almost un-
changed the other properties. An excellent improvement of the
pressures of 10, 100 and 1000 bar, for temperatures varying be- supercritical properties is also achieved when considering H2.
tween 150 and 750 K with a step of 50 K. For comparison purpose,  Finally, regression on solely supercritical properties is not
regression of the a-function parameters was also performed on (i) desirable as errors on saturated pressures and vaporization
only subcritical data and (ii) only supercritical data. Results of these enthalpies skyrocket while the benefit on the prediction accu-
calculations are reported in Fig. 18 and Fig. 19 and numerical values racy of supercritical properties is not significant. More generally,
of the MAPE are reported in Table 6. one can observe that on the temperature range for which
Results shown in Figs. 18 and 19 are in accordance with previous experimental data were generated and used in the regression
literature studies as almost all a-functions are increasing functions procedure, the a-function is a strictly decreasing, convex func-
of temperature in the low temperature range. The single exception tion of temperature. A graphical representation of the third
is the a-curve of H2 with parameters regressed on supercritical data derivative would show that it is negative. This confirms that
which is a monotonically decreasing, convex function of tempera- molecules which behave as “quantic fluids” in the low temper-
ture. This result is highly reassuring as at temperatures above 150 K ature range do not keep this unusual tendency at higher tem-
(the smallest temperature of the experimental supercritical data) perature and that the consistency test, defined by Eq. (19), is
we expect the fluid not to behave as a “quantic fluid” anymore, i.e. satisfied.
to respect the constraints defined by Eq. (19). The following con-
clusions can be drawn from Table 6, Figs. 18 and 19: In conclusion, a-functions of quantic fluids:
526 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

Table 6
MAPE calculated for He and H2 with respect to the type of data used in the regression procedure of the Twu91 a-function parameters (subcritical, supercritical or both).

Type of data used in the regression procedure MAPE (%)

P sat Dvap H Supercritical hres Supercritical cres

v

He Subcritical 0.001 3.68 20.79 100.00

Supercritical 141.44 78.29 21.30 24.17
Both 1.75 7.00 18.09 42.77
H2 Subcritical 0.02 3.17 54.52 58.91
Supercritical 38.23 50.23 4.90 44.22
Both 2.16 4.66 6.08 56.27

- Exhibit a bell shape (increasing and concave functions) below decreasing functions of temperature; therefore, special caution
40e50 K. should be taken when those are applied to light supercritical
- Become decreasing, convex and with a negative third derivative compounds at reduced temperatures above the temperature at
at temperatures higher than 50 K (i.e. they satisfy Eq. (19)). which the minimum of these a-functions occurs.
- Should be regressed on both subcritical and supercritical data to 2 The higher the number of adjustable parameters in a-functions,
achieve satisfying representation accuracy of the thermody- the harder the optimization procedure will be since more con-
namic properties on the whole temperature range. straints on the parameters will have to be added, to ensure that
the function passes the consistency test. Two-parameter a-
After this debate on the particular behavior of quantic fluids, the functions appear as a good trade-off between complexity and
next section aims at discussing which a-functions, published in the flexibility of the a-function formulation.
open literature, pass the proposed consistency test (see Eq. (19)). 3 Piecewise a-functions should at least have 2 parameters in the
subcritical and supercritical temperature domains. Application
9. Do current a-functions available in the open literature pass of the continuity constraint of the first and second derivatives of
the proposed consistency test? an a-function permits to establish relations between subcritical
and supercritical parameters. As a consequence, thermody-
Mathematical conditions on the a-function and its successive namic properties, such as h and cv , would be kept continuous at
derivatives aimed at ensuring consistent and accurate calculations the critical temperature.
of thermodynamic properties in both subcritical and supercritical 4 It should be noticed that, among the two generalized a-func-
domains have been identified along the previous part of this paper tions studied in this paper (Soave and Twu 1995), none of them
(see Eq. (19)). This section evaluates which of the mostly used a- satisfies all the constraints of the proposed consistency test. For
functions of the literature fulfil all these conditions. Mathematical molecules with an acentric factor out of the range ½0 ; 1 or for
derivations are detailed in Appendices A to L and only main results the prediction of heat capacities just above the critical temper-
are reported in Table 7. ature, the Soave a-function should be preferred.
The twelve selected a-functions are those which are generally
available in process simulation software. We selected:
10. Conclusion
1- Soave 1972
2- Soave 1984  The main conclusion of this article is that a-functions must:
3- Gibbons-Laughton 1. Be of class C2, i.e., their first ðda=dTÞ and second ðd2 a=dT 2 Þ de-
4- Stryjek-Vera rivatives must exist and must be continuous for any tempera-
5- Mathias-Copeman ture value,
6- Mathias 2. Be positive for any temperature value ðaðTÞ  0Þ,
7- Twu 1995 3. Be monotonically decreasing for any temperature value
8- Boston-Mathias ðda=dT  0Þ,
9- Trebble-Bishnoi 4. Be convex for any temperature value ðd2 a=dT 2  0Þ,
10- Melhem-Saini-Goodwin 5. Also verify d3 a=dT 3  0 for any temperature value.
11- Twu 1988
12- Twu 1991 All these requirements are absolutely necessary to get accurate
and physically meaningful behaviors in both the subcritical and
As highlighted in Table 7, eight over twelve of these well-known supercritical domains. Other key conclusions follow.
a-functions (Soave 1972, Soave 1984, Gibbons-Laughton, Stryjek-
Vera, Boston-Mathias, Mathias-Copeman, Mathias, Twu 1995) fail  The role of the a-function is at least as important as the role of
the proposed consistency test whatever the values of the parame- mixing rules. As highlighted with the CO2-argon binary system,
ters are. The four remaining a-functions (Trebble-Bishnoi, Melhem- even with the most elaborated mixing rules, an inconsistent a-
Saini-Goodwin, Twu 1988, Twu 1991) pass the test but only if function will always lead to poor results if at least one compo-
constraints (sometimes drastic) are added to the parameters. nent is supercritical. The reverse is also true: a consistent a-
We can thus conclude that, presently, it is impossible to find an function coupled with unsuitable mixing rules cannot lead to
a-function which passes our consistency test for any values of the accurate results.
parameters.  When the parameters of component-specific a-functions are
From a practical point of view: fitted in order to pass our consistency test [Eq. (19)], pure
component properties in the supercritical region are highly
1 All a-functions which derive from the Soave a-function, the so- improved without deteriorating the accuracy of the calculated
called “polynomial a-functions”, are not monotonically subcritical properties.
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 527

Table 7
Mathematical analysis of the mostly used a-functions.

a-function Class a  0 da  0 d2 a 0 d3 a  0 Condition(s)

dTr dTr2 dTr3
C2a
8 pffiffiffiffiffi
Soave (1972) [3] < a ¼ ½1 þ mðuÞð1  Tr Þ2 mð1 þ mÞ  0
mRK ðuÞ ¼ 0:480 þ 1:574u  0:176u2
:
mPR ðuÞ ¼ 0:37464 þ 1:54226u  0:26992u2
 
Soave (1984) a0 b>0
a ¼ 1 þ ð1  Tr Þ a þ Tbr
[24]

a 0

pffiffiffiffiffi
Gibbons- a ¼ 1 þ XðTr  1Þ þ Yð Tr  1Þ X  0 pffiffiffiffi pffiffiffiffi
Laughton 2X  2 X  Y  minð0; 2X þ 2 X Þ
[25]
pffiffiffiffiffi
Stryjek-Vera a ¼ ½1 þ kð1  pTrffiffiffiffiÞffi 2 kð1 þ kÞ  0
[11] k ¼ k0 þ k1 ð1 þ Tr Þð0:7  Tr Þ
k0 ¼ 0:378893 þ 1:4897153u  0:17131848u2 þ 0:0196554u3
k1 ðT  0:7Tc Þs0
k1 ðT > 0:7Tc Þ ¼ 0
pffiffiffiffiffi pffiffiffiffiffi 2 pffiffiffiffiffi 3 2
Mathias- aSUB ¼ aðTr  1Þ ¼ ½1 þ c1 ð1  Tpr ffiffiffiffiÞþ
ffi c22 ð1  Tr Þ þ c3 ð1  Tr Þ 
Copeman aSUP ¼ aðTr > 1Þ ¼ ½1 þ mðuÞð1  Tr Þ
[10]
pffiffiffiffiffi
Mathias [28] aSUB ¼ aðTr  1Þ ¼ ½1 þ mðuÞð1  Tr Þ  pð1  Tr Þð0:7  Tr Þ2
8
>
> a ¼ aðTr > 1Þ ¼ ½exp½Lð1  Trg Þ2
> SUP
>
>
>
>
< m
L ¼ 1 þ þ 0:3p
> 2
>
>
>
>
>
> L1
: g¼
L
8
Twu (1995) >
> a ¼ að0Þ þ u½að1Þ  að0Þ 
u2½0; 1
[16,17] < ð0Þ ð0Þ
að0Þ ¼ TrN ½M 1 exp½Lð0Þ ð1  TrM N Þ
ð0Þ ð0Þ

>
>
: að1Þ ¼ T Nð1Þ ½Mð1Þ 1 exp½Lð1Þ ð1  T Mð1Þ Nð1Þ Þ
r r
(see coefficients in Appendix C)
pffiffiffiffiffi 2
Boston-Mathias aSUB ¼ aðTr  1Þ ¼ aSOAVE ¼ ½1 þ mðuÞð1  Tr Þ 0m3
[19]
aSUP ¼ aðTr  1Þ ¼ exp½Lð1  Trg Þ
8
> 2m
>
>L ¼ 1 þ m
<
>
>
:g ¼ 1 þ m
>
2
Trebble-Bishnoi a ¼ exp½Lð1  Tr Þ L0
[26]
pffiffiffiffiffi
Melhem-Saini- a ¼ exp½mð1  Tr Þ þ nð1  Tr Þ2  mn0
Goodwin
[27]
Twu (1988) [14] a ¼ T 2ðM1Þ exp½Lð1  T 2M Þ
r
LM  0
r
M  0:8909
8
Twu (1991) [15] a ¼ Trd exp½Lð1  Trg Þ <d  0
Lg  0
:
g1d
8
or
>
> d0
< Lg  0
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
>
> g  1  2d þ 2 dðd  1Þ
: 3
4Y þ 4ZX þ 27Z 2  18XYZ  X 2 Y 2  0
3

with
8
>
< X ¼ 3ðg þ d  1Þ
Y ¼ g2 þ 3dg  3g þ 3d2  6d þ 2
>
: Z ¼ dðd2  3d þ 2Þ

a
a function is said to be of class C2 if the function, its first derivative and its second derivative are all continuous.
528 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

 The parameter-optimization procedure, aimed at determining It is recalled that a given function is said to be of class C2 if its
the parameters of a given a-function by minimizing the de- first and second derivatives exist and are continuous.
viations between calculated and experimental vapor-liquid
equilibrium data, can generate many sets of parameters lead-  Is the a-function always positive? Yes.
ing to similar objective function values (meaning that the same  Is the first derivative of the a-function always negative?
behavior is observed when temperature is below the critical
temperature of the pure compound) but inducing totally The first derivative of this a-function is
different behaviors in the supercritical domain. The constraints
proposed in this paper make it possible to select the best set of 8
parameters. It is believed that similar constraints should be < da ðTr Þ ¼ aðTr Þ PðTr Þ
>
developed for SAFT-type EoS [41]. It is indeed well-known that dTr Tr (A.4)
>
:
the fitting of the three parameters m, s and ε on saturation PðTr Þ ¼ d þ LgTr
g
pressures and liquid densities leads to many triplets of solutions
The quantity dTda ðT Þ is positive provided PðT Þ is positive itself.
for which the objective function is the same. r
r r
 None of the a-functions currently published in the open litera- Study of the sign of PðTr Þ (reported in Table A.1): Tr is positive so
ture passes the proposed consistency test. Researches in this P is a monotonic function of Tr . Therefore, PðTr Þ is positive for any
area thus need to be intensified. value of Tr > 0 if the limits of P are both positive when Tr tends to
0 or to þ∞.

Table A.1
Limit of P at the infinite temperature limit with respect to the signs of L and g.
8 g 8 g
< lim Tr ¼ 0 < lim Tr ¼ þ∞
Tr /0 Tr /0
g > 00 g < 00
: lim Trg ¼ þ∞ : lim Trg ¼ 0
Tr /þ∞ Tr /þ∞

L>0 L<0 L>0 L<0

lim PðTr Þ d d ∞ þ∞
Tr /0
lim PðTr Þ þ∞ ∞ d d
Tr /þ∞
Acceptable solution? Yes, if d < 0 No No Yes, if d < 0

Acknowledgment Summary: PðTr Þ is positive for any positive value of Tr if:

The French Petroleum Company TOTAL and more particularly Lg > 0

 and Freddy (A.5)
Dr. Pierre DUCHETeSUCHAUX, Dr. Laurent AVAULLEE d<0
GARCIA (experts in thermodynamics) are gratefully acknowledged
for sponsoring this research.

 Is the second derivative of the a-function always positive?

Appendices
Following Eq. (A.5), it is now assumed that Lg > 0 and d < 0. The
Appendix A: Mathematical study of the Twu91 a-function [15] second derivative of the Twu91 a-function writes
8
Although classically expressed by introducing 3 adjustable pa- >
> d2 a aðTr Þ
>
> ðTr Þ ¼ Q ðqÞ
rameters denoted L, M and N (see Eq. (A.1)), the Twu91 a-function >
< dTr2 Tr2
can be alternatively written under a simpler form involving three (A.6)
>
> 2
Q ðqÞ ¼ q  ð2d þ g  1Þq þ dðd  1Þ
different adjustable parameters denoted L, d and g (this form ap- >
>
>
:
pears more convenient to perform a mathematical analysis of the q ¼ LgTrg
function and its derivatives):
h  i Tr and Lg being positive, the variable q is also positive. Polynomial
aðTr Þ ¼ TrNðM1Þ exp L 1  TrMN (A.1) Q ðqÞ is positive if [42]:
8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
   < ð2d þ g  1Þ >  2 dðd  1Þ
aðTr Þ ¼ Trd exp L 1  Trg (A.2)
:
and (A.7)
dðd  1Þ > 0
where d and g are related to M and N following:
Since d < 0, the second condition of Eq. (A.7) is already true and
8 can be left aside. Finally, Eq. (A.7) leads to a single inequality:
<N ¼ g  d
d ¼ NðM  1Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
⇔ g (A.3)
g ¼ NM :M ¼ g < 1  2d þ 2 dðd  1Þ (A.8)
gd

 Is the a-function of class C2? Yes.  Is the third derivative of the a-function always negative?
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 529

Inequalities (A.5) and (A.8) are now assumed to be true. The Conclusion:
third derivative of the Twu91 a-function is: Therefore, the Twu91 a-function is a positive decreasing convex
8 function of temperature with a negative third derivative if one of
>
> d3 a aðTr Þ the two following sets of constraints is verified.
>
> ðT Þ ¼  3 RðqÞ
>
> 3 r 8
>
> dTr Tr
8
>
> >
> d<0
>
>   <d<0 < Lg > 0
< pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
RðqÞ ¼ q3  3ðg þ d  1Þq2 þ g2 þ 3dg  3g þ 3d2  6d þ 2 q Lg > 0 or
> : > g < 1  2d þ 2 dðd  1Þ
>
>
>   g<1  d : 3
>
> 2 4Y þ 4ZX 3 þ 27Z 2  18XYZ  X 2 Y 2 > 0
>
>  d d  3d þ 2
>
> (A.13)
>
>
>
:
q ¼ LgTrg
Where the quantities X, Y and Z are defined by Eq. (A.10).
(A.9)

with q > 0. The quantities X, Y and Z are defined as Appendix B: mathematical study of the Twu88 a-function [14]
8
< X ¼ 3ðg þ d  1Þ
>
2 This a-function writes
Y ¼ g2 þ 3dg  3g þ
 3d  6d þ 2 (A.10)
> h  i
: Z ¼ d d2  3d þ 2 ¼ dðd  1Þðd  2Þ
aðTr Þ ¼ Tr2ðM1Þ exp L 1  Tr2M (B.1)

Polynomial RðqÞ is positive if [42]: This a-function is a particular case of the Twu91 a-function with
8 N set to 2. Therefore, results derived in Appendix A for the Twu91 a-
<X>0
Z >0 function can be used to determine the conditions on the L and M
Y > 0 or
: H ¼ 4Y 3 þ 4ZX 3 þ 27Z 2  18XYZ  X 2 Y 2 > 0 parameters that lead to a consistent Twu88 a-function.
Z>0
For this a-function, the relations between parameters M, N and
(A.11) d, g are
Since d < 0 (see Eq. (A.5)), the condition Z > 0 (where Z is defined 8
< d ¼ 2ðM  1Þ
by Eq. (A.10)) is always verified. According to Eqs. (A.9) and (A.11), it
g ¼ 2M (B.2)
can be concluded that dT d3 a ðT Þ is negative provided: :
3
r
r
N¼2
X>0
or H ¼ 4Y 3 þ 4ZX 3 þ 27Z 2  18XYZ  X 2 Y 2 > 0
Y >0
(A.12)  Is the a-function of class C2? Yes.
A graphical representation of the constraints Y < 0 (red area)  Is the a-function always positive? Yes.
and X < 0 (yellow area) is proposed in Figure A.1. Colored areas are  Is the first derivative of the a-function always negative?
forbidden to parameters g and d, only the white area is permitted.
The concavity constraint of the a-function is also represented (the By combining Eqs. (A.5) and (B.2), one obtains:
2
blue area is an area wherein the condition ddTa2 ðTr Þ  0 is not satis-
r Lg > 0 M<1
fied). An analysis of the graphic shows that if the constraint X > 0 is ⇔ (B.3)
satisfied then the convexity constraint and the Y > 0 constraint are
d<0 LM > 0
automatically satisfied.

 Is the second derivative of the a-function always positive?

It is now assumed that M < 1 and LM > 0. By combining Eqs. (A.8)

and (B.2), one obtains:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
g < 1  2d þ 2 dðd  1Þ ⇔ 6M  5 < 2 2ðM  1Þð2M  3Þ
(B.4)
A mathematical study of inequality (B.4) makes it possible to
demonstrate that the only acceptable values of M are such that:

1 1
M þ pffiffiffiz0:947 (B.5)
2 5
Note that the combination of inequalities M < 1 and (B.5) results
in the single inequality (B.5).

 Is the third derivative of the a-function always negative?

One assumes that LM > 0 and M  12 þ p1ffiffiffi. According to Eqs.

Fig. A.1. Representation of the set of constraints defined by Eq. (A.11). 5
(A.10), (A.12) and (B.2), two possible sets of conditions ensure the
3
d a
negativity of dT 3 ðTr Þ:
r
530 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

1. Either X > 0 and Y > 0: Table C.1

. Universal parameters regressed for the Twu95 a-function.
As for the Twu91 a-function, the constraint X > 0 entails the Tr  1 Tr > 1
convexity of the Twu88 a-function and the inequality Y > 0. (0) (1)
a a a(0) a(1)
X>0 ⇔ g<1  d ⇔ M < 0:75 (B.6) Redlich-Kwong L 0.141599 0.500315 0.441411 0.032580
M 0.919422 0.799457 6.500018 1.289098
N 2.496441 3.291790 0.200000 8.000000
Peng-Robinson L 0.125283 0.511614 0.401219 0.024955
M 0.911807 0.784054 4.963070 1.248089
2. Or H ¼ 4Y 3 þ 4ZX 3 þ 27Z 2  18XYZ  X 2 Y 2 > 0: N 1.948150 2.812520 0.200000 8.000000

In the case of the Twu88 a-function, this inequality writes:

16448M 6 þ 55296M5  64416M 4  29376M 3  4020M2

þ 216M  4  0
Different sets of parameters are used in the subcritical and su-
(B.7)
percritical domains. These parameters were regressed to obtain an
By solving numerically this inequality, one obtains: a-function of class C2.

M2½0:0448; 0:8909∪½1; 1:2902 (B.8)  Is the a-function of class C2? Yes, due to the addition of this
constraint during the regression procedure of the universal
Combining inequalities (B.5) and (B.8) leads to: parameters (see Figure C.1)

M2½0:0448; 0:8909 (B.9) The presence of salient points on the second derivative
curves of this a-function can however be noticed (which is the
As previous, it can be proved that choosing the M parameter in consequence of a non-continuous third derivative of the a-
this range ensures also the convexity of the Twu88 a-function. function).
Finally, according to Eqs. (B.6) and (B.9), one deduces that the
third derivative of the Twu88 a-function is negative provided:  Is the a-function always positive? No (see a-function curves
plotted for negative or high positive values of u in Figure C.1(a)).
M < 0:8909 (B.10)
 Is the first derivative of the a-function always negative? No
(see Figure C.1(b)).
 Is the second derivative of the a-function always positive? No
(see Figure C.1(c)).
Conclusion:  Is the third derivative of the a-function always negative? No
To summarize, the Twu88 a-function is a positive decreasing (see Figure C.1(d)).
convex function of temperature with a negative third derivative if
the following set of constraints is verified (resulting from the
combination of Eqs. (B.3) and (B.10)):
Conclusion:
LM > 0 To conclude, Figure C.1 shows that the constraints on the first
(B.11)
M < 0:8909 and second derivatives (see Eq. (19)) are passed by molecules
having an acentric factor ranging between 0 and 1. For mole-
cules having a higher acentric factor, the first and second de-
rivatives of the Twu95 a-function fail the proposed consistency
Appendix C: mathematical analysis of the Twu95 a-function [16,17] test (this is for instance the case of n-alkanes heavier than n-
C22).
This a-function is said generalized since it does not involve Let us note also that regardless of the u-value, the sign of the
adjustable parameters. The sole knowledge of an experimental third derivative of the Twu95 a-function is not constant on the
value of the acentric factor (u) makes it possible to evaluate aðTr Þ. whole temperature range (see Figure C.1(d)). Inevitably, poor pre-
dictions of heat capacities are obtained in the temperature domain
where the third derivative of the a-function is positive. For all
8 h i molecules, this temperature domain is located in the vicinity of the
>
> aðTr Þ ¼ að0Þ ðTr Þ þ u að1Þ ðTr Þ  að0Þ ðTr Þ
< h  i critical temperature, i.e., it can be reached in many practical cases.
N ð0Þ ½M ð0Þ 1 ð0Þ ð0Þ
að0Þ ðTr Þ ¼ Tr exp Lð0Þ 1  TrM N (C.1) Predictions of heat capacities with the Twu95 a-function could thus
>
> h  i be less accurate than those calculated with other simpler a-func-
: ð1Þ N ð1Þ ½M ð1Þ 1 ð1Þ ð1Þ
a ðTr Þ ¼ Tr exp Lð1Þ 1  TrM N tions.

The universal parameters (LðkÞ ,M ðkÞ ,N ðkÞ for k ¼ 0 and k ¼ 1)

involved in this a-function are reported in Table C.1.
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 531

Fig. C.1. The Twu95 a-function (a) and its first (b), second (c) and third derivatives (d) (case of the RK EoS).

Appendix D: mathematical study of the Trebble and Bishnoi

a-function [26] d2 a
ðTr Þ ¼ L2 ,aðTr Þ (D.3)
dTr2
This a-function contains only one adjustable parameter (L):

 Is the third derivative of the a-function always negative?

aðTr Þ ¼ exp½Lð1  Tr Þ (D.1)
This a-function can be seen as a particular case of the Twu91 a- The third derivative of the Trebble-Bishnoi a-function is:
function with N set to 1 and M set to 1.
d3 a
ðTr Þ ¼ L3 ,aðTr Þ (D.4)
 Is the a-function of class C ? Yes.
2
dTr3
 Is the a-function always positive? Yes.
 Is the first derivative of the a-function always negative?

The first derivative of this a-function writes It is assumed that L > 0. Doing so, a > 0 and d3 a ðT Þ < 0.
dTr3 r

da Conclusion:
ðTr Þ ¼ L,aðTr Þ (D.2)
dTr To conclude, the Trebble and Bishnoi a-function passes the
proposed consistency test provided L > 0.
As a consequence, da/dTr is negative provided L is positive.

 Is the second derivative of the a-function always positive? Appendix E: mathematical study of the Melhem-Saini-Goodwin a-
Yes. function [27]

The second derivative of this a-function indeed writes This a-function involves 2 adjustable parameters (m and n):
532 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

By noting that:
 pffiffiffiffiffi 2
aðTr Þ ¼ exp mð1  Tr Þ þ n 1  Tr (E.1)  
m > n > 006mn > 6n2 04n3 þ 6mn  6n2 > 0 (E.6)

it can be claimed that under assumption (E.3), all the coefficients of

pffiffiffiffiffi
the polynomial R are positive and thus Rð Tr Þ  0. Consequently,
 Is the a-function of class C ? Yes. 2
the third derivative of the Melhem-Saini-Goodwin a-function is
 Is the a-function always positive? Yes.
always negative.
 Is the first derivative of the a-function always negative?

The first derivative of the a-function is: Conclusion:

The Melhem-Saini-Goodwin a-function is consistent (i.e., posi-
8
>
< da ðTr Þ ¼ apðTr Þ tive, decreasing, convex with a negative third derivative) provided
ffiffiffiffiffi PðTr Þ
dTr Tr (E.2) m  n  0.
>
: pffiffiffiffiffi
PðTr Þ ¼ ðm  nÞ Tr þ n
Appendix F: mathematical study of the Soave72 a-function [24]
The quantity PðTr Þ is positive if:
The Soave72 a-function involved a generalized parameter (m)
mn0 correlated to the acentric factor (u):
⇔m  n  0 (E.3)
n0
8 h  pffiffiffiffiffi i2
< aðTr Þ ¼ 1 þ mðuÞ 1  Tr
: mðuÞ ¼ 0:480 þ 1:574u  0:176u for the RK EoS
2

 Is the second derivative of the a-function always positive? mðuÞ ¼ 0:37464 þ 1:54226u  0:26992u2 for the PR EoS
Yes. (F.1)

The second derivative of the Melhem-Saini-Goodwin a-function The Soave72 function is certainly the most employed a-function
is: in the oil and gas industries. It is in particular used in the PPR78
[43e45] and PR2SRK predictive models [46].
8
>
> aðTr Þ pffiffiffiffiffi   Is the a-function of class C2? Yes.
>
> d2 a
>
> ðTr Þ ¼ 3=2 Q Tr
 Is the a-function always positive? Yes.
>
> dTr2
2Tr
>
>  Is the first derivative of the a-function always negative?
>  ffi pffiffiffiffiffi 3 pffiffiffiffiffi 2
< Q pffiffiffiffi
>
Tr ¼ 2ðm  nÞ2 Tr þ 4nðm  nÞ Tr
|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} The first derivative of the Soave72 a-function is:
>
> >0 >0
>
> sffiffiffiffiffiffiffiffiffiffiffi
>
> pffiffiffiffiffi  mðuÞ h  pffiffiffiffiffi i aðTr Þ
>
> þ 2n2 Tr þ n
da
ðTr Þ ¼  pffiffiffiffiffi 1 þ mðuÞ 1  Tr ¼ mðuÞ
>
>
>
> |ffl{zffl} |{z} dTr Tr Tr
: >0 >0 (F.2)
(E.4) The limit of the first derivative of the a-function when the
temperature tends to þ∞ is:
Under assumption (E.3), all the coefficients of the polynomial Q
are positive, therefore, Q is always positive. da
lim ðTr Þ ¼ ½mðuÞ2 > 0 (F.3)
Tr /þ∞ dTr
 Is the third derivative of the a-function always negative?
It appears that the first derivative of the Soave72 a-function is
The third derivative of this a-function is: not always negative.
pffiffiffiffiffi
Due to its polynomial expression in Tr , the Soave72 a-function
can only exhibit two types of behavior, depending on the u-value:

8 3
>
> d a aðTr Þ pffiffiffiffiffi 
>
> ðTr Þ ¼  5=2 R Tr
>
> 3
dTr
>
> 4Tr
>
> pffiffiffiffiffi 5 pffiffiffiffiffi 4 pffiffiffiffiffi 3
>
> pffiffiffiffiffi 
>
< R Tr ¼ 4ðm  nÞ3 Tr þ 12nðm  nÞ2 Tr þ 12n2 ðm  nÞ Tr
|fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} (E.5)
>
>
>
> >0 >0 >0
>
> pffiffiffiffiffi 2 pffiffiffiffiffi 
>
>  3
>
> þ 4n þ 6mn  6n2 Tr þ 6n2 Tr þ 3n
>
>
>
: |ffl{zffl} |{z}
>0 >0
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 533

- either the function is decreasing then increasing for increasing (i.e., a large range of acentric-factor values is covered). It is observed
values of Tr , that in practice, Tmin is often higher than 1500 K, largely above
- or, the function is monotonically increasing. common temperatures encountered in practical applications.

To discriminate between these two scenarios, the limit of the

first derivative of the a-function when the temperature tends to 0 is Table F.1
now studied: Values of Tmin for some chemical compounds using either the PR or RK CEoS

Component Acentric factor u Tc/K RK CEoS PR CEoS

da
lim ðTr Þ ¼ sign½mð1 þ mÞ∞ (F.4) Tr;min Tmin =K Tr;min Tmin =K
Tr /0 dTr
Methane 0.012 190.564 9.045 1723.566 12.591 2399.325
Case 1: mð1 þ mÞ > 0 n-decane 0.492 617.700 3.330 2057.108 3.748 2314.913
If mð1 þ mÞ is positive (what is the case for more than 99% of the n-eicosane 0.907 768.000 2.456 1886.582 2.705 2077.277
da ðT Þ ¼ ∞ and the Soave72 a-func- n-triacontane 1.307 844.000 2.094 1767.355 2.305 1945.584
compounds), one has lim dT r
Tr /0 r Water 0.345 647.096 3.992 2583.530 4.595 2973.520
tion exhibits a minimum at a reduced temperature denoted Tr;min . CO2 0.224 304.21 4.905 1492.258 5.839 1776.253
For Tr 2½0 ; Tr;min , a is a decreasing function while it is an Argon 0.000 150.860 9.507 1434.218 13.463 2031.065
increasing function for Tr > Tr;min. Such a behavior is observed for Ethanol 0.644 514.000 2.904 1492.795 3.228 1659.020
the following values of the acentric factor: Ammonia 0.253 405.650 4.641 1882.512 5.469 2218.676

u21;0:8580∪0:2952;9:2384½∪½9:8012;þ∞½
(F.5) Case 2: mð1 þ mÞ < 0
for the RK EoS
This case is observed for the following values of the acentric
factor:
u21;0:7838∪0:2334;5:9472½∪½6:4976;þ∞½
(F.6)
for the PR EoS u2  0:8580; 0:2952∪9:2384; 9:8012½ for the RK EoS
(F.8)
To understand whether the irregular behavior of the Soave72 a-
function (i.e., increasing of a for Tr >Tr;min ) may affect practical
applications, the relation between Tr;min and u is now investigated. u2  0:7838; 0:2334∪5:9472; 6:4976½ for the PR EoS
Tr;min is obtained by solving: (F.9)

da   In practice, this instance arises for molecules having negative

T ¼0 acentric factors (e.g., helium, aluminum and arsenic). In such a case:
dTr r;min
  pffiffiffiffiffiffiffiffiffiffiffiffi  lim da ðTr Þ ¼ þ∞ and the Soave72 a-function is an increasing
⇔mðuÞ 1 þ mðuÞ 1  Tr;min ¼ 0 (F.7) Tr /0 dTr
function of temperature for any positive value of Tr . This behavior is
1 þ mðuÞ 2
not desirable and therefore, the Soave72 a-function should be used
⇔Tr;min ¼
mðuÞ only for molecules having an acentric factor u in the ranges defined
by Eq. (F.5) or Eq. (F.6).
Note that at Tr;min , the Soave72 a-function cancels: aðTr;min Þ ¼ 0.
The evolution of Tr;min with respect to the acentric factor is re-
 Is the second derivative of the a-function always positive?
ported in Figure F.1.
It is now assumed that mð1 þ mÞ > 0 (which is obtained for
acentric factor values lying in the ranges defined by Eq. (F.5) or Eq.
To complete the analysis, some values of Tr;min and Tmin are
(F.6)). The second derivative of the Soave72 a-function is:
reported in Table F.1 for a series of real compounds of various sizes
d2 a mðuÞ½mðuÞ þ 1
ðTr Þ ¼ (F.10)
dTr2 3=2
2Tr
Since mð1 þ mÞ > 0, the convexity constraint is verified.

 Is the third derivative of the a-function always negative?

The third derivative of the Soave72 a-function is:

d3 a 3 mðuÞ½mðuÞ þ 1
ðTr Þ ¼  (F.11)
dTr3 4 5=2
Tr
As previous, assuming mð1 þ mÞ > 0, the negativity of the third
derivative is ensured.

Conclusion:
When the quantity mð1 þ mÞ is positive (i.e., for acentric-factor
values lying in the ranges defined by Eq. (F.5) or Eq. (F.6), which is
verified in practice for more than 99% of real components), the
Soave a-function is consistent up to the temperature defined by Eq.
Fig. F.1. Tr;min as a function of u.
(F.7), which is generally higher than 1500 K (it can thus be claimed
534 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

that the Soave72 a-function is consistent on temperature ranges of

practical interest). d2 a 2b
ðTr Þ ¼ 3 (G.6)
Such an a-function can thus be considered as fully consistent in dTr2 Tr
most cases explaining why it has achieved great success in the d2 a ðT Þ
It appears that dTr2 r is positive provided b is positive.
industry.

 Is the third derivative of the a-function always negative?

The third derivative of the a-function writes

Appendix G: mathematical study of the Soave84 a-function [24]
d3 a 6b
This a-function involves two adjustable parameters a and b: ðTr Þ ¼  4 (G.7)
dTr3 Tr
  It appears that d3 a ðT Þ is negative provided b is positive.
b dTr3 r
aðTr Þ ¼ 1 þ ð1  Tr Þ a þ (G.1)
Tr
Conclusion:
The only possibility for satisfying both the positivity condition
for the Soave84 a-function and the negativity condition for its first
 Is the a-function of class C2? Yes. derivative would be to set:
 Is the a-function always positive?
a¼0 (G.8)
The Soave84 a-function can be written as a second order poly-
This would lead to modify the mathematical expression of the
nomial of the reduced temperature
Soave84 a-function. As a conclusion, the Soave84 a-function cannot
8 pass the proposed consistency test.
< aðTr Þ ¼ 1 FðTr Þ
>
Tr (G.2)
>
: Appendix H: mathematical analysis of the Stryjek-Vera a-function
FðTr Þ ¼ aTr2 þ Tr ð1 þ a  bÞ þ b [11]
Polynomial FðTr Þ is positive provided [42]:
The Stryjek-Vera a-function involves one adjustable parameter
8 (k1 ) and one universal parameter (k0 ðuÞ):
<a  0
b0 pffiffiffiffiffiffiffiffiffi (G.3) h  pffiffiffiffiffi i2
: aðTr Þ ¼ 1þ kðTr Þ, 1 Tr
1 þ a  b  2 ab
 pffiffiffiffiffi 
kðTr Þ ¼ k0 ðuÞþ k1 1þ Tr ð0:7Tr Þ
k0 ðuÞ ¼ 0:378893þ1:4897153u 0:17131848u2 þ0:0196554u3
 Is the first derivative of the a-function always negative? k1 ðTr  0:7Þs0
k1 ðTr >0:7Þ ¼ 0
The first derivative of the Soave84 a-function is: (H.1)
8 The Stryjek-Vera a-function is identical to the Soave72 a-func-
>
> da PðTr Þ tion for reduced temperatures above 0.7. Therefore, it suffers from
< ðTr Þ ¼  2
dTr Tr (G.4) all the identified defects of the Soave72 a-function described in
>
>
: 2 Appendix F. In addition, the Stryjek-Vera a-function exhibits a
PðTr Þ ¼ aTr þ b
discontinuous behavior in Tr ¼ 0:7 which is prejudicial to all derived
Polynomial PðTr Þ is positive provided: EoS-properties (enthalpy, entropy, heat capacity etc.).

8
<a  0 Appendix I: mathematical study of the Mathias-Copeman a-
b  0 pffiffiffiffiffiffi (G.5) function [10]
:
0  2 ab
This a-function uses two different expressions depending on the
Note that the last condition is automatically satisfied. It appears
temperature domain considered (sub- or super-critical) and in-
that constraints (G.3) and (G.5) are contradictory and cannot be
volves three adjustable parameters (c1 , c2 , c3 )
satisfied simultaneously.
Therefore, if a is positive: 8 "  pffiffiffiffiffi   pffiffiffiffiffi 2 #2
>
< if Tr  1 : a ðTr Þ ¼ 1 þ c1 1  Tr þ c2 1  Tr
SUB pffiffiffiffiffi 3
- aðTr Þ is negative, at least, on a certain temperature range aðTr Þ ¼ þc3 1  Tr
da ðT Þ is negative for any (positive) temperature value
- dT r
>
: h  pffiffiffiffiffi i2
r if Tr > 1 : aSUP ðTr Þ ¼ 1 þ c1 1  Tr
On the other hand, if a is negative: (I.1)

- aðTr Þ is positive for any (positive) temperature value

da ðT Þ is positive, at least, on a certain temperature range
- dT
r
r
 Is the second derivative of the a-function always positive?  Is the a-function of class C2?

The second derivative of the Soave84 a-function is: The first derivatives of the Mathias-Copeman a-function are:
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 535

8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi"  pffiffiffiffiffi 
c1 þ 2c2 1  Tr #
>
> da a ðT Þ
>
> if Tr  1 : SUB
ðTr Þ ¼  SUB r
>
>  pffiffiffiffiffi 
< dTr Tr
da þ3c3 1 þ Tr  2 Tr
ðTr Þ ¼ (I.2)
dTr >
> sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
>
>
>
> da aSUP ðTr Þ
: if Tr > 1 : SUP
ðTr Þ ¼ c1
dTr Tr

lim aðTr Þ ¼ lim aSUP ðTr Þ ¼ þ∞ (I.6)

Tr /þ∞ Tr /þ∞
It appears that:

daSUB da
ðTr ¼ 1Þ ¼ SUP ðTr ¼ 1Þ ¼ c1 (I.3)
dTr dTr
The Mathias-Copeman a-function is thus at least of class C1.
The second derivatives of the Mathias-Copeman a-function are: Appendix J: mathematical study of the Boston-Mathias a-function
[19]
8 2
d aSUB PðTr Þ
>
> if Tr  1 : dT 2 ðTr Þ ¼ 2T 3=2
>
> The Boston-Mathias a-function contains one universal param-
>
> r r
>
> eter mðuÞ:
>
>
>
>
5=2
PðTr Þ ¼ 12c3 Tr  15c3 ðc2 þ 3c3 ÞTr2
2
>
>
>
>   8 h  pffiffiffiffiffi i2
>
> 3=2
>
> þ4Tr
> 2c1 c3 þ 10c2 c3 þ c22 þ 15c23 > if Tr  1 : aSUB ðTr Þ ¼ aSOAVE ðTr Þ ¼ 1þmðuÞ, 1 Tr
>
>
< >
>
d2 a   <  
ðTr Þ ¼ 3T 4c c þ c c þ c þ 2c2 þ 10c2 þ 10c c g 
dTr2 >
> r 1 3 1 2 3 2 3 2 3 aðTr Þ ¼ if Tr >1 : aSUP ðTr Þ ¼ exp L 1Tr
>
> >
>
> >
>
>
>
>
> þ3c23 þ 2c22 þ c21 þ 3c3 þ c1 þ 3c1 c2 þ 5c2 c3 : with : L ¼ mðuÞ and g ¼ 1þ mðuÞ
>
>
> g 2
>
>
>
> þ4c1 c3 þ 2c2
>
> (J.1)
>
>
>
> d2 aSUP c ð1 þ c Þ
: if Tr > 1 : ðTr Þ ¼ 1 3=2 1 The relations between parameters L, g and mðuÞ were derived
2
dTr 2Tr by Boston and Mathias assuming that
(I.4) pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi
- The first derivatives of aSUB and aSUP are equal at the critical
Therefore, one has: temperature,
pffiffiffiffiffiffiffiffiffiffi
8 2 - The second derivative of aSUP has to be equal to zero at the
>
> d aSUB c ð1 þ c1 Þ
>
> ðTr ¼ 1Þ ¼ 1 critical temperature.
< dT 2 2
r
(I.5)
>
> d2 aSUP c ð1 þ c1 Þ þ 2c2
A priori, these assumptions do not guaranty for the function to
>
>
: ðTr ¼ 1Þ ¼ 1 be of class C2 and more generally, to pass the proposed consistency
dTr2 2
test for a-functions.
As a consequence, the Mathias-Copeman a-function is of class
C2 provided c2 is set equal to 0. Doing so, this a-function would  Is the a-function of class C2?
loose one of its three adjustable parameters.
The first derivatives of the Boston-Mathias a-function are:
8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Conclusion: >
> daSUB aSUB ðTr Þ
>
> if Tr  1 : ðTr Þ ¼ mðuÞ
At this step, there is no need to pursue the study: it is obvious da < dT T
r r
that the Mathias-Copeman a-function cannot pass the consistency ðTr Þ ¼ (J.2)
dTr >
>
test since the C2-class condition requires to annihilate coefficient c2 >
> da
: if Tr > 1 : SUP ðTr Þ ¼ LgTrg1 aSUP ðTr Þ
and consequently to modify the a-function expression. dTr
In addition, it is remarkable that the Mathias-Copeman a-
function is similar to the Soave72 a-function for Tr > 1 except that At Tr ¼ 1, one has:
the universal mðuÞ parameter of the Soave function is replaced by 8
an adjustable parameter c1. >
> daSUB
If, as often seen in practice, the parameter c1 is set to mðuÞ, the < dT ðTr ¼ 1Þ ¼ mðuÞ
>
r
Mathias-Copeman a-function would loose a second adjustable (J.3)
>
> da mðuÞ
parameter (c1 ). It would result that the Mathias-Copeman a-func- >
: SUP
ðTr ¼ 1Þ ¼ Lg ¼  g ¼ mðuÞ
dTr g
tion would only contain one adjustable parameter (c3 ).
Furthermore, similarly to the Soave72 a-function, the Mathias- This proves that the Boston-Mathias a-function is at least of
Copeman a-function exhibits an unphysical behavior at infinite class C1.
temperature: The second derivatives of the Boston-Mathias a-function are:
536 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

The first derivatives of the Mathias a-function are:

8
>
> d2 aSUB 1 mðuÞ½mðuÞ þ 1 8
>
> if Tr  1 : ðTr Þ ¼ > daSUB pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mðuÞ
< dT 2 2 3=2 >
>if Tr 1: ðTr Þ¼2 aSUB ðTr Þ pffiffiffiffiffi þpð2Tr 1:7Þ
d2 a r T r da < dT 2 Tr
ðTr Þ ¼ ðTr Þ¼
r
dTr2 >
> 2   >
>
> d aSUP g 2 g dTr >
> da
: if Tr > 1 : ðTr Þ ¼ LgTr LgTr  g þ 1 aSUP ðTr Þ :if Tr >1: SUP ðTr Þ¼2LgTrg1 aSUP ðTr Þ
dTr2 dTr
(J.4) (K.2)
At Tr ¼ 1, one has:

8 2
>
> d aSUB mðuÞ½mðuÞ þ 1
>
> ðTr ¼ 1Þ ¼
< dT 2 2
r
(J.5)
>
> d2 aSUP mðuÞ ½mðuÞ2
>
>
: ðTr ¼ 1Þ ¼ Lg ½L g  g þ 1 ¼ mð uÞ mðu Þ  1  þ 1 ¼
dTr2 2 2

2 2
aSUB ðT ¼ 1Þ
Since d dT s d aSUP
ðTr ¼ 1Þ, the Boston-Mathias At Tr ¼ 1:
2 r dTr2
r
a-function is thus not of class C2. 8
>
> daSUB mðuÞ mðuÞ
>
< ðTr ¼ 1Þ ¼ 2 þ pð1:7  2Þ ¼ 2 þ 0:3p
Conclusion: dTr 2 2
At this step, there is no need to pursue the study. Not being of >
> daSUP
> mðuÞ
class C2, the Boston-Mathias a-function does not pass the proposed : ðTr ¼ 1Þ ¼ 2Lg ¼ 2 þ 0:3p
dTr 2
consistency test.
(K.3)
Appendix K: mathematical study of the Mathias a-function [28]

The Mathias a-function contains one universal parameter

This proves that the Mathias a-function is at least of class C1.
(mðuÞ) and one (component-dependent) adjustable parameter (p).
The second derivatives of the Mathias a-function are:
8 h  pffiffiffiffiffi i2
>if Tr 1:aSUB ðTr Þ¼ 1þmðuÞ 1 Tr pð1Tr Þð0:7Tr Þ 8
>
> >
> d2 aSUB mðuÞ 2
>
< >
> if T  1 : ðT Þ ¼ 2 p ffiffiffiffi
ffi þpð2T 1:7Þ
    >
>
r
dTr2
r
2 Tr
r
aðTr Þ¼ if Tr >1:aSUP ðTr Þ¼ exp L 1Trg 2 >
>
>
>
> >
> !
>
> < pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mðuÞ
: mðuÞ L1 d2 a
with: L¼1þ þ0:3p and g¼ ðTr Þ ¼ þ2 a SUB r ðT Þ 2p
2 L dTr2 >
> 3=2
4Tr
>
>
(K.1) >
>
>
>
>
> d2 aSUP g2  g 
The expressions of parameters L and g were derived by Mathias >
: if Tr >1 : ðTr Þ ¼ 2LgTr 2LgTr  g þ1 aSUP ðTr Þ
2
dTr
by assuming that:
(K.4)
pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi
- The first derivatives of aSUB and aSUP are equal at the critical
temperature,
pffiffiffiffiffiffiffiffiffiffi
- The second derivative of aSUP is equal to zero at the critical
temperature.
At Tr ¼ 1, one has:
A priori, these assumptions do not guaranty for the function to
be of class C2 and more generally, to pass the proposed consistency
test for a-functions.

 Is the a-function of class C2?

8
>
> d2 aSUB
>
> ðTr ¼ 1Þ ¼ 0:5m2 þ 0:6mp þ 0:18p2 þ 0:5m  4p
>
> dT 2
>
> r
>
>
< 2
d aSUP
ðTr ¼ 1Þ ¼ 2Lg½2Lg  g þ 1 (K.5)
>
> 2
> dTr
>
>
>
>
>
>
> 10m2 þ 12mp þ 5m3 þ 9m2 p þ 5:4mp2 þ 10m þ 3:6p2 þ 1:08p3 þ 6p
:¼
10 þ 5m þ 3p
 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538 537

d2 aSUB d2 aSUP
Since dTr2
ðTr ¼ 1Þ s dTr2
ðTr ¼ 1Þ, the Boston-Mathias a-
da
function is thus not of class C . 2 lim ðTr Þ ¼ X (L.8)
Tr /þ∞ dTr
Therefore, X has to be negative. This condition is contradictory
Conclusion: with condition (L.3). The reason of this contradiction is that the
At this step, there is no need to pursue the study. Not being of Gibbons and Laughton a-function exhibits a minimum when
class C2, the Mathias a-function does not pass the proposed con- X0
(this behavior is similar to the one of the Soave72 a-
sistency test. Y0
function). As a conclusion, the Gibbons and Laughton a-function
can only be decreasing on a limited temperature range ½0 ; Tr;min .
Appendix L: mathematical study of the Gibbons and Laughton a- As illustrated in Figure L.1, the minimum of the a-function is found
function [25] at

The a-function involved two adjustable parameters denoted X qffiffiffiffiffiffiffiffiffiffiffiffi Y

Tr;min ¼ (L.9)
and Y: 2X
pffiffiffiffiffi 
aðTr Þ ¼ 1 þ XðTr  1Þ þ Y Tr  1 (L.1)

 Is the a-function of class C2? Yes.

 Is the a-function always positive?

Eq. (L.1) can be alternatively written:

pffiffiffiffiffi 2 pffiffiffiffiffi
aðTr Þ ¼ X Tr þ Y Tr þ ð1  X  YÞ (L.2)

Which is a second-order polynomial of the square root of the

reduced temperature. Therefore, this polynomial is positive if [42]
8
<X  0
Y  1 p
Xffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (L.3)
:
Y  2 Xð1  X  YÞ

As will be shown below, the parameter Y is necessarily negative

to ensure that the a-function is decreasing at low temperature. The
third inequality of Eq. (L.3) can be squared, both sides of the
equation being negative, in order to solve this inequality with Fig. L.1. Changes of Tr,min with respect to the change of Y for different values of X. Note
respect to Y. The resolution of the inequality provides the following that the value of the a-function at Tr;min is not equal to zero.
existence interval for parameter Y   Y 2
a Tr;min ¼ þ1XY (L.10)
pffiffiffiffi pffiffiffiffi 4X
2X  2 X  Y  2X þ 2 X (L.4)
pffiffiffiffi
One can demonstrate that 2X þ 2 X  1  X. Therefore, pos-  Is the second derivative of the a-function always positive?
itivity condition of the a-function is, reminding that Y has to be
negative: The second derivative of the Gibbons and Laughton a-function
is:
X0 
pffiffiffiffi pffiffiffiffi  (L.5)
2X  2 X  Y  min 0; 2X þ 2 X d2 a 1 Y
ðTr Þ ¼  (L.11)
dTr2 4 T 3=2
r

As a consequence, d2 a ðT Þ is positive provided Y is negative.

dTr2 r

 Is the first derivative of the a-function always negative?

 Is the third derivative of the a-function always negative?

The first derivative of the Gibbons and Laughton a-function is:

The third derivative of the Gibbons and Laughton a-function is:
da Y
ðTr Þ ¼ X þ pffiffiffiffiffi (L.6) d3 a 3 Y
dTr 2 Tr ðTr Þ ¼ (L.12)
dTr3 8 T 5=2
r
When Tr tends to 0, the limit of Eq. (L.6) is:
As a consequence, d3 a ðT Þ is negative provided Y is negative.
dTr3 r
da ∞ if Y < 0
lim ðTr Þ ¼ (L.7)
Tr /0 dTr þ∞ if Y > 0
Conclusion:
Therefore, Y has to be negative The Gibbons and Laughton a-function behaves as a consistent a-
When Tr tends to þ∞, the limit of Eq. (L.6) is: function, up to a certain reduced temperature Tr,min defined by Eq.
538 Y. Le Guennec et al. / Fluid Phase Equilibria 427 (2016) 513e538

(L.9) if j.supflu.2016.05.012 (in press).

[24] G. Soave, Improvement of the van der Waals equation of state, Chem. Eng. Sci.
39 (1984) 357e369, http://dx.doi.org/10.1016/0009-2509(84)80034-2.
X0 
pffiffiffiffi pffiffiffiffi  (L.13)
[25] R.M. Gibbons, A.P. Laughton, An equation of state for polar and non-polar
2X  2 X  Y  min 0; 2X þ 2 X substances and mixtures, J. Chem. Soc. Faraday Trans. 2 (80) (1984) 1019,
http://dx.doi.org/10.1039/f29848001019.
[26] M.A. Trebble, P.R. Bishnoi, Development of a new four-parameter cubic
equation of state, Fluid Phase Equilib. 35 (1987) 1e18, http://dx.doi.org/
10.1016/0378-3812(87)80001-8.
References [27] G.A. Melhem, R. Saini, B.M. Goodwin, A modified Peng-Robinson equation of
state, Fluid Phase Equilib. 47 (1989) 189e237, http://dx.doi.org/10.1016/
0378-3812(89)80176-1.
[1] J.O. Valderrama, The state of the cubic equations of state, Ind. Eng. Chem. Res. [28] P.M. Mathias, A versatile phase equilibrium equation of state, Ind. Eng. Chem.
42 (2003) 1603e1618, http://dx.doi.org/10.1021/ie020447b. Process Des. Dev. 22 (1983) 385e391, http://dx.doi.org/10.1021/i200022a008.
[2] G.M. Kontogeorgis, E.C. Voutsas, I.V. Yakoumis, D.P. Tassios, An equation of [29] J.L. Daridon, B. Lagourette, H. Saint-Guirons, P. Xans, A cubic equation of state
state for associating fluids, Ind. Eng. Chem. Res. 35 (1996) 4310e4318. model for phase equilibrium calculation of alkane þ carbon dioxide þ water
[3] G. Soave, Equilibrium constants from a modified Redlich-Kwong equation of using a group contribution kij, Fluid Phase Equilib. 91 (1993) 31e54, http://
state, Chem. Eng. Sci. 27 (1972) 1197e1203. dx.doi.org/10.1016/0378-3812(93)85077-Y.
[4] M.S. Wertheim, Fluids with highly directional attractive forces. I. Statistical [30] U.K. Deiters, Comments on the modeling of hydrogen and hydrogen-
thermodynamics, J. Stat. Phys. 35 (1984) 19e34. containing mixtures with cubic equations of state, Fluid Phase Equilib. 352
[5] J.D. Van der Waals, Over de continuiteit van den gas- en vloeistoftoestand, (2013) 93e96, http://dx.doi.org/10.1016/j.fluid.2013.05.032.
1873. [31] C. Colina, J. Santos, C. Olivera-Fuentes, High-temperature Behaviour of the
[6] J.V. Sengers, R.F. Kayser, C.J. Peters, H.J. White, Equations of State for Fluids and Cohesion Parameter of Cubic Equations of State 29, High Temp.-High Press,
Fluid Mixtures, first ed., Elsevier, Amsterdam ; New York, 2000. 1997, pp. 525e532, http://dx.doi.org/10.1068/htec462.
[7] R. Clausius, Über das Verhalten der Kohlens€ aure in Bezug auf Druck, Volumen [32] M.M. Abbott, J.M. Prausnitz, Generalized van der waals theory: a classical
und Temperatur, Ann. Phys. 245 (1880) 337e357. perspective, Fluid Phase Equilib. 37 (1987) 29e62, http://dx.doi.org/10.1016/
[8] O. Redlich, J.N.S. Kwong, On the thermodynamics of solutions; an equation of 0378-3812(87)80042-0.
state; fugacities of gaseous solutions, Chem. Rev. 44 (1949) 233e244. [33] S.I. Sandler, The generalized Van der Waals partition function as a basis for
[9] V. Kalikhman, D. Kost, I. Polishuk, About the physical validity of attaching the excess free energy models, J. Supercrit. Fluids 55 (2010) 496e502, http://
repulsive terms of analytical EOS models by temperature dependencies, Fluid dx.doi.org/10.1016/j.supflu.2010.10.014.
Phase Equilib. 293 (2010) 164e167, http://dx.doi.org/10.1016/ [34] R. Privat, J.-N. Jaubert, Classification of global fluid-phase equilibrium be-
j.fluid.2010.03.003. haviours in binary systems, Chem. Eng. Res. Des. 91 (2013) 1807e1839, http://
[10] P.M. Mathias, T.W. Copeman, Extension of the Peng-Robinson equation of dx.doi.org/10.1016/j.cherd.2013.06.026.
state to complex mixtures: evaluation of the various forms of the local [35] S. Lasala, P. Chiesa, R. Privat, J.-N. Jaubert, VLE properties of CO2 e based bi-
composition concept, Fluid Phase Equilib. 13 (1983) 91e108, http:// nary systems containing N2, O2 and Ar: experimental measurements and
dx.doi.org/10.1016/0378-3812(83)80084-3. modelling results with advanced cubic equations of state, Fluid Phase Equilib.
[11] R. Stryjek, J.H. Vera, PRSV: an improved PengdRobinson equation of state for (2016), http://dx.doi.org/10.1016/j.fluid.2016.05.015 (in press).
pure compounds and mixtures, Can. J. Chem. Eng. 64 (1986) 323e333, http:// [36] M.-J. Huron, J. Vidal, New mixing rules in simple equations of state for rep-
dx.doi.org/10.1002/cjce.5450640224. resenting vapour-liquid equilibria of strongly non-ideal mixtures, Fluid Phase
[12] B.E. Poling, J.M. Prausnitz, J.P. O'Connell, The Properties of Gases and Liquids, Equilib. 3 (1979) 255e271, http://dx.doi.org/10.1016/0378-3812(79)80001-1.
McGraw-Hill, New York, 2001 (accessed July 22, 2015), http:// [37] G.M. Wilson, Vapor-liquid equilibrium. XI. A new expression for the excess
accessengineeringlibrary.com/browse/properties-of-gases-and-liquids-fifth- free energy of mixing, J. Am. Chem. Soc. 86 (1964) 127e130, http://dx.doi.org/
edition. 10.1021/ja01056a002.
[13] G. Heyen, A cubic equation of state with extended range of application, Chem. [38] R. Span, W. Wagner, A new equation of state for carbon dioxide covering the
Eng. Thermodyn., SA Newman1981, Proceedings of the 2nd world congress of fluid region from the triple-point temperature to 1100 K at Pressures up to
chemical engineering, Montreal, Canada, pages 4e9. 800 MPa, J. Phys. Chem. Ref. Data 25 (1996) 1509, http://dx.doi.org/10.1063/
[14] C.H. Twu, A modified redlich-kwong equation of state for highly polar, su- 1.555991.
percritical systems, Proc. Int. Symp. Thermodyn. Chem. Eng. Ind. (1988) 148. [39] C. Tegeler, R. Span, W. Wagner, A new equation of state for argon covering the
[15] C.H. Twu, D. Bluck, J.R. Cunningham, J.E. Coon, A cubic equation of state with a fluid region for temperatures from the melting line to 700 K at pressures up to
new alpha function and a new mixing rule, Fluid Phase Equilib. 69 (1991) 1000 MPa, J. Phys. Chem. Ref. Data 28 (1999) 779, http://dx.doi.org/10.1063/
33e50, http://dx.doi.org/10.1016/0378-3812(91)90024-2. 1.556037.
[16] C.H. Twu, J.E. Coon, J.R. Cunningham, A new generalized alpha function for a [40] K.A.M. Gasem, W. Gao, Z. Pan, R.L. Robinson, A modified temperature
cubic equation of state Part 1. Peng-Robinson equation, Fluid Phase Equilib. dependence for the PengeRobinson equation of state, Fluid Phase Equilib. 181
105 (1995) 49e59, http://dx.doi.org/10.1016/0378-3812(94)02601-V. (2001) 113e125, http://dx.doi.org/10.1016/S0378-3812(01)00488-5.
[17] C.H. Twu, J.E. Coon, J.R. Cunningham, A new generalized alpha function for a [41] Y. Chen, F. Mutelet, J.-N. Jaubert, Modeling the solubility of carbon dioxide in
cubic equation of state Part 2. Redlich-Kwong equation, Fluid Phase Equilib. imidazolium-based ionic liquids with the PC-SAFT equation of state, J. Phys.
105 (1995) 61e69, http://dx.doi.org/10.1016/0378-3812(94)02602-W. Chem. B 116 (2012) 14375e14388, http://dx.doi.org/10.1021/jp309944t.
[18] C. Coquelet, A. Chapoy, D. Richon, Development of a new alpha function for [42] J.W. Schmidt, W. Heß, Positivity of cubic polynomials on intervals and positive
the PengeRobinson equation of state: comparative study of alpha function spline interpolation, BIT 28 (1988) 340e352, http://dx.doi.org/10.1007/
models for pure gases (natural gas components) and water-gas systems, Int. J. BF01934097.
Thermophys. 25 (2004) 133e158, http://dx.doi.org/10.1023/B: [43] J.-N. Jaubert, F. Mutelet, VLE predictions with the PengeRobinson equation of
JOT.0000022331.46865.2f. state and temperature dependent kij calculated through a group contribution
[19] J.F. Boston, P.M. Mathias, Phase equilibria in a third-generation process sim- method, Fluid Phase Equilib. 224 (2004) 285e304, http://dx.doi.org/10.1016/
ulato, in: Proc. 2nd Int. Conf. Phase Equilib. Fluid Prop. Chem. Process Ind., j.fluid.2004.06.059.
West Berlin, 1980, pp. 17e21. [44] S. Vitu, R. Privat, J.-N. Jaubert, F. Mutelet, Predicting the phase equilibria of
[20] E. Neau, O. Herna ndez-Garduza, J. Escandell, C. Nicolas, I. Raspo, The Soave,
CO2þhydrocarbon systems with the PPR78 model (PR EOS and kij calculated
Twu and BostoneMathias alpha functions in cubic equations of state, Fluid through a group contribution method), J. Supercrit. Fluids 45 (2008) 1e26,
Phase Equilib. 276 (2009) 87e93, http://dx.doi.org/10.1016/ http://dx.doi.org/10.1016/j.supflu.2007.11.015.
j.fluid.2008.09.023. [45] J.-N. Jaubert, R. Privat, F. Mutelet, Predicting the phase equilibria of synthetic
[21] E. Neau, I. Raspo, J. Escandell, C. Nicolas, O. Hern andez-Garduza, The Soave, petroleum fluids with the PPR78 approach, AIChE J. 56 (2010) 3225e3235,
Twu and BostoneMathias alpha functions in cubic equations of state. Part II. http://dx.doi.org/10.1002/aic.12232.
Modeling of thermodynamic properties of pure compounds, Fluid Phase [46] J.-N. Jaubert, R. Privat, Relationship between the binary interaction parame-
Equilib. 276 (2009) 156e164, http://dx.doi.org/10.1016/j.fluid.2008.10.010. ters (kij) of the PengeRobinson and those of the SoaveeRedlicheKwong
[22] P. Mahmoodi, M. Sedigh, Soave alpha function at supercritical temperatures, equations of state: application to the definition of the PR2SRK model, Fluid
J. Supercrit. Fluids 112 (2016) 22e36, http://dx.doi.org/10.1016/ Phase Equilib. 295 (2010) 26e37, http://dx.doi.org/10.1016/
j.supflu.2016.01.004. j.fluid.2010.03.037.
[23] P. Mahmoodi, M. Sedigh, Second derivative of alpha functions in cubic
equations of state, J. Supercrit. Fluids (2016), http://dx.doi.org/10.1016/