586
Ind. Eng. Chem. Process Des. Dev. 1984,
1436611, Unlverslty of Oklahoma, Norman, OK, 1980. Swanson, A. C.; Chueh, C. F. “Estimation of Liquid Heat Capacity by the Halcon Physlcal Data System,“ Proc. 85th Annual A l C M Meeting, New York, Nov 1972. Taylor, R. D.:Johnson, B. H.; Kilpatrick, J. E. J. Cbem. phvs. 1856. 23, 1225. Vick, G. K.;Epperly, W. R. Sclence 1962, 217(4557), 311.
23,586-589
Watson, K. M.; Nelson. E. F. Ind. Eng. Chem. 1833, 25, 880. Wilson. G. M.; Johnston, R. H.; Hwang, S. C.; Tsonopoulos, C. I n d . Eng. Che” Process Des. D e v . 1861, 20,94.
Receiued for reuiew April 27, 1983 Accepted October 7, 1983
Prediction of Multicomponent Vapor-Liquid Equilibrium with the Wilson Equation: Effect of the Mlnimization Function and of the Quality of Binary Data Norman Silverman and Pimltrios T a d o s New Jersey Institute of Tschnology, Newark, New Jersey 07102
Ten objective functions are used in the regresskm of the vapor-liquid equillbrium WLE) data for 247 binary systems with the Wilson equation and the obtained parameters are used in the predlctlon of the VLE behavior of 73 mulucMnponent systems. All models except model 1 give good results with model 5 being the best one. The resutts also suggest that the accuracy of prediction is not very sensitive to the quality at the binary data involved.
Introduction The prediction of multicomponent vapor-liquid equilibrium behavior is very important in industrial applications, for while such applications call in the typical case for multicomponent data most experimental work has dealt with binary systems. Therefore, several expressions have been presented through the years that attempt to utilize the available binary experimental VLE data to predict multicomponent behavior. Greatest success has been achieved with the more recent expression: Wilson (1964); NRTL (Renon and Prausnitz, 1968);LEMF (Marina and Tassios, 1973), and UNIQUAC (Abrams and Prausnitz, 1975). Among these models, the Wilson equation has found wide application, both in the correlation of binary data and in the prediction of multicomponent VLE behavior. Its popularity is based upon ita simplicity of expression coupled with accuracy comparable to that obtained by the other more complex models. Prediction of multicomponent VLE behavior with the Wilson equation, as well as with the others listed above, involves two basic steps: first, regression of the available binary experimental VLE data against some objective (minimization) function to establish the values of the binary parameters; second, prediction of the multicomponent VLE behavior using these binary parameter values. Several objective functions have been proposed through the years for the regression of binary data (Prausnitz et al., 1967; Nagahama et ai., 1971; etc.). W e the fact that the use of different objective functions leads to different values of the binary correlating parameters has been recognized, their impact on the accuracy of the resultant multicomponent predictions has attracted limited attention (Brinkman et al., 1974). It is the first objective of this study to examine the effect of the binary parameter values obtained from ten different objective functions (Table I) on the accuracy of the predicted multicomponent VLE behavior. For this purpose, the Wilson equation is used along with VLE data for 247 binary and 73 multicomponent systems (Silver”, 1978). Most of the multicomponent systems involve three com-
ponenta, and the remaining ones involve four and five. In most prediction studies, especially where different models are compared, the binary data are screened so that only those of acceptable quality are used. Yet the design engineer may have an occasion to use binary data of poorer quality. It is the second objective of this study to examine the effect of the quality of the binary data on the accuracy of the predicted multicomponent VLE behavior. The Wilson Equation For a binary system, the Wilson equation expresses the activity coefficients y1 and y2 in terms of two adjustable parameters, (Alz - A& and (Azl - All), which are considered temperature independent over a modest temperature range
and
where
and u k L is the liquid molar volume, component k, cm3/gmol. For a system containing m components, the activity coefficient of constituent i is given by
The VLE behavior, hence, of a multicomponent system can be expressed in terms of binary parameters, obtained form binary data alone, without the need for multicom-
Q196-43Q5/84l1123-Q586$Q1.5QlQ 0 1984 American Chemical Society
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984
587
Table I. Objective Functions
model 3 model 4
model 5
model 10
ponent interaction parameters. Evaluation of The Experimental Data Activity coefficients were calculated by the method of Prausnitz et al. (1967). The quality of the binary data was examined by using the area test. For this purpose the calculated activity coefficients were fitted to the following polynomial In (y1/y2) = A
+ Bxl + CrI2
(4)
and a consistency index C.I. C.I. =
lApl- lAnl lApl+ lAnl
where A, = positive area, above the x axis, and A,, = negative area, below the x axis, was calculated. It should be noticed that for some very complex systems more parameters may be needed than the three involved in eq 4 (Prausnitz, 1969). Evaluation of a larger number of parameters, however, requires a large and very accurate set of experimental measurements. Such data are not available in the typical case. Quation 4 was used therefore for the systems in this study. For isothermal data the C.I. value depends of the volume of mixing which is small. Hence Prausnitz (1969) suggest that for good quality data, C.I. 5 0.02. For isobamic data, on the other hand, the C.I. value depends on the heat of mixing and it can be substantial. For these systems the Herington (1951) test was applied. According to this test, the factor J is calculated 1.5lAq J.:(6) Tmin where AT is the difference between the maximum and minimum temperatures observed acrose the concentration range and Tmhis the minimum temperature in K. Data are considered of acceptable quality if ((3.1.- J) C 0.1. It should be emphasized, however, that this criterion is empirical and that both tests represent necessary but not
Table 11. Classification of Binary Data According t o Area Thermodynamic Consistency Test isothermal, isobaric, quality C.I.x 100 (C.1.-J ) x 100 no. good 95 4 10 147 10-20 47 fair 5-15 >15 > 20 37 poor 16 unacceptable area test not possible total 247
sufficient conditions. The results of these tests are presented by Silverman (1978) and are summarized in terms of their quality in Table 11. Very few isothermal systems had a value of C.I. < 0.02, but 59% of all binary data used can be considered of good quality: isothermal, C.I. < 0.05; isobamic, (C.I. - J) < 0.10; and an additional 20% of fair quality. On the other hand, the results for 15% were poor while for the remaining 6% the scattering of the data was such that it was not possible to identify a positive and a negative area. They are classified as unacceptable. All binary systems, however, were used for multicomponent predictions. No attempt was made to evaluate the quality of the multicomponent data. Quation 4 was also used to evaluate the infinite dilution activity coefficients In ylm = A In y2- = -(A
+ B + C)
(74 (7b)
Evaluation of the Binary Parameters In theory the pair of Wilson equation parameters for a binary system can be evaluated from a single experimental point. Hudson and Van Winkle (1969), for example, showed that parameter values obtained from one experimental point around the middle of the concentration range generally lead to successful correlation of the binary data and reliable prediction of multicomponent VLE behavior. In practice, however, the parameters are estimated by regression all available binary data points against the Wilson equation until some objective function is mini-
588
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984
Table 111. Typical Variation in Wilson Parameters acetone-chloroform acetone-acetonitrile model no. AI2 - A22 AI, - A21 All - All A12 -- A , , -374.41 124.13 -34.11 -86.42 -424.4 8 24.98 31.94 43.67 -394.61 -12.91 55.41 1.59 -34.41 -376.3 5 40.55 21.42 58.35 -448.15 -118.24 274.60 74.45 -133.20 -4 60 .O 0 302.91 -456.54 31.28 -140.89 315.67 -464.3 6 -92.97 87.15 227.50
mized. The objective function is a statement of the residual between the experimental and calculated values of some variable of interest. Starting with an estimate of the parameter values, the regression subroutine varies them until a minimum value of the objective function is obtained. Both Nagahama et al. (1971) and Verhoeye (1970) have shown that the parameter estimates are insensitive to the manner in which they are varied in the regression routine. They are, however, very sensitive to the form of the objective function. Ten objective functions are considered in this study and are presented in Table I. The first involves use of the infinite dilution activity coefficients calculated through eq 7a and 7b. Shreiber and Eckert (1971) and Hankinson et al. (1972) have shown that reasonably good prediction of multicomponent VLE behavior can be achieved by using this model. Model 2 involves the Gibbs excess energy and has been used, among others, by Nagahama et al. (1971) and Verhoye (1970). Models 3 and 4 involve the difference and relative difference between experimental and calculated activity coefficients. Several authors (e.g., Nagahama et al. (1971), Nankinson et ai. (1972), Marina and Tassios (1973)) have considered these models. Model 5 considers the residuals in both the vapor phase composition and total pressure. The relative error in pressure in used to obtain values that are of the same order of magnitude as those in the composition term. Renon and Prausnitz (1968) have used this model. A bubble point pressure calculation, using the corresponding experimental temperature and liquid phase composition values for both isothermal and isobaric data, was involved in the calculation. Convergence is assumed when the change between two consecutive iterations, in both the relative pressure and the sum of the vapor mole fractions, are each less than Models 6 and 7 are structurally the same in that they both measure the residuals in the calculated vapor phase composition. Model 6 uses the aforementioned bubble point pressure subroutine and has been considered, among others, by Nagahama et al. (1971). Model 7, however, represents a variation of model 6 in which the experimental value of the pressure, in addition to temperature and liquid composition, is used instead of the value calculated from the bubble point pressure subroutine. Convergence of this latter routine is reached when two consecutive iterations give values of each vapor mole fraction that do not differ by more than While this approach has the serious disadvantage that the converged value of (yl + yz) may not be equal to unity, it utilizes more experimental information than model 6. It can normally be expected that without the normalization step in the subroutine the sum of the mole fractions will not be unity. Such cases resulted in parameter values that different significantly from those obtained by model 6. Similarly, models 8 and 9 are structurally the same in that they both measure the relative error in total pressure. While model 8 uses the bubble point pressure subroutine to calculate the total pressure from the experimental T-x,
hexane-benzene - All 221.94 162.80 243.56 253.80 242.44 297.03 338.67 97.59
A12
-6001
heptane-toluene
- A22 131.40 194.06 131.66 125.96 135.68 104.68 75.78 232.53
A,,
A,, - A l l
126.12 -19.35 66.92 51.02 -21.01 -37.33 -37.93 8.23
135.35 255.12 185.15 194.14 257.80 271.32 272.15 233.68
A12 -
A12
' " l l ' " ' ' ' ~ ' ' T - - T - - -
-803 -KO -400 -200 .O
-200 400
600
1
803 000
A/? - A / /
Figure 1. Standard deviation in activity coefficients for various Wilson parameters: methyl acetate (1)-benzene (2) (Nagata, 1962).
data, model 9 uses a modification of the subroutine that utilizes, in addition, the experimental yi data. Convergence is assumed when two consecutive iterations give relative errors in total pressure that do not differ by more than lo5. Although the calculational procedure for models 8 and 9 differ somewhat, very close parameter estimates were obtained in all cases. Note that model 8 has found extensive use for it used in the monograph of Prausnitz et al. (1967). Finally, model 10 is similar to model 3 in that it uses the difference between the experimental and calculated activity coefficients. The error terms, however, are weighted by a pair of factors reflecting the deviation of the experiment activity coefficients from unity. Unless the weighting factor term is the same for both components, use of this model has the effect of skewing the search in the three-space defined by the Wilson parameters and the objective function relative to that of model three. This method was used by Hudson and Van Winkle (1969). The values of the binary parameters for all binary systems and models are given by Silverman (1978). Typical results are presented in Table I11 and indicate that distinctly different parameter values may be obtained from different objective functions. For systems with negative deviations from Raoult's law, where as many as three pairs of parameters are possible, the smallest in absolute value set was used for, as shown by Silverman and Tassios (1977),it generally gives the best results in the correlation of binary data and the prediction of multicomponent VLE behavior. The variation in the parameter values of Table I11 can be explained with the help of the contour map for the system methyl acetate-benzene presented in Figure 1 (Larson et al., 1971). The curves represent constant values of the objective function Q, model 4. It can be seen that substantially different parameter pairs can give equally good fit of the experimental data. For example, the pairs (590, -350) and (-170, 350) give Q values of 0.05. In essence, Figure 1 demonstrates that the parameters are correlatable and such behavior is common to all models
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984
Table V. Effect of the Quality of the Binary Data on the Prediction of Multicomponent Vapor Liquid Eauilibrium (Model 5 ) group no. quality of binary data overall Ay x 1000 10.2 A 13 all good 9.1 B 11 mixture of good and fair C 23 mixture of good and fair 9.2 and one poor D 26 one in unacceptable 12.5
Table IV. Overall Average Absolute Deviation in Vapor Phase Composition (Ay ) and Total Pressure (AP) model ay x 1000 u,mmHg 1 2 3 4 5 6 7 8 9 10
19.2 12.3 11.4 10.9 10.6 11.6 12.5 10.7 10.7 11.5
33.9 16.5 16.8 15.3 14.4 21.0 18.1 14.5 14.5 17.2
where the parameers reflect energy differences between pairs of molecules. Parameter evaluation could also be carried out by the principle of maximum likelihood (Prausnitz et al., 1980) but the appropriate subroutines were not available at the time of this study. It has been reported recently (Skjold-Jorgensen, 1983), however, that such an approach should not be expected to yield better multicomponent predictions than conventional methods. Results a n d Discussion Vapor phase concentrations and total pressures were calculated for the multicomponent systems of this study from the experimental values of liquid phase concentrations and temperatures following the procedure of Prausnitz et al. (1967) and using the binary parameter values and computer programs given by Silverman (1978). Average absolute deviations between calculated and experimental vapor phase concentrations (Ay) and total pressure (AP)for each of the 73 multicomponent systems and 10 models were calculated as m m
m
AP = {ClpJcalcd) - Pi (exptl)l)/n
589
(9)
I
where m is the number of components and n is the number of experimental data points for the multicomponent system. Detailed results for all systems and models are given by Silverman (1978) and are summarized in Table IV. The poor results of model 1 reflect the uncertainty involved in determining infinite dilution activity coefficients. The remaining models provide reasonably good predictions of multicomponent VLE behavior. Best resulta are obtained with models 5, 8, and 9 which require a computer time consuming bubble point calculation. This can be avoided by using model 4 where the relative difference between calculated and experimental activity coefficients is minimized. The use of experimental data, in excess to those required from the bubble point calculations -models 7 and 9-while tempting, does not provide better results. Obviously the relaxation of the requirement that the sum of the calculated y i s is equal to one negates any positive effect that these additional information has on the binary parameters. The effect of the quality of the binary data on the accuracy of the predicted vapor phase concentrations is presented in Table V and suggests no negative impact until one of the binaries involved in the prediction is of unacceptable quality. The overall error in vapor phase concentration for groups A, B, and C involving 47 multicomponent systems is 0.009, which can be considered very good, especially if the uncertainty in the multicomponent data is taken into account. It appears that the interdependency between the Wilson equation parameters, as
demonstrated by its one-parameter form, tends to remove some of the uncertainty in the binary data. Even for group D, the overall error is 0.013 and average errors over 0.025 were observed only twice, with the worst case of 0.032. Conclusions Good prediction of multicomponent vapor liquid equilibria from binary data is obtained from the Wilson equation with several minimization functions in the regression routine. Best results are obtained with models 5, 8, and 9 that require a computer time consuming calculation. This can be avoided, however, by using model 4, which provides practically equally good performance. The accuracy of prediction, probably because of the interdependence of the two Wilson parameters, does not appear to be very sensitive to the quality of the binary data. The effect becomes substantial only in the presence of a binary system where the scattering of the data is such that it is impossible to determine a positive and negative area in the consistency test. For 47 multicomponent systems that do not involve such binaries (64% of the total), the overall error in y is 0.009 and in total pressure 12" Hg (model 5). Nomenclature C.I. = thermodynamic consistency index m = number of components n = number of experimental points P = total pressure, mmHg R = gas constant, 1.987 cal/g-mol K T = temperature, K vkL= liquid molar volume of species k, cm3/g-mol x k = liquid phase mole fraction, component k Yk = vapor phase mole fraction, component k Y k = liquid phase activity coefficient, component k Ykw = infinite dilution activity coefficient, component k Literature Cited Abrams, D. S.; Prausnltz, J. AIChEJ. 1975, 2 1 , 116. Brinkman, N. D.; Tao, L. C.; Weber, J. H. Ind. Eng. Chem. Fundem. 1974,
13, 156. Hankinson, R. W.; Longfm, 6.; Tassios, D. Can. J . Chem. Eng. 1972, 50,
511. Herrington, E. F. G. T . Inst. Pepo/. 1951, 37, 457. Hudson, J. W.; Van Wlnkle, M. J . Chem. Eng. Date 1969, 1 4 , 310. Larson, D.; Sinkoff, S.; Tassios, D. l6lst Meeting of the American Chemical Society, Los Angeles, CA. March 1971. Marina, J.; Tassios, D. P. I d . €ng. Chem. Process D e s . D e v . 1973, 12,
67. Nagahama. K., Juzuki, I.; Hirata, U. J . Chem. Eng. Jpn. 1971, 4 , 1. Preusnitz, J.; Eckert, C. A.; Orye. R. V.; O'Connell, J. P. "Computer Calcuiatlons for Muhcomponent Vaporliquid Equilibria"; RenticeSlall; Englewood Cliffs, NJ, 1967. Prausnitz, J. "Molecular Thermodynamics of hM-Phase Equilibria"; RenticeHall; Englewood Cllffs, NJ, 1969. Prausnitz, J.; Anderson, T.; Grens, E.; Eckert. C.; Hsieh, R.; O'Connell "Computer Calculations for Multicomponent Vapor-Liquid Equilibria"; Prentice-Hall; Englewood Cliffs, NJ, 1980. Renon, H.; Prausnltz, J. M. A I & € J. 1968, 1 4 , 135. Shrelber, L. B.; Eckert C. A. Ind. Eng. Chem. Process Des. Dev. 1971, 70,
572. Silverman, N.; Tassios. D. Ind. Eng. Chem. Process Des. D e v . 1977, 73,
13. Silverman, N. M.S. Thesis, New Jersey Institute of Technology, Newark, NJ,
1976. Verhoeye, L. A. Chem. Eng. Scl. 1870, 25, 1903. Wilson, G. M. J . Am. Chem. SOC.1964, 86, 127.
Received for review June 3, 1982 Accepted October 24,1983
