- Email: [email protected]
Matching equation of state mixing rules to activity coefficient model expressions
[LIIIDPHAS[ EQUILIBRIA ELSEVIER
Fluid Phase Equilibria 121 (1996) 15 26
Matching equation of state mixing rules to activity coefficient model expressions Michael L. M i c h e l s e n lnstitut for Kemiteknik, Bygning 229, DTU, DK 2800 Lyngby, Denmark Received 28 September 1995; accepted t6 January 1996
Abstract
Procedures for deriving equation of state mixing rules based on matching the zero pressure excess Gibbs energy to that of a previously parameterized activity coefficient model are discussed. I have investigated the conditions required for accurate reproduction of the original model and provide a simple graphical tool for evaluating the error of the equation of state mixing rule. An extension of the MHV-2 mixing rule to incorporate mixtures containing components with very low volatility is suggested.
1. I n t r o d u c t i o n
Equation of state mixing rules capable of incorporation excess Gibbs energy model expressions without reparameterization have found widespread use (Heideman and Khokal, 1990; Michelsen, 1990a,b; Holderbaum and Gmehling, 1991; Soave, 1992; Wong and Sandler, 1992; Boukouvalas et al., 1994). Recently, the ability of such mixing rules to reproduce the behaviour of the underlying gE model has been questioned (Kalospiros et al., 1995). In the present article we analyze the mixing rules derived from Michelsen's exact model (Michelsen, 1990a). The aim is to clarify which mixtures they are suitable for, and to provide a simple tool for estimating their accuracy. Finally, a simple modification to the MHV-2 mixing rule (Dahl and Michelsen, 1990), which extends its range of applicability, is suggested.
2. Matching gE at z e r o pressure For the Redlich-Kwong equation of state p=
RT
a
v-b
0378-3812/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved. PH S 0 3 7 8 - 3 8 1 2 ( 9 6 ) 0 3 0 2 0 - 8
(l)
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26
16
Michelsen (1990b) derived the following expression for the liquid phase excess Gibbs energy at zero pressure: gE)
b + '~'~zi I n - - = q(c~) Eos.
P=0
Eziq(o¢-ii)
i
(2)
i
Here, the mixture parameter c~ = a / ( b R T ) , and the function q is given by u0+ 1 q ( a ) = - 1 - ln(u 0 - 1) - o~ l n - -
(3)
/A0
where u 0 is the reduced liquid phase volume at zero pressure
u o = ( v / b ) e = o = ~'( c ~ - l - ( o L 2 - 6 o ~ + 1 ) '/ 2)
(4)
The exact method of Michelsen (1990a) utilizes Eq. (2) as a density independent mixing rule for the mixture parameter o~. Replacing the excess Gibbs energy by the term calculated from an activity coefficient model results in
b(g )
q(o~) = ~ z i q ( e ~ i i ) + ] ~ z i l n - i
i
(5)
+
bii
"~
AM
and the use of this mixing rule for o~ results in an exact match between the excess Gibbs energy of the equation of state at zero pressure and that of the underlying activity coefficient model, regardless of the mixing rule used for the b parameter. The excess volume using the conventional linear mixing rule for b is very modest, and pressure effects on the EOS excess Gibbs energy therefore are of no importance at low to moderate pressures. Application of the exact mixing rule requires that the value of o~ for the pure components as well as for the mixture exceeds 5.83. To enable extrapolation to elevated temperatures Michelsen (1990a) suggested to select a limiting value of ~, o~ = a lim, where c~lim is chosen well above 5.83, and to replace q by a relation of the form -- Sli m =
g l ( q -- q(ot lira)) -+- g z ( q
-- q(~
lim)) 2
(6)
when c~ is less than the limiting value. The constants K~ and K 2 are selected such that the resulting q(o~) has continuous first and second derivatives at C~l~m. A variety of alternatives have later been proposed (Michelsen, 1990b; Holderbaum and Gmehling, 1991; Soave, 1992). They all replace q ( a ) by an approximating function q * ( a ) which is subsequently used for all values of e~. Michelsen (1990b) proposed the linear approximation MHV-I: q* (o~) = q0 + q , ~ , q, = - 0 . 5 9 3
(7)
and the quadratic approximation, MHV-2: q* (c~) = qo + q,c~ + q2oL2, q, = - 0 . 4 7 8 , q2 = - 0 . 0 0 4 7
(8)
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
17
Holderbaum and Gmehling used MHV-1 with a different value of ql ( - 0 . 6 4 3 ) and Soave used a somewhat more complex functional form: OL *(or) = 3.365 - - - l n 2 2
2 1 n 2 ) 2 + 17.2511/2
(9)
With all these mixing rules the EOS fugacity coefficients are determined from ln¢i = In p ( v Z b )
+
(
'
v-b
°
" + b J/--ffC
Jw
,j + I n
v
1--~ni ]
wn
(10)
where the composition derivative of e~ is calculated from
(dq)(tno) ) -~
Oni
c~ = q * ( et ii ) -- q * ( ot ) + l n "y A M + l n b i---~i+ -~- - - 1
(11)
All these variants require virtually identical computational effort. The exact model does not permit explicit solution of Eq. (5) for or, but numerical solution by Newton's method is straightforward as q ( a ) is monotonically decreasing, with derivatives given by: dq --
dc~
u0 + 1 -
I n - - ,
u0
d2q de~2
1 -
e~(2Uo + 1 - e ~ )
(12)
and the additional computational effort in this step is negligible.
3. R a n g e
of applicability
The EOS excess Gibbs energy calculated with the approximate methods all deviate to some extent from that of the activity coefficient model, the errors associated with each method depending on the e~ range in which they are used. As the purpose of these methods is to utilize g E correlations for data measured at conditions, where such correlations have been used for representing vapour-liquid equilibrium, i.e, data measured around atmospheric pressure or lower, we shall first investigate the relation between the et parameter and the saturation pressure. The pure component liquid fugacity at zero pressure is given by In RT)
= q(e~)
(13)
where, for simplicity, we have omitted the component index. At saturation, the fugacities in the liquid and the vapour phase are identical, i.e. fl(Psat) = L ( P s a t ) =
Psat~v(Psat)
(14)
18
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
10 -L
~
\ 10
-2
(73 (/3 &) L10
-3
\
RD ©
~
1 0 -'-
\
(1) Ct2 10
-5
\
5
110
1'5 0<
2~D
\ 2~5
Fig. 1. Relationship between the pure component o~ parameter and the reduced saturation pressure for a component with an acentric factor of 0.4. If we neglect the pressure dependence of the liquid phase fugacity and assume the vapour phase to be ideal, Eq. (14) yields, and Eq. (13) can be written
RTl =ln - - ~ ] =lnlBb +lnP~.sat--lnT~=q(c~ )
In
(15)
where ~-~b = 0.08664. In this expression the relation between Tr and o~ is component specific, but the variation of In T~ with o~ is modest, the reduced temperature typically decreasing from about 0.7 to about 0.4 when OL increases from 10 to 20. An approximate relation between T~ and o~ can be obtained from Soave's expression (Soave, 1972) for the temperature dependence of the a parameter a - [1 + m(co)(1 - T°5)] 2, m ( o ~ ) = 0 . 4 8 + 1 . 5 7 4 o ) - 0 . 1 7 4 ( o 2 (16) acrit or
O~
(1 + m ( c o ) ( 1 -
Tr°'5)) 2
(17) (1 crit
Tr
where o~cr~t= ~ a / ~ ' ~ b = 4.73. Substitution of this relation between the reduced temperature and o~ into Eq. (15) yields the desired relation between c~ and the reduced saturation pressure, as shown in Fig. 1 for a representative eo value of 0.40. A saturation pressure of 1 bar corresponds to reduced pressures in the range 0.005-0.04, the lower limit applying for water. From Fig. 1 w e notice that the corresponding oL range is 10-13, and the coefficients used in MHV-1 and MHV-2 are selected to provide particularly accurate approximations in this range. The reduced pressure decreases very rapidly with increasing values of OL and is 10 4 at o~= 19.
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26
19
4. Error estimate for approximate mixing rules
The error in the EOS Gibbs free energy, i.e. the difference between the g Z calculated from the equation of state and that of the activity coefficient model, can be determined by subtracting Eq. (5) (with q * replacing q) from Eq. (2): AgE
(gE)
R-T =
~
where e ( a ) = q(a and write
(gE) EOS-- ~
AM = e ( c ¢ ) -
~zie(c~ii) i
(18)
- q * (¢x). This expression can be simplified as follows: We define fi = ~i zi a ii
]
~--~ = e ( N ) - z i e ( ~ i i
) +[e(a)-e(N)]=e,+e
2
(19)
The first term, e~, is not related to excess properties but depends only on the pure component oL parameters. It is identically zero when these are equal. The second term can be approximated by a Taylor series expansion from e~ = (de/ (a-~)-e(~)=( e(c¢) - e ( ~ ) = e ( ~ ) + I dcx ]~
de] (cx-~) I dot ]u
(20)
For the present purpose an adequate estimate of a = ~ can be obtained from the MHV-1 approximation, a=
1 1 Y'~ziaii+--r(z)=~+--r(z) i q= ql
(21)
where qj = - 0 . 5 9 3 and
r(z):
~
+ ~ziln AM
i
bii
We thus arrive at l(de)
e2=--
r(z)
(22)
and this term is therefore proportional to the activity coefficient model excess Gibbs energy, modified by the b term. Fig. 2 shows e ( a ) and d e / d R for the MHV-2 approximation, q* = q0 + qta + q2 a2. The value used for qo is of no importance for the analysis, and it is here chosen to match q at a = 12. The e~ contribution can be determined by a simple geometrical construction. Consider two components with pure component oL values of 16 (point B) and 11 (point A), respectively. The line connecting A and B corresponds to z~ e(al~) + z2e(a22) and e~ is equal to the vertical distance from the line AB to the e(e~) curve, as marked by the arrow at zl = 0.7 (~ = 14.5) It is evident that e~ is particularly small when the a values for the pure components are close to the value at the inflection point (~ = 12), and a very modest error contribution results when both
20
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26 0.03
e end de/dc~ for MHV-2
/ / /
0.02
\~
\\
---
0
L_ L
/ /
de/do~
\
0.01
b_J
e
///
//////
0.00
-0.01
-0.02
i
8
I0
l~2
14
IJ6
18
2~0
G
Fig. 2. Magnitude of e and d e / d ~ for the MHV-2 mixing rule.
values are located in the interval 10 to 15, corresponding roughly to pure component partial pressures in the range 0.1-1 bar. Inside this range the magnitude of e 2 is also very modest, the factor of r(z) being below 0.01 in the entire range 9.5 < a < 16. Acceptable results from MHV-2 can thus be expected for mixtures with pure component partial pressures in the range form about 0.01 bar to a few bar. Fig. 2 also illustrates the dramatic increase in the magnitude of e~ for mixtures containing a light ((~ = 10-12) and a very heavy (c~ > 18) component, and in practice the use of MHV-2 is not advisable outside the range 8-18. Activity coefficients at 298.15 K and 1 bar for the mixture hexane-hexadecane, investigated earlier by Kalospiros et al. (1995), are shown in Fig. 3 using UNIFAC as the excess Gibbs energy model. The pure component c~ values are 12 and 28, respectively, and the behaviour of MHV2 is clearly unsatisfactory. Fig. 4 shows the e and d e / d ( ~ for Soave's approximation, Eq. (9). This expression has the advantage that the correct value for the derivative, dq/da = - I n 2, is obtained in the limit a --* oo. We notice that the Soave expression provides a remarkably good match in the entire range 9 < e~ < 20, being clearly superior to MHV-2 at large a values. Among the mixing rules investigated here the Soave approximation appears to be the superior with respect to accuracy as well as c~ range. When the mixture contains components with small (10-12) as well as large (above 20) a values, the magnitude of e~ can however become quite substantial. The actual difference between the EOS excess Gibbs energy and that of the activity coefficient model is compared with that predicted by Eq. (19) for the acetone/water mixture at 340 K, 1 bar. The UNIQUAC equation with parameters estimated in the Dechema series (Vol I, 1, p.237) from the data of Kojima et al. are used for the excess Gibbs energy calculation. The predicted error is in close agreement with the calculation (Fig. 5), and the the RMS deviation between the reduced excess Gibbs energy calculated from the equation of state and the activity coefficient is 0.0022. For the Soave approximation the RMS error in reproducing the UNIQUAC excess Gibbs energy has an RMS value
21
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26 1.0
Hexone-Hexodecane at 298 K 0.5
---
©
UNIFAC MHV-2
0.0 C:3 O> c- - O 5
/
l
\
\
/ / -I
/
.0
/
\
/
\ \
/
\ \ \
-1.5
0.'2 0.'4 0.'6 0 '8 1.'0 H e x o n e mole froction Fig. 3. Activity coefficients in the mixture hexane-hexadecane at 298 K, 1 bar, calculated with MHV-2 and with UNIFAC. 0.0
of only 0.0015. The difference in the pure component c~ values (10.1 and 15.2, respectively) for this mixture is quite substantial. The error curves for the linear MHV-1 and Holderbaum-Gmehling ( H / G ) approximations are shown in Fig. 6. It is readily shown that e~ is only related to the curvature of the error, and e~ will therefore be the same for the two methods. It is also clear that e~ is always positive, and increasing in magnitude with increasing difference between the pure component ct values. In the o~ range 10-15,
0.0.3 e
and
002
de/dG
for
--
e
- - -
l~'e'~ ~
Soove
001 k_ 0 L k_ I,I 0.00
-0.01
-0.02
8
110
Ii2
114
I~6
118
2~0
O~
Fig. 4. Magnitude of e and de/dot for the Soave approximation.
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
22
0.008 - -
Actual Estim.
error error
0006 Od
/
LLI 0.004 C7~
~_ 0
f
--
~ \
/
\
0 002
tm
0.000
-0.002
O. 0
0.~2 0.~4 0.'6 0.'8 Acetone mole fraction
1 .'0
Fig. 5. Comparison of error estimates and actual errors for the MHV-2 mixing rule applied to a mixture of acetone and water at 340 K, 1 bar.
the error in e 2 for MHV-1 is substantially smaller than that obtained with the Holderbaum-Gmehling coefficient, as illustrated in Fig. 7 for the acetone/water mixture. The factor of r(x) in the Holderbaum-Gmehling approximation is however negative, and as the excess Gibbs energy for most mixtures is positive the error contributions from e I and e 2 in this approach will partly cancel. As a result the RMS error for the Holderbaum-Gmehling approximation (0.014) is actually lower than that
0.15 e and
de/dc<
for
0.10
MHV-1
and
H/O
e "'""
0.05
\
""-...
---
-...
- ....
de/dc~,
de/dcL
MHV-1
1'4 O~
1'6
1'8
H/O
\
k_
0.00 k._
LJ
-0.05
-0.10
-0.15
8
1~0
1'2
210
Fig. 6. Magnitude of e and d e/do~ for the MHV-1 and the Holderbaum-Gmehling mixing rules.
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26 0 050
23
Acetone-Woter, 540 K using MHV-1 and H/G -
-
el
-- - 0025
. . . . . . .
e 2, M H V - 1 e2, H / G
Od DJ C~
0.000
© k_ LJ
- 0 025
-0050
- ~
09
02
0.'4
0.16
0.8
10
Acetone mole froction Fig. 7. Error contributions for the MHV-I and the Holderbaum-Gmehling mixing rules applied to a mixture of acetone and water at 340 K, 1 bar.
of MHV-1 (0.022). As expected, the errors in both approximations are substantially larger than that of MHV-2 and the Soave approximation.
5. Extending the range The majority of the data published in the Dechema series (Gmehling and Onken, 1977) correspond to pure component a parameters in the range 9 - 1 7 , the most important exception being data at a very low pressures. An ability to be able to incorporate components with high et values could also be of importance in connection with representation of liquid-liquid equilibria or for use with polymer mixtures (Kalospiros et al., 1995). Finally very large a values may occur when a generalized correlation such as UNIFAC is used as the excess Gibbs energy model (Kalospiros et al., 1995), but such cases are probably less important since the excess Gibbs energy model is likely to be inaccurate when an extreme volatility range has to be covered. The hexane-hexadecane mixture analyzed by Kalospiros and Tassios, 1995 is e.g. characterized by a difference of 5 orders of magnitude in the pure component saturation pressures. In order to enable the MHV-2 approximation to handle such mixtures our suggestion is to combine MHV-2 with the exact model of Michelsen. Thus we use for the approximating function q* (et) = q ( a ) q*(o~) =qo+qlot+q2ot
for ot _> O~lim 2
(23)
for ~ < ~li m
At tx = a lira we require continuity in q * and its first two derivatives. Finally, the limiting value is chosen such that the coefficients q~ and q2 match those of MHV-2 as well as possible. A suitable choice is et lira 12.3, where we obtain q'(a) -- - 0.5986 and q"(et) = - 0.00918, the corresponding derivatives using the conventional MHV-2 constants being - 0 . 5 9 4 6 and - 0 . 0 0 9 4 , respectively. The =
M.L. Micbelsen / Fluid Phase Equilibria 121 (1996) 15-26
24
Hexane-Hexadecane
at 550 K
0.2
-
<3
0.0
-
UNIFAC
- -
Modif.
MHV-2
E E o cr~
-0.2
c
-O.4
-0.6
O.0
0.12 0.14 0 J6 0.18 Hexane mole fraction
1.10
Fig. 8. Activity coefficients in the mixture hexane-hexadecane at 350 K, 1 bar, calculated with the modified MHV-2 mixing rule and with UNIFAC. close match implies that correlations based on MHV-2, such as the gas solubility tables of Dalai et al. (1991), can be used with unchanged parameters. When the mixture ct parameter exceeds oLlira in the entire composition range the zero pressure excess Gibbs energy is equal to that of the corresponding activity coefficient model. A situation of particular interest is that where one component (1) has an c~ value lower than ct lira while the c~ value for the mixture and for the other component (2) is larger than ct lim" For such a mixture we have gE A--
= e(ct)
RT
-
Ezie(oLii
) =
-
Zle(Otll )
(24)
i
since e(o~2) = e(o~) = 0. The error in the activity coefficient for component 1 therefore is 0 Aln~, 1 ---- -ff~-(--ne(ct it)) = --e(t~ll )
(25)
IJft i.
and is thus given by the error in the pure component q value. For the hexane-hexadecane mixture at 298 K, the ct value of the light component is very close to the limiting value, and essentially exact reproduction of the activity coefficient model is obtained. At 350 K, the pure component a values are 9.4 and 21.8, respectively, and Fig. 8 shows that the activity coefficient of the light component is now well represented at all composition. The infinite dilution activity coefficient of hexadecane in hexane is somewhat in error, but this deviation is of no importance for the calculation of vapour-liquid equilibrium since the vapour pressure of the heavy component is extremely small. For the acetone-water mixture the accuracy of the modified MHV-2 mixing rule, with an RMS error of 0.001 l, is better than that of the Soave approximation. It should, however, be kept in mind that the accuracies of MHV-2, the Soave approximation and the modified MHV-2 mixing rule in
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
25
reproducing the activity coefficient model are all well below the accuracy of the activity coefficient model in representing the original data. Improvements of any of the 3 methods for this mixture therefore serves no practical purpose.
6. Conclusion A simple and illustrative graphical tool to estimate the error associated with a number of widely used approximations to the exact method of Michelsen for matching the zero pressure excess Gibbs energy of an equation of state to that of an activity coefficient model has been developed. In addition we have provided a modification of the MHV-2 mixing rules that permits accurate handling of mixtures containing components of widely differing volatilities.
7. List of symbols Mixture parameter, Redlich-Kwong equation Mixture parameter, Redlich-Kwong Equation Difference between approximate and exact q function e Contributions to e, Eq. (19) e I , e2 Fugacity f Excess Gibbs energy gE Constants in extrapolation method Kl, K2 Function of acentric factor in SRK equation, Eq. (16) m Pressure (bar) P Reduced pressure q Function defined by Eq. (3) q* Approximation to q qo, ql, q2 Coefficients in polynomial approximation to q Gas constant R Temperature (K) T Reduced temperature bl Reduced molar volume, u = v / b Molar volume U Mixture mole fraction Z a
b
7.1. Subscripts 0 AM EOS lim i ii
Zero pressure value Activity coefficient model property Equation of state property Limit for eL value in extrapolation method Component index Pure component index
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
26 7.2. G r e e k
(l)
f~A, f~B q~
letters
Mixture parameter, a = a/(bRT) Acentric factor Constants of Redlich-Kwong equation Fugacity coefficient
References Boukouvalas, C., Spilitois, N., Coutsikos, Ph., Tzouvaras, N. and l'assios, D., 1994. Prediction of vapor-liquid equilibrium with the LCVM model: a linear combination of the Vidal and Michelsen mixing rules coupled with the original UNIFAC and the t-mPR equation of state. Fluid Phase Equilibria, 92: 75. Dahl, S. and Michelsen, M.L., 1990. High-pressure vapor-liquid equilibrium with a UNIFAC-based equation of state. AIChE. J., 36: 1829. Dahl, S., Fredenslund, Aa, Rasmussen, P., 1991. The MHV-2 model: a UNIFAC-based equation of state model for prediction of gas solubility and vapor-liquid equilibrium at low and high pressures, Ind. Eng. Chem. Res., 30: 1936. Gmehling, J., Onken, U., 1977. Vapor-liquid Equilibrium Data Collection. Chemistry Data Series, Vol !. Part 1. DECHEMA, Frankfurt. Heideman, R.A. and Khokal, S.L., 1990. Combined excess free energy models and equations of state. Fluid Phase Equilibria, 56: 17. Holderbaum, T, Gmehling, J., 1991. A group contribution equation of state based on UNIFAC, Fluid Phase Equilibria, 70: 251. Kalospiros, N.S., Tzouvaras, N., Coutsikos, P.. Tassios, D.P., 1995. Analysis of zero-reference-pressure EoS/G E models, AIChE.J., 41: 928. Kalospiros, N.S., Tassios, D., 1995. Prediction of vapor-liquid equilibria in polymer solutions using an equation of state/excess Gibbs free energy model, Ind. Eng. Chem. Res., 34:2117. Michelsen, M.L., 1990a. A method for incorporating excess Gibbs energy models in equations of state, Fluid Phase Equilibria, 60: 47. Michelsen, M.L., 1990b. A modified Huron-Vidal mixing rule for cubic equations of state, Fluid Phase Equilibria, 60: 213. Soave, G., 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci., 27: 1197. Soave, G., 1992. A new expression of q(a) for the modified Huron-Vidal methods. Fluid Phase Equilibria, 72: 325. Wong, D.S.H. and Sandler, S.I., 1992. A theoretically correct mixing rule for cubic equations of state. AIChE. J., 38: 671.
Fluid Phase Equilibria 121 (1996) 15 26
Matching equation of state mixing rules to activity coefficient model expressions Michael L. M i c h e l s e n lnstitut for Kemiteknik, Bygning 229, DTU, DK 2800 Lyngby, Denmark Received 28 September 1995; accepted t6 January 1996
Abstract
Procedures for deriving equation of state mixing rules based on matching the zero pressure excess Gibbs energy to that of a previously parameterized activity coefficient model are discussed. I have investigated the conditions required for accurate reproduction of the original model and provide a simple graphical tool for evaluating the error of the equation of state mixing rule. An extension of the MHV-2 mixing rule to incorporate mixtures containing components with very low volatility is suggested.
1. I n t r o d u c t i o n
Equation of state mixing rules capable of incorporation excess Gibbs energy model expressions without reparameterization have found widespread use (Heideman and Khokal, 1990; Michelsen, 1990a,b; Holderbaum and Gmehling, 1991; Soave, 1992; Wong and Sandler, 1992; Boukouvalas et al., 1994). Recently, the ability of such mixing rules to reproduce the behaviour of the underlying gE model has been questioned (Kalospiros et al., 1995). In the present article we analyze the mixing rules derived from Michelsen's exact model (Michelsen, 1990a). The aim is to clarify which mixtures they are suitable for, and to provide a simple tool for estimating their accuracy. Finally, a simple modification to the MHV-2 mixing rule (Dahl and Michelsen, 1990), which extends its range of applicability, is suggested.
2. Matching gE at z e r o pressure For the Redlich-Kwong equation of state p=
RT
a
v-b
0378-3812/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved. PH S 0 3 7 8 - 3 8 1 2 ( 9 6 ) 0 3 0 2 0 - 8
(l)
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26
16
Michelsen (1990b) derived the following expression for the liquid phase excess Gibbs energy at zero pressure: gE)
b + '~'~zi I n - - = q(c~) Eos.
P=0
Eziq(o¢-ii)
i
(2)
i
Here, the mixture parameter c~ = a / ( b R T ) , and the function q is given by u0+ 1 q ( a ) = - 1 - ln(u 0 - 1) - o~ l n - -
(3)
/A0
where u 0 is the reduced liquid phase volume at zero pressure
u o = ( v / b ) e = o = ~'( c ~ - l - ( o L 2 - 6 o ~ + 1 ) '/ 2)
(4)
The exact method of Michelsen (1990a) utilizes Eq. (2) as a density independent mixing rule for the mixture parameter o~. Replacing the excess Gibbs energy by the term calculated from an activity coefficient model results in
b(g )
q(o~) = ~ z i q ( e ~ i i ) + ] ~ z i l n - i
i
(5)
+
bii
"~
AM
and the use of this mixing rule for o~ results in an exact match between the excess Gibbs energy of the equation of state at zero pressure and that of the underlying activity coefficient model, regardless of the mixing rule used for the b parameter. The excess volume using the conventional linear mixing rule for b is very modest, and pressure effects on the EOS excess Gibbs energy therefore are of no importance at low to moderate pressures. Application of the exact mixing rule requires that the value of o~ for the pure components as well as for the mixture exceeds 5.83. To enable extrapolation to elevated temperatures Michelsen (1990a) suggested to select a limiting value of ~, o~ = a lim, where c~lim is chosen well above 5.83, and to replace q by a relation of the form -- Sli m =
g l ( q -- q(ot lira)) -+- g z ( q
-- q(~
lim)) 2
(6)
when c~ is less than the limiting value. The constants K~ and K 2 are selected such that the resulting q(o~) has continuous first and second derivatives at C~l~m. A variety of alternatives have later been proposed (Michelsen, 1990b; Holderbaum and Gmehling, 1991; Soave, 1992). They all replace q ( a ) by an approximating function q * ( a ) which is subsequently used for all values of e~. Michelsen (1990b) proposed the linear approximation MHV-I: q* (o~) = q0 + q , ~ , q, = - 0 . 5 9 3
(7)
and the quadratic approximation, MHV-2: q* (c~) = qo + q,c~ + q2oL2, q, = - 0 . 4 7 8 , q2 = - 0 . 0 0 4 7
(8)
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
17
Holderbaum and Gmehling used MHV-1 with a different value of ql ( - 0 . 6 4 3 ) and Soave used a somewhat more complex functional form: OL *(or) = 3.365 - - - l n 2 2
2 1 n 2 ) 2 + 17.2511/2
(9)
With all these mixing rules the EOS fugacity coefficients are determined from ln¢i = In p ( v Z b )
+
(
'
v-b
°
" + b J/--ffC
Jw
,j + I n
v
1--~ni ]
wn
(10)
where the composition derivative of e~ is calculated from
(dq)(tno) ) -~
Oni
c~ = q * ( et ii ) -- q * ( ot ) + l n "y A M + l n b i---~i+ -~- - - 1
(11)
All these variants require virtually identical computational effort. The exact model does not permit explicit solution of Eq. (5) for or, but numerical solution by Newton's method is straightforward as q ( a ) is monotonically decreasing, with derivatives given by: dq --
dc~
u0 + 1 -
I n - - ,
u0
d2q de~2
1 -
e~(2Uo + 1 - e ~ )
(12)
and the additional computational effort in this step is negligible.
3. R a n g e
of applicability
The EOS excess Gibbs energy calculated with the approximate methods all deviate to some extent from that of the activity coefficient model, the errors associated with each method depending on the e~ range in which they are used. As the purpose of these methods is to utilize g E correlations for data measured at conditions, where such correlations have been used for representing vapour-liquid equilibrium, i.e, data measured around atmospheric pressure or lower, we shall first investigate the relation between the et parameter and the saturation pressure. The pure component liquid fugacity at zero pressure is given by In RT)
= q(e~)
(13)
where, for simplicity, we have omitted the component index. At saturation, the fugacities in the liquid and the vapour phase are identical, i.e. fl(Psat) = L ( P s a t ) =
Psat~v(Psat)
(14)
18
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
10 -L
~
\ 10
-2
(73 (/3 &) L10
-3
\
RD ©
~
1 0 -'-
\
(1) Ct2 10
-5
\
5
110
1'5 0<
2~D
\ 2~5
Fig. 1. Relationship between the pure component o~ parameter and the reduced saturation pressure for a component with an acentric factor of 0.4. If we neglect the pressure dependence of the liquid phase fugacity and assume the vapour phase to be ideal, Eq. (14) yields, and Eq. (13) can be written
RTl =ln - - ~ ] =lnlBb +lnP~.sat--lnT~=q(c~ )
In
(15)
where ~-~b = 0.08664. In this expression the relation between Tr and o~ is component specific, but the variation of In T~ with o~ is modest, the reduced temperature typically decreasing from about 0.7 to about 0.4 when OL increases from 10 to 20. An approximate relation between T~ and o~ can be obtained from Soave's expression (Soave, 1972) for the temperature dependence of the a parameter a - [1 + m(co)(1 - T°5)] 2, m ( o ~ ) = 0 . 4 8 + 1 . 5 7 4 o ) - 0 . 1 7 4 ( o 2 (16) acrit or
O~
(1 + m ( c o ) ( 1 -
Tr°'5)) 2
(17) (1 crit
Tr
where o~cr~t= ~ a / ~ ' ~ b = 4.73. Substitution of this relation between the reduced temperature and o~ into Eq. (15) yields the desired relation between c~ and the reduced saturation pressure, as shown in Fig. 1 for a representative eo value of 0.40. A saturation pressure of 1 bar corresponds to reduced pressures in the range 0.005-0.04, the lower limit applying for water. From Fig. 1 w e notice that the corresponding oL range is 10-13, and the coefficients used in MHV-1 and MHV-2 are selected to provide particularly accurate approximations in this range. The reduced pressure decreases very rapidly with increasing values of OL and is 10 4 at o~= 19.
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26
19
4. Error estimate for approximate mixing rules
The error in the EOS Gibbs free energy, i.e. the difference between the g Z calculated from the equation of state and that of the activity coefficient model, can be determined by subtracting Eq. (5) (with q * replacing q) from Eq. (2): AgE
(gE)
R-T =
~
where e ( a ) = q(a and write
(gE) EOS-- ~
AM = e ( c ¢ ) -
~zie(c~ii) i
(18)
- q * (¢x). This expression can be simplified as follows: We define fi = ~i zi a ii
]
~--~ = e ( N ) - z i e ( ~ i i
) +[e(a)-e(N)]=e,+e
2
(19)
The first term, e~, is not related to excess properties but depends only on the pure component oL parameters. It is identically zero when these are equal. The second term can be approximated by a Taylor series expansion from e~ = (de/ (a-~)-e(~)=( e(c¢) - e ( ~ ) = e ( ~ ) + I dcx ]~
de] (cx-~) I dot ]u
(20)
For the present purpose an adequate estimate of a = ~ can be obtained from the MHV-1 approximation, a=
1 1 Y'~ziaii+--r(z)=~+--r(z) i q= ql
(21)
where qj = - 0 . 5 9 3 and
r(z):
~
+ ~ziln AM
i
bii
We thus arrive at l(de)
e2=--
r(z)
(22)
and this term is therefore proportional to the activity coefficient model excess Gibbs energy, modified by the b term. Fig. 2 shows e ( a ) and d e / d R for the MHV-2 approximation, q* = q0 + qta + q2 a2. The value used for qo is of no importance for the analysis, and it is here chosen to match q at a = 12. The e~ contribution can be determined by a simple geometrical construction. Consider two components with pure component oL values of 16 (point B) and 11 (point A), respectively. The line connecting A and B corresponds to z~ e(al~) + z2e(a22) and e~ is equal to the vertical distance from the line AB to the e(e~) curve, as marked by the arrow at zl = 0.7 (~ = 14.5) It is evident that e~ is particularly small when the a values for the pure components are close to the value at the inflection point (~ = 12), and a very modest error contribution results when both
20
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26 0.03
e end de/dc~ for MHV-2
/ / /
0.02
\~
\\
---
0
L_ L
/ /
de/do~
\
0.01
b_J
e
///
//////
0.00
-0.01
-0.02
i
8
I0
l~2
14
IJ6
18
2~0
G
Fig. 2. Magnitude of e and d e / d ~ for the MHV-2 mixing rule.
values are located in the interval 10 to 15, corresponding roughly to pure component partial pressures in the range 0.1-1 bar. Inside this range the magnitude of e 2 is also very modest, the factor of r(z) being below 0.01 in the entire range 9.5 < a < 16. Acceptable results from MHV-2 can thus be expected for mixtures with pure component partial pressures in the range form about 0.01 bar to a few bar. Fig. 2 also illustrates the dramatic increase in the magnitude of e~ for mixtures containing a light ((~ = 10-12) and a very heavy (c~ > 18) component, and in practice the use of MHV-2 is not advisable outside the range 8-18. Activity coefficients at 298.15 K and 1 bar for the mixture hexane-hexadecane, investigated earlier by Kalospiros et al. (1995), are shown in Fig. 3 using UNIFAC as the excess Gibbs energy model. The pure component c~ values are 12 and 28, respectively, and the behaviour of MHV2 is clearly unsatisfactory. Fig. 4 shows the e and d e / d ( ~ for Soave's approximation, Eq. (9). This expression has the advantage that the correct value for the derivative, dq/da = - I n 2, is obtained in the limit a --* oo. We notice that the Soave expression provides a remarkably good match in the entire range 9 < e~ < 20, being clearly superior to MHV-2 at large a values. Among the mixing rules investigated here the Soave approximation appears to be the superior with respect to accuracy as well as c~ range. When the mixture contains components with small (10-12) as well as large (above 20) a values, the magnitude of e~ can however become quite substantial. The actual difference between the EOS excess Gibbs energy and that of the activity coefficient model is compared with that predicted by Eq. (19) for the acetone/water mixture at 340 K, 1 bar. The UNIQUAC equation with parameters estimated in the Dechema series (Vol I, 1, p.237) from the data of Kojima et al. are used for the excess Gibbs energy calculation. The predicted error is in close agreement with the calculation (Fig. 5), and the the RMS deviation between the reduced excess Gibbs energy calculated from the equation of state and the activity coefficient is 0.0022. For the Soave approximation the RMS error in reproducing the UNIQUAC excess Gibbs energy has an RMS value
21
M.L. Michelsen /Fluid Phase Equilibria 121 (1996) 15-26 1.0
Hexone-Hexodecane at 298 K 0.5
---
©
UNIFAC MHV-2
0.0 C:3 O> c- - O 5
/
l
\
\
/ / -I
/
.0
/
\
/
\ \
/
\ \ \
-1.5
0.'2 0.'4 0.'6 0 '8 1.'0 H e x o n e mole froction Fig. 3. Activity coefficients in the mixture hexane-hexadecane at 298 K, 1 bar, calculated with MHV-2 and with UNIFAC. 0.0
of only 0.0015. The difference in the pure component c~ values (10.1 and 15.2, respectively) for this mixture is quite substantial. The error curves for the linear MHV-1 and Holderbaum-Gmehling ( H / G ) approximations are shown in Fig. 6. It is readily shown that e~ is only related to the curvature of the error, and e~ will therefore be the same for the two methods. It is also clear that e~ is always positive, and increasing in magnitude with increasing difference between the pure component ct values. In the o~ range 10-15,
0.0.3 e
and
002
de/dG
for
--
e
- - -
l~'e'~ ~
Soove
001 k_ 0 L k_ I,I 0.00
-0.01
-0.02
8
110
Ii2
114
I~6
118
2~0
O~
Fig. 4. Magnitude of e and de/dot for the Soave approximation.
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
22
0.008 - -
Actual Estim.
error error
0006 Od
/
LLI 0.004 C7~
~_ 0
f
--
~ \
/
\
0 002
tm
0.000
-0.002
O. 0
0.~2 0.~4 0.'6 0.'8 Acetone mole fraction
1 .'0
Fig. 5. Comparison of error estimates and actual errors for the MHV-2 mixing rule applied to a mixture of acetone and water at 340 K, 1 bar.
the error in e 2 for MHV-1 is substantially smaller than that obtained with the Holderbaum-Gmehling coefficient, as illustrated in Fig. 7 for the acetone/water mixture. The factor of r(x) in the Holderbaum-Gmehling approximation is however negative, and as the excess Gibbs energy for most mixtures is positive the error contributions from e I and e 2 in this approach will partly cancel. As a result the RMS error for the Holderbaum-Gmehling approximation (0.014) is actually lower than that
0.15 e and
de/dc<
for
0.10
MHV-1
and
H/O
e "'""
0.05
\
""-...
---
-...
- ....
de/dc~,
de/dcL
MHV-1
1'4 O~
1'6
1'8
H/O
\
k_
0.00 k._
LJ
-0.05
-0.10
-0.15
8
1~0
1'2
210
Fig. 6. Magnitude of e and d e/do~ for the MHV-1 and the Holderbaum-Gmehling mixing rules.
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26 0 050
23
Acetone-Woter, 540 K using MHV-1 and H/G -
-
el
-- - 0025
. . . . . . .
e 2, M H V - 1 e2, H / G
Od DJ C~
0.000
© k_ LJ
- 0 025
-0050
- ~
09
02
0.'4
0.16
0.8
10
Acetone mole froction Fig. 7. Error contributions for the MHV-I and the Holderbaum-Gmehling mixing rules applied to a mixture of acetone and water at 340 K, 1 bar.
of MHV-1 (0.022). As expected, the errors in both approximations are substantially larger than that of MHV-2 and the Soave approximation.
5. Extending the range The majority of the data published in the Dechema series (Gmehling and Onken, 1977) correspond to pure component a parameters in the range 9 - 1 7 , the most important exception being data at a very low pressures. An ability to be able to incorporate components with high et values could also be of importance in connection with representation of liquid-liquid equilibria or for use with polymer mixtures (Kalospiros et al., 1995). Finally very large a values may occur when a generalized correlation such as UNIFAC is used as the excess Gibbs energy model (Kalospiros et al., 1995), but such cases are probably less important since the excess Gibbs energy model is likely to be inaccurate when an extreme volatility range has to be covered. The hexane-hexadecane mixture analyzed by Kalospiros and Tassios, 1995 is e.g. characterized by a difference of 5 orders of magnitude in the pure component saturation pressures. In order to enable the MHV-2 approximation to handle such mixtures our suggestion is to combine MHV-2 with the exact model of Michelsen. Thus we use for the approximating function q* (et) = q ( a ) q*(o~) =qo+qlot+q2ot
for ot _> O~lim 2
(23)
for ~ < ~li m
At tx = a lira we require continuity in q * and its first two derivatives. Finally, the limiting value is chosen such that the coefficients q~ and q2 match those of MHV-2 as well as possible. A suitable choice is et lira 12.3, where we obtain q'(a) -- - 0.5986 and q"(et) = - 0.00918, the corresponding derivatives using the conventional MHV-2 constants being - 0 . 5 9 4 6 and - 0 . 0 0 9 4 , respectively. The =
M.L. Micbelsen / Fluid Phase Equilibria 121 (1996) 15-26
24
Hexane-Hexadecane
at 550 K
0.2
-
<3
0.0
-
UNIFAC
- -
Modif.
MHV-2
E E o cr~
-0.2
c
-O.4
-0.6
O.0
0.12 0.14 0 J6 0.18 Hexane mole fraction
1.10
Fig. 8. Activity coefficients in the mixture hexane-hexadecane at 350 K, 1 bar, calculated with the modified MHV-2 mixing rule and with UNIFAC. close match implies that correlations based on MHV-2, such as the gas solubility tables of Dalai et al. (1991), can be used with unchanged parameters. When the mixture ct parameter exceeds oLlira in the entire composition range the zero pressure excess Gibbs energy is equal to that of the corresponding activity coefficient model. A situation of particular interest is that where one component (1) has an c~ value lower than ct lira while the c~ value for the mixture and for the other component (2) is larger than ct lim" For such a mixture we have gE A--
= e(ct)
RT
-
Ezie(oLii
) =
-
Zle(Otll )
(24)
i
since e(o~2) = e(o~) = 0. The error in the activity coefficient for component 1 therefore is 0 Aln~, 1 ---- -ff~-(--ne(ct it)) = --e(t~ll )
(25)
IJft i.
and is thus given by the error in the pure component q value. For the hexane-hexadecane mixture at 298 K, the ct value of the light component is very close to the limiting value, and essentially exact reproduction of the activity coefficient model is obtained. At 350 K, the pure component a values are 9.4 and 21.8, respectively, and Fig. 8 shows that the activity coefficient of the light component is now well represented at all composition. The infinite dilution activity coefficient of hexadecane in hexane is somewhat in error, but this deviation is of no importance for the calculation of vapour-liquid equilibrium since the vapour pressure of the heavy component is extremely small. For the acetone-water mixture the accuracy of the modified MHV-2 mixing rule, with an RMS error of 0.001 l, is better than that of the Soave approximation. It should, however, be kept in mind that the accuracies of MHV-2, the Soave approximation and the modified MHV-2 mixing rule in
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
25
reproducing the activity coefficient model are all well below the accuracy of the activity coefficient model in representing the original data. Improvements of any of the 3 methods for this mixture therefore serves no practical purpose.
6. Conclusion A simple and illustrative graphical tool to estimate the error associated with a number of widely used approximations to the exact method of Michelsen for matching the zero pressure excess Gibbs energy of an equation of state to that of an activity coefficient model has been developed. In addition we have provided a modification of the MHV-2 mixing rules that permits accurate handling of mixtures containing components of widely differing volatilities.
7. List of symbols Mixture parameter, Redlich-Kwong equation Mixture parameter, Redlich-Kwong Equation Difference between approximate and exact q function e Contributions to e, Eq. (19) e I , e2 Fugacity f Excess Gibbs energy gE Constants in extrapolation method Kl, K2 Function of acentric factor in SRK equation, Eq. (16) m Pressure (bar) P Reduced pressure q Function defined by Eq. (3) q* Approximation to q qo, ql, q2 Coefficients in polynomial approximation to q Gas constant R Temperature (K) T Reduced temperature bl Reduced molar volume, u = v / b Molar volume U Mixture mole fraction Z a
b
7.1. Subscripts 0 AM EOS lim i ii
Zero pressure value Activity coefficient model property Equation of state property Limit for eL value in extrapolation method Component index Pure component index
M.L. Michelsen / Fluid Phase Equilibria 121 (1996) 15-26
26 7.2. G r e e k
(l)
f~A, f~B q~
letters
Mixture parameter, a = a/(bRT) Acentric factor Constants of Redlich-Kwong equation Fugacity coefficient
References Boukouvalas, C., Spilitois, N., Coutsikos, Ph., Tzouvaras, N. and l'assios, D., 1994. Prediction of vapor-liquid equilibrium with the LCVM model: a linear combination of the Vidal and Michelsen mixing rules coupled with the original UNIFAC and the t-mPR equation of state. Fluid Phase Equilibria, 92: 75. Dahl, S. and Michelsen, M.L., 1990. High-pressure vapor-liquid equilibrium with a UNIFAC-based equation of state. AIChE. J., 36: 1829. Dahl, S., Fredenslund, Aa, Rasmussen, P., 1991. The MHV-2 model: a UNIFAC-based equation of state model for prediction of gas solubility and vapor-liquid equilibrium at low and high pressures, Ind. Eng. Chem. Res., 30: 1936. Gmehling, J., Onken, U., 1977. Vapor-liquid Equilibrium Data Collection. Chemistry Data Series, Vol !. Part 1. DECHEMA, Frankfurt. Heideman, R.A. and Khokal, S.L., 1990. Combined excess free energy models and equations of state. Fluid Phase Equilibria, 56: 17. Holderbaum, T, Gmehling, J., 1991. A group contribution equation of state based on UNIFAC, Fluid Phase Equilibria, 70: 251. Kalospiros, N.S., Tzouvaras, N., Coutsikos, P.. Tassios, D.P., 1995. Analysis of zero-reference-pressure EoS/G E models, AIChE.J., 41: 928. Kalospiros, N.S., Tassios, D., 1995. Prediction of vapor-liquid equilibria in polymer solutions using an equation of state/excess Gibbs free energy model, Ind. Eng. Chem. Res., 34:2117. Michelsen, M.L., 1990a. A method for incorporating excess Gibbs energy models in equations of state, Fluid Phase Equilibria, 60: 47. Michelsen, M.L., 1990b. A modified Huron-Vidal mixing rule for cubic equations of state, Fluid Phase Equilibria, 60: 213. Soave, G., 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci., 27: 1197. Soave, G., 1992. A new expression of q(a) for the modified Huron-Vidal methods. Fluid Phase Equilibria, 72: 325. Wong, D.S.H. and Sandler, S.I., 1992. A theoretically correct mixing rule for cubic equations of state. AIChE. J., 38: 671.