COMPUTATION OF VAPOR-LIQUID EQUILIBRIUM DATA FROM SOLUTION VAPOR PRESSURE MEASUREMENTS F
. 0 . h1 I X 0 N ,
Research Triangle Institute, Durham, N . C.
B0 G DA N GU M0 W S
K I AND B
. H. CA RP ENT E R ,
Union Carbide Gorp., South Charleston, W. Va.
A new method for computing vapor-liquid equilibrium data from solution vapor pressure measurements for systems of any number of components involves determination of the excess Gibbs free energy by successive approximatioris to the solution vapor pressure surface, with activity coefficients obtained b y differentiation of the free energy function. Results are given for the chloroform-ethanol and acetone-chloroform-methanol systems.
HE usual technique for obtaining vapor-liquid equilibrium T d a t a is by direct measurement-Le., equilibrium is established, and phases are Isampled and analyzed. Normally in the operation of equilibrium stills, the experimental technique must be rather delicate in order to ensure meaningful results. This fact, coupled with the necessity for much analytical work, tends to enhance interrst in methods for the determination of equilibrium data that do not involve sampling and analysis of the vapor phase. I t has long been realized, in the case of binary solutions a t least, that the behavior of one component is coupled to that of the other through the Gibbs-Duhem equation and that the gross behavior of the solution (vapor pressure us. composition or equilibrium temperature us. composition) is thus uniquely dependent on the behavior of only one of its components. This dependency then suggests the possibility of inferring, within the framework of classical thermodynamics, the behavior of one component of a binary solution from a knowledge of the behavior of the solution itself. More specifically, it is possible to determine vapor-liquid equilibrium data from experimental measurements of solution vapor pressure over the composition range. Several methods have been employed for the calculation of component behavior frclm gross solution behavior. Ljunglin and Van Ness (4) have suggested the classification of these methods into two categories, direct and indirect methods. The direct methods involve calculation of vapor compositions by integration of the coexistence equation, a first-order differential equation derived from the Gibbs-Duhem equation relating phase compositions a t equilibrium. Hala, Pick, Fried, and Vilim (3)have given a detailed discussion of the basic direct method, and Ljunglin and Van Ness (4) have discussed techniques for handling nonconstant temperature or pressure conditions as well as nonideal vapor phase behavior. The indirect methods involve first the calculation, by some appropriate means, of thLe liquid phase activity coefficients and subsequent calculation of vapor compositions therefrom. These methods usually involve ascertaining which of selected solution equations to the Gibbs-Duhem equation lead to the best fit to the experimental vapor pressure data, and of the determination of the pal-ametric values producing the best fit. Barker (7), for example, has developed a procedure based on the assumption that the excess free energy can be represented as a polynomial function of composition. This assumption is, of course, equivalent to assuming that the Margules equations relate the activity coefficients to composition.
There is a basic difference in the degree of rigor associated with the direct method and the indirect method of Barker. In the former, one makes no assumptions about the solution behavior-Le., the nature of molecular interactions. Solution behavior is determined directly from the experimental vapor pressure data. The Barker method necessitates the assumption of a particular model and the estimation of its parameters, This deficiency in the method of Barker has been recognized by Tao (7), who has presented another indirect method in which the necessity for the a priori assumption of a particular functional form for the excess free energy has been removed. Tao’s procedure involves calculation of the activity coefficienu essentially by integration of an equation resembling the coexistence equation. His procedure, though indirect, retains the rigor usually associated with the direct method. Virtually all of the methods for calculating vapor phase compositions from solution vapor pressure measurements have been applied only to binary systems. Apart from numerical complexities, however, there is no reason why such calculations cannot be made for ternary and higher ordered systems. If one considers the application of the direct method to ternary systems, one finds that the problem increases in complexity from the solution of one first-order ordinary differential equation to the simultaneous solution of two first-order partial differential equations. Not only are two independent variables involved (the two independent liquid phase compositions) but also two dependent variables (the two independent vapor phase compositions). I t is easy to see that generalization of the direct method to higher ordered systems will be limited by the availability of computing machinery capable of handling the required dimensionality. An indirect method such as that of Barker is readily and easily generalized to ternary and higher ordered systems, but this method retains the disadvantage of lacking rigor as compared with the direct method. The method of Tao appears specific to binary systems and does not seem to be easily generalized. This paper presents an indirect method that retains the same degree of rigor that characterizes Tao’s method and integration of the coexistence equations. The a priori assumption of a functional form for the excess free energy is not necessary, and the method is comparatively easy to generalize to ternary and higher ordered systems. Basically, an expression is written for the solution vapor pressure in terms of the excess free energy and its composition derivatives. The expression is “inverted” to VOL. 4
NO. 4
NOVEMBER 1 9 6 5
455
give the excess free energy function that produces the observed vapor pressure behavior, and activity coefficients are then deduced from the excess free energy function. Thermodynamic Formulation
Assuming ideal vapor phase behavior, the total pressure for a system of N components is given by 1v-
where the correction terms, g(n)( E ) , are determined so that P c n + l ) ( ~will ) be a closer approximation to P ( E ) than P ( % ) ( E ) . If the excess free energy values are incremented according to Equation 6, the corresponding pressure change will be given to within the locally linear approximation by the Taylor expansion of Equation 5 to linear terms in g ( n ) ( ~ about ) G(')(?i). That is,
1
(7)
k=1
Equation 1 can be approximated by its finite difference representation in the following manner. The 1%' - 1 independent composition coordinates are regarded as orthogonal coordinates in N - 1 dimensional space and further regarded as being discretized with lattice interval A , so that xi = aiA for integral ai
(2)
With this representation, both the pressure and the excess free energy have values for all lattice points and can be regarded as entries in A: - 1 dimensional arrays, designated as P(a1, a2, . . ., as-l) and G(a1, (YZ, . . ., aAv--l). If, furthermore, i k is the N - 1 dimensional elementary vector with unity in the kth position and zeros elsewhere, i k = (0,0, . . ., 1, 0,. . .)
G(E
+ i k ) - G ( E - Zk)
bxl,
Direct calculation of the partial derivatives from Equation 5 and substitution into Equation 9 yield, after some rearrangement,
(3)
then the finite difference representation for the partial derivatives can be written as
bG _
where the partial derivatives are calculated from Equation 5 and are evaluated at G(')(E) and where the sum is to range over all admissible values of 8. One would like to determine the correction terms, g ( n ) ( E ) , from Equation 7 in such a ) identical with the experimanner that P ( n + l ) ( ~becomes mental pressure, P ( E ) . Substitution of this requirement into Equation 7 results in
(4)
2A
Then the finite difference representation of Equation 1 becomes
Equation 9, repeated for all lattice points, provides a linear set of equations in the correction terms g(E) that can be solved by relaxation techniques-Le., by either block or point relaxation, as preferred. I n this study, block relaxation techniques were' used for binary calculations, and point relaxation was used for the ternary calculation. Binary Calculations, Block Relaxation
N-1
c
i-1
G(E a4
+ pi)
-
G ( E - ii)
2A
For a binary system, Equation 9 becomes
I+
p1W - p l i d 2A
Equation 5 shows the explicit dependence of the pressure, P , a t the lattice point represented by E on the free energy a t lattice point 8 and the neighboring lattice points. The computational problem is that of finding values of the excess free energy function a t all lattice points that will reproduce the experimentally observed pressures according to Equation 5 .
g(.
+ 1) = P(a) - P'"'(a)
(10)
If Equation 10 is then repeated for all lattice points, there results, in matrix-vector form,
Computational Procedure
The computational procedure for solving Equation 5 is a combination of Newton's method with either block or point relaxation techniques. Assume that one has available the nth estimates of the free energy values, W ) ( E ) , a t the lattice points. Corresponding to these estimates of the free energy will be estimated values for the solution vapor pressure, P ( n )( E ) , where P ( n ) ( ~is) calculated from Equation 5 with the nth estimates of the free energy appropriately substituted. In general, (a)will differ from P(E),the experimental solution vapor pressure, and G ( " ) ( E ) must be incremented to G(%+l)(~)-i.e.,
G(n+l)(a)= G ( ' ~ ( E+ ) g(")(&) 456
I&EC FUNDAMENTALS
(6)
P(1) P(2) P(3)
&l
- P'"'(4) - P'"'(2) - P'"'(3)
.AM):
(11)
Table 1. Calculated Results for Chloroform-Ethanol System at 55" c.
The set of equations represented by Equation 11 can be solved very efficiently utilizing a method described by Bruce et al. ( 2 ) . The procedure is as follows: C(1)
A.
Set w(1) = B ( 1 ) , b(1) = -,f(l)
B.
Compute
=
41)
W(.)
=
B(.)
- A(a)b(a -
1)
(1 6)
AGE/RT y, YZ 0 1.0000 279.69 321.73 0.0290 1.8205 0.9990 365.95 0,0596 1.8572 0,9975 409.58 0.0908 1.8589 0.9974 450.64 0.1216 1.8334 1 , 0 0 0 4 487.83 0.1511 1.7876 I 0079 520.44 0.1786 i17274 1.0212 548.24 0.2035 1.6577 1.0417 571.38 0.2249 1.5831 1.0710 590.29 0.2423 1.5069 1.1109 605.58 0.2552 1.4322 1.1632 617.95 0.2630 1.3613 1 2305 628.05 0.2652 1.2958 1.3155 636.43 0.2614 1.2369 1.4218 643.39 0.2511 1.1846 1.5553 648.91 0.2339 1.1385 1.7275 652.54 0.2090 1.0960 1.9706 653.31 0.1734 1.0434 2.4976 649.60 0.1247 1.0979 2.9663 639.05 0.0667 1.0043 3.5001 1.0000 1 .o 618.49 0 a Pressure calculated f r o m following polynominal: P = 279.69 f 796.10~1 1751.98~1~ - 5506.50~1~ 6603.88~1~ - 2706.66~1~. y~
x1
P,"Mrn.Hg
0 0.1750 0.3139 0.4211 0.5033 0.5666 0.6158 0.6546 0.6854 0.7105 0.7314 0,7494 0.7657 0.7813 0.7972 0.8139 0.8311 0.8396 0.8723 0.9234
0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1 .o
~
+
+
Results
Starting with the initial approximation
G(a) = 0
(21)
iteration of the above procedure to a value of n of around 6 to 10 yields the excess free energy function producing the desired solution vapor pressure. Activity coefficients are then calculated as indicated by Equation l , and vapor phase compositions by the usual equations.
The proposed method has been tested on several binary systems and on one ternary system. I n all cases, polynomials, fitted to the experimental pressure data by least squares, have been used for smoothing the experimental data and for interpolation. For the chloroform-ethanol system at 55" C., the solution vapor pressure data of Scatchard and Raymond ( 5 ) Ivere fitted by the polynomial
P
= 279.69
Ternary Calculations, Point Relaxation
Because of difficulties in programming block relaxation calculations over the triangular composition diagram, point relaxation was used in preference to the more efficient block relaxation methods. Clne simply computes a new value of g a t each lattice point from the old values ofg a t the neighboring lattice points as indicated by Equation 9. When convergence is obtained on g, the free energy values are incremented according to Equation 6. New coefficients for the g values in Equation 9 are calculated, and the calculation is repeated until convergence to the experimental pressure is obtained. Again, activity coefficients are calculated as indicated by Equation 1, and vapor phase compositions by the usual equations.
1.0
+ 7 9 6 . 1 0 ~ 1+ 1 1 5 1 . 9 8 ~ 1-~ 5 5 0 6 . 5 0 ~ 1+ ~ 6 6 0 3 . 8 8 ~ 1~ 2 7 0 6 . 6 6 ~ 1 ~(22)
This equation fits the experimental data to within u limits of 2.18 mm. of Hg. Vapor phase compositions for this system have been calculated both by the present method and by integration of the coexistence equation as described by Hala et al. ( 3 ) . Both methods produce identical results; these results are summarized in Table I and are shown graphically in Figures 1 and L.
The calculated values of the vapor composition are higher by approximately 0.005 mole fraction of chloroform than the reported values over the entire composition range. This syste-
I
7001
bLCULbTED RESUL
0
'
Figure 1 .
'
bNO RbYHOND
02
' '
C'8 a , . MOL FRACTION CHLOROFORM I N LICIUID
0.4
0.6
1.0
2ooi
,
,
I
1
0.2 0.4 0.6 0.8 I!O x , - M O L FRACTION CHLOROFORM IN LIQUID
Chloroform-ethanol at 55" C.
A. Vopor-liquid equilibrium curve 6. Solution vapor pressure-liquid composition curve
X I - M O L FRACTION CHLOROFORM IN LIQUID
Figure 2. A.
6.
Chloroform-ethanol at 55" C. Excess free energy os a function of liquid composition Activity coefficients as o function of liquid composition
VOL. 4
NO. 4
NOVEMBER 1 9 6 5
457
Experimental (6) and Calculated Results for Acetone-Chloroform-Methanol System at 50' C.
Table II.
Liquid Composition, Mole &fao
Xl
Vapor Composition, Mole Y1 Y2
x3
E 0,051 0.904 C 0,050 0.900 E 0.200 0.590 C 0,200 0.600 E 0.453 0,047 C 0.450 0.050 E 0.045 0.053 C 0.050 0.050 E 0.195 0.188 C 0,200 0.200 E 0.896 0.052 C 0.900 0.050 E 0.401 0.140 C 0.400 0.150 a E. Exjerimental data. C.
%
X2
0.045 0.050 0.211 0.200 0.500 0.500 0.902 0.900 0.622 0.600 0.052 0.050 0.459 0.450
0.023 0.023 0.119 0.126 0.550 0.529 0.081 0.098 0.218 0.216 0.914 0.914 0.422 0.437
0 837 0.834 0.571 0.548 0.035 0.041 0.119 0.125 0.252 0.253 0.026 0.024 0.120 0.137
PI3
=
P23
=
+ XZPZO +
x1Plo
x3P3O x3P3O
+ +
x1x3[a31x1
+ +
a l m
X ~ X ~ [ ~ ~ Z aX Zm
- ~ U X I X ~ ] (24) - d 2 3 x m I (25)
For the ternary system : =
XlPlO
+ + + +
+
x2P20
xm[a21x1 xlx3[@1x1
+
x3P30
+ + + + + am2
~IZXIXZ]
am?
d l m x ~ ]
X ~ X ~ [ ~ ~ Z a23x3 X Z
xmxz[a21
013
d 2 3 ~ ~ ~ 3 1
023
-c
m - czxz
- ~2x31
(26)
This particular representation of the pressure data was chosen for the following reasons : Equation 26 reduces to 23, 24, or 25 for x3, x 2 , or x1 equal to zero, respectively. The binary parameters can be estimated from vapor pressure data on the binary systems alone, leaving only the c values to be estimated from ternary pressure data. For the system tested, this particular form fits the experimental pressure data to within an average deviation of about 1% of the total pressure. Vapor phase compositions have been calculated for the acetone-chloroform-methanol system using the proposed method. Computed results a t selected points are compared with the experimental data in Table 11. Though there is some 458
Y3
0.140 0.143 0.311 0,327 0.415 0.429 0.800 0,777 0.529 0.531 0.061 0.062 0.408 0.426
Activity Coefficients
Y1
YZ
Y3
0.424 0.414 0.571 0.594 1.170 1.117 1.400 1.543 1.010 0.984 1.010 1.002 1.110 1.011
1.010 1.001 1.080 1.032 0.814 0.940 1.990 2.387 1.400 1.374 0.568 0.574 0.935 1.010
4.190 3.818 2.010 2.260 1.140 1.208 0.997 1.010 1.100 1.192 1.640 1.790 1.190 1.298
Pressure, Mm. H g
577 555 586 580 588 583 47 3 485 552 558 604 605 575 568
Calculated datn.
matic discrepancy is the result of a slight overestimation of the activity coefficient of chloroform, as shown in Figure 2B. The computed values of the excess free energy are somewhat lower than the results derived from Scatchard and Raymond's data. This systematic discrepancy is apparently the result of failure to correct for vapor phase nonideality. Torgeson has reported a better fit to the experimental data when provisions for vapor phase nonideality are included in a procedure based on numerical integration of a coexistence equation (8). The fact that the present method and the coexistence equation produce identical results for the ideal vapor approximation suggests that inclusion of vapor phase nonideality into the present method would remove the systematic discrepancy between the experimental data and the calculated results. For the acetone-chloroform-methanol system, the solution vapor pressure data of Severns et al. (6) were fitted by the following form : For the binary pairs:
P
%
l&EC FUNDAMENTALS
discrepancy, the agreement is satisfactory. Part of the discrepancy between calculated and experimental results can be attributed to slight differences in the liquid compositions a t the chosen points. Also, Severns et al. used compressibility factors to account for the nonideality in the vapor phase. As with binaries, no corrections for vapor phase nonideality were made in the present work. Discussion
The present method for computing vapor compositions from solution vapor pressure measurements has several advantages over previous methods. For example, integration of some form of the coexistence equation for binaries, as pointed out by Hala et al. ( 3 ) ,requires special attention to the regions at the ends of the concentration range and at the azeotrope, if one exists. The coexistence equation reduces to an indeterminate form at these points, and an alternate formulation must be applied. The present method requires no such special attention at any point along the entire composition range. The present method is successful and calculations proceed smoothly for some cases in which the coexistence equation method meets with difficulty. For example, we have been unable to obtain a realistic numerical integration of this equation for the acetone-chloroform system, which exhibits a strong No difficulty was negative departure from ideality. experienced with the method described herein for this or any other system tested. The difficulty in integrating the coexistence equation seems to result from the fact that the experimental vapor pressure data must somehow be differentiated as part of the procedure. Our method for doing this was to obtain a least squares polynomial fit to the data and then to differentiate the polynomial. However, even though such polynomials can fit the vapor pressure data very well, there can also be local oscillations of the polynomial function slightly above and below the data. If such is the case, errors resulting from differentiation of the smoothed data can apparently cause difficulty in subsequent integration of the coexistence equation. No doubt such difficulties can be avoided by special treatment of the raw data to ensure sufficient smoothness in the derivatives, but the computational procedure presented in this paper completely avoids the problem. Differentiation of the data is not required, and the results are much more insensitive than those from the coexistence equation method to the particular technique chosen for smoothing the data.
If one considers the extension of the coexistence equation method to ternary systems, one is faced with two partial differential equations involving two independent variables (liquid phase compositions) aind two dependent variables (vapor compositions or partial pressures). I n the present formulation, one has a single partial differential equation in a single dependent variable (the excess free energy). The latter is somewhat more convenient from a computational point of view, and certainly becomes more attractive as one considers solutions with more components. Of the methods reported in the past for calculation of vapor compositions from solution vapor pressure data, only that of Barker ( 7 ) seems easily generalized to ternary and higher ordered systems. Barker's method involves assuming a polynomial form for the excess free energy and then using statistical techniques to estimate the parameters therein so as to get the best fit to the experimental pressure data. The present method has an advantage over Barker's method in that it is not necessary a t the outset to select a model. Though we used polynomials for data smoothing and interpolation, this is not necessary. I n principle, the method can be applied to graphically smoothed data or, if sufficient data are available, to the raw data themselves. Thus, the present method is not limited to the a priori choice of a particular excess free energy functionality. If deviations from ideality in the vapor phase are significant, these deviations can be taken into account rather easily. One first computes vapor phase compositions under the assumption that ideality holds, then uses these compositions to estimate vapor phase activity coefficients by the usual methods. With Equation 1 replaced by
the entire procedure is repeated for the final relaxation. The effect will be to change the coefficients in Equation 9 slightly by the introduction of the gas law correction terms.
Nomenclature
A
quantity defined by Equation 12 parameter in vapor pressure function, Equations 23 to 26 b = quantity defined by Equation 17 B = quantity defined by Equation 13 c = parameter in vapor pressure function C = quantity defined by Equation 14 d = parameter in vapor pressure function f = quantity defined by Equation 18 g @ ) ( n ) = nth estimate of correction to G ( ~ ) ( E ) G = GE/RT GE = excess Gibbs free energy of mixing G ( " ) ( E )= nth estimate of value of G at lattice point designated by E = number of comDonents in solution = nth estimate of partial pressure of kth component, = =
a
XkYi;(")Pk0
= = =
= W
=
xi
Yk
= = =
A
=
ak
=
6k
=
E
partial pressure of kth component in an ideal solution, xkPko solution vapor pressure vapor pressure of kth component nth estimate of solution vapor pressure at lattice point designated by E quantity defined by Equation 16 mole fraction of ith component in liquid vector of lattice indices liquid phase activity coefficient of kth component lattice interval elementary vector defined by Equation 3 vapor phase activity coefficient of kth component
literature Cited (1) Barker, J. A , Australian J . Chem. 6 , 207 (1953). (2) Bruce, G. H., et al., Trans. Am. Inst. Mining, Met. Petrol. Engrr. 198. 79 (1953). (3) Hala, 'E., Pick, J., Fried, V., Vilim, O., "Vapor-Liquid Equilibrium," Pergamon, New York, 1958. (4) Ljunglin, J. J., Van Ness, H. C., Chem. Eng. Sci. 17, 531 (1962). (5) Scatchard, G., Raymond, C. L., J . Am. Chem. Soc. 60, 1278
(1938). (6)' Sevkrns, \Y. H., ef ai., '4.I.Ch.E. J . 1, 401 (1955). (7) Tao, L. C.. Ind. Eng. Chem. 53, 307 (1961). (8) Torgeson, R. L., "Calculation of Vapor-Liquid Equilibria at Constant Temperature," M.S. thesis, bniver&ty of -Delaware, 1965. RECEIVED for review November 9, 1964 ACCEPTED July 20, 1965
GROUP CONTRIBUTIONS T O ACTIVITY COEFFICIIENTS The CH, and OH Groups W I LL IA M A
.
SC H ELL
E R , The Unnrersity o f x e b r a s k a , Lincoln, A'eb.
The group contribution work reported by Wilson and Deal i s modified and extended. The theory i s developed in terms of thermodynamic instead of empirical quantities. Extensive data for the CH2 and O H group contribultions to the residual partial molal energy of mixing are reported over a CH2 concentration range of 3.0 tlo 99.9%. A method for calculating the standard state group contributions needed for applying the method is presented. Additional information for using the method with diols and triols is provided. Logarithms of activity coefficients for eight binary systems near 50" C. were calculated from group contributions with an absolute average deviation from the experimental values of less than 10%.
FOR
process evaluation and design calculations it is highly desirable to account for solution nonidealities because of their effect on equilibrhm composition and equipment requirements. However, it is not always feasible to determine activity coefficients experimentally for this purpose. As a result attempts have been made to develop generalized methods
for predicting activity coefficients, especially for solutions of nonelectrolytes. Hildebrand's method (5) for regular solutions requires the use of solubility parameters. These can be calculated from heats of vaporization. In the method described by Pierotti, Deal, and Derr ( 7 7 ) , activity coefficients at infinite dilution VOL. 4
NO. 4
NOVEMBER 1 9 6 5
459
