- Email: [email protected]
Prediction of the Joule–Thomson inversion curve of air from cubic equations of state
Cryogenics 38 (1998) 721–728 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0011-2275/98/$—see front matter
PII: S0011-2275(98)00036-8
Prediction of the Joule–Thomson inversion curve of air from cubic equations of state Coray M. Colina and Claudio Olivera-Fuentes* TADiP Group, Thermodynamics and Transport Phenomena Department, Simo´n Bolı´var University, AP 89000, Caracas 1086-A, Venezuela
Accepted 4 February 1998 A modified van der Waals equation of state recommended in the literature for improved prediction of the inversion curve of air is shown to be thermodynamically inconsistent, giving large errors in the critical and two-phase regions. An alternative procedure is presented by means of which the cohesion function of any cubic equation of state can be adjusted to give arbitrarily accurate representation of an experimental inversion curve. New versions of the van der Waals, Redlich–Kwong and Peng– Robinson equations of state are developed based on experimental inversion data of air, and are shown to give better inversion predictions than more complex, multiparameter noncubic equations of state. 1998 Elsevier Science Ltd. All rights reserved Keywords: E. Joule–Thomson effect; equation of state; liquid air
Nomenclature
vR
a
VLE Z
attraction parameter in equation of state (kPa.m6.kg−2 ) AAD average absolute deviation b excluded volume parameter in equation of state (m3.kg−1 ) specific heat capacity at constant pressure cP (kJ.kg−1.K−1 ) EOS equation of state h specific enthalpy (kJ.kg−1 ) k1, k2 dimensionless constants in equation of state m1, …, m4 adjustable coefficients in new cohesion function MBWR modified Benedict–Webb–Rubin equation of state MVDW modified van der Waals equation of state P pressure (kPa) gas constant (kJ.kg−1.K−1 ) R s integration stepsize T temperature (K) specific volume (m3.kg−1 ) v
The Joule–Thomson effect is of considerable importance in refrigeration and gas liquefaction processes. Under appropriate conditions, adiabatic throttling of a gas will cause it to lower its temperature, an endothermic effect that *To whom correspondence should be addressed. Tel.: + 58-2906-37-56; Fax: + 58-2-906-37-43; e-mail: [email protected]
dimensionless volume, defined by Equation (9) vapour–liquid equilibrium compressibility factor, Pv/(RT)
Subscripts c r
critical value reduced, relative to critical value
Superscripts exp
experimental
Greek letters ␣
cohesion parameter in equation of state, a/ac Joule–Thomson coefficient (K.kPa−1 ) acentric factor dimensionless critical attraction parameter, acPc/(RTc )2 dimensionless excluded volume parameter, bPc/(RTc )
⍀ac ⍀b
may be strong enough even to induce condensation. For this cooling to take place, the thermodynamic state of the fluid must lie in the region bounded by the inversion curve or locus of points at which the Joule–Thomson coefficient is zero,
⬅
冉冊 ∂T ∂P
= 0.
(1)
h
Cryogenics 1998 Volume 38, Number 7 721
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes As usually represented in pressure–temperature co-ordinates, the Joule–Thomson inversion curve extends from a minimum temperature corresponding to a saturated state on the vapour pressure line, to a maximum temperature corresponding to the ideal-gas limit at zero density and pressure. The curve is parabolic in shape and exhibits a maximum or peak inversion pressure at an intermediate temperature. Knowledge of the inversion curve is obviously essential for good design and operation of Joule–Thomson cooling processes. Direct measurement of inversion points is difficult and unreliable. At near-inversion conditions, the vanishing Joule–Thomson coefficients imply that even large pressure changes will produce very small temperature differences, detectable only by extremely accurate measurement. It is generally preferable to make use of thermodynamic relations such as
=
冋冉 冊 册
∂v 1 T cP ∂T
−v =
P
冉冊
RT 2 ∂Z cPP ∂T
(2)
P
to derive Joule–Thomson coefficients and inversion points from experimental P–v–T data either by direct numerical treatment of the raw data, or more usually by first fitting the data by an equation of state (EOS). Even then, the inversion curves of the heavier fluids cannot be completely determined, because they extend into regions of high temperatures and pressures not attainable in experimental study. Efforts have been made to obtain generalised inversion curves from the known behaviour of light fluids. Gunn et al.1 calculated inversion points from volumetric data of argon, carbon monoxide, ethylene, methane, nitrogen and xenon. Since the data did not cover the upper temperature branch of the inversion curve, they computed additional theoretical points for argon from a truncated virial EOS with second and third virial coefficients based on the Kihara intermolecular potential. The entire set of points was then correlated by the empirical equation Pr = − 36.275 + 71.598Tr − 41.567T 2r + 11.826T 3r
(3)
− 1.6721T + 0.091167T . 4 r
5 r
An alternative, somewhat simpler equation based on published inversion data for argon, ammonia, carbon dioxide, carbon monoxide, ethylene, methane, nitrogen and propane was proposed by Miller2, Pr = 24.21 −
18.54 − 0.825T 2r. Tr
(4)
Equations (3) and (4) are two-parameter corresponding states correlations that give comparable results, since they derive from data for the same or very similar fluids, described1 as gases that have small, nearly spherical molecules with acentric factor close to zero, and2 as mostly gases that have a critical compressibility factor Zc of 0.29. The computation of Joule–Thomson inversion curves is one of the most demanding tests that can be performed on a set of experimental data1 or on an analytical EOS2. For pressure-explicit EOS’s, an equivalent formulation of the thermodynamic inversion criterion of Equations (1) and (2),
722
Cryogenics 1998 Volume 38, Number 7
冉冊 冉冊
T
∂P ∂T
+v
v
∂P ∂v
=0
(5)
T
shows that the EOS must accurately predict not only the P(T,v) function but also its first partial derivatives with respect to temperature and volume (or density), over wide temperature and pressure ranges that can exceed reduced values Tr = 5 and Pr = 12. Even state-of-the-art multiparameter EOS’s may fail to generate physically correct inversion curves when extrapolated to such extreme conditions3,4. Several authors have therefore used the prediction of inversion curves as a criterion to assess and rank the predictive capabilities of available EOS’s. Miller2 studied the van der Waals, Dieterici–Porter, Redlich–Kwong and Clausius–Martin EOS’s. Using Equation (4) as a basis of comparison, he concluded that the original Redlich– Kwong EOS was the best of the simple equations even though it underpredicted the peak inversion pressure by about 10%. He also suggested that analysis of the predicted inversion curve should become a standard test of any newly constructed EOS. Juris and Wenzel5 apparently unaware of Miller’s work, compared Equation (3) with predictions from the van der Waals, Dieterici, Berthelot, Redlich–Kwong, truncated virial (with second and third virial coefficients taken respectively from the Pitzer and Chueh–Prausnitz correlations), Su–Beattie–Bridgeman, several reduced forms of the Benedict–Webb–Rubin, and Martin–Hou EOS’s. They concluded that the van der Waals, Dieterici and Berthelot EOS’s had as might be expected very limited capability of predicting inversion behaviour. Among the more complex models the virial, Beattie–Bridgeman and Benedict– Webb–Rubin EOS’s gave poor predictions; only the Martin–Hou EOS yielded reasonable inversion values, although it could incorrectly predict two peaks of inversion pressure. Remarkably good results were obtained with the Redlich–Kwong EOS, which despite its relative simplicity predicted the inversion locus more accurately than any of the other cubic and noncubic models, even in the liquid region at more than twice the critical density. Dilay and Heidemann6 studied the Soave–Redlich– Kwong, Peng–Robinson, perturbed hard-chain and Lee– Kesler EOS’s, carrying out detailed comparisons for fluids with acentric factors ranging from 0 to 0.25. For simple fluids ( = 0) and nitrogen ( = 0.04) they observed that all four EOS’s gave good predictions of the low-temperature branch of the inversion curve up to a reduced pressure of about 8.0. The Soave–Redlich–Kwong EOS gave the best match of the peak inversion pressure, but failed to reproduce the high-temperature inversion region where only the Lee–Kesler EOS gave reasonably good results. Conclusions for carbon dioxide ( = 0.225) were similar, although less clear because of the greater scatter of the experimental high-temperature data. These authors remarked on the quite different results obtained from the Soave–Redlich–Kwong and Peng–Robinson EOS’s, which they found somewhat surprising given the similarities between these equations. Colazo et al.7 compared Equations (3) and (4) with predictions from the van der Waals, Soave–van der Waals, Redlich–Kwong, Peng–Robinson, Clausius–Martin, Schmidt–Wenzel, Valderrama–Cisternas, Trebble–Bishnoi and Trebble–Bishnoi–Salim EOS’s. They showed that for a generic cubic EOS of the form
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes P=
ac␣(Tr ) RT − 2 v − b v + k1bv + k2b2
(6)
the cohesion function ␣(Tr ) in the numerator of the second term is in principle responsible for the predicted lower and upper inversion temperatures, and the quadratic volume function in the denominator of the same term similarly dominates the predicted peak inversion pressure, which increases as the predicted critical compressibility factor Zc decreases. Cubic EOS’s having substantially identical cohesion functions but different predicted Zc will generate inversion loci that are also different. Since this is the case for the Soave–Redlich–Kwong (Zc = 1/3) and Peng– Robinson (Zc = 0.307) EOS’s, the simple explanation for the “somewhat surprising” differences in performance remarked upon earlier6 is that while these equations may indeed be similar for vapour–liquid equilibrium (VLE) computations, the similarity does not extend to inversion calculations because Joule–Thomson coefficients (unlike vapour pressures) are not translation-invariant properties8,9. The above results contrast sharply with conclusions reached by Geana and Feroiu10, who compared inversion curves predicted from a new four-parameter cubic EOS developed by themselves against values from the Redlich– Kwong, Soave–Redlich–Kwong and Peng–Robinson EOS’s. These authors stated that the translation in volume of a cubic EOS does not change the Pr − Tr inversion curve of the original equation, leading only to improvement in volume values that are displaced by a constant amount. This statement has subsequently been refuted11. It was also shown by Colazo et al.7 that cubic EOS’s which employ cohesion functions derived from pure-fluid (subcritical) saturation data cannot give good prediction when extended to the high (supercritical) temperatures typical of the upper branch of the inversion curve. This is partly due to extrapolation uncertainties, but the main factor is that the inversion behaviour of fluids at these conditions is dominated by the lower-order virial coefficients. This is why the original Redlich–Kwong and the Clausius–Martin EOS’s, both of which were developed to give improved representation of second virial coefficients, predict the upper inversion curve considerably better than the VLEorientated Soave–Redlich–Kwong and Peng–Robinson EOS’s. These authors suggested that new cohesion functions be developed from experimental data by solving Equations (5) and (6) simultaneously, for which they proposed a numerical procedure that will be described in greater detail in the next section. A different approach was followed by Najjar et al.12, who developed a two-constant cubic EOS identical in shape to the original van der Waals EOS but with new constants obtained by least-squares fit of experimental inversion data for air. Their modified van der Waals (MVDW) EOS can be written in SI units (kPa, m3.kg−1, K) as P=
ac RT − 2 , ac = 0.075228, b = 0.000762 v−b v
where vR is the dimensionless volume vR ⬅
Pcv . RTc
(9)
These authors found that the MVDW EOS predicted the inversion curve of air much better than other simple cubic EOS’s such as the original van der Waals, Redlich–Kwong, Soave–Redlich–Kwong and Peng–Robinson EOS’s, and reported an average standard deviation of 0.18 for the inversion curve. Inspection of their tabulated results shows an average absolute error of 5.87% and a maximum error of 21.78% in reduced inversion pressure. Unfortunately, by using modified constant parameters rather than introducing a temperature dependence for them as is usual nowadays, Najjar et al.12 made their MVDW thermodynamically inconsistent. In order to see this it suffices to compute the predicted critical point from this EOS, by setting the first and second partial derivatives of pressure with respect to volume equal to zero, or equivalently by requiring the EOS to give three equal volume roots at the critical point13. Doing this with Equation (7) yields Zc = 0.375 as expected for the van der Waals model, but with critical constants (in SI units) Pc =
ac 8ac = 101.9, = 4798.5, Tc = 2 27b 27bR
vc = 3b = 0.002286 which deviate substantially from the accepted values, e.g. Pc = 3770, Tc = 132.5, vc = 0.002913 used by Reynolds14. Repeating the calculations with Equation (8) gives the reduced co-ordinates of the MVDW “critical” point as Pr =
⍀ac 8⍀ac = 1.273, Tr = = 0.7693, 27⍀2b 27⍀b
vR = 3⍀b = 0.2266 instead of the correct values Pr = Tr = 1, vR = 0.375. The situation is illustrated in Figure 1, which shows three iso-
(7)
(the value actually reported by them, b = 0.00762 is obviously a misprint), and in reduced form as Pr =
Tr ⍀ac − , ⍀ac = 0.1961205, ⍀b = 0.07553501 vR − ⍀b v2R (8)
Figure 1 Isotherms of the MVDW EOS12 showing incorrect position of the critical point. (a) Tr = 1.0; (b) Tr = 0.7693; (c) Tr = 0.5; C.P. = predicted critical point (coordinates Pr = 1.273, vr = 0.2266)
Cryogenics 1998 Volume 38, Number 7 723
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes therms computed with Equation (8). The Tr = 1.0 isotherm does not exhibit the critical (inflexion) point required by classical thermodynamics; this point appears instead on the Tr = 0.7693 isotherm as explained above. Only below this temperature do the subcritical isotherms exhibit the van der Waals loop required for prediction of coexisting liquid and vapour phases; a third isotherm Tr = 0.5 is drawn as an example. It follows that although the MVDW EOS was intended for improved representation of air properties for cryogenic applications, it will fail precisely at low temperatures, predicting a phase coexistence region that extends over a too narrow temperature range (Tr ⱕ 0.7693) and a too wide pressure range (Pr ⱕ 1.273). It is of course basically incorrect to use a pure-fluid EOS such as the MVDW EOS to predict VLE of air which is a mixture, even if its bubble and dew points differ by only a few kelvins, but if equilibrium pressures and temperatures are taken from a suitable reference14 and Equation (7) is used just to compute the corresponding phase volumes, the results appear as in Figure 2. A liquid branch of the phase envelope is predicted for pressures ⬍ 180 kPa, i.e. at bubble temperatures under 85K only. Although this temperature lies well below the critical temperature predicted by Equation (7), the fact that the MVDW EOS overpredicts the critical pressure (and therefore all saturation pressures) results in only one real, vapour-like volume root at higher temperatures, so that the “liquid” and vapour branches become nearly identical. Since the inversion curve terminates on the liquid side of the saturation envelope, the MVDW EOS can be expected to give particularly poor prediction of air properties at the lower inversion temperatures. It is the purpose of the present work to show how cubic EOS’s can be developed that give an arbitrarily accurate prediction of the inversion curve, better in fact than that obtained from considerably more complex noncubic models, without sacrificing thermodynamic consistency. Examples of three such EOS’s are given, corresponding to modifications of the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s based on inversion data of air.
Development of equations of state from inversion curves The generic cubic EOS, Equation (6) is written in reduced terms as Pr =
Tr ⍀ac␣(Tr ) − 2 vR − ⍀b vR + k1⍀bvR + k2⍀2b
(10)
where the cohesion function ␣(Tr ) expresses the temperature dependence of the attraction term. On applying the inversion condition of Equation (5), the following differential equation results, also in reduced variables7, d␣ ⍀acTr 2 v + k1⍀bvR + k2⍀b dTr 2 R
−
(11)
⍀bTr ⍀acvR(2vR + k1⍀b ) ␣+ = 0. (v2R + k1⍀bvR + k2⍀2b )2 (vR − ⍀b )2
For an already existing cubic EOS, i.e. for known ␣(Tr ), the predicted inversion curve will be given by simultaneous solution of Equations (10) and (11). For the original van der Waals EOS (see Table 1), for example, on setting ␣(Tr ) ⬅ 1 we obtain Pr =
27 Tr − vR − 1/8 64v2R
(12)
and −
Tr 27 + = 0. 2 4vR (vR − 1/8)2
(13)
Equation (13) gives the inversion curve directly in Tr − vR co-ordinates. It may be combined with Equation (12) to give an equivalent Pr − vR version and also the well-known Pr − Tr formula15
冋
Tr = 3 1 ±
册
2 1 √9 − Pr . 6
(14)
Explicit or algorithmic solutions for other EOS’s have been presented by a number of authors, e.g. Gunn et al.1 for the truncated virial EOS, Miller2 for the Dieterici, Redlich– Kwong and Martin EOS’s, Dilay and Heidemann6 for the Soave–Redlich–Kwong and Peng–Robinson EOS’s, Colazo et al.7 and Olivera-Fuentes16 for the general cubic EOS of Equation (10), and Najjar et al.12 again for the van der Waals EOS. Conversely, given an inversion curve Pr(Tr ) the cohesion function may be computed by treating Equation (11) as a first order ordinary differential equation for ␣(Tr ), with vR an implicit function of Tr given by Equation (10). Since by definition ␣(Tr ) = 1 at Tr = 1, this is an initial-value problem that may solved by standard numerical methods. Colazo et al.7 used fourth-order Runge–Kutta integration with extrapolation to zero stepsize; their algorithm is also used in the present work and may be summarised as follows: Figure 2 Saturated liquid (䊊) and vapour (䊏) volumes predicted by the MVDW EOS12, compared with reference values (———) from Reynolds14
724
Cryogenics 1998 Volume 38, Number 7
1. Start at Tr = 1 with ␣ = 1. Set integration interval ⌬Tr, e.g. ⌬Tr = 0.1; 2. Obtain ␣(1) at Tr + ⌬Tr by Runge–Kutta integration
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes Table 1 Equations of state developed in this work EOS
EOS parameters
van der Waals Redlich– Kwong Peng– Robinson
Coefficients of new cohesion function
k1
k2
⍀ac
⍀b
0
0
27/64
1/8
− 14.1760
12.02581
1
0
0.42748…
0.08664…
− 5.72678
0.45723…
0.07779…
− 5.13662
−1
2
m1
with stepsize s(1), e.g. s(1) = 0.001. At each intermediate temperature, (a) Compute Pr from the known inversion curve Pr(Tr ), (b) Compute vR(Pr, Tr ) from the EOS Equation (10), (c) Compute d␣/dTr from the inversion condition Equation (11);
m4
R2
− 4.59465
0.514251
0.999985
4.404533
− 1.68415
0.235894
0.999941
4.501821
− 1.96490
0.345253
0.999421
m2
m3
comparison with the MVDW EOS proposed by those authors. Because the integration procedure requires a continuous representation of the entire inversion curve, the data were fitted by a fifth-degree polynomial similar to Equation (3), Pr = − 13.285 + 26.738Tr − 9.3407T 2r + 1.1469T 3r − 0.016527T − 0.0065313T 4 r
3.
Repeat the above computations to obtain a second estimate ␣(2) using stepsize s(2), e.g. s(2) = s(1)/2; Obtain ␣(Tr + ⌬Tr ) by linear extrapolation to zero stepsize,
4.
␣=
␣(1)s(2) − ␣(2)s(1) s(2) − s(1)
5.
Reset Tr to Tr + ⌬Tr and repeat from step (2).
(15)
The ad hoc extrapolation procedure of Equation (15), based on the observation that integration error approaches zero as a linear function of stepsize7, is essential to eliminate the accumulation of errors which would very rapidly invalidate the solution. The final accuracy of the computed cohesion function will in any case be assessed by predicting the entire inversion curve and comparing it with the original data. Using the above method, we have obtained new supercritical cohesion functions for air for use with the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s. Experimental inversion data in reduced variables Pr − Tr (see Table 2) were taken from Najjar et al.12 to allow direct
(16)
5 r
with regression coefficient R2 = 0.99991. Equation (16) represents the experimental data in Table 2 with an average absolute deviation (AAD) of 0.30%. A function similar to Equation (4) was also tried, but gave considerably worse results. Figure 3 shows the experimental data, the fitting equation, and the corresponding curve predicted from the generalised Equation (3), which is seen to deviate noticeably from the data, especially in the peak pressure and upper temperature regions. Since air (Zc = 0.289, = − 0.0037 based on bubble points14 ) may be expected to behave as a simple fluid for which Equation (3) should be applicable, these deviations might reflect some inaccuracy in the experimental data. We have been unable to refer to the original source cited by Najjar et al.12 for verification. The calculated cohesion parameters for the three cubic EOS’s are shown in Figure 4. For the van der Waals EOS, these are seen to become negative at reduced temperatures above 2.4. This is contrary to the traditional interpretation of the cohesion term as a measure of molecular attractions between fluid molecules, according to which cohesion parameters should be constrained to positive values only. As shown elsewhere17, negative cohesion parameters are in fact a direct consequence of the theoretical inadequacies of
Table 2 Comparison of inversion pressures predicted from several models Experimental12
Tr
Pr 1.010 1.145 1.290 1.474 1.728 1.971 2.268 2.494 2.824 3.158 3.482 3.817 4.211 4.535
5.376 6.730 8.068 9.349 10.701 11.550 11.850 11.750 11.350 10.445 9.217 7.680 5.450 3.331
MVDW12
Pr 6.207 7.452 8.539 9.609 10.626 11.194 11.454 11.386 10.941 10.137 9.069 7.708 5.819 4.060 AAD =
van der Waals (this work)
% error Pr 15.46 10.73 5.83 2.78 − 0.70 − 3.08 − 3.34 − 3.10 − 3.60 − 2.95 − 1.60 0.37 6.77 21.97 5.87
4.760 6.721 8.252 9.599 10.737 11.326 11.622 11.619 11.314 10.634 9.562 7.956 5.349 2.654 AAD =
Redlich–Kwong (this work)
% error Pr − 11.45 − 0.13 2.28 2.67 0.33 − 1.94 − 1.93 − 1.11 − 0.31 1.80 3.74 3.59 − 1.85 − 20.33 3.82
5.002 6.732 8.164 9.509 10.732 11.406 11.746 11.723 11.323 10.519 9.372 7.799 5.412 2.974 AAD =
% error Pr − 6.95 0.02 1.20 1.71 0.29 − 1.25 − 0.87 − 0.23 − 0.24 0.71 1.68 1.54 − 0.70 − 10.71 2.01
MBWR14
Peng–Robinson (this work)
5.033 6.743 8.159 9.491 10.715 11.403 11.766 11.753 11.349 10.515 9.325 7.719 5.374 3.115 AAD =
% error Pr − 6.38 0.20 1.12 1.52 0.13 − 1.27 − 0.71 0.02 − 0.01 0.67 1.17 0.51 − 1.39 − 6.48 1.54
4.328 6.286 7.924 9.443 10.761 11.397 11.570 11.361 10.632 9.476 8.029 6.263 3.879 1.677 AAD =
ALLPROPS20
% error Pr
% error
− 19.49 − 6.60 − 1.79 1.01 0.56 − 1.32 − 2.36 − 3.31 − 6.32 − 9.28 − 12.88 − 18.45 − 28.82 − 49.66 11.56
− 21.04 − 8.90 − 4.44 − 1.59 − 1.54 − 3.13 − 4.25 − 5.50 − 8.91 − 11.96 − 15.27 − 20.29 − 29.85 − 44.73 12.96
4.245 6.131 7.710 9.200 10.536 11.188 11.346 11.104 10.339 9.196 7.810 6.122 3.823 1.841 AAD =
Cryogenics 1998 Volume 38, Number 7 725
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes
Figure 3 Experimantal12 inversion data of air (䊊) and empirical crrelation by ( – – – ) Equation (3) and (———) Equation (16)
Figure 5 Inversion curve of air computed from the van der Waals EOS with new cohesion function of Equation (18). (䊊) Experimental data12
exceed the corresponding experimental value according to the formula ZEOS = 0.857Zexp + 0.084 c c
(17)
= 0.332 almost identical to the which for air gives ZEOS c Redlich–Kwong value. It should also be obvious that the mathematical functions traditionally employed to correlate cohesion parameters obtained from saturation pressures will lack the necessary flexibility to represent the values of Figure 4, among other reasons because they do not allow for negative cohesion parameters. We have tentatively correlated the present results by the following simple extension of the well-known Soave19 formula:
␣ = 1 + m1(T 1/2 − 1) + m2(Tr − 1) + m3(T 3/2 − 1) r r
(18)
+ m4(T − 1) 2 r
Figure 4 New cohesion functions computed from inversion data of air for the (a) van der Waals, (b) Redlich–Kwong and (c) Peng-Robinson EOS’s
the van der Waals repulsion term RT/(v − b), and consistently arise in studies of high-temperature fluid properties such as extrapolated vapour pressures, second virial coefficients, critical isochores, and inversion curves as in the present case. The computed cohesion function for the Peng–Robinson EOS is also non-traditional in that it exhibits a minimum at about Tr = 3.1 and positive derivatives at higher temperatures. It seems clear that in spite of their apparent similarities, the three cubic EOS’s studied here behave quite differently as regards the representation of inversion conditions and therefore require dissimilar supercritical cohesion functions to match the same experimental data. The most “normal” behaviour is that of the Redlich– Kwong EOS, which has monotonically decreasing but positive cohesion parameters over the entire range of inversion temperatures. This may be taken as a further indication that this EOS is the natural choice for a near-simple fluid such as air, as observed in previous studies2,5,7 and put on a quantitative basis by Martin18, who suggested that the critical compressibility factor predicted by a cubic EOS should
726
Cryogenics 1998 Volume 38, Number 7
with coefficients given in Table 1. By construction, Equation (18) satisfies the boundary condition ␣ = 1 at Tr = 1. Figures 5–7 show the inversion curves respectively pre-
Figure 6 Inversion curve of air computed from the Redlich– Kwong EOS with new cohesion function of Equation (18). (䊊) Experimental data12
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes
Figure 7 Inversion curve of air computed from the Peng– Robinson EOS with new cohesion function of Equation (18). (䊊) Experimental data12
dicted for air by the van der Waals, Redlich–Kwong and Peng–Robinson with the new cohesion functions. Some deviation from the experimental data is apparent in these Figures; this must be attributed to the incremental errors introduced at each stage of the computational process, i.e. in substituting the experimental data by the empirical Equation (16), performing the numerical integration over the range of reduced temperatures from 1.0 to 5.0, and correlating the resulting cohesion parameters by Equation (18). Table 2 shows the predicted inversion pressures and their percent deviation from the experimental data. The AAD is 3.82% for the van der Waals, 2.01% for the Redlich– Kwong and 1.54% for the Peng–Robinson EOS. The integration procedure has been carefully tested7 to ensure that computed cohesion parameters are accurate to at least four significant digits. For the van der Waals EOS, for example, when the raw values [i.e. from Figure 4 rather than from Equation (18)] of ␣(Tr ) and d␣/dTr generated by integration are used directly to predict the inversion curve, an AAD of just 0.18% is obtained. It follows that the largest contributing factor to overall deviation is the adoption of Equation (18). Further research into appropriate mathematical models for supercritical cohesion functions is clearly required. All three cubic EOS’s considered in this work predict the inversion curve of air more accurately than does the MVDW EOS. This is hardly surprising given the point-bypoint fitting procedure employed. What is possibly more important is that by following this procedure any simple cubic EOS can be made to perform much better than a more complex noncubic EOS. As an example, inversion pressures were predicted using the 34-constant modified Benedict–Webb–Rubin (MBWR) EOS developed specifically for air by Reynolds14. Since this model may have been superseded by more recent reference standards for air, calculations were also performed using the ALLPROPS v4.2 computer program from the Center for Applied Thermodynamics Studies at the University of Idaho20. Table 2 shows that these two noncubic models in fact yield very similar results. Their AAD exceeds that of the cubic EOS’s by one order of magnitude, with greater errors at the lower inversion pressures. It is natural to enquire whether the new cohesion func-
tions derived from Joule–Thomson data may lead to improved predictions of other closely related properties, such as residual enthalpies, specific heats etc. Any meaningful study in this respect should, for consistency, use the same source of standard reference data both for the inversion points from which the new cohesion functions will be developed, and for the derived thermodynamic properties against which the predictive abilities of the modified cubic EOS will be tested. This is not possible in the present case, because the experimental data12 cover only the inversion curve. As shown in Table 2, these data are in poor agreement with values obtained from the Reynolds and ALLPROPS models; comparison of other properties against these air standards therefore would be of limited significance. The impact of inversion-based cohesion parameters on the prediction of supercritical residual enthalpies has been analysed by Barre´ et al.8 taking as reference the Lee–Kesler corresponding states EOS, and by Castillo et al.21 on the basis of the Boublı´k–Alder–Chen–Kreglewski (BACK) augmented hard body EOS. Both studies show that the new cohesion functions can lead to better enthalpy predictions, with AAD lower than that of the original cubic EOS’s by as much as one order of magnitude, but for best results they should be combined with a volume translation of the EOS. Optimal inversion-based translation parameters and generalised supercritical cohesion functions for use with the Redlich–Kwong and Peng–Robinson EOS’s have been obtained21 for 52 fluids, including hydrocarbons, nonhydrocarbons and alcohols, whose BACK EOS parameters were available in the TRC Thermodynamic Tables22,23. Volume translation, however, has the undesirable effect of degrading the prediction of compressibility factors21, as might be expected from the well known limitations of the cubic models. Finally, we return to the point of thermodynamic consistency by presenting in Figure 8 the saturated liquid and vapour volumes predicted by the van der Waals EOS with the new cohesion function of Equation (18). These have again been computed at bubble and dew pressures and temperatures taken from the literature14, the EOS being used solely to calculate the corresponding phase volumes. Phase envelopes taken from Reynolds14, and also computed with
Figure 8 Saturated liquid (䊊) and vapour (䊏) volumes predicted by the van der Waals EOS, with new cohesion function of Equation (18), compared with reference values (———) from Reynolds14 and ( – – – ) from ALLPROPS
Cryogenics 1998 Volume 38, Number 7 727
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes ALLPROPS20 are included in Figure 8 for comparison purposes; very close agreement between these two reference sources is again evident. Since in computing our inversionorientated cohesion parameters the critical stability criteria have been kept valid by enforcing the boundary condition ␣ = 1 at Tr = 1, the EOS is now able to predict physically correct phases consistently up to critical point. The liquid volumes are of course still too high, but this is to be expected from the van der Waals model with its high predicted Zc. It should in particular be noticed that no subcritical inversion data were used in the development of Equation (18), i.e. that integration was performed upwards in temperature from the critical point. Although nothing prevents the method from being applied, with negative stepsizes, to the subcritical portion of the inversion curve, this would cover only a limited range of subcritical temperatures (down to the minimum reduced inversion temperature, which is typically around 0.85). Castillo et al.21 found that the new cohesion functions, when extrapolated below the critical temperature, very closely matched the traditional Soave functions developed from saturation data. It seems preferable therefore in this region to continue to rely on the VLEbased cohesion functions, moreover since these usually give satisfactory inversion predictions at low temperatures6. In Figure 8, nevertheless, Equation (18) has been extrapolated down to reduced temperatures as low as 0.577 without any adverse results, which gives further proof of the flexibility and consistency of the present approach.
Conclusions A modified van der Waals EOS previously proposed in the literature for accurate prediction of the Joule–Thomson inversion curve of air has been shown to be thermodynamically inconsistent, violating the stability criteria at the critical point and failing to predict the correct liquid phase volumes on the phase envelope at high saturation temperatures. We have presented instead a computational procedure that enables any cubic EOS to match the experimental inversion curve of a given fluid, and have applied it to fit the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s to experimental inversion data of air. The method involves the computation of cohesion parameters based on the thermodynamic criterion for inversion, and should be particularly useful as a means to develop new cohesion functions valid in the supercritical region.
References 1. Gunn, R.D., Chueh, P.L. and Prausnitz, J.M., Inversion temperatures and pressures for cryogenic gases and their mixtures. Cryogenics, 1966, 6, 324–329.
728
Cryogenics 1998 Volume 38, Number 7
2. Miller, D.G., Joule–Thomson inversion curve, corresponding states, and simpler equations of state. Ind. Eng. Chem. Fundam., 1970, 9(4), 585–589. 3. de Reuck, K.M. Fluorine. IUPAC International Thermodynamic Tables of the Fluid State—11. Blackwell, Oxford, 1990. 4. de Reuck, K.M. and Craven, R.J.B. Methanol. IUPAC International Thermodynamic Tables of the Fluid State—12. Blackwell, Oxford, 1993. 5. Juris, K. and Wenzel, L.A., A study of inversion curves. AIChE J., 1972, 18(4), 684–688. 6. Dilay, G.W. and Heidemann, R.A., Calculation of Joule–Thomson inversion curves from equations of state. Ind. Eng. Chem. Fundam., 1986, 25, 152–158. 7. Colazo, A.V., Da Silva, F.A., Mu¨ller, E.A. and Olivera-Fuentes, C., Joule–Thomson inversion curves and the supercritical parameters of cubic equations of state. Latin American Applied Research, 1992, 22, 135–147. 8. Barre´, Y., Colazo, A.V., Colina, C., Da Silva, F.A., Fuentes, J., Santos, J.W. and Olivera-Fuentes, C. An enthalpy-oriented criterion to determine the supercritical cohesion parameters of cubic equations of state. Paper presented at the 12th Symposium on Thermophysical Properties, Boulder, CO, 1994. 9. Leibovici, C.F., Variant and invariant properties from cubic equations of state. Fluid Phase Equilib., 1993, 84, 1–8. 10. Geana, D. and Feroiu, V., Calculation of Joule Thomson inversion curves from a general cubic equation of state. Fluid Phase Equilib., 1992, 77, 121–132. 11. Colina, C.M., Da Silva, F.A., Santos, J.W. and Olivera-Fuentes, C. Volume translation and Joule Thomson inversion curves predicted from cubic equations of state. Informacio´n Tecnolo´gica 9(1), in press (in Spanish) 1998. 12. Najjar, Y.S.H., Al-Beirutty, M.H. and Ismail, M.S., Two-constant equation of state for accurate prediction of the Joule Thomson inversion curve for air cryogenic applications. Cryogenics, 1993, 33(2), 169–174. 13. Martin, J.J. and Hou, Y.-C., Development of an equation of state for gases. AIChE J., 1955, 1, 142–151. 14. Reynolds, W.C. Thermodynamic Properties in SI Units. Department of Mechanical Engineering, Stanford University, Stanford, CA, 1979. 15. Partington, J.R. Advanced Treatise on Physical Chemistry, Vol. I. Longmans, Green, London, 1949. 16. Olivera-Fuentes, C. Thermodynamics and Transport Phenomena in Fluid Continua (in Spanish). Universidad Simo´n Bolı´var, Caracas, 1996. 17. Colina, C.M., Santos, J.W. and Olivera-Fuentes, C., High-temperature behaviour of the cohesion parameter of cubic equations of state. High Temperatures—High Pressures, 1997, 29(5), 525–532. 18. Martin, J.J., Cubic equations of state—which?. Ind. Eng. Chem. Fundam., 1979, 18(2), 81–97. 19. Soave, G., Equilibrium constants from a modified Redlich–Kwong equation of state. Chem. Engng. Sci., 1972, 27, 1197–1203. 20. Lemmon, E.W., Jacobsen, R.T., Penoncello, S.G. and Beyerlein, S.W. Computer programs for calculating thermodynamic properties of fluids of engineering interest. Center for Applied Thermodynamic Studies CATS Report 97-1. University of Idaho, Moscow, ID, 1997. 21. Castillo, M.G., Colina, C.M., Dubuc, J. and Olivera-Fuentes, C.G. Prediction of supercritical fluid properties from cubic equations of state using inversion-based cohesion functions. Submitted for publication, 1998. 22. TRC Thermodynamic Tables—Hydrocarbons, Thermodynamics Research Center. Texas A&M University, College Station, Texas, 1986. 23. TRC Thermodynamic Tables—Non-Hydrocarbons, Thermodynamics Research Center. Texas A&M University, College Station, Texas, 1986.
PII: S0011-2275(98)00036-8
Prediction of the Joule–Thomson inversion curve of air from cubic equations of state Coray M. Colina and Claudio Olivera-Fuentes* TADiP Group, Thermodynamics and Transport Phenomena Department, Simo´n Bolı´var University, AP 89000, Caracas 1086-A, Venezuela
Accepted 4 February 1998 A modified van der Waals equation of state recommended in the literature for improved prediction of the inversion curve of air is shown to be thermodynamically inconsistent, giving large errors in the critical and two-phase regions. An alternative procedure is presented by means of which the cohesion function of any cubic equation of state can be adjusted to give arbitrarily accurate representation of an experimental inversion curve. New versions of the van der Waals, Redlich–Kwong and Peng– Robinson equations of state are developed based on experimental inversion data of air, and are shown to give better inversion predictions than more complex, multiparameter noncubic equations of state. 1998 Elsevier Science Ltd. All rights reserved Keywords: E. Joule–Thomson effect; equation of state; liquid air
Nomenclature
vR
a
VLE Z
attraction parameter in equation of state (kPa.m6.kg−2 ) AAD average absolute deviation b excluded volume parameter in equation of state (m3.kg−1 ) specific heat capacity at constant pressure cP (kJ.kg−1.K−1 ) EOS equation of state h specific enthalpy (kJ.kg−1 ) k1, k2 dimensionless constants in equation of state m1, …, m4 adjustable coefficients in new cohesion function MBWR modified Benedict–Webb–Rubin equation of state MVDW modified van der Waals equation of state P pressure (kPa) gas constant (kJ.kg−1.K−1 ) R s integration stepsize T temperature (K) specific volume (m3.kg−1 ) v
The Joule–Thomson effect is of considerable importance in refrigeration and gas liquefaction processes. Under appropriate conditions, adiabatic throttling of a gas will cause it to lower its temperature, an endothermic effect that *To whom correspondence should be addressed. Tel.: + 58-2906-37-56; Fax: + 58-2-906-37-43; e-mail: [email protected]
dimensionless volume, defined by Equation (9) vapour–liquid equilibrium compressibility factor, Pv/(RT)
Subscripts c r
critical value reduced, relative to critical value
Superscripts exp
experimental
Greek letters ␣
cohesion parameter in equation of state, a/ac Joule–Thomson coefficient (K.kPa−1 ) acentric factor dimensionless critical attraction parameter, acPc/(RTc )2 dimensionless excluded volume parameter, bPc/(RTc )
⍀ac ⍀b
may be strong enough even to induce condensation. For this cooling to take place, the thermodynamic state of the fluid must lie in the region bounded by the inversion curve or locus of points at which the Joule–Thomson coefficient is zero,
⬅
冉冊 ∂T ∂P
= 0.
(1)
h
Cryogenics 1998 Volume 38, Number 7 721
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes As usually represented in pressure–temperature co-ordinates, the Joule–Thomson inversion curve extends from a minimum temperature corresponding to a saturated state on the vapour pressure line, to a maximum temperature corresponding to the ideal-gas limit at zero density and pressure. The curve is parabolic in shape and exhibits a maximum or peak inversion pressure at an intermediate temperature. Knowledge of the inversion curve is obviously essential for good design and operation of Joule–Thomson cooling processes. Direct measurement of inversion points is difficult and unreliable. At near-inversion conditions, the vanishing Joule–Thomson coefficients imply that even large pressure changes will produce very small temperature differences, detectable only by extremely accurate measurement. It is generally preferable to make use of thermodynamic relations such as
=
冋冉 冊 册
∂v 1 T cP ∂T
−v =
P
冉冊
RT 2 ∂Z cPP ∂T
(2)
P
to derive Joule–Thomson coefficients and inversion points from experimental P–v–T data either by direct numerical treatment of the raw data, or more usually by first fitting the data by an equation of state (EOS). Even then, the inversion curves of the heavier fluids cannot be completely determined, because they extend into regions of high temperatures and pressures not attainable in experimental study. Efforts have been made to obtain generalised inversion curves from the known behaviour of light fluids. Gunn et al.1 calculated inversion points from volumetric data of argon, carbon monoxide, ethylene, methane, nitrogen and xenon. Since the data did not cover the upper temperature branch of the inversion curve, they computed additional theoretical points for argon from a truncated virial EOS with second and third virial coefficients based on the Kihara intermolecular potential. The entire set of points was then correlated by the empirical equation Pr = − 36.275 + 71.598Tr − 41.567T 2r + 11.826T 3r
(3)
− 1.6721T + 0.091167T . 4 r
5 r
An alternative, somewhat simpler equation based on published inversion data for argon, ammonia, carbon dioxide, carbon monoxide, ethylene, methane, nitrogen and propane was proposed by Miller2, Pr = 24.21 −
18.54 − 0.825T 2r. Tr
(4)
Equations (3) and (4) are two-parameter corresponding states correlations that give comparable results, since they derive from data for the same or very similar fluids, described1 as gases that have small, nearly spherical molecules with acentric factor close to zero, and2 as mostly gases that have a critical compressibility factor Zc of 0.29. The computation of Joule–Thomson inversion curves is one of the most demanding tests that can be performed on a set of experimental data1 or on an analytical EOS2. For pressure-explicit EOS’s, an equivalent formulation of the thermodynamic inversion criterion of Equations (1) and (2),
722
Cryogenics 1998 Volume 38, Number 7
冉冊 冉冊
T
∂P ∂T
+v
v
∂P ∂v
=0
(5)
T
shows that the EOS must accurately predict not only the P(T,v) function but also its first partial derivatives with respect to temperature and volume (or density), over wide temperature and pressure ranges that can exceed reduced values Tr = 5 and Pr = 12. Even state-of-the-art multiparameter EOS’s may fail to generate physically correct inversion curves when extrapolated to such extreme conditions3,4. Several authors have therefore used the prediction of inversion curves as a criterion to assess and rank the predictive capabilities of available EOS’s. Miller2 studied the van der Waals, Dieterici–Porter, Redlich–Kwong and Clausius–Martin EOS’s. Using Equation (4) as a basis of comparison, he concluded that the original Redlich– Kwong EOS was the best of the simple equations even though it underpredicted the peak inversion pressure by about 10%. He also suggested that analysis of the predicted inversion curve should become a standard test of any newly constructed EOS. Juris and Wenzel5 apparently unaware of Miller’s work, compared Equation (3) with predictions from the van der Waals, Dieterici, Berthelot, Redlich–Kwong, truncated virial (with second and third virial coefficients taken respectively from the Pitzer and Chueh–Prausnitz correlations), Su–Beattie–Bridgeman, several reduced forms of the Benedict–Webb–Rubin, and Martin–Hou EOS’s. They concluded that the van der Waals, Dieterici and Berthelot EOS’s had as might be expected very limited capability of predicting inversion behaviour. Among the more complex models the virial, Beattie–Bridgeman and Benedict– Webb–Rubin EOS’s gave poor predictions; only the Martin–Hou EOS yielded reasonable inversion values, although it could incorrectly predict two peaks of inversion pressure. Remarkably good results were obtained with the Redlich–Kwong EOS, which despite its relative simplicity predicted the inversion locus more accurately than any of the other cubic and noncubic models, even in the liquid region at more than twice the critical density. Dilay and Heidemann6 studied the Soave–Redlich– Kwong, Peng–Robinson, perturbed hard-chain and Lee– Kesler EOS’s, carrying out detailed comparisons for fluids with acentric factors ranging from 0 to 0.25. For simple fluids ( = 0) and nitrogen ( = 0.04) they observed that all four EOS’s gave good predictions of the low-temperature branch of the inversion curve up to a reduced pressure of about 8.0. The Soave–Redlich–Kwong EOS gave the best match of the peak inversion pressure, but failed to reproduce the high-temperature inversion region where only the Lee–Kesler EOS gave reasonably good results. Conclusions for carbon dioxide ( = 0.225) were similar, although less clear because of the greater scatter of the experimental high-temperature data. These authors remarked on the quite different results obtained from the Soave–Redlich–Kwong and Peng–Robinson EOS’s, which they found somewhat surprising given the similarities between these equations. Colazo et al.7 compared Equations (3) and (4) with predictions from the van der Waals, Soave–van der Waals, Redlich–Kwong, Peng–Robinson, Clausius–Martin, Schmidt–Wenzel, Valderrama–Cisternas, Trebble–Bishnoi and Trebble–Bishnoi–Salim EOS’s. They showed that for a generic cubic EOS of the form
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes P=
ac␣(Tr ) RT − 2 v − b v + k1bv + k2b2
(6)
the cohesion function ␣(Tr ) in the numerator of the second term is in principle responsible for the predicted lower and upper inversion temperatures, and the quadratic volume function in the denominator of the same term similarly dominates the predicted peak inversion pressure, which increases as the predicted critical compressibility factor Zc decreases. Cubic EOS’s having substantially identical cohesion functions but different predicted Zc will generate inversion loci that are also different. Since this is the case for the Soave–Redlich–Kwong (Zc = 1/3) and Peng– Robinson (Zc = 0.307) EOS’s, the simple explanation for the “somewhat surprising” differences in performance remarked upon earlier6 is that while these equations may indeed be similar for vapour–liquid equilibrium (VLE) computations, the similarity does not extend to inversion calculations because Joule–Thomson coefficients (unlike vapour pressures) are not translation-invariant properties8,9. The above results contrast sharply with conclusions reached by Geana and Feroiu10, who compared inversion curves predicted from a new four-parameter cubic EOS developed by themselves against values from the Redlich– Kwong, Soave–Redlich–Kwong and Peng–Robinson EOS’s. These authors stated that the translation in volume of a cubic EOS does not change the Pr − Tr inversion curve of the original equation, leading only to improvement in volume values that are displaced by a constant amount. This statement has subsequently been refuted11. It was also shown by Colazo et al.7 that cubic EOS’s which employ cohesion functions derived from pure-fluid (subcritical) saturation data cannot give good prediction when extended to the high (supercritical) temperatures typical of the upper branch of the inversion curve. This is partly due to extrapolation uncertainties, but the main factor is that the inversion behaviour of fluids at these conditions is dominated by the lower-order virial coefficients. This is why the original Redlich–Kwong and the Clausius–Martin EOS’s, both of which were developed to give improved representation of second virial coefficients, predict the upper inversion curve considerably better than the VLEorientated Soave–Redlich–Kwong and Peng–Robinson EOS’s. These authors suggested that new cohesion functions be developed from experimental data by solving Equations (5) and (6) simultaneously, for which they proposed a numerical procedure that will be described in greater detail in the next section. A different approach was followed by Najjar et al.12, who developed a two-constant cubic EOS identical in shape to the original van der Waals EOS but with new constants obtained by least-squares fit of experimental inversion data for air. Their modified van der Waals (MVDW) EOS can be written in SI units (kPa, m3.kg−1, K) as P=
ac RT − 2 , ac = 0.075228, b = 0.000762 v−b v
where vR is the dimensionless volume vR ⬅
Pcv . RTc
(9)
These authors found that the MVDW EOS predicted the inversion curve of air much better than other simple cubic EOS’s such as the original van der Waals, Redlich–Kwong, Soave–Redlich–Kwong and Peng–Robinson EOS’s, and reported an average standard deviation of 0.18 for the inversion curve. Inspection of their tabulated results shows an average absolute error of 5.87% and a maximum error of 21.78% in reduced inversion pressure. Unfortunately, by using modified constant parameters rather than introducing a temperature dependence for them as is usual nowadays, Najjar et al.12 made their MVDW thermodynamically inconsistent. In order to see this it suffices to compute the predicted critical point from this EOS, by setting the first and second partial derivatives of pressure with respect to volume equal to zero, or equivalently by requiring the EOS to give three equal volume roots at the critical point13. Doing this with Equation (7) yields Zc = 0.375 as expected for the van der Waals model, but with critical constants (in SI units) Pc =
ac 8ac = 101.9, = 4798.5, Tc = 2 27b 27bR
vc = 3b = 0.002286 which deviate substantially from the accepted values, e.g. Pc = 3770, Tc = 132.5, vc = 0.002913 used by Reynolds14. Repeating the calculations with Equation (8) gives the reduced co-ordinates of the MVDW “critical” point as Pr =
⍀ac 8⍀ac = 1.273, Tr = = 0.7693, 27⍀2b 27⍀b
vR = 3⍀b = 0.2266 instead of the correct values Pr = Tr = 1, vR = 0.375. The situation is illustrated in Figure 1, which shows three iso-
(7)
(the value actually reported by them, b = 0.00762 is obviously a misprint), and in reduced form as Pr =
Tr ⍀ac − , ⍀ac = 0.1961205, ⍀b = 0.07553501 vR − ⍀b v2R (8)
Figure 1 Isotherms of the MVDW EOS12 showing incorrect position of the critical point. (a) Tr = 1.0; (b) Tr = 0.7693; (c) Tr = 0.5; C.P. = predicted critical point (coordinates Pr = 1.273, vr = 0.2266)
Cryogenics 1998 Volume 38, Number 7 723
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes therms computed with Equation (8). The Tr = 1.0 isotherm does not exhibit the critical (inflexion) point required by classical thermodynamics; this point appears instead on the Tr = 0.7693 isotherm as explained above. Only below this temperature do the subcritical isotherms exhibit the van der Waals loop required for prediction of coexisting liquid and vapour phases; a third isotherm Tr = 0.5 is drawn as an example. It follows that although the MVDW EOS was intended for improved representation of air properties for cryogenic applications, it will fail precisely at low temperatures, predicting a phase coexistence region that extends over a too narrow temperature range (Tr ⱕ 0.7693) and a too wide pressure range (Pr ⱕ 1.273). It is of course basically incorrect to use a pure-fluid EOS such as the MVDW EOS to predict VLE of air which is a mixture, even if its bubble and dew points differ by only a few kelvins, but if equilibrium pressures and temperatures are taken from a suitable reference14 and Equation (7) is used just to compute the corresponding phase volumes, the results appear as in Figure 2. A liquid branch of the phase envelope is predicted for pressures ⬍ 180 kPa, i.e. at bubble temperatures under 85K only. Although this temperature lies well below the critical temperature predicted by Equation (7), the fact that the MVDW EOS overpredicts the critical pressure (and therefore all saturation pressures) results in only one real, vapour-like volume root at higher temperatures, so that the “liquid” and vapour branches become nearly identical. Since the inversion curve terminates on the liquid side of the saturation envelope, the MVDW EOS can be expected to give particularly poor prediction of air properties at the lower inversion temperatures. It is the purpose of the present work to show how cubic EOS’s can be developed that give an arbitrarily accurate prediction of the inversion curve, better in fact than that obtained from considerably more complex noncubic models, without sacrificing thermodynamic consistency. Examples of three such EOS’s are given, corresponding to modifications of the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s based on inversion data of air.
Development of equations of state from inversion curves The generic cubic EOS, Equation (6) is written in reduced terms as Pr =
Tr ⍀ac␣(Tr ) − 2 vR − ⍀b vR + k1⍀bvR + k2⍀2b
(10)
where the cohesion function ␣(Tr ) expresses the temperature dependence of the attraction term. On applying the inversion condition of Equation (5), the following differential equation results, also in reduced variables7, d␣ ⍀acTr 2 v + k1⍀bvR + k2⍀b dTr 2 R
−
(11)
⍀bTr ⍀acvR(2vR + k1⍀b ) ␣+ = 0. (v2R + k1⍀bvR + k2⍀2b )2 (vR − ⍀b )2
For an already existing cubic EOS, i.e. for known ␣(Tr ), the predicted inversion curve will be given by simultaneous solution of Equations (10) and (11). For the original van der Waals EOS (see Table 1), for example, on setting ␣(Tr ) ⬅ 1 we obtain Pr =
27 Tr − vR − 1/8 64v2R
(12)
and −
Tr 27 + = 0. 2 4vR (vR − 1/8)2
(13)
Equation (13) gives the inversion curve directly in Tr − vR co-ordinates. It may be combined with Equation (12) to give an equivalent Pr − vR version and also the well-known Pr − Tr formula15
冋
Tr = 3 1 ±
册
2 1 √9 − Pr . 6
(14)
Explicit or algorithmic solutions for other EOS’s have been presented by a number of authors, e.g. Gunn et al.1 for the truncated virial EOS, Miller2 for the Dieterici, Redlich– Kwong and Martin EOS’s, Dilay and Heidemann6 for the Soave–Redlich–Kwong and Peng–Robinson EOS’s, Colazo et al.7 and Olivera-Fuentes16 for the general cubic EOS of Equation (10), and Najjar et al.12 again for the van der Waals EOS. Conversely, given an inversion curve Pr(Tr ) the cohesion function may be computed by treating Equation (11) as a first order ordinary differential equation for ␣(Tr ), with vR an implicit function of Tr given by Equation (10). Since by definition ␣(Tr ) = 1 at Tr = 1, this is an initial-value problem that may solved by standard numerical methods. Colazo et al.7 used fourth-order Runge–Kutta integration with extrapolation to zero stepsize; their algorithm is also used in the present work and may be summarised as follows: Figure 2 Saturated liquid (䊊) and vapour (䊏) volumes predicted by the MVDW EOS12, compared with reference values (———) from Reynolds14
724
Cryogenics 1998 Volume 38, Number 7
1. Start at Tr = 1 with ␣ = 1. Set integration interval ⌬Tr, e.g. ⌬Tr = 0.1; 2. Obtain ␣(1) at Tr + ⌬Tr by Runge–Kutta integration
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes Table 1 Equations of state developed in this work EOS
EOS parameters
van der Waals Redlich– Kwong Peng– Robinson
Coefficients of new cohesion function
k1
k2
⍀ac
⍀b
0
0
27/64
1/8
− 14.1760
12.02581
1
0
0.42748…
0.08664…
− 5.72678
0.45723…
0.07779…
− 5.13662
−1
2
m1
with stepsize s(1), e.g. s(1) = 0.001. At each intermediate temperature, (a) Compute Pr from the known inversion curve Pr(Tr ), (b) Compute vR(Pr, Tr ) from the EOS Equation (10), (c) Compute d␣/dTr from the inversion condition Equation (11);
m4
R2
− 4.59465
0.514251
0.999985
4.404533
− 1.68415
0.235894
0.999941
4.501821
− 1.96490
0.345253
0.999421
m2
m3
comparison with the MVDW EOS proposed by those authors. Because the integration procedure requires a continuous representation of the entire inversion curve, the data were fitted by a fifth-degree polynomial similar to Equation (3), Pr = − 13.285 + 26.738Tr − 9.3407T 2r + 1.1469T 3r − 0.016527T − 0.0065313T 4 r
3.
Repeat the above computations to obtain a second estimate ␣(2) using stepsize s(2), e.g. s(2) = s(1)/2; Obtain ␣(Tr + ⌬Tr ) by linear extrapolation to zero stepsize,
4.
␣=
␣(1)s(2) − ␣(2)s(1) s(2) − s(1)
5.
Reset Tr to Tr + ⌬Tr and repeat from step (2).
(15)
The ad hoc extrapolation procedure of Equation (15), based on the observation that integration error approaches zero as a linear function of stepsize7, is essential to eliminate the accumulation of errors which would very rapidly invalidate the solution. The final accuracy of the computed cohesion function will in any case be assessed by predicting the entire inversion curve and comparing it with the original data. Using the above method, we have obtained new supercritical cohesion functions for air for use with the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s. Experimental inversion data in reduced variables Pr − Tr (see Table 2) were taken from Najjar et al.12 to allow direct
(16)
5 r
with regression coefficient R2 = 0.99991. Equation (16) represents the experimental data in Table 2 with an average absolute deviation (AAD) of 0.30%. A function similar to Equation (4) was also tried, but gave considerably worse results. Figure 3 shows the experimental data, the fitting equation, and the corresponding curve predicted from the generalised Equation (3), which is seen to deviate noticeably from the data, especially in the peak pressure and upper temperature regions. Since air (Zc = 0.289, = − 0.0037 based on bubble points14 ) may be expected to behave as a simple fluid for which Equation (3) should be applicable, these deviations might reflect some inaccuracy in the experimental data. We have been unable to refer to the original source cited by Najjar et al.12 for verification. The calculated cohesion parameters for the three cubic EOS’s are shown in Figure 4. For the van der Waals EOS, these are seen to become negative at reduced temperatures above 2.4. This is contrary to the traditional interpretation of the cohesion term as a measure of molecular attractions between fluid molecules, according to which cohesion parameters should be constrained to positive values only. As shown elsewhere17, negative cohesion parameters are in fact a direct consequence of the theoretical inadequacies of
Table 2 Comparison of inversion pressures predicted from several models Experimental12
Tr
Pr 1.010 1.145 1.290 1.474 1.728 1.971 2.268 2.494 2.824 3.158 3.482 3.817 4.211 4.535
5.376 6.730 8.068 9.349 10.701 11.550 11.850 11.750 11.350 10.445 9.217 7.680 5.450 3.331
MVDW12
Pr 6.207 7.452 8.539 9.609 10.626 11.194 11.454 11.386 10.941 10.137 9.069 7.708 5.819 4.060 AAD =
van der Waals (this work)
% error Pr 15.46 10.73 5.83 2.78 − 0.70 − 3.08 − 3.34 − 3.10 − 3.60 − 2.95 − 1.60 0.37 6.77 21.97 5.87
4.760 6.721 8.252 9.599 10.737 11.326 11.622 11.619 11.314 10.634 9.562 7.956 5.349 2.654 AAD =
Redlich–Kwong (this work)
% error Pr − 11.45 − 0.13 2.28 2.67 0.33 − 1.94 − 1.93 − 1.11 − 0.31 1.80 3.74 3.59 − 1.85 − 20.33 3.82
5.002 6.732 8.164 9.509 10.732 11.406 11.746 11.723 11.323 10.519 9.372 7.799 5.412 2.974 AAD =
% error Pr − 6.95 0.02 1.20 1.71 0.29 − 1.25 − 0.87 − 0.23 − 0.24 0.71 1.68 1.54 − 0.70 − 10.71 2.01
MBWR14
Peng–Robinson (this work)
5.033 6.743 8.159 9.491 10.715 11.403 11.766 11.753 11.349 10.515 9.325 7.719 5.374 3.115 AAD =
% error Pr − 6.38 0.20 1.12 1.52 0.13 − 1.27 − 0.71 0.02 − 0.01 0.67 1.17 0.51 − 1.39 − 6.48 1.54
4.328 6.286 7.924 9.443 10.761 11.397 11.570 11.361 10.632 9.476 8.029 6.263 3.879 1.677 AAD =
ALLPROPS20
% error Pr
% error
− 19.49 − 6.60 − 1.79 1.01 0.56 − 1.32 − 2.36 − 3.31 − 6.32 − 9.28 − 12.88 − 18.45 − 28.82 − 49.66 11.56
− 21.04 − 8.90 − 4.44 − 1.59 − 1.54 − 3.13 − 4.25 − 5.50 − 8.91 − 11.96 − 15.27 − 20.29 − 29.85 − 44.73 12.96
4.245 6.131 7.710 9.200 10.536 11.188 11.346 11.104 10.339 9.196 7.810 6.122 3.823 1.841 AAD =
Cryogenics 1998 Volume 38, Number 7 725
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes
Figure 3 Experimantal12 inversion data of air (䊊) and empirical crrelation by ( – – – ) Equation (3) and (———) Equation (16)
Figure 5 Inversion curve of air computed from the van der Waals EOS with new cohesion function of Equation (18). (䊊) Experimental data12
exceed the corresponding experimental value according to the formula ZEOS = 0.857Zexp + 0.084 c c
(17)
= 0.332 almost identical to the which for air gives ZEOS c Redlich–Kwong value. It should also be obvious that the mathematical functions traditionally employed to correlate cohesion parameters obtained from saturation pressures will lack the necessary flexibility to represent the values of Figure 4, among other reasons because they do not allow for negative cohesion parameters. We have tentatively correlated the present results by the following simple extension of the well-known Soave19 formula:
␣ = 1 + m1(T 1/2 − 1) + m2(Tr − 1) + m3(T 3/2 − 1) r r
(18)
+ m4(T − 1) 2 r
Figure 4 New cohesion functions computed from inversion data of air for the (a) van der Waals, (b) Redlich–Kwong and (c) Peng-Robinson EOS’s
the van der Waals repulsion term RT/(v − b), and consistently arise in studies of high-temperature fluid properties such as extrapolated vapour pressures, second virial coefficients, critical isochores, and inversion curves as in the present case. The computed cohesion function for the Peng–Robinson EOS is also non-traditional in that it exhibits a minimum at about Tr = 3.1 and positive derivatives at higher temperatures. It seems clear that in spite of their apparent similarities, the three cubic EOS’s studied here behave quite differently as regards the representation of inversion conditions and therefore require dissimilar supercritical cohesion functions to match the same experimental data. The most “normal” behaviour is that of the Redlich– Kwong EOS, which has monotonically decreasing but positive cohesion parameters over the entire range of inversion temperatures. This may be taken as a further indication that this EOS is the natural choice for a near-simple fluid such as air, as observed in previous studies2,5,7 and put on a quantitative basis by Martin18, who suggested that the critical compressibility factor predicted by a cubic EOS should
726
Cryogenics 1998 Volume 38, Number 7
with coefficients given in Table 1. By construction, Equation (18) satisfies the boundary condition ␣ = 1 at Tr = 1. Figures 5–7 show the inversion curves respectively pre-
Figure 6 Inversion curve of air computed from the Redlich– Kwong EOS with new cohesion function of Equation (18). (䊊) Experimental data12
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes
Figure 7 Inversion curve of air computed from the Peng– Robinson EOS with new cohesion function of Equation (18). (䊊) Experimental data12
dicted for air by the van der Waals, Redlich–Kwong and Peng–Robinson with the new cohesion functions. Some deviation from the experimental data is apparent in these Figures; this must be attributed to the incremental errors introduced at each stage of the computational process, i.e. in substituting the experimental data by the empirical Equation (16), performing the numerical integration over the range of reduced temperatures from 1.0 to 5.0, and correlating the resulting cohesion parameters by Equation (18). Table 2 shows the predicted inversion pressures and their percent deviation from the experimental data. The AAD is 3.82% for the van der Waals, 2.01% for the Redlich– Kwong and 1.54% for the Peng–Robinson EOS. The integration procedure has been carefully tested7 to ensure that computed cohesion parameters are accurate to at least four significant digits. For the van der Waals EOS, for example, when the raw values [i.e. from Figure 4 rather than from Equation (18)] of ␣(Tr ) and d␣/dTr generated by integration are used directly to predict the inversion curve, an AAD of just 0.18% is obtained. It follows that the largest contributing factor to overall deviation is the adoption of Equation (18). Further research into appropriate mathematical models for supercritical cohesion functions is clearly required. All three cubic EOS’s considered in this work predict the inversion curve of air more accurately than does the MVDW EOS. This is hardly surprising given the point-bypoint fitting procedure employed. What is possibly more important is that by following this procedure any simple cubic EOS can be made to perform much better than a more complex noncubic EOS. As an example, inversion pressures were predicted using the 34-constant modified Benedict–Webb–Rubin (MBWR) EOS developed specifically for air by Reynolds14. Since this model may have been superseded by more recent reference standards for air, calculations were also performed using the ALLPROPS v4.2 computer program from the Center for Applied Thermodynamics Studies at the University of Idaho20. Table 2 shows that these two noncubic models in fact yield very similar results. Their AAD exceeds that of the cubic EOS’s by one order of magnitude, with greater errors at the lower inversion pressures. It is natural to enquire whether the new cohesion func-
tions derived from Joule–Thomson data may lead to improved predictions of other closely related properties, such as residual enthalpies, specific heats etc. Any meaningful study in this respect should, for consistency, use the same source of standard reference data both for the inversion points from which the new cohesion functions will be developed, and for the derived thermodynamic properties against which the predictive abilities of the modified cubic EOS will be tested. This is not possible in the present case, because the experimental data12 cover only the inversion curve. As shown in Table 2, these data are in poor agreement with values obtained from the Reynolds and ALLPROPS models; comparison of other properties against these air standards therefore would be of limited significance. The impact of inversion-based cohesion parameters on the prediction of supercritical residual enthalpies has been analysed by Barre´ et al.8 taking as reference the Lee–Kesler corresponding states EOS, and by Castillo et al.21 on the basis of the Boublı´k–Alder–Chen–Kreglewski (BACK) augmented hard body EOS. Both studies show that the new cohesion functions can lead to better enthalpy predictions, with AAD lower than that of the original cubic EOS’s by as much as one order of magnitude, but for best results they should be combined with a volume translation of the EOS. Optimal inversion-based translation parameters and generalised supercritical cohesion functions for use with the Redlich–Kwong and Peng–Robinson EOS’s have been obtained21 for 52 fluids, including hydrocarbons, nonhydrocarbons and alcohols, whose BACK EOS parameters were available in the TRC Thermodynamic Tables22,23. Volume translation, however, has the undesirable effect of degrading the prediction of compressibility factors21, as might be expected from the well known limitations of the cubic models. Finally, we return to the point of thermodynamic consistency by presenting in Figure 8 the saturated liquid and vapour volumes predicted by the van der Waals EOS with the new cohesion function of Equation (18). These have again been computed at bubble and dew pressures and temperatures taken from the literature14, the EOS being used solely to calculate the corresponding phase volumes. Phase envelopes taken from Reynolds14, and also computed with
Figure 8 Saturated liquid (䊊) and vapour (䊏) volumes predicted by the van der Waals EOS, with new cohesion function of Equation (18), compared with reference values (———) from Reynolds14 and ( – – – ) from ALLPROPS
Cryogenics 1998 Volume 38, Number 7 727
Prediction of the Joule–Thomson inversion curve of air: Coray M. Colina and Claudio Olivera-Fuentes ALLPROPS20 are included in Figure 8 for comparison purposes; very close agreement between these two reference sources is again evident. Since in computing our inversionorientated cohesion parameters the critical stability criteria have been kept valid by enforcing the boundary condition ␣ = 1 at Tr = 1, the EOS is now able to predict physically correct phases consistently up to critical point. The liquid volumes are of course still too high, but this is to be expected from the van der Waals model with its high predicted Zc. It should in particular be noticed that no subcritical inversion data were used in the development of Equation (18), i.e. that integration was performed upwards in temperature from the critical point. Although nothing prevents the method from being applied, with negative stepsizes, to the subcritical portion of the inversion curve, this would cover only a limited range of subcritical temperatures (down to the minimum reduced inversion temperature, which is typically around 0.85). Castillo et al.21 found that the new cohesion functions, when extrapolated below the critical temperature, very closely matched the traditional Soave functions developed from saturation data. It seems preferable therefore in this region to continue to rely on the VLEbased cohesion functions, moreover since these usually give satisfactory inversion predictions at low temperatures6. In Figure 8, nevertheless, Equation (18) has been extrapolated down to reduced temperatures as low as 0.577 without any adverse results, which gives further proof of the flexibility and consistency of the present approach.
Conclusions A modified van der Waals EOS previously proposed in the literature for accurate prediction of the Joule–Thomson inversion curve of air has been shown to be thermodynamically inconsistent, violating the stability criteria at the critical point and failing to predict the correct liquid phase volumes on the phase envelope at high saturation temperatures. We have presented instead a computational procedure that enables any cubic EOS to match the experimental inversion curve of a given fluid, and have applied it to fit the van der Waals, Redlich–Kwong and Peng–Robinson EOS’s to experimental inversion data of air. The method involves the computation of cohesion parameters based on the thermodynamic criterion for inversion, and should be particularly useful as a means to develop new cohesion functions valid in the supercritical region.
References 1. Gunn, R.D., Chueh, P.L. and Prausnitz, J.M., Inversion temperatures and pressures for cryogenic gases and their mixtures. Cryogenics, 1966, 6, 324–329.
728
Cryogenics 1998 Volume 38, Number 7
2. Miller, D.G., Joule–Thomson inversion curve, corresponding states, and simpler equations of state. Ind. Eng. Chem. Fundam., 1970, 9(4), 585–589. 3. de Reuck, K.M. Fluorine. IUPAC International Thermodynamic Tables of the Fluid State—11. Blackwell, Oxford, 1990. 4. de Reuck, K.M. and Craven, R.J.B. Methanol. IUPAC International Thermodynamic Tables of the Fluid State—12. Blackwell, Oxford, 1993. 5. Juris, K. and Wenzel, L.A., A study of inversion curves. AIChE J., 1972, 18(4), 684–688. 6. Dilay, G.W. and Heidemann, R.A., Calculation of Joule–Thomson inversion curves from equations of state. Ind. Eng. Chem. Fundam., 1986, 25, 152–158. 7. Colazo, A.V., Da Silva, F.A., Mu¨ller, E.A. and Olivera-Fuentes, C., Joule–Thomson inversion curves and the supercritical parameters of cubic equations of state. Latin American Applied Research, 1992, 22, 135–147. 8. Barre´, Y., Colazo, A.V., Colina, C., Da Silva, F.A., Fuentes, J., Santos, J.W. and Olivera-Fuentes, C. An enthalpy-oriented criterion to determine the supercritical cohesion parameters of cubic equations of state. Paper presented at the 12th Symposium on Thermophysical Properties, Boulder, CO, 1994. 9. Leibovici, C.F., Variant and invariant properties from cubic equations of state. Fluid Phase Equilib., 1993, 84, 1–8. 10. Geana, D. and Feroiu, V., Calculation of Joule Thomson inversion curves from a general cubic equation of state. Fluid Phase Equilib., 1992, 77, 121–132. 11. Colina, C.M., Da Silva, F.A., Santos, J.W. and Olivera-Fuentes, C. Volume translation and Joule Thomson inversion curves predicted from cubic equations of state. Informacio´n Tecnolo´gica 9(1), in press (in Spanish) 1998. 12. Najjar, Y.S.H., Al-Beirutty, M.H. and Ismail, M.S., Two-constant equation of state for accurate prediction of the Joule Thomson inversion curve for air cryogenic applications. Cryogenics, 1993, 33(2), 169–174. 13. Martin, J.J. and Hou, Y.-C., Development of an equation of state for gases. AIChE J., 1955, 1, 142–151. 14. Reynolds, W.C. Thermodynamic Properties in SI Units. Department of Mechanical Engineering, Stanford University, Stanford, CA, 1979. 15. Partington, J.R. Advanced Treatise on Physical Chemistry, Vol. I. Longmans, Green, London, 1949. 16. Olivera-Fuentes, C. Thermodynamics and Transport Phenomena in Fluid Continua (in Spanish). Universidad Simo´n Bolı´var, Caracas, 1996. 17. Colina, C.M., Santos, J.W. and Olivera-Fuentes, C., High-temperature behaviour of the cohesion parameter of cubic equations of state. High Temperatures—High Pressures, 1997, 29(5), 525–532. 18. Martin, J.J., Cubic equations of state—which?. Ind. Eng. Chem. Fundam., 1979, 18(2), 81–97. 19. Soave, G., Equilibrium constants from a modified Redlich–Kwong equation of state. Chem. Engng. Sci., 1972, 27, 1197–1203. 20. Lemmon, E.W., Jacobsen, R.T., Penoncello, S.G. and Beyerlein, S.W. Computer programs for calculating thermodynamic properties of fluids of engineering interest. Center for Applied Thermodynamic Studies CATS Report 97-1. University of Idaho, Moscow, ID, 1997. 21. Castillo, M.G., Colina, C.M., Dubuc, J. and Olivera-Fuentes, C.G. Prediction of supercritical fluid properties from cubic equations of state using inversion-based cohesion functions. Submitted for publication, 1998. 22. TRC Thermodynamic Tables—Hydrocarbons, Thermodynamics Research Center. Texas A&M University, College Station, Texas, 1986. 23. TRC Thermodynamic Tables—Non-Hydrocarbons, Thermodynamics Research Center. Texas A&M University, College Station, Texas, 1986.