Znd. Eng. Chem. Res. 1991,30, 524-532
524
it remembers either the initial bracketing interval or some updated smaller bracket for the root. The full algorithm follows, in a "pseudoprogram" format. The subroutine "brent" uses its arguments to bracket the root on the first call and thereafter updates its first argument to provide a new estimate of the root.
START: comment Find R, such that 3(R,,RR,)= 0 comment Root bracket R, E (Rl,l) Re := R1 lOOP call brent (R,,3(Re,R,),l,3(l,l)) until 3(R,,RR,)E 0 comment Find yield surfaces such that 3(R-,R+) = 0, S(R-8,) = 0 comment Outer root bracket is R- E (Rl,RR,-t)
R- := R1 loop comment Inner bracket is R+ E (RJ) R+ := Rloop call brent(R+,S(R-,R+),1,3(R-,l)) until 3(R-,R+)= 0 call brent(R-,O(R-,R+),R. - c,O(RI-e,RR.)) until $'(R-,R+) = 0
END. Literature Cited Anshus, B. E. Bingham plastic flow in annuli. Ind. Eng. Chem. Fundam. 1974,13,2,162-4. Bird, B. R.; Dai, G. C.; Yarusso, B. J. The rheology and flow of viscoplastic materials. Rev. Chem. Eng. 1983,1 (l),1-70. Brent, R. Algorithms for Minimisation without Derivatives; Prentice-Hall: Englewood Cliffs, NJ, 1973.
Casson, N. A flow equation for pigment oil suspensions of the printing ink type. In Rheology of Disperse Systems; Mill, C. C., Ed.; Pergamon: Oxford, UK, 1969. Dzuy, N. Q.; Boger, D. V. Yield stress measurements with the vane method. J. Rheol. 1985,29 (31,335-47. Frederickson, A. G.; Bird, R. B. Non-Newtonian flow in annuli. Ind. Eng. Chem. 1958,50(31,347-352. Froishteter, G. B.; Vinogadov, G. V. The laminar flow of plastic disperse systems in circular tubes. Rheol. Acta 1980,19,239-50. Grinchik, I. P.; Kim, A. Kh. Axial flow of a non-linear viscoplastic fluid through cylindrical pipes. J. Eng. Phys. 1974,23,1039-41 (published in Russian 1972). Hanks, R. W. The axial laminar flow of yield-pseudoplasticfluids in a concentric annulus. Ind. Eng. Chem. Process Des. Dev. 1979, 18 (3), 488-93. Herschel, W. H.; Bulkley, R. Konsistenzmessungen von GummiBenzollosungen. Kolloid 2. 1926,39,291-300. Jaisinghani, R. Annular flow of a Casson fluid. M.S. thesis, University of Wisconsin-Milwaukee, 1974. Kooijman, J. M.; van Zanten, D. C. The flow of a Cassonian fluid through parallel-plate channels and through cylindrical tubes. Chem. Eng. J. 1972,4,185-8. Lih, M. M A . Transport Phenomena in Medicine and Biology; Wiley: New York, 1975;Chapter 9. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vettering, W. T. Root-finding and non-linear sets of equations. In Numerical Recipes: the art of scientific computing; Cambridge University Press: Cambridge, UK, 1986,Chapter 5,see especially pp 251-4. Robertson, R. E.; Stiff, H. A. An improved mathematical model for relating shear stress to shear rate in drilling fluids and cement slurries. SOC.Pet. Eng. J. 1976,Feb, 31-36. Shah, V. L. Blood flow. Adu. Tramp. Processes 1980,1, 1-57. Shul'man, Z. P. Calculation of a laminar axial flow of a non-linear viscoplastic medium in an annular channel. J. Eng. Phys. 1973, 19, 1283-9 (published in Russian 1970).
Received for review February 21, 1990 Revised manuscript receioed August 24, 1990 Accepted September 12, 1990
Applications of a Generalized Equation of State for Associating Mixtures S. Jayaraman Suresh and J. Richard Elliott, Jr.* Chemical Engineering Department, University of Akron, Akron, Ohio 44325-3906
A generalized equation of state that explicitly recognizes self-association and cross-association has been evaluated for systems with two associating species a t high pressures as well as low pressures. Forty-five systems were evaluated for components ranging from water to propanol. The Soave equation of state and the UNIFAC model were used for comparison. Percent average absolute deviations in bubble-point pressure were 3.54 for the Soave equation, 2.53 for UNIFAC, and 1.62 for the association model. Generalized correlations of the pure-component parameters and a single binary interaction parameter were used for all models. Some limitations are discussed, but the results show that the reaction equilibrium approach is a clear improvement over the Soave equation of state and competitive with UNIFAC. Introduction For many years, one of the most challenging problems in thermodynamics hae been prediction of phase equilibria in hydrogen-bonding systems. These systems are commonly encountered in the chemical industry, and the ability to make quantitative predictions would facilitate optimization of many commercial processes. Previous successful treatments of these mixtures have used activity coefficient methods. The Van Laar equation, the Wilson (1964) equation, and UNIFAC (Fredenslund et al., 1975) are examples of significant developments by this approach. The Van Laar and the Wilson equations are usually considered to be correlative models whereas
the UNIFAC model is considered to be a predictive model. Nevertheless, we have implemented UNIFAC to evaluate the activity coefficients in the liquid phase only and used the Soave equation (1972) to evaluate the fugacity coefficients in the vapor phase. The Soave equation for the vapor phase requires an adjustable binary interaction coefficient. Hence, our implementation of UNIFAC in conjunction with the Soave equation becomes more a correlation with one adjustable parameter. This treatment is typical of the way these equations are used in process simulators like ASPEN. Many sets of UNIFAC parameters have been compiled; therefore the specific set of UNIFAC parameters was taken from Reid et al. (1987) for
088&5885/ 91/2630-0624$02.50/ 0 0 1991 American Chemical Society
Ind. Eng. Chem. Res., Vol. 30, No. 3,1991 525
our study. Even though these activity coefficient models have been quite successful, they have had some limitations. For example, application of thew methods to high-pressure vapol-liquid equilibria has led to large errors (Gupte and Daubert, 1986). Furthermore, application of these methods to supercritical components is often difficult. One reaponse to the difficuties encountered with activity coefficient methods has been to use an equation of state with sophisticated mixing rules. Then the fugacities of both phases are determined from the equation of state. Examples of contributions to this approach are papers by Whiting and Prausnitz (1982), Gupte and Daubert (1986), Schwartzentruber and Renon (1989), and Lee and Chao (1989). These models attribute anomalies in polar mixtures to peculiarities in the local composition behavior. Some of these models have provided reasonable representation of phase behavior in polar as well as nonpolar mixtures, but none has been identified as being the single best general method at the current time. Since no clearly superior method is available, we have implemented the original Soave (1972) equation as our basis for comparison to equation of state methods. Although many modifications of the Soave equation have been proposed, especially for polar systems, most modifications are for specific cases and may not provide a basis for comparison that a broad audience can appreciate. We describe the impact that typical modifications might have in the discussion below. As an alternative to the above approaches, we have recently been studying equations of state that describe polar mixtures by an association model. Previous works in this area include papers by Heidemann and Prausnitz (19761, Gmehling et al. (1979), Ikonomou and Donohue (19881, Anderko (19891, and our previous paper (Elliott et al., 1990). The physical basis for these association models can be illustrated by considering the marked rise in volatility when mixing alcohols into hydrocarbons. In the context of an association model, the hydrocarbons tend to break the hydrogen-bonding network of the alcohol, leaving the alcohol molecules relatively uninhibited in their escaping tendency. This breaking of the association network appears directly in the model. A model based solely on nonassociative interactions, like the Soave (1972) equation, cannot represent this behavior accurately. The emphasis of our work has been to develop a genr n a l i z e d equation of state with broad applicability and few adjustable parameters. The initial development of the equation of state for pure associating compounds and mixtures containing one associating species with a nonassociating diluent was performed by Elliott et al. (1990). That work showed that our theory compares favorably with available molecular simulation results in addition to describing real fluids accurately. This paper deals with extension of the previous work to mixtures containing two associating compounds over a range of higher pressures than in the previous work. Results are presented to show how the association model performs in comparison to UNIFAC and a standard nonassociating model, the Soave equation (1972). Generalized methods are used for purecomponent properties and mixtures are analyzed by obtaining a single binary interaction coefficient to represent the unlike pair interactions.
adding a nonassociating diluent tends to break up the association network in such a way that the volatility of the mixture is difficult to describe via a model based solely on physical interactions. In this work, we examine extension of this approach to more than one associating species. Unfortunately, the complexity of the mixtures formed when two species are mutually associating causes a truly rigorous analysis to be beyond our current capabilities. For example, consider the formation of pentamers that are three of one species and two of another species simultaneously with formation of decamers that are nine of one species and one of the other. Then, some of these species may form cyclic structures. Clearly, this is a complex process that will require considerable fundamental study to obtain a completely rigorous theory. Nevertheless, some approximate results have been developed by considering limiting cases that can be solved analytically (Ikonomou and Donohue, 1988; Anderko, 1989). Despite some limitations, the basic concept of association reducing the fugacity is still explicitly recognized, and this represents an improvement over a model based on entirely physical interactions. In our analysis, we have evaluated the approximate solutions proposed in both of the above references. The equation of state that was applied to pure associating species and mixtures of one associating species with a nonassociating diluent (Elliott et al., 1990) is
Theory When molecules associate to form dimers, trimers, etc., it is clear that the net number of moles ie reduced. The most significant effect of this association is the suppression of the fugacity of associating species relative to nonassociating species with equivalent size, shape, and attractive dispersion forces. In our previous work, we ehowed how
( v Y ) ~ (qtlY)o/Cxiqi
where c = a parameter accounting for shape effect on repulsive contributions q
E
1.0 + 1.90476(~- 1)
where q is a parameter that accounts for shape effects on attractive contributions NT true number of molecules after accounting for association No 1 number of molecules if there had not been any association 00
(No/V)Cxiui*
where u* is the volume of a monomer molecule ( cv ) o
(qqY)o Yij
where
cij
(No/ V) C C (cu*)ijxixj (No/V)CC(qu*)xixjYij
= exp(8tij) - 1.0617
is the potential well depth (CV*)ij = (CiUj* + Cj4*)/2 (qv*)ij = (qivj*
+ QjUi*)/2
where kij is the binary interaction parameter The subscript "0"implies that these quantities are calculated on the basis of the monomer parameters. The generalized values of c, elk, and u* are given for nonassociating species by c = 1.0 + 3.535w + a533w2 (2)
526 Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991
+ 0.945(~- 1) + 0.134(~- 1)' k Tcl.023 + 2.225(c - 1) + 0.478(c - 1)2 (3) kTc 0.0312 + 0.087(~- 1) + 0.008(~-1)2 u* = (4) fc' 1.OOOO + 2.455(~- 1) + 0.732(~- 1)' -=
1.OOO
where k is Boltzmann's constant. For nonassociating components, (1)is a cubic equation and hence is easy to solve for the roots of 2. For associating species, the quantity that needs to be determined from the reaction equilibrium method is NT/No,the ratio of the true number of moles to the apparent number of moles. For two associating species, we can assume three general equations for the reaction equilibria.
+ iAl = Ai+l
(5)
B1 + jB1 = Bj+1
(6)
AI Ai
+ Bj
AiBj
(7)
We assume that i and j vary from one to infinity. We also assume that the equilibrium constants of reactions of a particular type, reaction 5 for example, are all equal. Hence, K,,K2,and KI2are the equilibrium constants for reactions 5,6, and 7, respectively. K1 and K 2 are determined from the pure-component data and are not adjustable. In general, KI2must be determined from mixture data and that would increase the number of adjustable parameters. We have avoided the introduction of adjustable parameters, however, by implementing several approximate methods that do not require explicit knowledge of K12. These methods are described below. Following the approximate method outlined by Ikonomou and Donohue (1988), we obtain the relations
Table I. Valuer of KO*, a,c, t/k (K),and v* (cm8/molecule) for Some Arrociating Compoundr componenta KO* H c elk u* (1029) 1.727 methanol 22.2477 -3.25 2.3491 197.01 ethanol 11.2131 -3.35 2.5000 206.75 2.620 1-propanol 9.7953 -1.00 3.9412 195.30 2.989 2-propanol 13.5469 -1.70 3.7972 182.55 2.891 2-propanone 14.3122 -1.60 2.2794 204.44 3.384 0.5090 -2.80 2.7035 264.66 1-butanol 4.992 water 18.8730 -2.16 2.1322 258.86 0.865 ammonia 15.5071 -1.60 2.0765 166.84 1.135
was obtained by using the 2,criteria, namely, Z//ZCNA = Zcws/Z,HoMoas described by Elliott et al. (1990). Table I lists all the equation of state parameters for the various compounds used in this study. Note that we have assumed a12to be given by (13) and this approximation eliminates the need for any additional cross-association parameter. While future work should provide better approximations, (13) provides a basis for initial consideration of multiple associations. Given the pure-component parameters, the three equations ((8),(9), and (10)) need to be solved for the three unknowns ( Wl, W2, and NT/No). Solution of the above system is quite complex. In particular, the double summation must generally be carried to 30 X 30 terms and this summation must be performed at each iteration. Even so, the effort might be worthwhile for preliminary developments if the results were highly accurate. Unfortunately, the above system of equations leads to the conclusion that a dilute associating compound does not associate, not even with other species. This conclusion is unrealistic and leads to high errors in applications to dilute associating mixtures. Fortunately, Ikonomou and Donohue (1988) developed alternative expressions based on the assumption that K12= (K1K2)1/2 and compared results of the numerical solution to several analytical forms. The analytid forms they developed were accurate at most conditions of interest for systems considered in this evaluation. Although these approximations are based on the numerical solution, they do not follow the same trend as the numerical solution in the dilute regions. Therefore, the analytical equations ((16)-( 18) below)
w, = (1 + (14x1 + 4a1)1/2)2
where
w, =
aI2 = Kci,M*
= KciRTci/(u*)O
( u * ) , = zui*xi
(13) (14) (15)
W1, W2 = ratio of the total number of moramers of the respective components to the total apparent number of moles in the mixture
4x2 (1 + (1 + 4a2)1/2)2
(16) (17)
provide an accurate description of phase behavior even in the dilute regions. This fortuitous result is probably due to cancellation of errors. One additional problem in the analysis of multiple associations is that (16)-(18) do not satisfy the limiting condition at which association of one of the components (for example, component 2) approaches zero. The exact solution at this limiting condition is
X I , x2
apparent mole fractions of the respective components Cij = (i + j ) ! / i ! j ! It is to be noted that the factor 1- l.stloin the denominator of a,, a2, and a12is a result of not assuming the ci aqua1 to ic,, which is different from the theory of Ikonomou and Donohue (1988). The value of Kc* for pure components
In order to correct the inconsistency of the analytical solution at the limit when one of the components does not
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 527 0.111
Table 11. Errors in Bubble-Point Pressure Using the Approximations of Ikonomou and Donohue (1988) and Anderko (1989)
%AAD system methanol l-butanol acetone + 2-propanol 2-propanol + 1-propanol ethanol + water methanol + ethanol methanol + water 2-propanol water ethanol 1-propanol methanol + 2-propanol methanol + l-propanol acetone + ethanol acetone + methanol methanol + benzeneo grand average
+
+
+
T,K 422.75-491.45 328.15 298.15 423.15-523.15 298.15-433.15 312.90-523.15 423.15-548.15 323.15-353.15 328.15 303.15-333.17 305.15-400.0 298.15-328.15 373.15-493.15
Ikonomou's Anderko's approxn approxn 0.67 1.48 0.75 1.60 0.89 1.01 1.92 2.13 0.88 0.82 2.20 1.71 1.93 1.66 0.93 1.50 0.89 1.01 1.59 1.03 2.14 3.48 1.50 10.22 2.10 2.10 1.62 2.67
a Since benzene is a diluent, both methods are equivalent for this system.
associate, Anderko (1989) formulated an alternate solution which corrects the trend at this limit ((21)-(23)). It can
- -NT
No
1
4x1 (1 + (1 + 4(aixi
+ a1f12))'/2)2
(21)
4x2 w, = (1+ (1 + 4(azxz + a12X1))'/2)2
(22)
w1=
2x1 + (1+ 4(alx1 + u ~ , x Z ) ) ~ / ' + 2x2
1 + (1 + 4(a,x,
+ al~X1))'/~(23)
be argued, however, that accurate description of this limit is not of much practical interest since an association site either exists or does not. Therefore, we have evaluated both methods and compared the results. All of the above formulas for multiple association represent approximate analyses of an extremely complex system. In some sense, it is ironic that the theoretical f wlysis of these systems is so complex, because the phase oehavior of two associating species like ethanol + methanol is often quite simple (Prausnitz et al., 1986). The presence of nonassociating species is what causes complicated phase behavior, and the results of previous work (Elliott et al., 1990) indicate that the effects of these diluents are well represented by this approach. These considerations indicate that even rough approximations can be expected to describe the behavior of a system like ethanol + methanol while the complex effects of diluents are already well represented. From this perspective, the practical significance of the above extensions to multiple associations is to make it possible to treat a broad class of multicomponent mixtures by the same theory. One might note, however, that future work will likely reveal mixtures for which more sophisticated descriptions of the multiple associations will be necessary. Results and Discussion Several binary systems were analyzed by use of the approximatione proposed by Ikonomou and Donohue (1988) and Anderko (1989) to evaluate which set of approximations was more accurate. The pure-component properties
0.0
0.2
0.4
0 . 6
0.8
1.0
X( 1 ) , Y ( 1 ) Figure 1. Pressure (MPa) vs mole fraction for ethanol (1) + 1propanol (2) at 353.15 K. ( 0 )Experimental. (-) Association model using Ikonomou and Donohue's (1988) approximation. (- - -) Association model using Anderko's (1989) approximation.
for this analysis are listed in Table I. The only adjustable parameter in the case of mixtures was the binary interaction coefficient, which was fit to the experimental data by minimizing the percent average absolute deviation ( % AAD) in the bubble-point pressures. Table I1 lists the errors in bubble-point pressure using both of these approximations. It can be seen that, for most of the systems, both the approximations are equally accurate. However, for certain systems, Ikonomou and Donohue's approximation is more accurate than Anderko's approximation. Moreover, it was found that, in most cases, the mole fractions in the vapor phase are more accurately represented by Ikonomou and Donohue's approximation than Anderko's. Figure 1 shows a comparison of these two approximations for the ethanol + l-propanol system when applied to the association model. Thus, our overall evaluation showed that the approximations proposed by Ikonomou and Donohue (1988) ((16)-( 18)) provided a better representation of phase behavior than the approximations proposed by Anderko (1989) ((21)-(23)). Consequently, the approximations proposed by Ikonomou and Donohue were adopted for our association model to develop results for various binary mixtures. A sample calculation to evaluate the fugacity coefficients for a binary mixture (methanol + ethanol) is provided in the Appendix. Resulta of the association model were compared to those of the Soave equation of state and the UNIFAC model. The Soave equation was implemented in its original form using binary interaction coefficients. The UNIFAC equation was applied to the liquid phase only and the Soave equation was applied to the vapor phase. The following equation applies for phase equilibrium via UNIFAC: y$yjCiL@,iaat
=
IVP x Y
@.
(24)
where yi = activity coefficient of the component in the liquid phase as calculated by UNIFAC; aiv= fugacity coefficient of the component in the vapor phase as calculated by the Soave equation; and aimt= fugacity coefficient of the pure saturated liquid, also calculated by the Soave equation. Thus the UNIFAC analysis included a binary interaction coefficient (which was constrained to be less
528 Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 3.0
0.85 2 . 5
6
n
0.75
4
4
n I2
a I
. 0
W
v0.65
a
a
1.5
0.55
1
.o 0.0
Figure 2. Preesure ( m a ) vs mole fraction of 2-propanol (1) + water (2) at 423.16 K. ( 0 )Experimental. (-) This work. (---) Soave equation, constrained to give a stable single phase. (---) UNIFAC.
0 . 2
0 . 4
0 . 6
0 . 8
1.0
X( 1) ,Y( 1 ) Figure 4. Pressure (MPa) vs mole fraction for methanol (1) + benzene (2) at 453.15 K. ( 0 )Experimental. (-) This work. (- - -) Soave equation, constrained to give a stable, single phase. (---I UNIFAC.
0.020
0.015 n
4
a I v
a 0,010
0.005
0.0
0 . 2
0.4
0 . 6
0.8
1.0
X ( 1) , Y ( 1 ) Figure 3. Pressure (MPa) vs mole fraction for methanol (1) + ethanol (2) at 298.15 K. ( 0 )Experimental. (-) This work. (- -) Soave equation, constrained to give a stable, single phase. (-'-) UNIFAC.
-
than unity) via the vapor phase. Including the binary interaction coefficient in this way was found to be significant for several systems. The errors in bubble-point pressure are listed in Table 111. The overall %AAD (percent average absolute deviation) for 45 systems was found to be 1.62 for the association model, 3.54 for the Soave equation of state, and 2.53 for the UNIFAC model. For all the systems, a single value of kij was obtained at both high and low pressures. In Figure 2, which represents the water + 2-propanol system, it is evident that the Soave equation is inaccurate near the azeotrope. Figures 3 and 4 and P-XY plots of methanol + ethanol and methanol + benzene, respectively. It is evident that the Soave equation performs reasonably well in these systems. The error in the Soave equation is mainly due to the error in the vapor pressures as shown in Figure 3, but the vapor pressure error in some systems is small and the Soave equation is still inaccurate, as shown in Figure 2. The
typical impact of moat modifications of the Soave equation is to reduce the errors in vapor pressure. In general, the phase behavior is relatively simple in the mixtures that the Soave equation describes accurately. Simple phase behavior is possible in some mixtures where cross-association is very similar to self-association. However, it should be noted that the Soave equation is unreliable for some mixtures (e.g., water + 2-propanol) and it is difficult to anticipate when it will be unreliable. Furthermore, phase behavior is more nonideal in mixtures containing a single associating species and an inert diluent, wherein the association model was clearly superior to the Soave equation (Elliott et al., 1990). This means that the association model performs consistently well in mixtures of single associating species as well as mixtures of more than one associating species whereas the Soave equation can occasionally be unreliable. As for the UNIFAC method, it tended to be more accurate than the association model at low pressures but it became less reliable at high pressures. The values of the binary interaction coefficient tended to be quite large in many of the UNIFAC evaluations. We limited their magnitude to unity, but the evaluations were insensitive to the precise value of the binary interaction coefficient in many cases. This is because the binary interaction coefficient has little impact at low pressures where the vapor phase behaves like an ideal gas. In other words, this approach attempts to describe shortcomings in the UNIFAC equation with a vapor-phase parameter. The binary interaction coefficient played a significant role for highpressure systems only. The association model was very consistent in its reliability over the entire range of conditions. Furthermore, the extension of UNIFAC to supercritical components, like H a in an oil recovery process, is difficult. Extension of the association model to such a system would be very straightforward. Even though many of these mixtures appear very simple, the value of the theory should not be underestimated since it consistently explains (1)the different shapes of the vapor pressure curves in polar pure fluids, (2)the anomalies of mixtures of associating species with diluenta, and (3)the "simplicity" of the multiple associating mixtures. This means that the same phenomenological approach appears
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 529 Table 111. Errors in Bubble-Point Pressure system methanol + n-butanol acetone + 2-propanol 2-propanol + 1-propanol ethanol + water methanol + ethanol
methanol + water
2-propanol + water
ethanol
+ 1-propanol
+
methanol 2-propanol methanol + 1-propanol
acetone
acetone
+ ethanol
+ methanol
methanol
+ benzene
grand average
T,K 422.75-491.45 328.15 298.15 423.15 473.15 523.15 298.15 373.15 393.15 413.15 433.15 312.90 373.15 423.15 473.15 523.15 318.20 423.15 473.15 523.15 548.15 323.15 333.15 343.15 353.15 328.15 303.15 313.15 323.15 333.15 333.17 305.15 313.15 321.15 328.15 375.00 298.15 308.15 318.15 323.15 328.15 373.15 413.15 453.15 493.15
association model k, %AAD 0.062 0.67 0.040 0.75 -0.005 2.05 0.035 2.01 1.94 1.81 0.002 0.96 1.11 0.68 0.95 0.68 -0.003 1.53 1.75 2.03 2.18 2.91 0.023 1.62 1.63 2.55 1.48 2.07 0.016 1.56 0.83 0.80 0.53 0.003 0.89 0.009 2.21 1.31 1.86 1.40 1.17 0.028 1.97 1.38 1.09 0.96 5.32 0.035 2.44 1.29 1.18 1.32 1.29 0.082 1.21 2.50 3.17 1.53 1.62
Soave eauation UNIFAC kij %AAD kij(Soave) %AAD 0.007 2.33 0.60 0.56 0.78 0.00 1.09 0.043 -0.017 1.92 0.00 3.00 -0.061 2.23 0.65 1.22 2.24 1.37 2.31 4.50 -0.010 13.98 0.74 4.07 2.04 1.63 1.47 2.41 0.44 2.04 0.32 1.81 0.20 -0.080 10.79 0.80 2.85 2.19 1.78 1.03 0.92 1.25 1.38 0.40 -0.190 15.02 0.80 3.22 6.40 2.33 1.81 5.44 3.06 4.70 1.25 4.05 0.015 -1.00 2.20 4.96 1.46 4.94 1.27 5.63 0.61 4.39 -0.013 1.94 3.64 1.00 1.20 1.00 0.002 8.37 1.17 5.48 4.31 0.69 2.75 2.60 2.86 2.75 2.63 0.31 0.026 1.00 1.79 1.61 1.41 1.71 1.63 2.84 6.40 7.17 3.54 0.012 1.64 -0.30 2.06 0.75 0.97 0.94 1.89 0.25 1.31 0.39 2.07 5.91 0.117 1.00 6.50 5.81 5.57 7.63 4.57 6.29 3.54 2.53
to be valuable and reliable in explaining all these types of mixtures. Conclusions
A generalized equation of state has been extended to mixtures of two associating components and compared to the resulta of the Soave equation and the UNIFAC model. Phase equilibria of several systems have been studied using a single binary interaction parameter in all models. 1. The association model is more reliable than the Soave equation for mixtures containing one or more associating species. For binary mixtures containing two associating species, the Soave equation performs reasonably well owing to the ‘simplicity” of most mixtures of this type but is occasionally unreliable even for some of these mixtures. 2. The association model is more accurate than UNIFAC at high pressures, but it is slightly less accurate than UNIFAC at low pressures. Whereas the UNIFAC model is unreliable at high pressures, the association model is consistently reliable over the entire range of conditions. 3. In general, the association model is consistently accurate for mixtures containing a whole range of compo-
reference Kav and Donaham (1955) Frishwater and Pike (1967) Polak et al. (1970) Barr-David and Dodge (1959) Hall et al. (1979) Butcher and Robinson (1966)
Griswold and Wong (1952)
Sada and Morisue (1975) Barr-David and Dodge (1959)
Udovenko and Frid (1948)
Freshwater and Pike (1967) Schmidt (1926)
Gordon and Hines (1946) Vinichenko and Susarev (1966) Gomez-Nieto (1977) Tamir et al. (1981) Marinichev and Susarev (1965) Abbott and Vanness (1975) Marinichev and Susarev (1965) Butcher and Medani (1968)
nents, polar as well as nonpolar, supercritical as well as subcritical. Future work will focus on improving the theory and evaluating it for enthalpies of mixing and vapor-liquidliquid behavior. While the Ikonomou and Donohue method was most accurate for the mixtures studied in this evaluation, other mixtures may require more sophisticated analysis of the multiple associations. Techniques for developing improvements will include molecular simulations of various chain lengths with varying numbers of association sites. Clearly, more effort needs to be directed at developing fundamental insights about these complex associations. Some detailed information is already available from spectroscopic analysis (cf. Prausnitz et al., 1986) and molecular simulations (Jackson et al., 1988),but more of this information is needed. Nevertheless, the success of the approximate theory presented here establishes the viability of a kind of “middle ground”. For engineering calculationsof phase equilibria, it is important to recognize roughly how the fugacity of a component is affected by associations and mixing. As detailed in the Appendix, knowing the precise shape and nature of the oligomers is not always necessary. In this context, further
530 Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991
detailed insighte into the chemistry of these complex associations can go hand in hand with engineering applications.
Acknowledgment
calculation was performed for the liquid root (2,= PV/ N&T = 0.00796) at the following conditione: T = 433.15 K; X1 = 0.942. The expression for the apparent fugacity coefficient of the monomer of component i in a binary mixture is given as follows: r
Financial support of this research by BP America is gratefully acknowledged.
Nomenclature a (for pure component) = association parameter = [qo/(l 1.9v0)1Kc*exp[H(l - l/TJI C- = a constant = (i + j)!/I!j! A% O = standard heat capacity change of hydrogen bond Pormation c = quantity representing shape effect on repulsive term of
z
AHo = standard free energy change of hydrogen bond for-
mation H = AHo/RT - ACpa/R = parameter characterizing temperature dependence of association K = reaction equilibrium constant of association reaction Kc* (for pure component) = K,RTc/u* k = Boltzmann's constant kij = binary interaction parameter k3 = a constant = 1.90476 No = number of moles if there had not been any association NT = number of moles if there was association P = pressure q = 1+ k3(c - 1)= quantity representing the shape effect on the attractive term of 2 R = gas constant T = temperature V = total volume u* = characteristic size parameter u = specific volume W = ratio of the total number of monomers of a particular component to the total apparent number of moles in the mixture x = mole fraction Y = exp(c/kT) - 1.0617 = attractive energy parameter 2 = compressibility factor Greek Symbols c
= potential energy well depth
where
@ = fugacity coefficient y = activity coefficient 7 = reduced density o = acentric factor u
(A-2)
= effective molecular diameter
(A-3)
Subscripts
c = critical property i = ith associating species M = mixture property r = reduced property T = quantity calculated after accounting for association 0 = quantity independent of association Superscripts
L = liquid phase V = vapor phase sat = saturated EOS = equation of state HOMO = homomorph A = associating component NA = nonassociating component
Appendix This sample calculation is provided to illustrate the application of the equation of state to evaluate the fugacity coefficienta of the methanol (1)+ ethanol (2) system. The
(A-5)
- PV = - + - - NT NokT No
4(c7)0 1 - 1.970
z,,,(qsy)o
1 + ki(7Y)o
(A-6)
The constants in the equation of state are k, = 1.7745; k, = 1.0617 k3 = 1.90476; 2, 9.49 The pure-component parameters given by Elliott et al. (1990) are ~1 = 2.349; C, = 2.500 all/k = 197.01 K; ~ 2 2 / k = 206.75 K
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 631 ul* = 1.73 X
cm3/molecule; u2* = 2.62 X cm3/molecule
Note: ZcEosis the value of critical compressibility factor assumed by the equation of state. It is obtained from the pure component analysis of the equation of state (Elliott et al., 1990). (1 - k i j ) (197.01 X 206.75)'/'(1 - 0.002)
201.42 K
( u * ) ~= &*xi = 0.942(1.73 X + = 1.78 X lbZ3 cm3/molecule 0.058(2.62 X qi = 1
+ k3(ci - 11,
which gives q1 = 3.5695; q 2 = 3.858 qiuj
+ qjui 2
, which then provides
(qu*)11 =
6.17 X cm3/molecule; (qu*)22 =10.1 X lo-" cm3/molecule; ( Q U * ) ~=~ 8.01 X cm3/molecule
Yij = exp(flcij) - k2, hence Y11 = 0.5888; Y22 = 0.6301; Y12 = 0.6075 (klYu*)o= klCYiiui*xi= 1.87 X
cm3/molecule
(ZmqYu*)0= Z m C C ( q U * ) i j Y i j z i z j = 6.75 X Po =
P
No
- 1.699 X 'iT = WT 7=
P O ( U * ) ~=
1-
2)
u2*
= cm3/molecule
Finally, on substituting these values into the fugacity relation (A-1)
ZclEOs= 0.271 31; ZczEos= 0.30506
=
2kS
4.8 X
Hi = -3.25; H2 = -3.35
(qu*)ij
"""(
= 22.2478; Kc2,~* 11.2131
Kcl,M*
-e12k -- (f11e22)1'2 k
+
C(CU*)~,)Xj
cm3/molecule
Note: The expressions for fugacity coefficient (A-1) are the working equations. These equations have been arranged such that their mode of application is clear. Arranging them in this way sacrifices the clear relation of the equations to the underlying theory, but it helps to clarify the end result. The theoretical derivation of the fugacity coefficient takes into account the physical interactions between molecules of various oligomers via the equation of state and relates the chemical interactions (association) to an infiiite series of reactions. Hence, the binary mixture is actually considered to be a mixture of oligomers of various chain lengths. The theoretical derivation was summarized for dimerization (the first element of the series) by Elliott et al. (1990), and for the infinite series, the derivation of Heidemann and Prausnitz (1976) shows how the extension is made. One way of perceiving the end result is to consider the difference between the fugacity coefficient that comes directly from the equation of state for physical interactions relative to the "apparent" fugacity Coefficient which is used in the calculations. The fugacity coefficient from the physical interactions is based on the true mole fraction of species, mole fraction of the monomer, for instance. The apparent fugacity coefficient is based on the apparent mole fraction of a component given by the mass of the component divided by the molecular weight of the monomer. The relationship between the two is
molecules/cm3
0.303 molecule
This relationship shows very clearly that the apparent fugacity coefficient decreases steadily as monomer reacts away to form higher oligomers.
Literature Cited al = 28.73; a2 = 23.38
It should be noted that K 1 and K2 are included in al and a2 but that K12does not appear explicitly in the working
equations. From eqs A-2, A-3, and A-4
W1 0.027; W2 0.002; - - 0.171 NO NT
Z(CU*)lj,@j
(!"
+
2k3
1-
a)ul* =
3.78 X
cm3/molecule
Abbott, M. M.; Vanness, H. C. (1975)In VapoPLiquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA Frankfurt, Germany, 1982; Vol. 1/2c, p 66. Anderko, A. Extension of the AEOS model to Systems Containing Any Number of Associating and Inert Components. Fluid Phase Equilib. 1989,50, 21-52. Barr-David, F.; Dodge, B. F. Vapor Liquid Equilibrium at High Pressure. The Systems Ethanol-Water and 2-Propanol-Water. J. Chem. Eng. Data 1989,4, 107. Butcher, K. L.; Robinson, W. I. Apparatus for Determining High Pressure Liquid-Vapor Equilibrium Data. I. The Methanol Ethanol System. J. Appl. Chem. (London) 1966,16, 289. Butcher, K. L.; Medani, M. S. Thermodynamic Properties of Methanol-Benzene Mixtures A t Elevated Temperatures. J . Appl. Chem. 1968,18, 100-107. Elliott, J. R., Jr.; Suresh, S. J.; Donohue, M. D. A Simple Equation of State for Nonspherical and Associating Molecules. Znd. Eng. Chem. Res. 1990,29, 1476-1485. Fredenslund, Aa.; Jones, R. L.; Prauenitz, J. M. Group-Contribution Estimation of Activity Coefficients in Nonideal Liquid Mixtures. AZChE J. 1975,21 (61, 1086. Freshwater, D. C.; Pike, K. A. (1967) In Vapor-Liquid Equilibrium Dota Collection; Gmehling, J., Onken, U., Me.; DECHEMA Frankfurt, Germany, 1977; Vol. 1/2a, pp 76,126,626; Vol. 1/2b, P 45; Gmehling, J.; Liu, D. D.; Prausnitz, J. M. High Pressure VaporLiquid Equilibrium for Mixtures Containing One or More Polar
532
Znd. Eng. Chem. Res. 1991,30, 532-536
Components. Chem. Eng. Sci. 1979,34,951. Gomez-Nieto,M. A. High Pressure Vapor Liquid Equilibrium: The Propane-Ethanol-Acetone System. Ph.D. Dissertation, Northwestern University, Evanston, IL, 1977. Gordon, A. R.; Hines, W. G. (1946) In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA: Frankfurt, Germany, 1977; Vol. 1/2a, pp 323-325. Griswold, J.; Wong, S. Y. Phase Equilibria with Acetone-Methanol-Water System from 100°C into the Critical Region. Chem. Eng. Prog. Symp. Ser. 1952, 48 (3), 18. Gupte, P. A.; Daubert, T. E. Extension of UNIFAC to High Presswe Vapor-Liauid Equilibrium Using Vidal Mixing Rules. Fluid Phhse Eqiilib. 1986,28, 155-170 Hall, D. J.; Mash, C. J.; Pemberton, R. C. NPL Rep. Chem. 95, January 1979. Heidemarh, R. A.; Prausnitz, J. M. A Van der Waals-type Equation of State for Fluids with Associating Molecules. Proc. Natl. Acad. Sci. 1976, 73, 1773. Ikonomou, G. D.; Donohue, M. D. Extension of the Associating Perturbed Chain Theory to Mixtures with more than One Associating Component. Fluid Phase Equilib. 1988,39, 129. Jackson, G.; Chapman, W.G.; Gubbins, K. E. Phase Equilibria of Associating Fluids: Spherical Molecules with Multiple Bonding Sites. Mol. Phys. 1988,65, 1. Kay, W. B.; Donaham, W. E. Liquid-Vapor Equilibria in the is0 Butanol-n Butanol, Methanol-n Butanol and Diethyl Ether-n Butanol Systems. Chem. Eng. Sci. 1965,4, 1. Lee, R. J.; Chao, K. C. Local Composition Embedded Equation of State for Strongly Nonideal Fluid Mixturee. Znd.Eng. Chem. Res. 1989,28,1251-1261. Marinichev, A. N.; Suearev, M. P. (1965) In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA: Frankfurt, Germany, 1977; Vol. 1/2a, pp 79-81,629. Polak, J.; Murakami, S.; Benson, G. C.; Pflug, H. D. (1970) In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Ma.; DECHEMA: Frankfurt, Germany, 1982; Vol. 1/2c, p 488.
Prausnitz, J. M.; Lichtenthaler, R. N.; de Amvedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria. In Fugacities in Liquid Mixtures: Excess Functions; Neal, R. Amudson, Ed.; Prentice-Hall: Englewood Cliffs, NJ, 1986. Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1987. Sada, E.; Morisue, T. Isothermal Vapor-Liquid Equilibrium Data of IsoPropanol-Water System. J. Chem.Eng. Jpn. 1975, 8, 191. Schmidt, G. C. (1926) In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA Frankfurt, Germany, 1977; Vol. 1/2a, pp 116-122. Schwartzentruber, J.; Renon, H. Extension of UNIFAC to High Pressures and Temperatures by the Use of a Cubic Equation of State. Ind. Eng. Chem. Res. 1989,28,1049-1055. Soave, G. Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chem. Eng. Sci. 1972,27, 1197. Tamir, A.; Appelblat, A.; Wagner, M. Fluid Phase Equilib. 1981,6, 113. In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA: Frankfurt, Germany, 1982, Vol. 1/2c, p 77. Udovenko, V. V.; Frid, I. B. (1948) In Vapor-Liquid Equilibrium Data Collection; Gmehling, J., Onken, U., Eds.; DECHEMA: Frankfurt, Germany, 1977; Vol. 1/2a, pp 337-340. Vinichenko, I. G.; Susarev, M. P. (1966) In Vapor-Liquid Equilibrium Data Collection;Gmehling, J., Onken, U., Eds.;DECHEMA: Frankfurt, Germany, 1977, Vol. 1/2a, p 327. Whiting, W. B.; Prausnitz, J. M. Equations of State for Strongly Nonideal Fluid Mixtures: Application of Local Composition Toward Density-Dependent Mixing Rules. Fluid Phase Equilib. 1982, 9, 119. Wilson, G. M. Excess Free Energy of Mixing. J. Am. Chem. SOC. 1964, 86, 127.
Receiued for reuiew March 19, 1990 Reuised manuscript receiued August 2, 1990 Accepted August 31, 1990
Solubility of Carbon Dioxide in Tar Sand Bitumen: Experimental Determination and Modeling Milind D. Deo,* Chia J. Wang, and Francis V. Hanson Department of Fuels Engineering, Uniuersity of Utah, Salt Lake City, Utah 84112
An understanding of the solubility of carbon dioxide (CO,)in tar sand bitumen is essential for the development of in situ processes in the recovery of bitumen from tar sand deposita. The solubility of COz in the Tar Sand Triangle (Utah), the PR Spring Rainhow I (Utah),and the Athabasca (Canada) tar sand bitumens was determined with the use of a high-pressure microbalance a t temperatures of 358.2 and 393.2K and pressures up to 6.2 MPa As expected, the solubilities increased with pressure a t a given temperature and decreased with increases in temperature. The’PengRobinson and the Schmidt-Wenzel equations of state were used to match the experimentally observed solubilities. Correlations for the interaction parameters between C02and the bitumen were developed for both equations of state, wherein the interaction parameter could be obtained by using specific gravity and the UOP K factor for the bitumen. The correlations were developed with the optimum interaction parameters obtained for each of the samples at each temperature.
Introduction Depleting petroleum reserves and increasing energy demands will eventually require better exploitation of unconventional petroleum reserves like oil shale and tar sands. Tar sands are a mixture of sand and bitumen. The organic constituent, bitumen, cannot be recovered by ordinary petroleum production methods due to its high viscosity and the lack of reservoir energy in most deposita. The world’s known tar sand reservoirs contain more than 2100 billion barrels of bitumen in place, with about 30 billion barrels of bitumen in place in the United States. The state of Utah has about 96% of the total United States tar sand resource (excluding the state of Alaska). This
represents a potentially significant domestic energy resource, compared with the crude oil reserves of the continental United States. Information regarding the Utah tar sand and oil deposita is available in reporta by the Utah Geological and Mineralogical Survey (Ritzma, 1979; Wood and Ritzma, 1972). Much research effort has been spent on oil recovery mechanisms and other aspects of the recovery process when C02is used aa a recovery agent (Holm and Josendal, 1974; Mungan, 1981). As a result of these studies, a number of ways in which C 0 2 becomes effective in removing oil from porous rock were recognized. These findings were summarized by Holm and Josendal, 1974:
088S-SS85/91/2630-0532$02.5~~~ 0 1991 American Chemical Society
