- Email: [email protected]
An intensified analytic solution for finding the roots of a cubic equation of state in low temperature region
Journal of Molecular Liquids 206 (2015) 165–169
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
An intensified analytic solution for finding the roots of a cubic equation of state in low temperature region Behrooz M. Ziapour Department of Mechanical Engineering, University of Mohaghegh Ardabili, Ardabil, Iran
a r t i c l e
i n f o
Article history: Received 19 December 2014 Received in revised form 3 February 2015 Accepted 16 February 2015 Available online 19 February 2015 Keywords: Analytic Cubic Equation of state Low temperature Organic
a b s t r a c t Because of the effect of the enlarging process of errors in low temperature region, the analytical solutions of a cubic equation of state (CES) lead to illogical results. To overcome this problem, we recommend an enhanced analytic method, in this paper. Up to now, there is not any analytical method for this problem. In the present study, the general cubic equation x3 + Bx2 + Cx + D = 0 is transformed into simple equation γ3 − βγ + β = 0. For β ≥ 27/4, then the three roots are real. To compare, as an example, for 1-butene at T = 112.3 (K) and P = 3.79 × 10−9(Pa), the empirical value of the molar volume is 72.00 (m3mol−1). By using Patel–Teja (PT) CES, the molar volume value is obtained as −342.3 from traditional analytic method and 73.88 from both the new intensified analytic method and the iterative Newton–Raphson method. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The CES models are useful formulas for analyzing and correlating the thermodynamic equilibrium state of simple fluids and their mixtures, with agreeable approach [1–5]. Generally, for a given substance, the CES is an empirical relation of the pressure P, the temperature T and the molar volume v consist of a repulsive term and an attractive term as follows [6]: P ðT; vÞ ¼
RT α ðT Þ − v−b ðv−cÞðv−dÞ
ð1Þ
where R is the gas constant. b, c and d are constant parameters which are obtained based on the kind of equation. In Eq. (1), α(T) is a function and it has many empirical definitions in the literature. By expressing Eq. (1) in terms of molar volume, the following equation is obtained: RT 2 RT RT α ðT Þ 3 v þ bc þ bd þ cd þ cþ dþ v v − bþcþdþ P P P P RT α ðt Þ − bcd þ cd þ b ¼ 0: P P
ð2Þ
Eq. (2) is referred to as cubic in molar volume. Additionally we can substitute the definition of compressibility factor (i.e., Z = Pv/RT) into Eq. (2) and obtain different cubic equations in Z. The first CES was recommended by van der Waals in 1873 as P(T, v) = RT/(v − b) − α(T)/v2. It has been proven that the original van der Waals E-mail address: [email protected].
http://dx.doi.org/10.1016/j.molliq.2015.02.026 0167-7322/© 2015 Elsevier B.V. All rights reserved.
CES cannot provide simultaneous accurate predictions for all the properties of pure fluids and their mixtures. But, the most CES models (e.g., Peng-Robison CES in 1976 as P(T, v) = RT/(v − b) − α(T)/(v2 + 2bv − b2) being used widely today for practical design purposes have been derived from van der Waals CES [7]. For finding the roots of Eq. (2), the analytical methods can solve their three roots simultaneously. One of the relevant math formulas is the well-known formula due to Ferro and Tartaglia, communicated by Gerolamo Cardano in 1545 [8,9]. Cardano's method involves division. Thus it can encounter near “0/0”, resulting insignificant numerical errors. Thus it may encounter a numerically unstable case. Lagrange in 1770 gave the new method for solving the three roots of a cubic equation, with 18 possible interpretations [10]. The Lagrange formula, does not require division and thus avoids the “0/0” case. In the Lagrange formula some interpretations are correct (yielding the solutions), But the others are not. Zhao et al. [11] proposed a convection which provided correct interpretations of the Lagrange formula for all cubic equations with real coefficients. Zhi and Lee [12] revealed by several examples that, in low temperature region, the analytical solutions of CES lead to irrational results, while the iterative solutions of the CES, using Newton–Raphson method produced valid results. They appeared that errors caused by the limitation of significant figures of the computer languages are revealed, and a magnification of errors is defined which a main factor is bringing out the irrational results of the analytical solution of CES. Salim [13] showed that, in the case of low temperature region, the coefficient D in cubic equation x3 + Bx2 + Cx + D = 0 is very small. Therefore the main root is also small and becomes so close to zero. Thus, by eliminating the term of x3 in x3 + Bx2 + Cx + D = 0, then the main root can be
166
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
found via the approximated equation of the form Bx2 + Cx + D = 0. Loperena [14] proposed an iterative refinement of the solution obtained by the analytical method. This method allowed one to take advantage of the calculations fulfilled in the application of the analytical solution. His numerical results show that the proposed iterative refinement procedure is an attractive and numerical inexpensive option and can be easily incorporated into any process simulation program, that use Cardano's method, without significant modifications. Up to now, as mentioned above discussions, there is not any completely analytical method for this problem. To notice to complexity and overhead of analytical methods, it has been recommended to use the numeric iterative methods. In this paper, we describe an enhanced completely analytic method for solving the general cubic equation x3 + Bx2 + Cx + D = 0 consists of real coefficients. We verify this new method with the miscellaneous examples.
formula as: 8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 3 > 2 3 > β β β β β2 β3 A 27 > > @ A @ > − − γ ¼ − þ − ; for β ≤ þ − < 4 27 4 27 2 2 4 ! rffiffiffiffi > > > β 1 9 27 −1 > > pffiffiffiffiffiffiffiffiffi ; for β ≥ cosθ; θ ¼ cos : : γ ¼ −2 3 3 4 12β ð5Þ
2. New analytic solution and the miscellaneous examples
It is clear that (see Fig. 1), when β → − ∞, then γ → 1. Also, when β → + ∞, then γ → − ∞. To obtain the two formulas in Eq. (5), an analytic method was completely used. This method has been described in Appendix A. Therefore, through an intensified analytic solution, the real root of the cubic Eq. (3) is captured. The other roots of Eq. (3) i.e., x2, x3 can be calculated analytically from the following quadratic equation as:
An arbitrary cubic equation which consists of real coefficients B, C and D is as follows:
x þ ðB þ x1 Þx−
3
2
x þ Bx þ Cx þ D ¼ 0:
ð3Þ C3
In Eq. (3), one of the three roots is real. To find it, let β ¼ − D12 , 1
where C1 = 9C − 3B2 and D1 = 2B3 − 9BC + 27D. Make the substi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tution x ¼ − 13 3 ð1−γÞD1 þ B or x ¼ − 13 DC 11 γ þ B . By this substitution, Eq. (3) is transformed as the following simple equation for γ: 3
γ −βγ þ β ¼ 0:
ð4Þ
γ To find γ, we arrange Eq. (4) into the functional form as βðγ Þ ¼ γ−1 . 3
Fig. 1 shows the sketch of β(γ). It is seen that for β ≥ 27 4 , then all three roots are real. Therefore, Eq. (3) has three real roots (i.e., x1, x2, x3). Also, for β b 27 4 , Eq. (3) has one real root (i.e., x1). In this case, the others are conjugate imaginary roots (i.e., x2, x3). Therefore, like a new analytic method, for finding the one real root of pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a cubic equation, let x1 ¼ − 13 3 ð1−γ ÞD1 þ B or x1 ¼ − 13 DC 11 γ þ B be one real root of Eq. (3), where γ is related with the left branch of the curves β in Fig. 1 (i.e., γ ≤ 1); and it is obtained from the following
2
D ¼ 0: x1
ð6Þ
Note that, in order to reach the high precision of the solutions, especially in low temperature region, the two other roots may be obtained from Eq. (6). In this paper, all solutions were conducted by Microsoft Excel, 2003. Example 1. In equation: x3 − 3x2 + 4x − 2 = 0, the coefficients are: B = − 3, C = 4 and D = − 2. Then, C1 = 9 and D1 = 0. Therefore, β = − ∞. Thus, from Fig. 1: γ = 1. Therefore, x1 = − B/3 = 1. Finally from Eq. (6): x2 − 2x + 2 = 0. Then, x2 = 1 + i and x3 = 1 − i. Example 2. In equation: x3 + 6x2 + 3x − 10 = 0, the coefficients are: B = 6, C = 3, D = −10. Then C1 = −81 and D1 = 0. Therefore, β = + ∞. Thus, from Fig. 1: γ = − ∞. Therefore, x1 = − B/3 = −2. Finally from Eq. (6): x2 + 4x − 5 = 0. Then, x2 = 1 and x3 = −5. Example 3. In equation: x3 − 2x2 − 5x + 6 = 0, the coefficients are: B = − 2, C = − 5, D = 6. Then C1 = − 57 and D1 = 56. Then, β = 59.0539. Thus, from Eq. (5):γ = −8.1429. Therefore, x1 = −2.0000. Finally from Eq. (6):x2 − 4x + 3 = 0. Then, x2 = 1.0000 and x3 = 3.0000. Example 4. The cubic equation: (x − 1)2x = x3 − 2x2 + x = 0. The coefficients are: B = −2, C = 1, D = 0. Then C1 = −3 and D1 = 2. Then, β = 6.75. Thus, from Eq. (5):γ = −3. Therefore, x1 = 0. Finally from solution of the equation:x2 − 2x + 1 = 0, we have:x2 = x3 = 1.0000. Example 5. The molar volume of carbon dioxide (CO 2 ) at P = 10 atm and T = 300 K, by using the van der Waals CES (P(T, v) = RT/(v − b) − α(T)/v 2 ), is the cubic equation as: v 3 − 2.5046v 2 + 0.3598v − 0.0155 = 0. This equation and another equation in the form of v 3 − 7.8693v 2 + 13.3771v − 6.5354 = 0
Table 1 Results for example 5. Case I: v3 − 2.5046v2 + 0.3598v − 0.0155 = 0 Newton (numeric)
Ostrowski [15] (numeric)
New (analytic) Fig. 1. The curve of β(γ).
2:4616 2:3633 Iterations→ 2:3546 2:3545→v1 2:4616 Iterations→ 2:3545→v1 β ¼ 6:7165 γ ¼ −2:9934 v1 ¼ 2:3545
Case II: v3 − 7.8693v2 + 13.3771v − 6.5354 = 0 7:3807 6:2992 Iterations→ 5:8355 5:7397 5:7357→v1 7:3807 Iterations→ 5:8187 5:7357→v1 β ¼ 6:7396 γ ¼ −2:9977 v1 ¼ 5:7357
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
were solved by the some numerical iteration methods [15]. The initial guess for starting the iteration was selected from the ideal gas law. We collected both these numerical solutions and our new analytical calculations, in Table 1. As shown from the new solutions, the vales of β are less than 6.75. Thus, the selected cubic equations have only one real root (i.e., v1 ≡ x1). The other imaginary roots are unusable. 3. Results and discussions For substances, it is clear that, in the low temperature and pressure ranges, their CES will have three real roots (i.e., β ≥ 6.75). In order to compare our new method with the traditional methods, taken from Zhi and Lee [12], we assume that CES be formed versus the compressibility factor z as: z3 + c1z2 + c2z + c3 = 0. The smallest root (i.e., the compressibility factor of liquid phase:z2) is obtained from the traditional method as follows: pffiffiffiffiffiffiffiffiffiffi 0 1 arccos r= −q3 pffiffiffiffiffiffiffiffi c ∘ þ 120 A− 1 z2 ¼ 2 −q cos@ 3 3
167
Table 3 Comparison among experimental, old and new analytical methods calculated liquid volume for 1-butene (TC = 419.5 K, PC = 4020 kPa, ω = 0.194). P (KPa) P (MPa)
87.8 91.3 94.8 98.3 101.8 105.3 108.8 112.3 115.8 119.3 122.8 126.3 161.3 196.3 231.3 266.3
3.56E−13 1.87E−12 8.58E−12 3.51E−11 1.29E−10 4.32E−10 1.33E−09 3.79E−09 1.01E−08 2.51E−08 5.92E−08 1.32E−07 4.98E−05 1.80E−03 1.90E−02 9.86E−02
VExp
70.09 70.35 70.62 70.89 71.16 71.44 71.71 72.00 72.28 72.57 72.86 73.16 76.34 80.00 84.31 89.52
PR
PT
VNew
VOld
VNewton VNew
VOld
VNewton
70.32 70.47 70.62 70.78 70.94 71.11 71.28 71.45 71.62 71.80 71.99 72.17 74.30 76.98 80.45 85.10
1685.9 1603.6 – −92,975 −25,216 1332.8 – 0.460 116.28 78.07 72.18 72.91 74.30 76.98 80.45 85.10
70.32 70.47 70.62 70.78 70.94 71.11 71.28 71.45 71.62 71.80 71.99 72.17 74.30 76.98 80.45 85.10
– 1621.3 1544.6 – 1408.3 1347.3 1290.5 −342.3 129.1 72.17 73.20 74.45 76.89 79.73 83.39 88.27
72.68 72.84 73.01 73.17 73.34 73.52 73.69 73.88 74.06 74.25 74.45 74.64 76.89 79.73 83.39 88.27
72.68 72.84 73.01 73.17 73.34 73.52 73.69 73.88 74.06 74.25 74.45 74.64 76.89 79.73 83.39 88.27
ð7Þ
1 3 c1 . Then the molar volwhere, q ¼ 19 3c2 −c21 and r ¼ 16 c1 c2 − 12 c3 − 27 3 −1 ume of the liquid phase is given as: V Old ¼ z2 RT ). Suppose P (m mol that, VExp, VNew, VOld and VNewton be denoted to the molar volume of the liquid phase obtained from the experimental, the new method, the old method and the iterative Newton–Rephson method, respectively. These values for three pure substances (propylene, 1-butene and 1pentene), have been listed in Tables 2–4, by using both the PR and the PT CES models. In these tables, to find the VNewton, we consider Newton 3 2 2 scheme as: zkþ1 ¼ zk − zk þ c1 zk þ c2 zk þ c3 = 3zk þ 2c1 zk þ c2 .
The recorded results for VNewton in Tables 2–4, were obtained due to the three iterations as k = 3. In this iterative process, the initial guess for zk was as z1 = 0. As shown in Tables 2–4, the results indicate that both the present intensified analytic method and the iterative Newton–Rephson method have identical accuracies. But, because of the effect of the enlarging process of errors, especially in lower temperature region, the traditional analytic solutions of CES lead to illogical results. The magnification of error was analyzed for traditional method (i.e., Eq. (7)) by Zhi and Lee [12]. Through this analysis, it was found that the pffiffiffiffiffiffiffiffiffiffiffiffi critical term in the magnification of error is 1= 1−s2 , where, the value pffiffiffiffiffiffiffiffiffiffi of s is as: s ¼ r= −q3. They showed that in low temperature region, the values of s do tend to 1, so the error magnification (M) will become
Table 2 Comparison among experimental, old and new analytical methods calculated liquid volume for propylene (TC = 364.8 K, PC = 4613 kPa, ω = 0.142). P (KPa)
P (MPa)
VExp
87.9 89.4 90.9 92.4 93.9 95.4 96.9 98.4 99.9 101.4 102.9 127.9 152.9 177.9 202.9 227.9
9.18E−10 1.60E−09 2.74E−09 4.59E−09 7.56E−09 1.22E−08 1.95E−08 3.06E−08 4.73E−08 7.21E−08 1.08E−07 2.08E−05 6.03E−04 6.09E−03 3.24E−02 1.14E−01
55.07 55.18 55.29 55.39 55.50 55.61 55.72 55.84 55.95 56.06 56.17 58.17 60.39 62.91 65.79 69.16
PR
PT
VNew
VOld
VNewton VNew
VOld
VNewton
53.91 53.97 54.04 54.10 54.17 54.24 54.30 54.37 54.44 54.51 54.58 55.87 57.42 59.31 61.64 64.58
1002.3 980.71 −158.8 108.41 82.69 67.72 53.97 61.62 55.69 54.45 54.92 55.87 57.42 59.31 61.64 64.58
53.91 53.97 54.04 54.10 54.17 54.24 54.30 54.37 54.44 54.51 54.58 55.87 57.42 59.31 61.64 64.58
– 995.94 −395.28 275.55 96.98 81.71 42.68 56.99 60.93 58.77 58.24 58.76 60.44 62.48 65.00 68.18
56.62 56.69 56.76 56.83 56.91 56.98 57.05 57.13 57.20 57.28 57.35 58.76 60.44 62.48 65.00 68.18
56.62 56.69 56.76 56.83 56.91 56.98 57.05 57.13 57.20 57.28 57.35 58.76 60.44 62.48 65.00 68.18
infinity. The magnification of error analysis for present intensified analytic method was done as the following process: pffiffiffiffiffiffiffiffiffi Let s ¼ 9= 12β, so the second formula of Eq. (5) can be arranged as: γ¼
−3 1 −1 cos cos ðsÞ : s 3
ð8Þ
When Δs tend to zero, the magnification of error (M) can expressed as [12]: M¼
dγ s ¼ ds γ
3 1 1 1 s −1 −1 p ffiffiffiffiffiffiffiffiffiffiffi ffi cos ð s Þ þ ð s Þ cos cos : ð9Þ sin 2 3 3 γ s2 s 1−s
In low temperature region, the values of β, γ and s do tend to 27/4, 3 and 1 respectively. When s = 1, then the second term in Eq. (9) gets p1ffiffiffiffiffiffiffi sin 1 cos−1 ð1Þ ¼ 0. Let 3 0
1 1−1
dγ s ds γ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 3 1 1 1 s −1 −1 cos ðsÞ þ pffiffiffiffiffiffiffiffiffiffiffiffi 1− cos2 ð cos ðsÞÞ : ð10Þ ¼ 2 cos 3 3 γ s s 1−s2
M¼
Or, by inserting Eq. (8), we can write as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 0 dγ s @ 3 1 1 s2 γ 2 A s −1 cos ¼ cos ðsÞ þ pffiffiffiffiffiffiffiffiffiffiffiffi 1− : M¼ 2 9 ds γ 3 γ s2 s 1−s
ð11Þ
As mentioned above, in low temperature region, we have as: γ → − 3. Finally, by inserting this value into Eq. (11), then the value of M is obtained as follows: M¼
dγ s ¼ ds γ
3 1 1 s −1 cos ð s Þ þ cos : 3 s γ s2
ð12Þ
Therefore, as seen from value of M in Eq. (12), in low temperature region, there is not any infinity term for value of M, which is the reason why this present method can give corrective calculation results. For a better presentation, we have listed the new solution procedure in Tables 5 and 6. Tables 5 and 6 have been listed based on the PR model for propylene and 1-butane, respectively. In the following example, more details of the new method calculations have been introduced.
168
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
Table 4 Comparison among experimental, old and new analytical methods calculated liquid volume for 1-pentene (TC = 464.4 K, PC = 3555 kPa, ω = 0.233). P (KPa)
107.93 111.43 114.93 118.43 121.93 125.43 128.93 123.43 135.93 139.43 142.93 146.43 149.93 184.93 219.93 254.93
P (MPa)
4.48E−12 1.67E−11 5.70E−11 1.79E−10 5.24E−10 1.43E−09 3.68E−09 8.92E−09 2.06E−08 4.52E−08 9.50E−08 1.92E−07 3.73E−07 5.93E−05 1.48E−03 1.32E−02
VExp
86.64 86.94 87.25 87.56 87.87 88.19 88.51 88.83 89.16 89.49 89.83 90.17 90.51 94.18 98.36 103.21
PR
PT
VNew
VOld
VNewton
VNew
VOld
VNewton
88.38 88.56 88.74 88.92 89.10 89.29 89.48 89.18 89.88 90.08 90.29 90.50 90.72 93.16 96.19 100.03
−810,425 1849.6 – 1706.4 −6205 1580.1 341.14 139.16 91.06 94.20 90.84 90.39 90.78 93.16 96.19 100.03
88.38 88.56 88.74 88.92 89.10 89.29 89.48 89.18 89.88 90.08 90.29 90.50 90.72 93.16 96.19 100.03
90.31 90.49 90.67 90.86 91.05 91.25 91.45 91.14 91.86 92.07 92.29 92.51 92.73 95.26 98.40 102.37
1944 1864.2 1789.6 1719.7 1654.1 1592.5 87.49 79.74 83.53 96.38 94.34 92.43 92.78 95.26 98.40 102.37
90.31 90.49 90.67 90.86 91.05 91.25 91.45 91.14 91.86 92.07 92.29 92.51 92.73 95.26 98.40 102.37
Example 6. From Table 5, for propylene using PR model, at T = 227.9(K) and P = 114(kPa), then the following CES is obtained:
from the following equation as: 3
z1 ¼ 3
2c1 −9c1 c2 þ 27c3 c1 − : 3 9c2 −3c1 2
ð13Þ
2
z −0:9969224047840840z þ 0:0303480538062392z−0:0001029289169739 ¼ 0:
Please note that, Eq. (13) has been adjusted based on the cubic equa D1 C 1 γ þ B ),
tion as: x3 + Bx2 + Cx + D = 0 and its one real root (x1 ¼ − 13 Therefore, C1 = − 2.7084303592252900 and D1 = − 1.7120783876879200. Thus: β = 6.7780650055531100. Then, θ = 0.0214638985371484. Hence from Eq. (5): = z1 = γ = − 3.0055377417229900. Then, x1 0.9656036987339030. From Eq. (6): x2 = z2 = 0.0038856563864182 and x3 = z3 = 0.0274330496637624. Finally we have: −1 m3 mol V 2 ¼ V New ¼ z2 RT . P ¼ 64:5849974569685000 According to the recorded results in Tables 5 and 6, one can see that, in lower temperature region, the value of γ is exactly closed to −3. A cheap analytic method may be introduced. Therefore, for finding the real root of the CES in lower temperature region, one can get as γ = −3; then the one root of z3 + c1z2 + c2z + c3 = 0 is easily obtained
where C1 = 9C − 3B2,D1 = 2B3 − 9BC + 27D and γ = −3. 4. Conclusions For substances, it is clear that, in the low temperature and pressure ranges, their CES will have three real roots. It is shown that, in the low temperature region, the traditional analytical solutions of the CES lead to irrational results, while the iterative solutions of the CES, produced valid results. In this work, as a new analytic method, by making the sub pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi stitution x ¼ − 13 3 ð1−γÞD1 þ B or x ¼ − 13 DC 11 γ þ B , the cubic equation of x3 + Bx2 + Cx + D = 0 was converted to equation γ3 − βγ + β = 0, where γ was related with the left branch of the curves β in Fig. 1 (i.e., γ ≤ 1); and it was obtained from Eq. (5). The magnification of error analysis for present intensified analytic method was done. Also,
Table 5 The details of the new method solutions, for propylene using PR model. T (K)
P (MPa)
β
γ
x2 or z2
V 2 ¼ V New ¼ z2 RT P
87.9 89.4 90.9 92.4 93.9 95.4 96.9 98.4 99.9 101.4 102.9 127.9 152.9 177.9 202.9 227.9
9.18E−10 1.60E−09 2.74E−09 4.59E−09 7.56E−09 1.22E−08 1.95E−08 3.06E−08 4.73E−08 7.21E−08 1.08E−07 2.08E−05 6.03E−04 6.09E−03 3.24E−02 1.14E−01
6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000100 6.7500000000000100 6.7500000000000300 6.7500000000000700 6.7500000000001700 6.7500000000003700 6.7500000000008100 6.7500000000016900 6.7500000212578000 6.7500070534817300 6.7503155582675100 6.7542936421364600 6.7780650055531100
−3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000100 −3.0000000000000100 −3.0000000000000300 −3.0000000000000700 −3.0000000000001600 −3.0000000000003300 −3.0000000041990700 −3.0000013932799600 −3.0000623317418200 −3.0008479870198000 −3.0055377417229900
0.0000000000677154 0.0000000001161804 0.0000000001959108 0.0000000003232501 0.0000000005245479 0.0000000008342117 0.0000000013143666 0.0000000020336576 0.0000000031002592 0.0000000046618211 0.0000000068901351 0.0000010928689822 0.0000272376250939 0.0002441832199728 0.0011837634526550 0.0038856563864182
53.9090441179269000 53.9732093075800000 54.0380360620892000 54.1035440340759000 54.1697359441803000 54.2366202312036000 54.3042045354741000 54.3724973467836000 54.4415067637341000 54.5112415115016000 54.5817100719311000 55.8731414016559000 57.4231736539261000 59.3064941722816000 61.6353326938832000 64.5849974569685000
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
169
Table 6 The details of the new method solutions, for 1-butene using PR model. T (K)
P (MPa)
β
γ
x2 or z2
V 2 ¼ V New ¼ z2 RT P
87.8 91.3 94.8 98.3 101.8 105.3 108.8 112.3 115.8 119.3 122.8 126.3 161.3 196.3 231.3 266.3
3.56E−13 1.87E−12 8.58E−12 3.51E−11 1.29E−10 4.32E−10 1.33E−09 3.79E−09 1.01E−08 2.51E−08 5.92E−08 1.32E−07 4.98E−05 1.80E−03 1.90E−02 9.86E−02
6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000200 6.7500000000001300 6.7500000000006300 6.7500000000027300 6.7500001118430200 6.7500505206741100 6.7522237938572300 6.7770982605060000
−3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000300 −3.0000000000001200 −3.0000000000005400 −3.0000000220924500 −3.0000099793730500 −3.0004392304094800 −3.0053471830713100
0.0000000000000343 0.0000000000001736 0.0000000000007688 0.0000000000030398 0.0000000000108125 0.0000000000350869 0.0000000001047942 0.0000000002900143 0.0000000007513435 0.0000000018169613 0.0000000041739155 0.0000000090723641 0.0000027589488605 0.0000849021935213 0.0007948730581291 0.0037895336524678
70.3240828714209000 70.4704699593276000 70.6244455337848000 70.7823792699400000 70.9432771425772000 71.1077416597013000 71.2758287868141000 71.4476066679774000 71.6231678087690000 71.8026029653191000 71.9860082326774000 72.1734841896668000 74.2978857333483000 76.9829386848662000 80.4539903478669000 85.0957361869379000
a new analytic method was tested for three pure substances as follows: propylene, 1-butene and 1-propylene. Then the results were listed in Tables 2–6. As shown in these tables, the appreciable solutions have been obtained for the abovementioned new analytic method. Also, it was found that, in lower temperature region, the value of γ is exactly closed to −3 so that the one root of z3 + c1z2 + c2z + c3 = 0 was easily obtained from Eq. (13). Appendix A. Analytic extraction of the two formulas in Eq. (5) We start with the following cubic equation as: 3
y þ Py þ Q ¼ 0:
ðA 1Þ
Now, replace y = u + v in Eq. (A-1); then we have: 3
3
u þ v þ ð3uv þ P Þðu þ vÞ þ Q ¼ 0:
ðA 2Þ
Impose a condition for the variables u and v as: 3uv + P = 0, then, we get as: u3 + v3 = − Q and u3v3 = − P3/27. Thus, u3 and v3 are the two roots of the following quadratic equation as:
If β ≤ 27/4, then the sign of the square root
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β2 =4−β3 =27 is posi-
tive; and the formula in Eq. (A-5) can be directly solved. If β N 27/4, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi then the sign of the square root β2 =4−β3 =27 is negative. In this case, we split Eq. (A-5) into the two complex terms as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 γ @ β β2 β3 A iθ − þ ib ¼ re ¼ − þ 4 27 2 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 γ @ β β2 β3 A −iθ − −ib ¼ re : ¼ − − 4 27 2 2
ðA 6Þ
pffiffiffiffiffiffiffiffiffi Thus, we find as: γ = r(eiθ + e− θ) = 2r cos θ, r ¼ − β=3 and θ ¼ pffiffiffiffiffiffiffiffiffi 1=3 cos−1 9= 12β . Therefore, the formula in Eq. (5) is extracted as follows: rffiffiffiffi β γ ¼ −2 cosθ; 3
! 1 9 −1 pffiffiffiffiffiffiffiffiffi ; θ ¼ cos 3 12β
for β ≥
27 : 4
ðA 7Þ
3
2
z þ Q z−
P ¼ 0: 27
ðA 3Þ
By solving this quadratic equation, we find as follows:u3 ¼ − Q2 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi Q2 Q Q2 P3 3 P3 4 þ 27 and v ¼ − 2 þ 4 þ 27. Therefore, we get as follows: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u u u 3 3 t Q Q 2 P3 t Q Q 2 P3 þ − − þ þ : y¼uþv¼ − þ 4 27 4 27 2 2
ðA 4Þ
Now, by comparing Eq. (4): γ3 − βγ + β = 0 with Eq. (A-1): y3 + Py + Q = 0, we have as: γ = y, β = − P = Q. Therefore, the first formula in Eq. (5) is extracted as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 3 3 2 3 β β β β β2 β3 A − A þ @− − − : γ ¼ @− þ 4 27 4 27 2 2 0
ðA 5Þ
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
S.A. Hong, J.D. Kim, J. Kim, J.W. Kang, I.J. Kang, J. Ind. Eng. Chem. 16 (2010) 859. S. Jiuxun, Fluid Phase Equilib. 193 (2002) 1. M.M. Papari, J. Moghadasi, M. Hosseini, F. Akbari, J. Mol. Liq. 158 (2011) 57. P.M. Mathias, H.C. Klotz, Chem. Eng. Prog. 6 (1994) 67. H. Orbey, S.I. Sandler, Modeling vapor–liquid equilibria, Cubic Equations of State and Their Mixing Rules, Cambridge University Press, Cambridge, U.K., 1998 F.J.A. Guevara-Rodriguez, Fluid Phase Equilib. 307 (2011) 190. D.Y. Peng, D.B.A. Robinson, Industrial and Engineering Chemistry: Fundamentals, 151976. 59. N. Jacobson, Basic Algebra, 2nd, Dover, New York, 2009. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes, The Art of Scientific Computing3rd., Cambridge University Press, Cambridge, U.K., 2007 D.E. Smith, History of Mathematics, Dover, New York, 2003. T. Zhao, D. Wang, H. Hong, J. Symb. Comput. 46 (2011) 904. Y. Zhi, H. Lee, Fluid Phase Equilib. 201 (2002) 287. P.H. Salim, Fluid Phase Equilib. 240 (2006) 224. R.M. Loperena, Ind. Eng. Chem. Res. 51 (2012) 6972. M.Ç. Koçak, J. Comput. Appl. Math. 235 (2011) 4736.
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
An intensified analytic solution for finding the roots of a cubic equation of state in low temperature region Behrooz M. Ziapour Department of Mechanical Engineering, University of Mohaghegh Ardabili, Ardabil, Iran
a r t i c l e
i n f o
Article history: Received 19 December 2014 Received in revised form 3 February 2015 Accepted 16 February 2015 Available online 19 February 2015 Keywords: Analytic Cubic Equation of state Low temperature Organic
a b s t r a c t Because of the effect of the enlarging process of errors in low temperature region, the analytical solutions of a cubic equation of state (CES) lead to illogical results. To overcome this problem, we recommend an enhanced analytic method, in this paper. Up to now, there is not any analytical method for this problem. In the present study, the general cubic equation x3 + Bx2 + Cx + D = 0 is transformed into simple equation γ3 − βγ + β = 0. For β ≥ 27/4, then the three roots are real. To compare, as an example, for 1-butene at T = 112.3 (K) and P = 3.79 × 10−9(Pa), the empirical value of the molar volume is 72.00 (m3mol−1). By using Patel–Teja (PT) CES, the molar volume value is obtained as −342.3 from traditional analytic method and 73.88 from both the new intensified analytic method and the iterative Newton–Raphson method. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The CES models are useful formulas for analyzing and correlating the thermodynamic equilibrium state of simple fluids and their mixtures, with agreeable approach [1–5]. Generally, for a given substance, the CES is an empirical relation of the pressure P, the temperature T and the molar volume v consist of a repulsive term and an attractive term as follows [6]: P ðT; vÞ ¼
RT α ðT Þ − v−b ðv−cÞðv−dÞ
ð1Þ
where R is the gas constant. b, c and d are constant parameters which are obtained based on the kind of equation. In Eq. (1), α(T) is a function and it has many empirical definitions in the literature. By expressing Eq. (1) in terms of molar volume, the following equation is obtained: RT 2 RT RT α ðT Þ 3 v þ bc þ bd þ cd þ cþ dþ v v − bþcþdþ P P P P RT α ðt Þ − bcd þ cd þ b ¼ 0: P P
ð2Þ
Eq. (2) is referred to as cubic in molar volume. Additionally we can substitute the definition of compressibility factor (i.e., Z = Pv/RT) into Eq. (2) and obtain different cubic equations in Z. The first CES was recommended by van der Waals in 1873 as P(T, v) = RT/(v − b) − α(T)/v2. It has been proven that the original van der Waals E-mail address: [email protected].
http://dx.doi.org/10.1016/j.molliq.2015.02.026 0167-7322/© 2015 Elsevier B.V. All rights reserved.
CES cannot provide simultaneous accurate predictions for all the properties of pure fluids and their mixtures. But, the most CES models (e.g., Peng-Robison CES in 1976 as P(T, v) = RT/(v − b) − α(T)/(v2 + 2bv − b2) being used widely today for practical design purposes have been derived from van der Waals CES [7]. For finding the roots of Eq. (2), the analytical methods can solve their three roots simultaneously. One of the relevant math formulas is the well-known formula due to Ferro and Tartaglia, communicated by Gerolamo Cardano in 1545 [8,9]. Cardano's method involves division. Thus it can encounter near “0/0”, resulting insignificant numerical errors. Thus it may encounter a numerically unstable case. Lagrange in 1770 gave the new method for solving the three roots of a cubic equation, with 18 possible interpretations [10]. The Lagrange formula, does not require division and thus avoids the “0/0” case. In the Lagrange formula some interpretations are correct (yielding the solutions), But the others are not. Zhao et al. [11] proposed a convection which provided correct interpretations of the Lagrange formula for all cubic equations with real coefficients. Zhi and Lee [12] revealed by several examples that, in low temperature region, the analytical solutions of CES lead to irrational results, while the iterative solutions of the CES, using Newton–Raphson method produced valid results. They appeared that errors caused by the limitation of significant figures of the computer languages are revealed, and a magnification of errors is defined which a main factor is bringing out the irrational results of the analytical solution of CES. Salim [13] showed that, in the case of low temperature region, the coefficient D in cubic equation x3 + Bx2 + Cx + D = 0 is very small. Therefore the main root is also small and becomes so close to zero. Thus, by eliminating the term of x3 in x3 + Bx2 + Cx + D = 0, then the main root can be
166
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
found via the approximated equation of the form Bx2 + Cx + D = 0. Loperena [14] proposed an iterative refinement of the solution obtained by the analytical method. This method allowed one to take advantage of the calculations fulfilled in the application of the analytical solution. His numerical results show that the proposed iterative refinement procedure is an attractive and numerical inexpensive option and can be easily incorporated into any process simulation program, that use Cardano's method, without significant modifications. Up to now, as mentioned above discussions, there is not any completely analytical method for this problem. To notice to complexity and overhead of analytical methods, it has been recommended to use the numeric iterative methods. In this paper, we describe an enhanced completely analytic method for solving the general cubic equation x3 + Bx2 + Cx + D = 0 consists of real coefficients. We verify this new method with the miscellaneous examples.
formula as: 8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 3 > 2 3 > β β β β β2 β3 A 27 > > @ A @ > − − γ ¼ − þ − ; for β ≤ þ − < 4 27 4 27 2 2 4 ! rffiffiffiffi > > > β 1 9 27 −1 > > pffiffiffiffiffiffiffiffiffi ; for β ≥ cosθ; θ ¼ cos : : γ ¼ −2 3 3 4 12β ð5Þ
2. New analytic solution and the miscellaneous examples
It is clear that (see Fig. 1), when β → − ∞, then γ → 1. Also, when β → + ∞, then γ → − ∞. To obtain the two formulas in Eq. (5), an analytic method was completely used. This method has been described in Appendix A. Therefore, through an intensified analytic solution, the real root of the cubic Eq. (3) is captured. The other roots of Eq. (3) i.e., x2, x3 can be calculated analytically from the following quadratic equation as:
An arbitrary cubic equation which consists of real coefficients B, C and D is as follows:
x þ ðB þ x1 Þx−
3
2
x þ Bx þ Cx þ D ¼ 0:
ð3Þ C3
In Eq. (3), one of the three roots is real. To find it, let β ¼ − D12 , 1
where C1 = 9C − 3B2 and D1 = 2B3 − 9BC + 27D. Make the substi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tution x ¼ − 13 3 ð1−γÞD1 þ B or x ¼ − 13 DC 11 γ þ B . By this substitution, Eq. (3) is transformed as the following simple equation for γ: 3
γ −βγ þ β ¼ 0:
ð4Þ
γ To find γ, we arrange Eq. (4) into the functional form as βðγ Þ ¼ γ−1 . 3
Fig. 1 shows the sketch of β(γ). It is seen that for β ≥ 27 4 , then all three roots are real. Therefore, Eq. (3) has three real roots (i.e., x1, x2, x3). Also, for β b 27 4 , Eq. (3) has one real root (i.e., x1). In this case, the others are conjugate imaginary roots (i.e., x2, x3). Therefore, like a new analytic method, for finding the one real root of pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a cubic equation, let x1 ¼ − 13 3 ð1−γ ÞD1 þ B or x1 ¼ − 13 DC 11 γ þ B be one real root of Eq. (3), where γ is related with the left branch of the curves β in Fig. 1 (i.e., γ ≤ 1); and it is obtained from the following
2
D ¼ 0: x1
ð6Þ
Note that, in order to reach the high precision of the solutions, especially in low temperature region, the two other roots may be obtained from Eq. (6). In this paper, all solutions were conducted by Microsoft Excel, 2003. Example 1. In equation: x3 − 3x2 + 4x − 2 = 0, the coefficients are: B = − 3, C = 4 and D = − 2. Then, C1 = 9 and D1 = 0. Therefore, β = − ∞. Thus, from Fig. 1: γ = 1. Therefore, x1 = − B/3 = 1. Finally from Eq. (6): x2 − 2x + 2 = 0. Then, x2 = 1 + i and x3 = 1 − i. Example 2. In equation: x3 + 6x2 + 3x − 10 = 0, the coefficients are: B = 6, C = 3, D = −10. Then C1 = −81 and D1 = 0. Therefore, β = + ∞. Thus, from Fig. 1: γ = − ∞. Therefore, x1 = − B/3 = −2. Finally from Eq. (6): x2 + 4x − 5 = 0. Then, x2 = 1 and x3 = −5. Example 3. In equation: x3 − 2x2 − 5x + 6 = 0, the coefficients are: B = − 2, C = − 5, D = 6. Then C1 = − 57 and D1 = 56. Then, β = 59.0539. Thus, from Eq. (5):γ = −8.1429. Therefore, x1 = −2.0000. Finally from Eq. (6):x2 − 4x + 3 = 0. Then, x2 = 1.0000 and x3 = 3.0000. Example 4. The cubic equation: (x − 1)2x = x3 − 2x2 + x = 0. The coefficients are: B = −2, C = 1, D = 0. Then C1 = −3 and D1 = 2. Then, β = 6.75. Thus, from Eq. (5):γ = −3. Therefore, x1 = 0. Finally from solution of the equation:x2 − 2x + 1 = 0, we have:x2 = x3 = 1.0000. Example 5. The molar volume of carbon dioxide (CO 2 ) at P = 10 atm and T = 300 K, by using the van der Waals CES (P(T, v) = RT/(v − b) − α(T)/v 2 ), is the cubic equation as: v 3 − 2.5046v 2 + 0.3598v − 0.0155 = 0. This equation and another equation in the form of v 3 − 7.8693v 2 + 13.3771v − 6.5354 = 0
Table 1 Results for example 5. Case I: v3 − 2.5046v2 + 0.3598v − 0.0155 = 0 Newton (numeric)
Ostrowski [15] (numeric)
New (analytic) Fig. 1. The curve of β(γ).
2:4616 2:3633 Iterations→ 2:3546 2:3545→v1 2:4616 Iterations→ 2:3545→v1 β ¼ 6:7165 γ ¼ −2:9934 v1 ¼ 2:3545
Case II: v3 − 7.8693v2 + 13.3771v − 6.5354 = 0 7:3807 6:2992 Iterations→ 5:8355 5:7397 5:7357→v1 7:3807 Iterations→ 5:8187 5:7357→v1 β ¼ 6:7396 γ ¼ −2:9977 v1 ¼ 5:7357
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
were solved by the some numerical iteration methods [15]. The initial guess for starting the iteration was selected from the ideal gas law. We collected both these numerical solutions and our new analytical calculations, in Table 1. As shown from the new solutions, the vales of β are less than 6.75. Thus, the selected cubic equations have only one real root (i.e., v1 ≡ x1). The other imaginary roots are unusable. 3. Results and discussions For substances, it is clear that, in the low temperature and pressure ranges, their CES will have three real roots (i.e., β ≥ 6.75). In order to compare our new method with the traditional methods, taken from Zhi and Lee [12], we assume that CES be formed versus the compressibility factor z as: z3 + c1z2 + c2z + c3 = 0. The smallest root (i.e., the compressibility factor of liquid phase:z2) is obtained from the traditional method as follows: pffiffiffiffiffiffiffiffiffiffi 0 1 arccos r= −q3 pffiffiffiffiffiffiffiffi c ∘ þ 120 A− 1 z2 ¼ 2 −q cos@ 3 3
167
Table 3 Comparison among experimental, old and new analytical methods calculated liquid volume for 1-butene (TC = 419.5 K, PC = 4020 kPa, ω = 0.194). P (KPa) P (MPa)
87.8 91.3 94.8 98.3 101.8 105.3 108.8 112.3 115.8 119.3 122.8 126.3 161.3 196.3 231.3 266.3
3.56E−13 1.87E−12 8.58E−12 3.51E−11 1.29E−10 4.32E−10 1.33E−09 3.79E−09 1.01E−08 2.51E−08 5.92E−08 1.32E−07 4.98E−05 1.80E−03 1.90E−02 9.86E−02
VExp
70.09 70.35 70.62 70.89 71.16 71.44 71.71 72.00 72.28 72.57 72.86 73.16 76.34 80.00 84.31 89.52
PR
PT
VNew
VOld
VNewton VNew
VOld
VNewton
70.32 70.47 70.62 70.78 70.94 71.11 71.28 71.45 71.62 71.80 71.99 72.17 74.30 76.98 80.45 85.10
1685.9 1603.6 – −92,975 −25,216 1332.8 – 0.460 116.28 78.07 72.18 72.91 74.30 76.98 80.45 85.10
70.32 70.47 70.62 70.78 70.94 71.11 71.28 71.45 71.62 71.80 71.99 72.17 74.30 76.98 80.45 85.10
– 1621.3 1544.6 – 1408.3 1347.3 1290.5 −342.3 129.1 72.17 73.20 74.45 76.89 79.73 83.39 88.27
72.68 72.84 73.01 73.17 73.34 73.52 73.69 73.88 74.06 74.25 74.45 74.64 76.89 79.73 83.39 88.27
72.68 72.84 73.01 73.17 73.34 73.52 73.69 73.88 74.06 74.25 74.45 74.64 76.89 79.73 83.39 88.27
ð7Þ
1 3 c1 . Then the molar volwhere, q ¼ 19 3c2 −c21 and r ¼ 16 c1 c2 − 12 c3 − 27 3 −1 ume of the liquid phase is given as: V Old ¼ z2 RT ). Suppose P (m mol that, VExp, VNew, VOld and VNewton be denoted to the molar volume of the liquid phase obtained from the experimental, the new method, the old method and the iterative Newton–Rephson method, respectively. These values for three pure substances (propylene, 1-butene and 1pentene), have been listed in Tables 2–4, by using both the PR and the PT CES models. In these tables, to find the VNewton, we consider Newton 3 2 2 scheme as: zkþ1 ¼ zk − zk þ c1 zk þ c2 zk þ c3 = 3zk þ 2c1 zk þ c2 .
The recorded results for VNewton in Tables 2–4, were obtained due to the three iterations as k = 3. In this iterative process, the initial guess for zk was as z1 = 0. As shown in Tables 2–4, the results indicate that both the present intensified analytic method and the iterative Newton–Rephson method have identical accuracies. But, because of the effect of the enlarging process of errors, especially in lower temperature region, the traditional analytic solutions of CES lead to illogical results. The magnification of error was analyzed for traditional method (i.e., Eq. (7)) by Zhi and Lee [12]. Through this analysis, it was found that the pffiffiffiffiffiffiffiffiffiffiffiffi critical term in the magnification of error is 1= 1−s2 , where, the value pffiffiffiffiffiffiffiffiffiffi of s is as: s ¼ r= −q3. They showed that in low temperature region, the values of s do tend to 1, so the error magnification (M) will become
Table 2 Comparison among experimental, old and new analytical methods calculated liquid volume for propylene (TC = 364.8 K, PC = 4613 kPa, ω = 0.142). P (KPa)
P (MPa)
VExp
87.9 89.4 90.9 92.4 93.9 95.4 96.9 98.4 99.9 101.4 102.9 127.9 152.9 177.9 202.9 227.9
9.18E−10 1.60E−09 2.74E−09 4.59E−09 7.56E−09 1.22E−08 1.95E−08 3.06E−08 4.73E−08 7.21E−08 1.08E−07 2.08E−05 6.03E−04 6.09E−03 3.24E−02 1.14E−01
55.07 55.18 55.29 55.39 55.50 55.61 55.72 55.84 55.95 56.06 56.17 58.17 60.39 62.91 65.79 69.16
PR
PT
VNew
VOld
VNewton VNew
VOld
VNewton
53.91 53.97 54.04 54.10 54.17 54.24 54.30 54.37 54.44 54.51 54.58 55.87 57.42 59.31 61.64 64.58
1002.3 980.71 −158.8 108.41 82.69 67.72 53.97 61.62 55.69 54.45 54.92 55.87 57.42 59.31 61.64 64.58
53.91 53.97 54.04 54.10 54.17 54.24 54.30 54.37 54.44 54.51 54.58 55.87 57.42 59.31 61.64 64.58
– 995.94 −395.28 275.55 96.98 81.71 42.68 56.99 60.93 58.77 58.24 58.76 60.44 62.48 65.00 68.18
56.62 56.69 56.76 56.83 56.91 56.98 57.05 57.13 57.20 57.28 57.35 58.76 60.44 62.48 65.00 68.18
56.62 56.69 56.76 56.83 56.91 56.98 57.05 57.13 57.20 57.28 57.35 58.76 60.44 62.48 65.00 68.18
infinity. The magnification of error analysis for present intensified analytic method was done as the following process: pffiffiffiffiffiffiffiffiffi Let s ¼ 9= 12β, so the second formula of Eq. (5) can be arranged as: γ¼
−3 1 −1 cos cos ðsÞ : s 3
ð8Þ
When Δs tend to zero, the magnification of error (M) can expressed as [12]: M¼
dγ s ¼ ds γ
3 1 1 1 s −1 −1 p ffiffiffiffiffiffiffiffiffiffiffi ffi cos ð s Þ þ ð s Þ cos cos : ð9Þ sin 2 3 3 γ s2 s 1−s
In low temperature region, the values of β, γ and s do tend to 27/4, 3 and 1 respectively. When s = 1, then the second term in Eq. (9) gets p1ffiffiffiffiffiffiffi sin 1 cos−1 ð1Þ ¼ 0. Let 3 0
1 1−1
dγ s ds γ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 3 1 1 1 s −1 −1 cos ðsÞ þ pffiffiffiffiffiffiffiffiffiffiffiffi 1− cos2 ð cos ðsÞÞ : ð10Þ ¼ 2 cos 3 3 γ s s 1−s2
M¼
Or, by inserting Eq. (8), we can write as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 0 dγ s @ 3 1 1 s2 γ 2 A s −1 cos ¼ cos ðsÞ þ pffiffiffiffiffiffiffiffiffiffiffiffi 1− : M¼ 2 9 ds γ 3 γ s2 s 1−s
ð11Þ
As mentioned above, in low temperature region, we have as: γ → − 3. Finally, by inserting this value into Eq. (11), then the value of M is obtained as follows: M¼
dγ s ¼ ds γ
3 1 1 s −1 cos ð s Þ þ cos : 3 s γ s2
ð12Þ
Therefore, as seen from value of M in Eq. (12), in low temperature region, there is not any infinity term for value of M, which is the reason why this present method can give corrective calculation results. For a better presentation, we have listed the new solution procedure in Tables 5 and 6. Tables 5 and 6 have been listed based on the PR model for propylene and 1-butane, respectively. In the following example, more details of the new method calculations have been introduced.
168
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
Table 4 Comparison among experimental, old and new analytical methods calculated liquid volume for 1-pentene (TC = 464.4 K, PC = 3555 kPa, ω = 0.233). P (KPa)
107.93 111.43 114.93 118.43 121.93 125.43 128.93 123.43 135.93 139.43 142.93 146.43 149.93 184.93 219.93 254.93
P (MPa)
4.48E−12 1.67E−11 5.70E−11 1.79E−10 5.24E−10 1.43E−09 3.68E−09 8.92E−09 2.06E−08 4.52E−08 9.50E−08 1.92E−07 3.73E−07 5.93E−05 1.48E−03 1.32E−02
VExp
86.64 86.94 87.25 87.56 87.87 88.19 88.51 88.83 89.16 89.49 89.83 90.17 90.51 94.18 98.36 103.21
PR
PT
VNew
VOld
VNewton
VNew
VOld
VNewton
88.38 88.56 88.74 88.92 89.10 89.29 89.48 89.18 89.88 90.08 90.29 90.50 90.72 93.16 96.19 100.03
−810,425 1849.6 – 1706.4 −6205 1580.1 341.14 139.16 91.06 94.20 90.84 90.39 90.78 93.16 96.19 100.03
88.38 88.56 88.74 88.92 89.10 89.29 89.48 89.18 89.88 90.08 90.29 90.50 90.72 93.16 96.19 100.03
90.31 90.49 90.67 90.86 91.05 91.25 91.45 91.14 91.86 92.07 92.29 92.51 92.73 95.26 98.40 102.37
1944 1864.2 1789.6 1719.7 1654.1 1592.5 87.49 79.74 83.53 96.38 94.34 92.43 92.78 95.26 98.40 102.37
90.31 90.49 90.67 90.86 91.05 91.25 91.45 91.14 91.86 92.07 92.29 92.51 92.73 95.26 98.40 102.37
Example 6. From Table 5, for propylene using PR model, at T = 227.9(K) and P = 114(kPa), then the following CES is obtained:
from the following equation as: 3
z1 ¼ 3
2c1 −9c1 c2 þ 27c3 c1 − : 3 9c2 −3c1 2
ð13Þ
2
z −0:9969224047840840z þ 0:0303480538062392z−0:0001029289169739 ¼ 0:
Please note that, Eq. (13) has been adjusted based on the cubic equa D1 C 1 γ þ B ),
tion as: x3 + Bx2 + Cx + D = 0 and its one real root (x1 ¼ − 13 Therefore, C1 = − 2.7084303592252900 and D1 = − 1.7120783876879200. Thus: β = 6.7780650055531100. Then, θ = 0.0214638985371484. Hence from Eq. (5): = z1 = γ = − 3.0055377417229900. Then, x1 0.9656036987339030. From Eq. (6): x2 = z2 = 0.0038856563864182 and x3 = z3 = 0.0274330496637624. Finally we have: −1 m3 mol V 2 ¼ V New ¼ z2 RT . P ¼ 64:5849974569685000 According to the recorded results in Tables 5 and 6, one can see that, in lower temperature region, the value of γ is exactly closed to −3. A cheap analytic method may be introduced. Therefore, for finding the real root of the CES in lower temperature region, one can get as γ = −3; then the one root of z3 + c1z2 + c2z + c3 = 0 is easily obtained
where C1 = 9C − 3B2,D1 = 2B3 − 9BC + 27D and γ = −3. 4. Conclusions For substances, it is clear that, in the low temperature and pressure ranges, their CES will have three real roots. It is shown that, in the low temperature region, the traditional analytical solutions of the CES lead to irrational results, while the iterative solutions of the CES, produced valid results. In this work, as a new analytic method, by making the sub pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi stitution x ¼ − 13 3 ð1−γÞD1 þ B or x ¼ − 13 DC 11 γ þ B , the cubic equation of x3 + Bx2 + Cx + D = 0 was converted to equation γ3 − βγ + β = 0, where γ was related with the left branch of the curves β in Fig. 1 (i.e., γ ≤ 1); and it was obtained from Eq. (5). The magnification of error analysis for present intensified analytic method was done. Also,
Table 5 The details of the new method solutions, for propylene using PR model. T (K)
P (MPa)
β
γ
x2 or z2
V 2 ¼ V New ¼ z2 RT P
87.9 89.4 90.9 92.4 93.9 95.4 96.9 98.4 99.9 101.4 102.9 127.9 152.9 177.9 202.9 227.9
9.18E−10 1.60E−09 2.74E−09 4.59E−09 7.56E−09 1.22E−08 1.95E−08 3.06E−08 4.73E−08 7.21E−08 1.08E−07 2.08E−05 6.03E−04 6.09E−03 3.24E−02 1.14E−01
6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000100 6.7500000000000100 6.7500000000000300 6.7500000000000700 6.7500000000001700 6.7500000000003700 6.7500000000008100 6.7500000000016900 6.7500000212578000 6.7500070534817300 6.7503155582675100 6.7542936421364600 6.7780650055531100
−3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000100 −3.0000000000000100 −3.0000000000000300 −3.0000000000000700 −3.0000000000001600 −3.0000000000003300 −3.0000000041990700 −3.0000013932799600 −3.0000623317418200 −3.0008479870198000 −3.0055377417229900
0.0000000000677154 0.0000000001161804 0.0000000001959108 0.0000000003232501 0.0000000005245479 0.0000000008342117 0.0000000013143666 0.0000000020336576 0.0000000031002592 0.0000000046618211 0.0000000068901351 0.0000010928689822 0.0000272376250939 0.0002441832199728 0.0011837634526550 0.0038856563864182
53.9090441179269000 53.9732093075800000 54.0380360620892000 54.1035440340759000 54.1697359441803000 54.2366202312036000 54.3042045354741000 54.3724973467836000 54.4415067637341000 54.5112415115016000 54.5817100719311000 55.8731414016559000 57.4231736539261000 59.3064941722816000 61.6353326938832000 64.5849974569685000
B.M. Ziapour / Journal of Molecular Liquids 206 (2015) 165–169
169
Table 6 The details of the new method solutions, for 1-butene using PR model. T (K)
P (MPa)
β
γ
x2 or z2
V 2 ¼ V New ¼ z2 RT P
87.8 91.3 94.8 98.3 101.8 105.3 108.8 112.3 115.8 119.3 122.8 126.3 161.3 196.3 231.3 266.3
3.56E−13 1.87E−12 8.58E−12 3.51E−11 1.29E−10 4.32E−10 1.33E−09 3.79E−09 1.01E−08 2.51E−08 5.92E−08 1.32E−07 4.98E−05 1.80E−03 1.90E−02 9.86E−02
6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000000 6.7500000000000200 6.7500000000001300 6.7500000000006300 6.7500000000027300 6.7500001118430200 6.7500505206741100 6.7522237938572300 6.7770982605060000
−3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000000 −3.0000000000000300 −3.0000000000001200 −3.0000000000005400 −3.0000000220924500 −3.0000099793730500 −3.0004392304094800 −3.0053471830713100
0.0000000000000343 0.0000000000001736 0.0000000000007688 0.0000000000030398 0.0000000000108125 0.0000000000350869 0.0000000001047942 0.0000000002900143 0.0000000007513435 0.0000000018169613 0.0000000041739155 0.0000000090723641 0.0000027589488605 0.0000849021935213 0.0007948730581291 0.0037895336524678
70.3240828714209000 70.4704699593276000 70.6244455337848000 70.7823792699400000 70.9432771425772000 71.1077416597013000 71.2758287868141000 71.4476066679774000 71.6231678087690000 71.8026029653191000 71.9860082326774000 72.1734841896668000 74.2978857333483000 76.9829386848662000 80.4539903478669000 85.0957361869379000
a new analytic method was tested for three pure substances as follows: propylene, 1-butene and 1-propylene. Then the results were listed in Tables 2–6. As shown in these tables, the appreciable solutions have been obtained for the abovementioned new analytic method. Also, it was found that, in lower temperature region, the value of γ is exactly closed to −3 so that the one root of z3 + c1z2 + c2z + c3 = 0 was easily obtained from Eq. (13). Appendix A. Analytic extraction of the two formulas in Eq. (5) We start with the following cubic equation as: 3
y þ Py þ Q ¼ 0:
ðA 1Þ
Now, replace y = u + v in Eq. (A-1); then we have: 3
3
u þ v þ ð3uv þ P Þðu þ vÞ þ Q ¼ 0:
ðA 2Þ
Impose a condition for the variables u and v as: 3uv + P = 0, then, we get as: u3 + v3 = − Q and u3v3 = − P3/27. Thus, u3 and v3 are the two roots of the following quadratic equation as:
If β ≤ 27/4, then the sign of the square root
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β2 =4−β3 =27 is posi-
tive; and the formula in Eq. (A-5) can be directly solved. If β N 27/4, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi then the sign of the square root β2 =4−β3 =27 is negative. In this case, we split Eq. (A-5) into the two complex terms as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 γ @ β β2 β3 A iθ − þ ib ¼ re ¼ − þ 4 27 2 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 3 γ @ β β2 β3 A −iθ − −ib ¼ re : ¼ − − 4 27 2 2
ðA 6Þ
pffiffiffiffiffiffiffiffiffi Thus, we find as: γ = r(eiθ + e− θ) = 2r cos θ, r ¼ − β=3 and θ ¼ pffiffiffiffiffiffiffiffiffi 1=3 cos−1 9= 12β . Therefore, the formula in Eq. (5) is extracted as follows: rffiffiffiffi β γ ¼ −2 cosθ; 3
! 1 9 −1 pffiffiffiffiffiffiffiffiffi ; θ ¼ cos 3 12β
for β ≥
27 : 4
ðA 7Þ
3
2
z þ Q z−
P ¼ 0: 27
ðA 3Þ
By solving this quadratic equation, we find as follows:u3 ¼ − Q2 þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi Q2 Q Q2 P3 3 P3 4 þ 27 and v ¼ − 2 þ 4 þ 27. Therefore, we get as follows: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u u u 3 3 t Q Q 2 P3 t Q Q 2 P3 þ − − þ þ : y¼uþv¼ − þ 4 27 4 27 2 2
ðA 4Þ
Now, by comparing Eq. (4): γ3 − βγ + β = 0 with Eq. (A-1): y3 + Py + Q = 0, we have as: γ = y, β = − P = Q. Therefore, the first formula in Eq. (5) is extracted as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 3 3 2 3 β β β β β2 β3 A − A þ @− − − : γ ¼ @− þ 4 27 4 27 2 2 0
ðA 5Þ
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
S.A. Hong, J.D. Kim, J. Kim, J.W. Kang, I.J. Kang, J. Ind. Eng. Chem. 16 (2010) 859. S. Jiuxun, Fluid Phase Equilib. 193 (2002) 1. M.M. Papari, J. Moghadasi, M. Hosseini, F. Akbari, J. Mol. Liq. 158 (2011) 57. P.M. Mathias, H.C. Klotz, Chem. Eng. Prog. 6 (1994) 67. H. Orbey, S.I. Sandler, Modeling vapor–liquid equilibria, Cubic Equations of State and Their Mixing Rules, Cambridge University Press, Cambridge, U.K., 1998 F.J.A. Guevara-Rodriguez, Fluid Phase Equilib. 307 (2011) 190. D.Y. Peng, D.B.A. Robinson, Industrial and Engineering Chemistry: Fundamentals, 151976. 59. N. Jacobson, Basic Algebra, 2nd, Dover, New York, 2009. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes, The Art of Scientific Computing3rd., Cambridge University Press, Cambridge, U.K., 2007 D.E. Smith, History of Mathematics, Dover, New York, 2003. T. Zhao, D. Wang, H. Hong, J. Symb. Comput. 46 (2011) 904. Y. Zhi, H. Lee, Fluid Phase Equilib. 201 (2002) 287. P.H. Salim, Fluid Phase Equilib. 240 (2006) 224. R.M. Loperena, Ind. Eng. Chem. Res. 51 (2012) 6972. M.Ç. Koçak, J. Comput. Appl. Math. 235 (2011) 4736.