Ind. Eng. Chem. Fundam. 1980, 19, 291-294
a = interfacial area
t = interphase contact time co = initial concentration in organic phase at t = 0 cd = initial concentration in aqueous phase at t = 0 D = aqueous to organic equilibrium distribution coefficient At = time change corresponding to concentration changes ct - co or C< - cd AH = heat of reaction, AH = H,- Hf [CTho] = molar concentration of thorium in organic phase [TBP] = free TBP concentration ([TBPiJ - 2[CTh,O] (excludes TBP complexed with HNOJ Ht = reverse (0-A) activation energy Hf= forward (A-0) activation energy
29 1
Literature Cited Bjerrum, J., Schwarzenbach, G., Slilen, L. G., "Stability Constants. Part 11. Inorganic Ligands", p 54, The Chemical Society, London, 1958. Flagg, J. E., "Chemical hocessing of Reactor Fuels", pp 227, 228, Academic Press, New York, 1981. Horner, D. E., Mallen, J. C., Thiel, S. W., Scott, T. C., Yates, R. G., Ind. f n g . Chern. Fundarn., 19, 103 (1980). Strenllng, C. Y., Scrlven, L. E., AIChEJ., 5(4), 514 (1959).
Receiued for review October 24, 1979 Accepted May 14, 1980
Research sponsored by the Division of Chemical Sciences, U.S. Department of Energy, under Contract W-7405-eng-26with the Union Carbide Corporation.
Three-Parameter Cubic Equation of State for Normal Substances A. Harmens' Petrocarbon Developments Ltd., Manchester M22 4TB, England
H. Knapp Technische Universitat Berlin, Instifut fur Thermodynamik und Anlagentechnik, 1000 Berlin 12, Germany
A cubic equation of state is proposed which possesses three adjustable parameters and therefore is more flexible than the well-known Redlich-Kwong and Peng-Robinson equations. The three parameters are adapted to give an optimal representation of the critical isotherm up to about pr = 5 and in addition are fitted to the vapor pressure from about T, = 0.3 to T, = 1 and to the volume along the critical isobar up to about T, = 2.5. The parameters were correlated in terms of T,, pc, and acentric factor w . The equation was then subjected to tests on pure substances, covering an extensive area of the p , Tfield. It was found in overall performance to be superior to four other generalized cubic equation of state procedures.
Introduction For process calculations particularly in petroleum and cryogenic technology, cubic equations of state have proved to be extremely useful. They are simple and can be solved with a straightforward algebraic procedure, so that they lead to robust computer programs for the prediction of thermodynamic data and t o relatively short computing times. They contain only a small number of adaptable parameters which can easily be related to the critical properties, so that the equations lend themselves well for generalization and application to mixtures. A t least one of the parameters must be treated as an empirical function of temperature, but such a function can also be generalized. Comparisons have shown (Oellrich et al., 1977) that thermodynamic computation procedures with generalized cubic equations of state are on the whole not noticeably inferior to the more expensive procedures using extended virial equations. However, the very fact that these equations are cubic in volume also imparts to them their main characteristic weakness: their failure to give a good account of the volume a t the critical point and at temperatures and pressures immediately above the critical. In that region the pressure changes not as a third-order function but as a fifth-order function of volume (Levelt Sengers, 1976) so that the cubic equation of state is fundamentally inadequate. Fortunately, for technical reasons of flow stability and process control it would be extremely unattractive to run a chemical process at or immediately above the critical 0196-4313l80/1019-0291$01.OO/O
point of one of its fluids: process designers avoid such conditions. The prediction of phase equilibrium can of course suffer from said inadequacy, but in many cases this leads to only minor errors since the calculated fugacity is much less seriously affected than the calculated volume itself (Harmens, 1975). So the shortcomings of the cubic equation are not always experienced as particularly troublesome. The various cubic equations of state also differ as to the extent and the severity of their supercritical unreliability. From the various cubic equations of state and associated procedures for parameter evaluation, two have acquired wide popularity: the Redlich-Kwong-Soave (Soave, 1972) and the Peng-Robinson (1976) procedure. They employ different two-parameter equations of state, but the treatment of the parameters is very similar. In both cases the repulsion parameter is kept constant at i h critical point value, while the attraction parameter contains a generalized function of temperature, fitted to the vapor pressure. For superheated vapors that subcritical temperature function is simply extrapolated above the critical temperature. The present study originated in the conviction that not all avenues had been explored which might lead to even better procedures. In particular, the use of a three-parameter equation and the actual fitting of the generalized temperature function in the region of superheated vapor seemed to deserve further investigation. Three different three-parameter equations were subjected to an explora-
0 1980 American Chemical Society
292
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
tory study of the practicability of various fitting procedures for the parameters. These three equations were the Clausius equation, as studied also by Freze and Chevalier (1978) and by Martin (19791, the equation studied lately by Usdin and McAuliffe (1976) which had previously been investigated by Himpan (1952) and by Fuller (1976), and finally the equation to be discussed below. The results of these preliminary investigations led us to rejecting the first two equations in favor of the third. The Three-Parameter Cubic Equation The preferred cubic equation is RT a (1) = - u2 + ucb - (c - 1) b2
u-b
Our study will largely be restricted to normal fluids, which will be characterized by their T,, pc,and w. The parameters a , b, and c will be expressed in terms of these three constants, and a will also be a function of reduced temperature. The shape of (1)was chosen for the simple reason that it can be regarded as the generating equation for several two-parameter equations with which we have favorable experience. For c = 1, (1)reduces to the Redlich-Kwong equation, for c = 2 to the Peng-Robinson equation, and for c = 3 to an equation used for phase equilibrium prediction in the N2-Ar-02 system (Harmens, 1977). Mollerup (1979) had already drawn attention to this relationship, and recently the equation was also studied by Schmidt and Wenzel (1980). The three adjustable parameters are evaluated by simultaneous solution of the equation of state itself and the two classical critical point conditions
at p = p,, T = T,, and u = u,. This leads to a cubic equation for p = b / u c in terms of the calculated critical compressibility factor ( = (pc/RTc)u,,cdc.The parameters a, b, and c are all functions of 0. For (1)the expression for the solution of this cubic equation is very involved and quite unsuitable for a routine calculation procedure, although the graph of p vs. ( shows only a very slight curvature. Therefore a simple polynomial was fitted to the calculated (3 values = 0.10770 + 0.76405(- 1.24282p + 0.96210$ (3) which over the range of 0.20 < ( < 0.34 represents the actual 0 values with an average absolute error of only 0.003 % . This leaves ( to be defined. For an equation with only two adjustable parameters, the simultaneous solution of the system referred to above yields a fixed value of { which is characteristic for the equation. So we have: RedlichKwong, { = 0.333333; Peng-Robinson, { = 0.307401; Harmens, { = 0.286186. In (l),the third parameter of course allows { to be chosen freely. Abbott (1973) has already warned against setting ( equal to the physical critical compressibility factor 2, = ( ~ , / R T , ) U ~ Although , , ~ ~ ~ . that makes the equation predict the correct critical volume, it compromises the calculated volumes at low and high pressures much more than can be tolerated. Our introductory study confirmed this. Better overall results are obtained with { > 2,. To define (, we have studied the critical isotherms of seven substances, from argon to n-decane, going from atmospheric pressure up to about p, = 5. For each substance a value of (was derived which made the average relative
error in predicted volume along the critical isotherm zero. So for each substance the {was chosen so as to give on average the best calculated critical isotherm. For these ( values it was then found that A{ = { - 2, was a quadratic function of w A( = 0.0301 + 0.0384~' Since for 2, we have after Pitzer 2, = 0.291 - 0 . 0 8 0 ~ this led to the relation ( = 0.3211 - 0 . 0 8 0 ~+ 0.0384~'
(4) With ( and 0thus defined, the simultaneous solution of (1)and (2) yields the following parameter values at the critical point 1 - 3{ R2T,2 RTC a = Ra -; b = ab -; c = 1 - (5) Pc Pc 0f in which R, = 1 - 3{+ 3 p + p((3 - 6{+ 0 {)
+
a b = Pf (6) The Temperature Dependence of P a r a m e t e r a A t temperatures other than the critical the parameter a will now be modified by a coefficient a which depends on T,, so that R2T,2 a = a(Tr)Qap~ (7)
The parameters b and c retain their critical point value throughout. Evidently a t the critical temperature a(T,) = a(1) = 1. Below the critical temperature a has been fitted to vapor pressures; above that temperature to volumes along the critical isobar. Vapor pressure data for 20 substances ranging from argon to n-decane (mostly hydrocarbons) were taken from reliable literature sources. For the normal alkanes the low-temperature data from Carruth (1970) were included, and for several substances Carruth's correlation was used for calculating liquid vapor pressures below the freezing point. So for almost all substances the data went from T, = 1 down to T , = 0.3 to 0.4. Values of a(T,) were computed from (7) by solving (1)plus its associated fugacity formula for that value of a which yielded the given vapor pressure at the corresponding T,. From the RedlichKwong-Soave and Peng-Robinson procedures we know that plots of al/' vs. T;I2 should be rectilinear. Our results confirmed this, but toward the lowest temperatures the plots tended slightly to turn away from the straight line. Therefore a two-constant correlation appeared to be indicated, and after several trials the formula 1
+ A(1- fi) -B
was chosen ( T , 51.0). The parameters A and B were correlated with w. For w > 0.2 the correlations were linear, but below that value parabolas were required. At w = 0.2 there is a smooth transition between the two sections. w I0.2 A = 0.50 + 0.27767~+ 2 . 1 7 2 2 5 ~ ~ B = -0.022 + 0 . 3 3 8 ~- 0.8450' w
> 0.2
A = 0.41311 + 1.14657~ B = 0.0118
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
For superheated vapor, volumes along the critical isobar were calculated for the same 20 substances with the LeeKesler procedure (1975). We used 24 volume values between TI = 1.02 and TI = 2.6. There is no point in extending the fitting beyond that range since the attractive term becomes very small compared with the repulsive term. Plots of calculated a values vs. TI showed that for TI > 1.03 the a was noticeably lower than the extrapolation of the subcritical a according to (8). Curves of a vs. In TI resembled parabolas, which led to their correlation by
+ 1 . 5 2 2 7 ~ In) TI + (0.1533 + 0.41w)(ln T J 2 (TI > 1.0) (10)
293
Table I. Average Absolute Errors in u along Critical Isotherm
Ar CH4 Na C3H8 nGH14 n-C.H.. n-c,,di, weighted average
RKS
PR
SW
HK
4.41% 6.63 4.39 7.85 10.83 13.66 16.22 9.13
3.95% 4.18 4.46 3.64 2.95 6.04 8.23 4.81
4.49% 6.39 3.71 4.80 3.64 4.90 5.11 4.78
2.84% 4.11 3.16 3.66 1.63 4.49 4.82 3.61
a(?',) = 1.0 - (0.6258
As a result of being defined in two parts, by (8) and (101, the function a ( T J has a discontinuous slope a t TI = 1. Going toward lower temperatures, the slope suddenly drops by about 28% at that point. This does not affect the fugacity, but the discontinuity will show up in the calculated enthalpies and entropies at the critical temperature. A computer program can be made to rectify this by interpolating locally between values calculated a short distance on either side of T,. Summarizing, the procedure for calculating the parameters a, b, and c consists in finding { and P with (4)and (3); finding fl, and flb with (6); finding b and c with (5); finding a(TI)with (8) and (9), or with (lo), and finally finding a with (7). Calculations of Volume of Superheated Vapor When the equation, with the use of the above procedure for TI > 1.0, is solved for u in the superheated vapor region, a singularity may occur which is unknown with the twoparameter equations. For T, > 1.0 one expects the cubic equation to yield only one real solution for u, as indeed the two-parameter equations do. Equation 1, however, in which both b and c are maintained at their critical point values throughout, appears unexpectedly to revert to the threefold solution a t a certain distance above the critical temperature. Above that point three volumes are calculated, of which only the largest should be used; the two other values are usually negative and thus meaningless. The largest volume is continuous with the volumes calculated a t lower temperatures as single solutions. The transition to the multiple root solution occurs at about TI = 2.6 for w = 0 and about TI = 1.4 for w = 0.5. If the equation is solved with an iterative process, the transition requires no extra provisions, but if the standard algebraic solution is used, the multiple root procedure must also be provided for in the superheated region, and the largest root must be selected. Tests and Comparisons with Other Procedures The new procedure (HK) was tested by comparing calculated results with literature data for pure substances. These comparisons were also run simultaneously for four other cubic equation of state procedures, namely Redlich-Kwong-Soave (RKS), Redlich-Kwong-GraboskiDaubert (1978) (RKGD), Peng-Robinson (PR), and Schmidt-Wenzel (SW). The Lee-Edmister equation (1973) does not enter these comparisons because it holds only for the vapor phase. Likewise, the procedure by Usdin and McAuliffe was left out because it requires four instead of three characterizing parameters. All five procedures mentioned above were subjected to five different tests. (1) The Critical Isotherm. Volumes were calculated along the critical isotherm from about atmospheric pressure up to pr = 4 to 5. Between 17 and 24 points were calculated on each isotherm, 145 points in total. Table I
Table 11. Average Absolute Errors in Boiling Point Temperatures PR
SW
HK
0.26% 0.24 0.21 0.05 0.35 0.35 0.50 0.30 0.49 0.11 0.29 0.31 0.42 0.33 0.311
0.06% 0.09 0.07 0.09 0.14 0.18 0.22 0.11 0.20 0.16 0.13 0.10 0.17 0.25 0.141
0.25% 0.26 0.17 0.13 0.14 0.13 0.10 0.10 0.24 0.22 0.10 0.18 0.24 0.21 0.175
RKS RKGD
Ar
0.26% 0.41 0.27 0 2 0.14 NZ CZH4 0.10 0.11 Has 0.10 C3H8 n-C4H10 0.13 C6H6 0.27 COZ 0.22 i-C5Hlz 0.16 n-C6H14 0.18 n-C,H,, 0.17 n-C;,,H';, 0.24 weighted average 0.197
CH4
0.34% 0.52 0.34 0.18 0.13 0.14 0.10 0.14 0.26 0.23 0.18 0.18 0.17 0.24 0.225
Table 111. Average Absolute Errors in Saturated Liquid Volumes RKS 4.62% 2.86 3.44 0 2 4.72 NZ Ca H4 6.35 6.35 HZS 10.20 C3H8 n-C4H10 10.01 C6H6 12.42 COZ 12.66 i-C5Hla 11.10 n-C6 13.74 n-C,H,, 16.67 n-C,oHz, 22.32 weightedaverage 9.82
Ar CH4
RKGD
PR
SW
HK
4.62% 2.85 3.44 4.71 6.34 6.35 10.20 10.01 12.42 12.66 11.10 13.74 16.67 22.32 9.82
9.69% 9.73 9.25 8.55 7.21 6.28 5.46 4.09 3.69 4.66 4.96 2.10 4.78 8.32 6.32
4.38% 2.47 3.27 4.22 3.21 2.64 4.42 2.51 4.18 4.04 3.34 2.03 3.43 3.08 3.35
5.68% 4.80 5.13 5.61 4.74 3.79 4.70 3.98 3.76 5.20 5.30 3.16 3.59 2.94 4.41
lists the average absolute errors. RKGD is not shown separately since at T = T , it is identical with RKS. The weighted average errors for the various procedures show that HK gives the best results, as was to be expected. (2) Boiling Point Temperatures. Table I1 lists the average absolute errors in calculated boiling point temperatures a t given vapor pressures, from TI = 0.3 to 0.4 upward. Between 8 and 17 points were calculated along each saturation curve, 206 points in total. SW shows the best performance, which stands to reason since it was specially made for the saturated state. HK is the second best, The results for PR are adversely influenced by relatively large errors at subatmospheric pressures. I t is surprising that RKGD is noticeably worse than RKS. (3) Saturated Liquid Volume. Table I11 lists the average absolute errors in saturated liquid volume calculated simultaneously with the above boiling point temperatures for the same 206 points. Again SW and HK come out noticeably better than the other procedures. (4) Volumes along the Critical Isobar. Table IV gives the average absolute errors in volumes calculated for
294
Ind. Eng.
Chem. Fundam., Vol. 19, No. 3,
1980
Table IV. Average Absolute Errors along Isobars pr = 1
RKS
Ar
0.24% 0.39 0 2 0.42 NZ 0.24 C,H, 0.70 1.66 H,S 0.82 C3H.S ‘bH6 1.61 CO, 1.45 n-C,H,, 1.17 weighted average 0.845
CH,
RKGD 0.26% 0.42 0.46 0.25 0.73 1.68 0.83 1.62 1.46 1.18 0.864
PR 2.03% 1.99 2.28 2.09 1.54 1.82 1.33 1.24 0.69 1.01 1.590
SW 1.21% 1.32 1.68 1.69 1.86 2.04 2.37 2.46 2.29 2.95 1.982
HK 0.12% 0.18 0.21 0.30 0.29 1.45 0.32 0.84 0.61 0.59 0.458
Acknowledgment A. H. wishes to thank the Management of Petrocarbon Developments Ltd. for permission to publish the results of this investigation. Appendix The equation of state (1)is associated with the following formulas for fugacity coefficient, enthalpy, and entropy departure. In these formulas we use the symbols
K=d1 L=-ln Kb
Table V. Average Absolute Errors along Isobars pr = 2
RKS Ar CH,
2.69% 3.22 3.07 2.88 4.01 4.04 4.37 5.97 5.19 5.85 weighted average 4.159
RKGD 2.80% 3.32 3.17 2.96 4.07 4.11 4.41 5.99 5.21 5.86 4.220
PR 2.05% 1.91 2.15 2.09 1.22 2.30 1.18 1.60 1.10 1.77 1.725
SW 2.51% 2.81 3.18 3.41 3.21 4.24 4.20 4.42 3.91 4.90 3.654
HK 2.90% 3.29 3.06 2.48 3.28 3.31 3.04 4.27 3.26 3.67 3.279
superheated vapor along isobars close to pr = 1, from T , just above 1.0 to T I about 2.5. Between 8 and 17 points were calculated on each isobar, 137 points in total. HK understandably gives the smallest deviations. (5) Volumes along the Isobar pr = 2. Table V is the sequel to Table IV for volumes calculated along isobars close to pr= 2, for the same 137 temperatures. All methods give rather large errors at temperatures just above the critical, as was to be expected. For HK errors of between 8 and 15% were found at around T I = 1.1. At lower T I the error diminishes. As a rule the error is well below 1% for T I > 2. PR gives the best results for this test, followed by HK. Tests 2 to 5 were also carried out for ammonia, a substance of considerable polarity. In the case of HK, the errors for test 2 were perfectly in keeping with those for the other substances, but for tests 3 to 5 the errors were about three times the weighted average. The other procedures gave similar results. Summarizing, eq 1 and four other generalized cubic equations of state were tested by comparing average absolute deviations between calculated values and literature data for boiling point temperatures, for volumes along the critical isotherm, the critical isobar and the isobar p r = 2, and for liquid volumes along the saturation line. These comparisons, shown in Tables I-V, indicate that the proposed equation represents the 831 test points with the lowest overall weighted average error (HK, 2.38%; SW, 2.63%; PR, 3.03%; RKS, 4.90% and RKGD, 4.92%). It is intended to apply the proposed equation to mixtures, and to compare calculated VLE results with empirical data.
In
cp
PV RT
= -- 1
Ah=
(11) 2u+b(c+K) 2u + b ( c - K )
(12)
RT a + In -p ( u - b) RT
(a - T -3L - p u + R T
RT da As = R In -- U - b dT Ah is the isothermal pressure correction for enthalpy, to be subtracted from ideal gas-state enthalpy. As relates in the same way to ideal gas-state entropy at unit pressure. Nomenclature a, b, c = adjustable parameters in (1) K = symbol defined by (11) L = symbol defined by (12) p = pressure R = universal gas constant T = absolute temperature u = molal volume 2 = physical value of compressibility factor cy = coefficient in (7) /3 = b / u , as defined by (3) .( = calculated value of compressibility factor cp = fugacity coefficient Qa, Qb = coefficients as defined by (6) w = acentric factor index c = “critical” index r = “reduced”: divided by critical Literature Cited Abbott, M. M., AIChE J., IS, 596 (1973). Carruth, G. F.. Thesis, Rice University, Houston, Texas, Nov. 1970, Part 2. Freze, R., Chevalier, J. L., Can. J . Chem. Eng., 58, 375 (1978). Fuller, G. G.. I d . Eng. Chem. Fundam., 15, 254 (1976). Graboski, M. S.,Daubeft, T. E., Ind. Eng. Chem. Process Des. Dev., 17, 443 (1976). Harmens, A., Cryogenics, 15, 217 (1975). Harmens, A., Cryogenics, 17, 519 (1977). Himpan, J., Z. phys. 131, 17, 130 (1951); 133, 647 (1952). Lee, E. I., Erbar, J. H., Edmister, W. C., AICM J., IS, 349 (1973). Lee, E. I., Kesler, M. G., AIChE J., 21, 510 (1975). Levell Sengers, J. M. H., physica, 82A, 319 (1976). Martin, J. J., I d . Eng. Chem. Fundam., 18, 81 (1979). Molierup, J., private communication, May 1979. Oeilrich, L., Pliicker, U., RausnRz, J. M., Knapp, H., Chem. Ing. Tech., 49, 955 (1977). Peng, D. Y.. Robinson, D. E., Ind. Eng. Chem. Fundem., 15, 59 (1976). Schmidt, G., Wenzei, H., Chem. Eng. Sci., in press, 1980. Soave, G., Chem. Eng. Sci., 27, 1197 (1972). Usdin, E., McAuliffe, J. C., Chem. Eng. Sci., 31, 1077 (1976).
Received for review November 12, 1979 Accepted April 24, 1980
