- Email: [email protected]
Fallibility of analytic roots of cubic equations of state in low temperature region
Fluid Phase Equilibria 201 (2002) 287–294
Fallibility of analytic roots of cubic equations of state in low temperature region Yun Zhi a,∗ , Huen Lee b a
b
Department of Chemical Engineering, Separation Engineering Research Center, Nanjing University of Chemical Technology, Nanjing 210009, PR China Department of Chemical Engineering, Korea Advanced Institute of Science and Technology, 373-1 Kusong-dong, Yusong-gu, Taejon 305-701, South Korea Received 11 March 2000; accepted 24 January 2002
Abstract It is showed by several examples that, in low temperature region, the analytical solutions of cubic equation of state (EOS) lead to irrational results, while the iterative solutions of cubic EOS by using the Newton–Raphson method produce valid results. Errors caused by the limitation of significant figures of the computer languages are revealed, and a magnification of errors is defined which is a main factor bringing out the irrational results of the analytical solution of cubic EOS. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Cubic equation of state; Analytic roots; Fallibility; Low temperature
1. Introduction Since van der Walls presented his famous equation of state (EOS) about 100 year ago, various EOS have been developed. Among varieties of type of EOS, the number of cubic EOS may be the largest one. As cubic EOS has a notable advantage over complex EOS in its mathematical simplicity, it has attracted numerous researchers and played an important role in engineering calculations. The complex EOS, such as the celebrated BWR EOS [1] etc., can not be solved analytically for compressibility factors or specific volumes with temperature and pressure specified. The trial or iterative methods are always needed to solve the complex EOS for getting their roots one by one, and after these roots being obtained, to discriminate them is often needed too. On the other hand, for the cubic EOS, analytical solutions can be easily got by using a relevant math formula, i.e. Cardan’s method to solve EOS for their three roots simultaneously, as recommended by some textbooks on thermodynamics ∗
Corresponding author. E-mail address: [email protected] (Y. Zhi). 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 2 ) 0 0 0 7 2 - 9
288
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
[2,3]. Generally speaking, with this merit, the cubic EOS make the programming and calculations of thermodynamics concerned easier and simpler than the complex EOS do. In recent years, some non-cubic EOS have been successfully developed, but no direct solution for them like a cubic EOS was considered as a drawback [4]. However, for phase equilibria calculation involving infeasible conditions, Mathias et al. [5] developed a heuristic and iterative procedure to calculate density roots using EOS, no matter cubic or non-cubic, and the procedure is as efficient as analytical approaches. Furthermore, a detailed study made by Mathias and Benson [6] showed that, for most cases of phase equilibria calculations, the computational load of cubic equations are approximately equivalent to more complicated equations. In the case of a pure component, the cubic equations are more computationally rapid than the more complicated equations, but the difference becomes small for systems with many components. The larger the number of components, the smaller the difference, as most of the CPU time is spent in calculating mixing-rule summations. But, it seems that the valid of analytical roots of cubic EOS is beyond suspicion, and has not been explored in literature yet. Even in the case of a pure component, there is a problem with analytic roots of cubic EOS which was found in low temperature region by authors. The aim of this work is to demonstrate several examples of the fallibility and discuss the reasons which cause the fallibility.
2. Calculations and results Three well known cubic EOSs, PR [7], PT [8], and CCOR EOSs [9] are adopted as examples in this work. In the literature [2,10,11], there are different expressions of the math formula needed for the analytic solution of cubic equations, but the essences of which are in fact the same, i.e. Cardan’s method. In this work, the formula recommended by Tester and Modell [2] is used. As in the temperature and pressure ranges concerned in this work, the cubic EOS will certainly have three real roots, for lucidity, the cases of less than three real roots are not taken into account. Suppose a cubic EOS being arranged as a polynomial form f (z) = z3 + c1 z2 + c2 z + c3 = 0
(1)
where z is compressibility factor, c1 , c2 and c3 the coefficients dependent on temperature and pressure. According to the approach illustrated by Tester and Modell [2], let q = 19 (3c2 − c12 )
(2)
r = 16 c1 c2 − 21 c3 −
1 3 c 27 1
(3)
Then, three real roots are given by arccos r/ −q 3 √ − c1 z1 = 2 −q cos 3 3
√ z2 = 2 −q cos
arccos r/ −q 3 3
(4)
c1 + 120 − 3 ◦
(5)
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
√ z3 = 2 −q cos
arccos r/ −q 3 3
289
c1 ◦ + 240 − 3
(6)
where z2 is the smallest root representing the compressibility factor of liquid phase, z3 is the largest root representing the compressibility factor of vapor phase, z1 has no physical significance. It is not the accuracy of the model we are concerned, but the validity of the compressibility factor calculation. For the facility of discussion, only liquid compressibility factors and liquid volumes are concerned in the following calculation. Liquid molar volume is given by V = z2
RT P
(7)
where R is the gas constant, T the temperature, and P the saturated vapor pressure which, in this work, is given by [12] q2 (8) + q3 ln T + q4 T q5 P = exp q1 + T where q1 , q2 , q3 , q4 and q5 are constants depended on substances. The liquids molar volumes of three pure substances (propylene, 1-butene and 1-pentene) are calculated by the above procedure, and the results are listed in the Tables 1–3, denoted as Vf .
Table 1 Comparison among experimental and calculated liquid volumes for propylene Number
T (K)
Tr
P (MPa)
Ve
PT
PR
CCOR
Vn
Vf
Vn
Vf
Vn
Vf
−2197.4 1008.7 −132.6 135.3
53.91 53.97 54.04 54.10
1002.6 981.03 −159.66 108.88
55.50 55.64 55.78 55.92
521.22 88.626 27.537 63.887
1 2 3 4
87.9 89.4 90.9 92.4
0.241 0.245 0.249 0.253
9.18E−10 1.60E−09 2.74E−09 4.59E−09
55.07 55.18 55.29 55.39
57.94 58.02 58.09 58.16
5 6 7 8
93.9 95.4 96.9 98.4
0.257 0.262 0.266 0.270
7.56E−09 1.22E−08 1.95E−08 3.06E−08
55.50 55.61 55.72 55.84
58.24 58.31 58.39 58.47
109.19 55.98 62.74 57.49
54.17 54.23 54.30 54.37
83.30 71.18 54.68 58.28
56.06 56.20 56.34 56.48
59.85 51.98 56.71 55.90
9 10 11 12
99.9 101.4 102.9 127.9
0.274 0.278 0.282 0.351
4.73E−08 7.21E−08 1.08E−07 2.08E−05
55.95 56.06 56.17 58.17
58.54 58.62 58.70 60.16
60.66 58.30 59.35 60.16
54.44 54.51 54.58 55.87
52.62 54.25 54.43 55.87
56.62 56.75 56.89 59.19
56.43 56.82 56.89 59.19
13 14 15 16
152.9 177.9 202.9 227.9
0.419 0.488 0.556 0.625
6.03E−04 6.09E−03 3.24E−02 1.14E−01
60.39 62.91 65.79 69.16
61.90 64.01 66.62 69.91
61.90 64.01 66.62 69.91
57.42 59.30 61.63 64.58
57.42 59.30 61.63 64.58
61.53 63.98 66.67 69.73
61.53 63.98 66.67 69.73
290
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
Table 2 Comparison among experimental and calculated liquid volumes for 1-butene Number
T (K)
Tr
P (MPa)
Ve
PT
PR
Vn
Vf
Vn
1628.5 1550.2 −370983 1411.1
70.37 70.50 70.63 70.77
CCOR Vf
Vn
Vf
1671.6 1601.8 1536.6 1475.6
69.63 69.95 70.27 70.59
−8310460 −2053680 6370.11 5940.21
1418.4 1364.8 −5679.9 1266.9
70.91 71.22 71.54 71.86
1 2 3 4
87.8 91.3 94.8 98.3
0.209 0.218 0.226 0.234
3.56E−13 1.87E−12 8.58E−12 3.51E−11
70.09 70.35 70.62 70.89
75.79 75.96 76.15 76.34
5 6 7 8
101.8 105.3 108.8 112.3
0.243 0.251 0.259 0.268
1.29E−10 4.32E−10 1.33E−09 3.79E−09
71.16 71.44 71.71 72.00
76.53 76.72 76.92 77.12
1349.1 1291.4 1237.7 −36.61
70.90 71.04 71.19 71.33
9 10 11 12
115.8 119.3 122.8 126.3
0.276 0.284 0.293 0.301
1.01E−08 2.51E−08 5.92E−08 1.32E−07
72.28 72.57 72.86 73.16
77.33 77.55 77.76 77.98
78.65 77.75 75.89 78.12
71.48 71.63 71.78 71.94
71.69 51.00 50.81 76.89
72.17 72.49 72.80 73.12
−91.96 92.06 71.29 72.85
13 14 15 16
161.3 196.3 231.3 266.3
0.384 0.468 0.551 0.635
4.98E−05 1.80E−03 1.90E−02 9.86E−02
76.34 80.00 84.31 89.52
80.49 83.63 87.66 93.01
80.49 83.63 87.66 93.01
72.10 74.20 76.86 80.28
71.94 74.20 76.86 80.28
73.43 77.11 80.94 85.15
73.07 77.11 80.94 85.15
5548.8 5191.4 −2129.9 202.2
Table 3 Comparison among experimental and calculated liquid volumes for 1-pentene Number
T (K)
Tr
P (MPa)
Ve
PT Vn
PR Vf
Vn
1 2 3 4
107.93 111.43 114.93 118.43
0.232 0.240 0.247 0.255
4.48E−12 1.67E−11 5.70E−11 1.79E−10
86.64 86.94 87.25 87.56
91.05 91.23 91.42 91.61
1961.51 1881.09 1805.88 1735.39
89.11 89.28 89.46 89.64
5 6 7 8
121.93 125.43 128.93 132.43
0.262 0.270 0.277 0.285
5.24E−10 1.43E−09 3.68E−09 8.92E−09
87.87 88.19 88.51 88.83
91.80 92.00 92.20 92.40
1669.2 1607.0 99.54 122.2
89.83 90.02 90.21 90.41
9 10 11 12
135.93 139.43 142.93 146.43
0.292 0.300 0.308 0.315
2.06E−08 4.52E−08 9.50E−08 1.92E−07
89.16 89.49 89.83 90.17
92.61 92.82 93.04 93.26
84.82 91.59 92.45 93.18
13 14 15 16
149.93 184.93 219.93 254.93
0.323 0.398 0.473 0.548
3.73E−07 5.93E−05 1.48E−03 1.32E−02
90.51 94.18 98.36 103.21
93.49 96.04 99.20 103.20
93.45 96.04 99.20 103.20
CCOR Vf
Vn
1946.1 1866.4 1791.8 −20549
Vf
86.44 86.86 87.27 87.69
8026.3 7468.2 6961.1 6499.1
1656.3 1594.5 87.57 22.16
88.10 88.51 88.92 89.33
−1771.1 161.2 46.3 81.1
90.61 90.82 91.03 91.24
82.88 91.30 92.39 91.27
89.75 90.16 90.57 90.98
85.69 90.18 90.38 90.92
91.46 93.92 96.97 100.83
91.41 93.92 96.97 100.83
91.40 95.59 100.0 104.84
91.40 95.59 100.0 104.84
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
291
For the purpose of comparisons, the Newton–Raphson method [11], i.e. zk+1 = zk −
f (zk ) f (zk )
(9)
is used iteratively for solving Eq. (1) too. The results denoted as Vn are also listed in the tables. The following equation [11] is adopted to produce liquid molar volumes for representing experimental data which are also listed in the tables for the comparisons. (1+(1−(T /S2 )))S3
Ve =
S1
S4
(10)
where S1 , S2 , S3 and S4 are constants depended on substances. From the Tables 1–3, it is easy to notice that in the low temperature region, i.e. in the region of about Tr < 0.3, the analytical solutions of three cubic EOS give unreasonable results. It should be pointed out that all calculations carried out in this work are based on double-precision arithmetic. If single-precision arithmetic is used, the results given by the analytical solutions will be worse. The tables also show that the results produced by Newton–Raphson iterative method are reasonable in comparison with experimental data, which means it is arithmetic, not cubic EOS themselves that are responsible for the errors. 3. Discussion 3.1. The limitation of signification figure For any algorithmic language used in computers, the numbers of significant figure are finite. In this work, the C language is used. Based on double-precision arithmetic, all variables used in the program made by the C language have 16 significance digits. It is easy to understand that when more than 16 significance digits is needed by a operation during a calculation, the calculated result will certainly more or less be in error. As examples, the following discussion will be focused on the calculation of liquid volume of 1-pentene by PT EOS at the condition of T = 121.93 K, P = 5.24 × 10−10 MPa which is produced by Eq. (8). For PT EOS, c3 in Eq. (1) is given by c3 = B 3 C + BC − AB
(11)
where A=
Pa R2T 2
Pb RT Pc C= RT B=
where a, b and c are the constants of PT EOS which depend on substances.
(12) (13) (14)
292
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
At the condition given above, B 3 C = 3.24620 × 10−42
(15)
BC = 1.60586 × 10−21
(16)
AB = 8.11908 × 10−21
(17)
Obviously, accurate algebraic operation among B3 C, BC and AB needs more than 16 significant figure. However, due to the limitations of computers and their language, the contribution made by B3 C to c3 in Eq. (11) is in fact nothing. Then, let Eq. (3) be checked. Under the same condition defined above, three terms of it have the following values, respectively 1 3 c 27 1
= −3.70370 × 10−2
1 cc 6 1 2 1 c 2 3
= −2.87523 × 10−10
= −3.97925 × 10−20
(18) (19) (20)
Once again, because of the limitation of significant figure, the contribution of c3 /2 to r in Eq. (3) is ignored. 3.2. The magnification of error Let
r s= −q 3
(21)
and the error of s is s. So the Eq. (5), which representing liquid compressibility factor can be expressed as √ c1 arccos(s) ◦ z2 = 2 −q cos + 120 − (22) 3 3 The magnification of error is defined as M=
z/z2 s/s
(23)
where z is the error of z2 . When s tend to zero, M can be expressed as
√ dz s 2 −q arccos(s) 1 s ◦ M= = (24) sin + 120 √ ds z2 3 3 1 − s 2 z2 √ It is easy to see that when the value of s equals 1, the term 1/ 1 − s 2 will turn to infinity, and at the same time, the other parts of the right side of the above Eq. (24) still remain finite, so the error magnification M will become infinity too. Unfortunately, in low temperature region, the values of s do tend to 1. For example, in the calculation of liquid volume of 1-pentene by PT EOS at the condition of temperature and pressure given above, there
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
293
Table 4 The values of variables in Eqs. (3), (11) and (24) at different temperatures 107.93 K
121.93 K
135.93 K
149.93 K
254.93 K
Eq. (11)
B3 C BC AB
2.84 × 10−50 1.50 × 10−25 8.86 × 10−24
3.25 × 10−42 1.61 × 10−21 8.12 × 10−20
4.99 × 10−36 1.99 × 10−18 8.75 × 10−17
3.64 × 10−31 5.38 × 10−16 2.08 × 10−14
6.81 × 10−14 2.33 × 10−7 4.31 × 10−6
Eq. (3)
(1/6)c1 c2 (1/2)c3 (1/27)c13
−3.27 × 10−12 −4.36 × 10−24 −3.70 × 10−2
−2.88 × 10−10 −3.98 × 10−20 −3.70 × 10−2
−8.74 × 10−9 −4.27 × 10−17 −3.70 × 10−2
−1.25 × 10−7 −1.01 × 10−14 −3.70 × 10−2
−1.17 × 10−3 −2.04 × 10−6 −3.70 × 10−2
Eq. (24)
s √ 1/ 1 − s 2
1 ∝
1 ∝
1−8.0 × 10−16 7.8 × 106
1−1.6 × 10−13 5.5 × 105
1−1.1 × 10−4 6.6 × 101
√ is s = 1, which leads 1/ 1 − s 2 = ∞. Here s may not be exactly equal to 1, it could be in the region of (1 − 1 × 10−16 ) < s < (1 + 1 × 10−16 ). But the computers can only output s = 1, for that the exact expression of it may beyond the ability of the computer’s language. In the Table 4, the values of variables in Eqs. (3), (11) and (24) for which PT EOS and 1-pentene are concerned at different temperatures are listed. The comparisons among the values show that, as temperature arises, the errors caused by the limitation of significant figures will become smaller, and at same time, the error magnifications also reduce obviously. From the Table 4, it can be noticed that even in low temperature region, the errors caused by the limitation of√significant figures are so small that they can be ignored in the most cases, but the values of the term 1/ 1 − s 2 are large √ enough to lead the errors being magnified tremendously. However, Eq. (9) has nothing to do with 1/ 1 − s 2 , which is the reason why the Newton–Raphson method can give corrective calculation results. Finally, it should be pointed out that if Eq. (1) is transformed into a polynomial in molar volume V instead compressibility factor z, the results of calculation are almost the same as that given above. 4. Conclusions Due to the limitation of significant figures of the computer languages, and especially the effect of the magnification of errors in low temperature region, the analytical solutions of cubic EOS lead to irrational results. In contrast, the iterative solutions of cubic EOS by using the Newton–Raphson method produce valid results. It is safer to use the iterative method than the analytical method for solving cubic EOS. Generally speaking, at least in the region of Tr < 0.3, it is suggested to adopt the iterative method instead of the analytical method. List of symbols a A b B c
parameter in PT EOS variable in Eq. (11) parameter in PT EOS variable in Eq. (11) parameter in PT EOS
294
c1 , c2 , c3 C CCOR f M P PR PT q q1 , q2 , q3 , q4 , q5 r R s S1 , S2 , S3 , S4 T V z Superscripts k
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
coefficients of Eq. (1) variable in Eq. (11) cubic chain-of-rotator EOS polynomial form of cubic EOS magnification of error saturated vapor pressure Peng–Robinson EOS Patel–Teja EOS variable in Eq. (2) constants in Eq. (8) variable in Eq. (3) gas constant variable in Eq. (21) variables in Eq. (10) temperature liquid molar volume compressibility factor
◦
kth iteration degree
Subscripts e f n r
experimental analytical solution Newton–Raphson method reduced property
References [1] M. Benedict, G.B. Webb, L.G. Rubin, J. Chem. Phys. 8 (1940) 334–345. [2] J.W. Tester, M. Modell, Thermodynamics and its Applications, 3rd Edition, Prentice-Hall, New Jersey, 1996, p. 925. [3] J.M. Prausnitz, R.N. Lichtenthaler, E.G.D. Azevedo, Molecular Thermodynamics of Fluid-Phase Equilibria, Prentice-Hall, New Jersey, 1999, p. 50. [4] G. Soave, Fluid Phase Equilib. 56 (1990) 39–57. [5] P.M. Mathias, J.F. Boston, S. Watanasiri, AIChE J. 30 (1984) 182–186. [6] P.M. Mathias, M.S. Benson, AIChE J. 32 (1986) 2087–2090. [7] D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 10 (1976) 59–64. [8] N.C. Petal, A.S. Teja, Chem. Eng. Sci. 37 (1982) 464–473. [9] H. Kim, H-M. Lin, K.C. Cao, Ind. Eng. Chem. Fundam. 25 (1986) 75–84. [10] R.H. Perry, D.W. Green, J.O. Maloney, Perry’s Chemical Engineers’ Handbook, 6th Edition, McGraw-Hill, New York, 1984, pp. 2–15. [11] D. Zwillinger et al. Standard Mathematical Tables and Formulae, CRC Press, Boca Raton, 1996, p. 82. [12] T.E. Daubert, Danner, Physical and Thermodynamic Properties of Pure Chemicals, Hemisphere, New York, 1989.
Fallibility of analytic roots of cubic equations of state in low temperature region Yun Zhi a,∗ , Huen Lee b a
b
Department of Chemical Engineering, Separation Engineering Research Center, Nanjing University of Chemical Technology, Nanjing 210009, PR China Department of Chemical Engineering, Korea Advanced Institute of Science and Technology, 373-1 Kusong-dong, Yusong-gu, Taejon 305-701, South Korea Received 11 March 2000; accepted 24 January 2002
Abstract It is showed by several examples that, in low temperature region, the analytical solutions of cubic equation of state (EOS) lead to irrational results, while the iterative solutions of cubic EOS by using the Newton–Raphson method produce valid results. Errors caused by the limitation of significant figures of the computer languages are revealed, and a magnification of errors is defined which is a main factor bringing out the irrational results of the analytical solution of cubic EOS. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Cubic equation of state; Analytic roots; Fallibility; Low temperature
1. Introduction Since van der Walls presented his famous equation of state (EOS) about 100 year ago, various EOS have been developed. Among varieties of type of EOS, the number of cubic EOS may be the largest one. As cubic EOS has a notable advantage over complex EOS in its mathematical simplicity, it has attracted numerous researchers and played an important role in engineering calculations. The complex EOS, such as the celebrated BWR EOS [1] etc., can not be solved analytically for compressibility factors or specific volumes with temperature and pressure specified. The trial or iterative methods are always needed to solve the complex EOS for getting their roots one by one, and after these roots being obtained, to discriminate them is often needed too. On the other hand, for the cubic EOS, analytical solutions can be easily got by using a relevant math formula, i.e. Cardan’s method to solve EOS for their three roots simultaneously, as recommended by some textbooks on thermodynamics ∗
Corresponding author. E-mail address: [email protected] (Y. Zhi). 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 2 ) 0 0 0 7 2 - 9
288
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
[2,3]. Generally speaking, with this merit, the cubic EOS make the programming and calculations of thermodynamics concerned easier and simpler than the complex EOS do. In recent years, some non-cubic EOS have been successfully developed, but no direct solution for them like a cubic EOS was considered as a drawback [4]. However, for phase equilibria calculation involving infeasible conditions, Mathias et al. [5] developed a heuristic and iterative procedure to calculate density roots using EOS, no matter cubic or non-cubic, and the procedure is as efficient as analytical approaches. Furthermore, a detailed study made by Mathias and Benson [6] showed that, for most cases of phase equilibria calculations, the computational load of cubic equations are approximately equivalent to more complicated equations. In the case of a pure component, the cubic equations are more computationally rapid than the more complicated equations, but the difference becomes small for systems with many components. The larger the number of components, the smaller the difference, as most of the CPU time is spent in calculating mixing-rule summations. But, it seems that the valid of analytical roots of cubic EOS is beyond suspicion, and has not been explored in literature yet. Even in the case of a pure component, there is a problem with analytic roots of cubic EOS which was found in low temperature region by authors. The aim of this work is to demonstrate several examples of the fallibility and discuss the reasons which cause the fallibility.
2. Calculations and results Three well known cubic EOSs, PR [7], PT [8], and CCOR EOSs [9] are adopted as examples in this work. In the literature [2,10,11], there are different expressions of the math formula needed for the analytic solution of cubic equations, but the essences of which are in fact the same, i.e. Cardan’s method. In this work, the formula recommended by Tester and Modell [2] is used. As in the temperature and pressure ranges concerned in this work, the cubic EOS will certainly have three real roots, for lucidity, the cases of less than three real roots are not taken into account. Suppose a cubic EOS being arranged as a polynomial form f (z) = z3 + c1 z2 + c2 z + c3 = 0
(1)
where z is compressibility factor, c1 , c2 and c3 the coefficients dependent on temperature and pressure. According to the approach illustrated by Tester and Modell [2], let q = 19 (3c2 − c12 )
(2)
r = 16 c1 c2 − 21 c3 −
1 3 c 27 1
(3)
Then, three real roots are given by arccos r/ −q 3 √ − c1 z1 = 2 −q cos 3 3
√ z2 = 2 −q cos
arccos r/ −q 3 3
(4)
c1 + 120 − 3 ◦
(5)
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
√ z3 = 2 −q cos
arccos r/ −q 3 3
289
c1 ◦ + 240 − 3
(6)
where z2 is the smallest root representing the compressibility factor of liquid phase, z3 is the largest root representing the compressibility factor of vapor phase, z1 has no physical significance. It is not the accuracy of the model we are concerned, but the validity of the compressibility factor calculation. For the facility of discussion, only liquid compressibility factors and liquid volumes are concerned in the following calculation. Liquid molar volume is given by V = z2
RT P
(7)
where R is the gas constant, T the temperature, and P the saturated vapor pressure which, in this work, is given by [12] q2 (8) + q3 ln T + q4 T q5 P = exp q1 + T where q1 , q2 , q3 , q4 and q5 are constants depended on substances. The liquids molar volumes of three pure substances (propylene, 1-butene and 1-pentene) are calculated by the above procedure, and the results are listed in the Tables 1–3, denoted as Vf .
Table 1 Comparison among experimental and calculated liquid volumes for propylene Number
T (K)
Tr
P (MPa)
Ve
PT
PR
CCOR
Vn
Vf
Vn
Vf
Vn
Vf
−2197.4 1008.7 −132.6 135.3
53.91 53.97 54.04 54.10
1002.6 981.03 −159.66 108.88
55.50 55.64 55.78 55.92
521.22 88.626 27.537 63.887
1 2 3 4
87.9 89.4 90.9 92.4
0.241 0.245 0.249 0.253
9.18E−10 1.60E−09 2.74E−09 4.59E−09
55.07 55.18 55.29 55.39
57.94 58.02 58.09 58.16
5 6 7 8
93.9 95.4 96.9 98.4
0.257 0.262 0.266 0.270
7.56E−09 1.22E−08 1.95E−08 3.06E−08
55.50 55.61 55.72 55.84
58.24 58.31 58.39 58.47
109.19 55.98 62.74 57.49
54.17 54.23 54.30 54.37
83.30 71.18 54.68 58.28
56.06 56.20 56.34 56.48
59.85 51.98 56.71 55.90
9 10 11 12
99.9 101.4 102.9 127.9
0.274 0.278 0.282 0.351
4.73E−08 7.21E−08 1.08E−07 2.08E−05
55.95 56.06 56.17 58.17
58.54 58.62 58.70 60.16
60.66 58.30 59.35 60.16
54.44 54.51 54.58 55.87
52.62 54.25 54.43 55.87
56.62 56.75 56.89 59.19
56.43 56.82 56.89 59.19
13 14 15 16
152.9 177.9 202.9 227.9
0.419 0.488 0.556 0.625
6.03E−04 6.09E−03 3.24E−02 1.14E−01
60.39 62.91 65.79 69.16
61.90 64.01 66.62 69.91
61.90 64.01 66.62 69.91
57.42 59.30 61.63 64.58
57.42 59.30 61.63 64.58
61.53 63.98 66.67 69.73
61.53 63.98 66.67 69.73
290
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
Table 2 Comparison among experimental and calculated liquid volumes for 1-butene Number
T (K)
Tr
P (MPa)
Ve
PT
PR
Vn
Vf
Vn
1628.5 1550.2 −370983 1411.1
70.37 70.50 70.63 70.77
CCOR Vf
Vn
Vf
1671.6 1601.8 1536.6 1475.6
69.63 69.95 70.27 70.59
−8310460 −2053680 6370.11 5940.21
1418.4 1364.8 −5679.9 1266.9
70.91 71.22 71.54 71.86
1 2 3 4
87.8 91.3 94.8 98.3
0.209 0.218 0.226 0.234
3.56E−13 1.87E−12 8.58E−12 3.51E−11
70.09 70.35 70.62 70.89
75.79 75.96 76.15 76.34
5 6 7 8
101.8 105.3 108.8 112.3
0.243 0.251 0.259 0.268
1.29E−10 4.32E−10 1.33E−09 3.79E−09
71.16 71.44 71.71 72.00
76.53 76.72 76.92 77.12
1349.1 1291.4 1237.7 −36.61
70.90 71.04 71.19 71.33
9 10 11 12
115.8 119.3 122.8 126.3
0.276 0.284 0.293 0.301
1.01E−08 2.51E−08 5.92E−08 1.32E−07
72.28 72.57 72.86 73.16
77.33 77.55 77.76 77.98
78.65 77.75 75.89 78.12
71.48 71.63 71.78 71.94
71.69 51.00 50.81 76.89
72.17 72.49 72.80 73.12
−91.96 92.06 71.29 72.85
13 14 15 16
161.3 196.3 231.3 266.3
0.384 0.468 0.551 0.635
4.98E−05 1.80E−03 1.90E−02 9.86E−02
76.34 80.00 84.31 89.52
80.49 83.63 87.66 93.01
80.49 83.63 87.66 93.01
72.10 74.20 76.86 80.28
71.94 74.20 76.86 80.28
73.43 77.11 80.94 85.15
73.07 77.11 80.94 85.15
5548.8 5191.4 −2129.9 202.2
Table 3 Comparison among experimental and calculated liquid volumes for 1-pentene Number
T (K)
Tr
P (MPa)
Ve
PT Vn
PR Vf
Vn
1 2 3 4
107.93 111.43 114.93 118.43
0.232 0.240 0.247 0.255
4.48E−12 1.67E−11 5.70E−11 1.79E−10
86.64 86.94 87.25 87.56
91.05 91.23 91.42 91.61
1961.51 1881.09 1805.88 1735.39
89.11 89.28 89.46 89.64
5 6 7 8
121.93 125.43 128.93 132.43
0.262 0.270 0.277 0.285
5.24E−10 1.43E−09 3.68E−09 8.92E−09
87.87 88.19 88.51 88.83
91.80 92.00 92.20 92.40
1669.2 1607.0 99.54 122.2
89.83 90.02 90.21 90.41
9 10 11 12
135.93 139.43 142.93 146.43
0.292 0.300 0.308 0.315
2.06E−08 4.52E−08 9.50E−08 1.92E−07
89.16 89.49 89.83 90.17
92.61 92.82 93.04 93.26
84.82 91.59 92.45 93.18
13 14 15 16
149.93 184.93 219.93 254.93
0.323 0.398 0.473 0.548
3.73E−07 5.93E−05 1.48E−03 1.32E−02
90.51 94.18 98.36 103.21
93.49 96.04 99.20 103.20
93.45 96.04 99.20 103.20
CCOR Vf
Vn
1946.1 1866.4 1791.8 −20549
Vf
86.44 86.86 87.27 87.69
8026.3 7468.2 6961.1 6499.1
1656.3 1594.5 87.57 22.16
88.10 88.51 88.92 89.33
−1771.1 161.2 46.3 81.1
90.61 90.82 91.03 91.24
82.88 91.30 92.39 91.27
89.75 90.16 90.57 90.98
85.69 90.18 90.38 90.92
91.46 93.92 96.97 100.83
91.41 93.92 96.97 100.83
91.40 95.59 100.0 104.84
91.40 95.59 100.0 104.84
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
291
For the purpose of comparisons, the Newton–Raphson method [11], i.e. zk+1 = zk −
f (zk ) f (zk )
(9)
is used iteratively for solving Eq. (1) too. The results denoted as Vn are also listed in the tables. The following equation [11] is adopted to produce liquid molar volumes for representing experimental data which are also listed in the tables for the comparisons. (1+(1−(T /S2 )))S3
Ve =
S1
S4
(10)
where S1 , S2 , S3 and S4 are constants depended on substances. From the Tables 1–3, it is easy to notice that in the low temperature region, i.e. in the region of about Tr < 0.3, the analytical solutions of three cubic EOS give unreasonable results. It should be pointed out that all calculations carried out in this work are based on double-precision arithmetic. If single-precision arithmetic is used, the results given by the analytical solutions will be worse. The tables also show that the results produced by Newton–Raphson iterative method are reasonable in comparison with experimental data, which means it is arithmetic, not cubic EOS themselves that are responsible for the errors. 3. Discussion 3.1. The limitation of signification figure For any algorithmic language used in computers, the numbers of significant figure are finite. In this work, the C language is used. Based on double-precision arithmetic, all variables used in the program made by the C language have 16 significance digits. It is easy to understand that when more than 16 significance digits is needed by a operation during a calculation, the calculated result will certainly more or less be in error. As examples, the following discussion will be focused on the calculation of liquid volume of 1-pentene by PT EOS at the condition of T = 121.93 K, P = 5.24 × 10−10 MPa which is produced by Eq. (8). For PT EOS, c3 in Eq. (1) is given by c3 = B 3 C + BC − AB
(11)
where A=
Pa R2T 2
Pb RT Pc C= RT B=
where a, b and c are the constants of PT EOS which depend on substances.
(12) (13) (14)
292
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
At the condition given above, B 3 C = 3.24620 × 10−42
(15)
BC = 1.60586 × 10−21
(16)
AB = 8.11908 × 10−21
(17)
Obviously, accurate algebraic operation among B3 C, BC and AB needs more than 16 significant figure. However, due to the limitations of computers and their language, the contribution made by B3 C to c3 in Eq. (11) is in fact nothing. Then, let Eq. (3) be checked. Under the same condition defined above, three terms of it have the following values, respectively 1 3 c 27 1
= −3.70370 × 10−2
1 cc 6 1 2 1 c 2 3
= −2.87523 × 10−10
= −3.97925 × 10−20
(18) (19) (20)
Once again, because of the limitation of significant figure, the contribution of c3 /2 to r in Eq. (3) is ignored. 3.2. The magnification of error Let
r s= −q 3
(21)
and the error of s is s. So the Eq. (5), which representing liquid compressibility factor can be expressed as √ c1 arccos(s) ◦ z2 = 2 −q cos + 120 − (22) 3 3 The magnification of error is defined as M=
z/z2 s/s
(23)
where z is the error of z2 . When s tend to zero, M can be expressed as
√ dz s 2 −q arccos(s) 1 s ◦ M= = (24) sin + 120 √ ds z2 3 3 1 − s 2 z2 √ It is easy to see that when the value of s equals 1, the term 1/ 1 − s 2 will turn to infinity, and at the same time, the other parts of the right side of the above Eq. (24) still remain finite, so the error magnification M will become infinity too. Unfortunately, in low temperature region, the values of s do tend to 1. For example, in the calculation of liquid volume of 1-pentene by PT EOS at the condition of temperature and pressure given above, there
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
293
Table 4 The values of variables in Eqs. (3), (11) and (24) at different temperatures 107.93 K
121.93 K
135.93 K
149.93 K
254.93 K
Eq. (11)
B3 C BC AB
2.84 × 10−50 1.50 × 10−25 8.86 × 10−24
3.25 × 10−42 1.61 × 10−21 8.12 × 10−20
4.99 × 10−36 1.99 × 10−18 8.75 × 10−17
3.64 × 10−31 5.38 × 10−16 2.08 × 10−14
6.81 × 10−14 2.33 × 10−7 4.31 × 10−6
Eq. (3)
(1/6)c1 c2 (1/2)c3 (1/27)c13
−3.27 × 10−12 −4.36 × 10−24 −3.70 × 10−2
−2.88 × 10−10 −3.98 × 10−20 −3.70 × 10−2
−8.74 × 10−9 −4.27 × 10−17 −3.70 × 10−2
−1.25 × 10−7 −1.01 × 10−14 −3.70 × 10−2
−1.17 × 10−3 −2.04 × 10−6 −3.70 × 10−2
Eq. (24)
s √ 1/ 1 − s 2
1 ∝
1 ∝
1−8.0 × 10−16 7.8 × 106
1−1.6 × 10−13 5.5 × 105
1−1.1 × 10−4 6.6 × 101
√ is s = 1, which leads 1/ 1 − s 2 = ∞. Here s may not be exactly equal to 1, it could be in the region of (1 − 1 × 10−16 ) < s < (1 + 1 × 10−16 ). But the computers can only output s = 1, for that the exact expression of it may beyond the ability of the computer’s language. In the Table 4, the values of variables in Eqs. (3), (11) and (24) for which PT EOS and 1-pentene are concerned at different temperatures are listed. The comparisons among the values show that, as temperature arises, the errors caused by the limitation of significant figures will become smaller, and at same time, the error magnifications also reduce obviously. From the Table 4, it can be noticed that even in low temperature region, the errors caused by the limitation of√significant figures are so small that they can be ignored in the most cases, but the values of the term 1/ 1 − s 2 are large √ enough to lead the errors being magnified tremendously. However, Eq. (9) has nothing to do with 1/ 1 − s 2 , which is the reason why the Newton–Raphson method can give corrective calculation results. Finally, it should be pointed out that if Eq. (1) is transformed into a polynomial in molar volume V instead compressibility factor z, the results of calculation are almost the same as that given above. 4. Conclusions Due to the limitation of significant figures of the computer languages, and especially the effect of the magnification of errors in low temperature region, the analytical solutions of cubic EOS lead to irrational results. In contrast, the iterative solutions of cubic EOS by using the Newton–Raphson method produce valid results. It is safer to use the iterative method than the analytical method for solving cubic EOS. Generally speaking, at least in the region of Tr < 0.3, it is suggested to adopt the iterative method instead of the analytical method. List of symbols a A b B c
parameter in PT EOS variable in Eq. (11) parameter in PT EOS variable in Eq. (11) parameter in PT EOS
294
c1 , c2 , c3 C CCOR f M P PR PT q q1 , q2 , q3 , q4 , q5 r R s S1 , S2 , S3 , S4 T V z Superscripts k
Y. Zhi, H. Lee / Fluid Phase Equilibria 201 (2002) 287–294
coefficients of Eq. (1) variable in Eq. (11) cubic chain-of-rotator EOS polynomial form of cubic EOS magnification of error saturated vapor pressure Peng–Robinson EOS Patel–Teja EOS variable in Eq. (2) constants in Eq. (8) variable in Eq. (3) gas constant variable in Eq. (21) variables in Eq. (10) temperature liquid molar volume compressibility factor
◦
kth iteration degree
Subscripts e f n r
experimental analytical solution Newton–Raphson method reduced property
References [1] M. Benedict, G.B. Webb, L.G. Rubin, J. Chem. Phys. 8 (1940) 334–345. [2] J.W. Tester, M. Modell, Thermodynamics and its Applications, 3rd Edition, Prentice-Hall, New Jersey, 1996, p. 925. [3] J.M. Prausnitz, R.N. Lichtenthaler, E.G.D. Azevedo, Molecular Thermodynamics of Fluid-Phase Equilibria, Prentice-Hall, New Jersey, 1999, p. 50. [4] G. Soave, Fluid Phase Equilib. 56 (1990) 39–57. [5] P.M. Mathias, J.F. Boston, S. Watanasiri, AIChE J. 30 (1984) 182–186. [6] P.M. Mathias, M.S. Benson, AIChE J. 32 (1986) 2087–2090. [7] D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 10 (1976) 59–64. [8] N.C. Petal, A.S. Teja, Chem. Eng. Sci. 37 (1982) 464–473. [9] H. Kim, H-M. Lin, K.C. Cao, Ind. Eng. Chem. Fundam. 25 (1986) 75–84. [10] R.H. Perry, D.W. Green, J.O. Maloney, Perry’s Chemical Engineers’ Handbook, 6th Edition, McGraw-Hill, New York, 1984, pp. 2–15. [11] D. Zwillinger et al. Standard Mathematical Tables and Formulae, CRC Press, Boca Raton, 1996, p. 82. [12] T.E. Daubert, Danner, Physical and Thermodynamic Properties of Pure Chemicals, Hemisphere, New York, 1989.