A Note on the Analytical Solution of Cubic Equations of State in

Research Note pubs.acs.org/IECR

A Note on the Analytical Solution of Cubic Equations of State in Process Simulation Rosendo Monroy-Loperena* ROMON Paseo de los Pirules 124, 04250 Distrito Federal, México Departamento de Energı ́a Universidad Autónoma Metropolitana-Azcapotzalco. Avenida San Pablo Xalpa 180, 02200, Distrito Federal, México ABSTRACT: This work address the problem of round-off errors in the analytical solution (Cardano’s Method) of cubic equations of state inside the simulation of chemical processes in the low-temperature region, as it is the case of cryogenic processes involving hydrates calculations. It is proposed as a strategy that can be taken as an iterative refinement of the solution obtained by the analytical method and allows one to take advantage of the calculations fulfilled in the application of the analytical solution (Cardano’s Method). Numerical experimentations are presented in order to show the application of the proposed strategy.

1. INTRODUCTION Cubic equations of state are widely used in process simulation.1 They have long been applied to the representation of properties for pure compounds and mixtures. Well-known cubic equations of state are those proposed by Soave2 and Peng−Robinson.3 Even when the theoretical basis of cubic equations of state is not strong, they can, in principle, accurately represent the relation among temperature, pressure, and phase compositions in binary and multicomponent systems. Further developments on the cubic equation on state have improved their predictive capabilities, particularly for polar mixtures. (See, for instance, the work of Orbey and Sandler.4) The cubic equations of state can be written in a general form as follows: a(T ) RT p= − v−b (v + δ1b)(v + δ2b)

a(T ) = a(Tc)f (t )

where f(t) is a characteristic function of the system under study. Values of δ1, δ2, Ωa, and Ωb are given in Table 1 for the Soave Table 1. Values of the Parameters of the Generalized Cubic Equation of State for the Soave−Redlich−Kwong and Peng− Robinson Equations of State δ1 Soave−Redlich−Kwong Peng−Robinson

δ2

1

1+

0

2

1−

2

Ωa

Ωb

0.42748 0.45724

0.08664 0.07780

and Peng−Robinson cubic equations of state. The above parameters can be replaced by the adimensional ones:

(1)

A=

where δ1 and δ2 are specific constants for each cubic equation of state, p is the pressure, R is the universal gas constant, T is the absolute temperature, v is the molar volume, b is the covolume (which takes into account the packing of spherical molecules of finite size,) and a is the cohesion parameter (usually called the attractive term, which represents attractive interactions between molecules). The values of a and b at the critical temperature are found by setting the first and second derivatives of pressure, with respect to volume to zero at the critical point; that is, ⎛ R2T 2 ⎞ c ⎟ a(Tc) = Ωa⎜⎜ ⎟ p ⎝ c ⎠

(4)

B=

⎛ p ⎞ = Ωa⎜ r2 ⎟f RT ⎝ Tr ⎠

(5)

⎛p⎞ bp = Ω b⎜ r ⎟ RT ⎝ Tr ⎠

(6)

ap

2 2

Ω ⎛f⎞ A a = = a⎜ ⎟ B RTb Ωb ⎝ Tr ⎠

(7)

where the subscript “r” denotes reduced conditions. Introduc(2)

⎛ RT ⎞ b = Ωb⎜⎜ c ⎟⎟ ⎝ pc ⎠

ing the compressibility factor, defined as Z = pv/(RT), eq 1 can be written in the cubic form,

(3)

Received: Revised: Accepted: Published:

where the subscript “c” denotes the critical point. The parameter b in eq 1 is assumed to be independent of temperature, and a is a function of temperature, written as © 2012 American Chemical Society

6972

October 6, 2011 February 7, 2012 April 11, 2012 April 11, 2012 dx.doi.org/10.1021/ie2023004 | Ind. Eng. Chem. Res. 2012, 51, 6972−6976

Industrial & Engineering Chemistry Research

Research Note

have reported that the analytical solution of the Peng− Robinson equation of state, via Cardano’s Method, for the calculation of liquid−solid equilibria in cryogenic process, specifically to predict CO2 freezing points needed for the design of cryogenic systems, can produce meaningless results, since the required calculations are sensitive to round-off errors. Because of the roundoff obtained when three real roots must be found, diverse solutions have been proposed to face this problem; the majority of them are based in the application of iterative Newton-like methods9−11 and others by means of eigenvalue-based polynomial solutions12 Another raised solution consists of finding a real root, by some iterative method, then perform a synthetic division and solve the resultant polynomial by means of the general formula of second-order polynomials.11 Even more, it is also suggested to consider the cubic equation as quadratic when d3 in eq 9 is very small.13 In order to correct the round-off errors of Cardano’s Method, because of several referrals to root or trigonometric functions, when three real roots are found, in the following section, a strategy for reducing the computational errors of the calculation is posed, followed by some numerical aspects and commentaries of this strategy; similarly, numerical experimentations appear to show the efficiency of the proposed strategy, and, finally, conclusions of the raised strategy are given.

Z3 − Z2[1 − B(δ1 + δ2 − 1)] ⎡ A⎤ − Z ⎢δ1 + δ2 − B(δ1δ2 − δ1 − δ2) − ⎥B ⎣ B⎦ − [A + B(1 + B)δ1δ2]B (8)

=0

Kamath et al.5 described the appropriate way of assigning the roots to the liquid and vapor phases, based on the first and second derivatives of the cubic equation. According to many authors (c.f. Tester and Model,6 Prausnitz et al.7), for cubic equations of state, using the form of eq 8, analytical solutions can be readily determined by the method proposed by Cardano8 to solve the cubic equation, which is outlined below. For the cubic equation x 3 + d1x 2 + d 2x + d3 = 0

(9)

with real coefficients d1, d2, and d3, first compute Q = − 3d2)/9 and R = (2d13 − 9d1d2 + 27d3)/54. If R2 ≤ Q3, then the cubic equation has three real roots. Find them by computing θ (d12

= arccos(R/ Q 3 ), in terms of which the three roots are ⎛θ⎞ d x1 = −2 Q cos⎜ ⎟ − 1 ⎝3⎠ 3

(10)

⎛ θ + 2π ⎞ d1 ⎟ − x 2 = −2 Q cos⎜ ⎝ 3 ⎠ 3

(11)

⎛ θ − 2π ⎞ d1 ⎟ − x3 = −2 Q cos⎜ ⎝ 3 ⎠ 3

(12)

2. A STRATEGY FOR REDUCING THE ROUND-OFF ERRORS IN CARDANO’S METHOD WHEN THE CUBIC EQUATION OF STATE POSSESSES THREE REAL ROOTS For the case of a cubic equation, if we assume that the real roots are α, β and γ; then, eq 9 is equivalent to

⎡ otherwise, compute A = −sgn(R)⎣ R +

⎤1/3 R2 − Q 3 ⎦ , where the positive square root is assumed. Next, compute

(x − α)(x − β)(x − γ ) = 0

by expanding eq 16 and comparing with eq 9, the following relations are found:

⎧Q (A ≠ 0) ⎪ B = ⎨A ⎪ ⎩ 0 (A = 0)

α + β + γ = −d1 αβ + αγ + βγ = d 2 αβγ = −d3

in terms of which the three roots are x1 = (A + B) −

d1 3

(16)

The system described by eq 17 can be seen as three nonlinear equations. The solution of the system described by eq 17 is equivalent to obtaining the roots of the polynomial described by eq 9; therefore, a robust, efficient, and reliable approach for solving system (17), together with its solution, will give, as a result, the real roots of the polynomial described by eq 9 without round-off errors. For this purpose, we suggest to cast system (17) as a fixed-point system, u = g(u), with the requirement for convergence being that the eigenvalue of G = [∂gi/∂uj] of the largest modulus must satisfy |λi| < 1.14 The system that fulfills this requirement has the form

(13)

d 1 3 x 2 = − (A + B ) − 1 + i (A − B ) 2 3 2

(14)

d 1 3 x 3 = − (A + B ) − 1 − i (A − B ) 2 3 2

(15)

(17)

the previous three equations are arranged both to minimize round-off error, and also to ensure that no choice of branch for the complex cube root can result in the spurious loss of a distinct root. Disadvantages of using Cardano’s Method when the cubic equation possesses three real roots (eqs 10−12 has been reported by several authors (c.f. Zhi and Lee,9 Gosset et al.,10 Deiters11), because of the fact that the solution implies several referrals to root or trigonometric functions. This can lead to a significant loss of numerical precision, even if double-precision numbers are used throughout a simulation computer program. In the context of process simulation, Eggeman and Chafin12

α = −d1 − (β + γ ) β =

d 2 − αγ α+γ

γ =

− d3 αβ

(18)

System (18) can be solved using the well-known successive substitution method. Going further, a possible improvement is 6973

dx.doi.org/10.1021/ie2023004 | Ind. Eng. Chem. Res. 2012, 51, 6972−6976

Industrial & Engineering Chemistry Research

Research Note

to compute the values of the real roots using the most recently calculated values, since they are supposedly better approximations to the actual solution; that is,

|| g(u) − g(u*)|| ≤ σ || u − u*||

for all u sufficiently close to u*; that is, ∥u − u*∥ < δ for some fixed δ > 0. Remark. The notion of a contraction depends on the underlying choice of matrix norm; it may be more convenient to adopt either the 1 or the ∞ norms rather than the standard Euclidean norm. Theorem 1. If u* = g(u*) is a f ixed point for the system (23) and g is a contraction at u*, then u* is an asymptotically stable f ixed point. Proof. Since ∥uk+1 −u*∥ = ∥g(uk) − g(u*)∥ ≤ σ ∥uk − u*∥, using the assumed estimate described by eq 24. Iterating this basic inequality immediately demonstrates that

α k + 1 = −d1 − (β k + γ k) β k+1 =

d 2 − α k + 1γ k

γ k+1 =

− d3 k+1 k+1

αk+1 + γ k α

β

(19)

where k is the iteration number index. The above system is suggested to be solved using Steffensen’s Method,15 because it gives near-quadratic convergence without evaluating the derivative, and usually a few iterations steps are sufficient.

|| u k − u*|| ≤ σ k || u 0 − u*||

(26)

Here ⎡ 0 ⎢ 2 ⎢ d2 + γ − 2 g′(u) = ⎢⎢ (α + γ ) ⎢ d3 ⎢ α 2β ⎣

and

that form an iterative system (22) 3

3

where k refers to the iteration number and g: →  is a real vector-valued function. A solution is a collection of points uk ∈ 3, in which the index k = 0, 1, 2, 3, ..., takes on non-negative integer values. Definition 1. The fixed point of system (22) is a vector u* ∈ 3, such that

g(u*) = u*

−1 0 d3 αβ 2

⎤ ⎥ d2 + α 2 ⎥ − (α + γ )2 ⎥ ⎥ ⎥ 0 ⎥ ⎦ −1

(27)

and denotes the 3 × 3 Jacobian matrix of the vector-valued function g, whose entries are the partial derivatives of its individual components. Since u* is fixed, then the right-hand side of eq 26 is a linear function of u. Moreover, u* remains a fixed point of the linear approximation. A linear function will converge to the fixed point if and only if its coefficient matrix namely, g′(u*)is a convergent matrix, meaning that its spectral radius ρ(g′(u*)) < 1. This observation motivates the following theorem and corollary. Theorem 2. Let u* be a f ixed point for the system described by eq 23. If the Jacobian matrix norm ∥g′(u*)∥ < 1, then g is a contraction at u*, and, hence, the f ixed point u* is asymptotically stable. Proof. The first-order Taylor expansion of g(u) at the fixed point u* takes the form g(u) = g(u*) + g′(u)(u−u*) + R(u − u*), where the remainder term satisfies limu→u* R(u − u*)/∥u − u*∥ = 0. Let ε > 0 be such that σ = ∥g′(u*)∥ + ε < 1. Choose 0 < δ < 1 such that R(u − u*) ≤ ε ∥u − u*∥ whenever ∥u − u*∥ < δ. For such u, we have, by the Triangle Inequality,

(21)

u k + 1 = g (u k )

(25)

g(u) ≈ g(u*) + g′(u)(u − u*) = u* + g′(u)(u − u*)

(20)

⎛−d1 − (β + γ )⎞ ⎜ ⎟ ⎜ d 2 − αγ ⎟ ⎜ ⎟ g=⎜ α+γ ⎟ ⎜ ⎟ − d3 ⎜ ⎟ αβ ⎝ ⎠

for k = 1, 2, 3 , ...

Since σ < 1, the right-hand side tends to 0 as k→∞, hence, uk→ u*. □ Since g is differentiable, it can be approximated by its firstorder Taylor series:

3. SOME ASPECTS ABOUT THE PROPOSED FIXED-POINT SYSTEM USED TO REDUCE THE ROUND-OFF ERRORS IN CARDANO’S METHOD WHEN THE CUBIC EQUATION OF STATE POSSESSES THREE REAL ROOTS PROVIDING A CONTRACTION MAPPING Let us define two vectors as ⎛α ⎞ ⎜ ⎟ u = ⎜ β⎟ ⎜ ⎟ ⎝γ ⎠

(24)

(23)

We easily see that every fixed point provides a constant solution to the system described by eq 23, namely, uk = u* for all k. Moreover, it is not hard to prove that any convergent solution necessarily converges to a fixed point. Proposition 1. If a solution to system (23) converges, limk→∞uk = u*, then the limit u* is a f ixed point. Proof. This is a simple consequence of the continuity of g. We have u* = limk→∞uk+1 = limk→∞g(uk) = g(limk→∞uk) = g(u*), the last two equalities following from the continuity of g. □ Definition 2. The function g:3 → 3 is a contraction at a point u* ∈ 3 if there exists a constant 0 ≤ σ < 1, such that

|| g(u) − g(u*)|| ≤ || g′(u*)(u − u*)|| + || R(u − u*)|| ≤ (|| g′(u*)|| + ε)|| u − u*|| = σ || u − u*||

which establishes the contraction inequality described by eq 25. □ Corollary 1. If the Jacobian matrix g′(u*) is a convergent matrix, meaning that its spectral radius satisf ies ρ(g′(u*)) < 1, then u* is an asymptotically stable f ixed point. 6974

dx.doi.org/10.1021/ie2023004 | Ind. Eng. Chem. Res. 2012, 51, 6972−6976

Industrial & Engineering Chemistry Research

Research Note

Table 2. Analytical Solution of the Peng−Robinson Equation of State for 1-Pentene Analytical Solution 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Analytical Solution with Iterative Refinement

T [K]

Tr

p [MPa]

ν1[m3/kg-mol]

ν2[m3/kg-mol]

ν3[m3/kg-mol]

ν1[m3/kg-mol]

ν2[m3/kg-mol]

ν3[m3/kg-mol]

107.93 111.43 114.93 118.43 121.93 125.43 128.93 132.43 135.93 139.43 142.93 146.43 149.93 184.93 219.93 254.93

0.232 0.240 0.247 0.255 0.262 0.270 0.277 0.285 0.293 0.300 0.308 0.315 0.323 0.398 0.473 0.549

4.48 × 10−12 1.67 × 10−11 5.70 × 10−11 1.79 × 10−10 5.24 × 10−10 1.43 × 10−9 3.68 × 10−9 8.92 × 10−9 2.06 × 10−8 4.52 × 10−8 9.50 × 10−8 1.92 × 10−7 3.73 × 10−7 5.93 × 10−5 1.48 × 10−3 1.32 × 10−2

2.0031 × 1011 5.5477 × 1010 1.6764 × 1010 5.5010 × 109 1.9347 × 109 7.2928 × 108 2.9130 × 108 1.2344 × 108 5.4863 × 107 2.5648 × 107 1.2509 × 107 6.3410 × 106 3.3420 × 106 2.5927 × 104 1.2342 × 103 1.5952 × 102

1.7132 × 100 1.6429 × 100 6.9567 × 101 1.5154 × 100 1.4575 × 100 4.3607 × 100 2.5331 × 100 2.5295 × 100 2.4349 × 100 2.3442 × 100 2.2684 × 100 2.1896 × 100 2.1153 × 100 1.5392 × 100 1.1597 × 100 8.9563 × 10−1

1.7133 × 100 1.6429 × 100 −6.6413 × 101 1.5154 × 100 1.4575 × 100 −1.5547 × 100 1.7031 × 10−1 7.7030 × 10−2 8.0189 × 10−2 8.4376 × 10−2 7.8130 × 10−2 7.9236 × 10−2 7.9612 × 10−2 8.1672 × 10−2 8.4308 × 10−2 8.7648 × 10−2

2.0031 × 1011 5.5477 × 1010 1.6764 × 1010 5.5010 × 109 1.9347 × 109 7.2928 × 108 2.9130 × 108 1.2344 × 108 5.4863 × 107 2.5648 × 107 1.2509 × 107 6.3410 × 06 3.3420 × 106 2.5927 × 104 1.2342 × 103 1.5952 × 102

3.3490 × 100 3.2081 × 100 3.0763 × 100 2.9527 × 100 2.8368 × 100 2.7277 × 100 2.6249 × 100 2.5279 × 100 2.4363 × 100 2.3496 × 100 2.2674 × 100 2.1894 × 100 2.1154 × 100 1.5392 × 100 1.1597 × 100 8.9563 × 10−1

7.7522 × 10−2 7.7673 × 10−2 7.7827 × 10−2 7.7984 × 10−2 7.8144 × 10−2 7.8308 × 10−2 7.8475 × 10−2 7.8645 × 10−2 7.8819 × 10−2 7.8997 × 10−2 7.9178 × 10−2 7.9363 × 10−2 7.9552 × 10−2 8.1672 × 10−2 8.4308 × 10−2 8.7648 × 10−2

Proof. By definition, g′(u*) is convergent if and only if ρ(g′(u*)) < 1. Choosing ε > 0 such that ρ(g′(u*)) + ε < 1. Any norm that satisfies ρ(g′(u*)) ≤ ∥g′(u*)∥ < ρ(g′(u*)) + ε has the desired property. □ Definition 3. The function g: 3 → 3 is a contraction mapping on a domain Ω ⊂ 3 if (a) it maps Ω to itself, so g(u) ∈ Ω whenever u ∈ Ω, and (b) there exists a constant 0 ≤ σ < 1 such that ∥g(u) − g(v)∥ ≤ σ∥u − v∥ for all u,v ∈ Ω. In other words, applying a contraction mapping reduces the mutual distance between points. Therefore, as its name indicates, a contraction mapping effectively shrinks the size of its domain. Theorem 3. If g is a contraction mapping on a closed bounded domain Ω ⊂ 3, then g admits a unique fixed point u* ∈ Ω. Moreover, starting with any initial point u0 ∈ Ω, the iterates uk+1=g(uk) necessarily converge to the f ixed point: uk→u*. Proof. The proof is evident. □ In particular, if ρ (g′(u)) < 1 for all u ∈ Ω, then the conclusions of the contraction mapping Theorem 3 hold. The Jacobian described by eq 27 has three real eigenvalues, namely,

4. COMMENTS ON THE PROPOSED APPROACH TO FACE THE ROUND-OFF ERRORS IN CARDANO’S METHOD WHEN THE CUBIC EQUATION OF STATE POSSESSES THREE REAL ROOTS The real advantage of the proposed approach to face the roundoff errors in Cardano’s Method due to several referrals to root or trigonometric functions when the cubic equation of state possesses three real roots can be summarized as follows: (i) The proposed approach results in an inexpensive iterative refinement procedure to face the round-off errors due to several referrals to root or trigonometric functions in the calculation of three real roots of the cubic equations of state using Cardano’s Method. (ii) Since system (18) is cast to obey the condition of having a contraction mapping, its convergence is guaranteed, starting from the solution obtained in Cardano’s Method, because it always ensures a direct-forward step toward the solution. (iii) The solution to system (19) can be seen as another step in Cardano’s Method; because of the use of Steffensen’s method,15 convergence is obtained as much in three iterations, when round-off errors exist, because it gives near-quadratic convergence without evaluating the derivative. (iv) The principal advantage in the suggested approach is that the point at which it is applied to test the round-off error obtained by Cardano’s Method in the solution of the three real roots appears to be the first step of the solution of system (19). In the case that round-off errors do not exist, the solution of the system described by eq 19 may be considered as converged and the calculation is fulfilled. (v) The implementation of the proposed approach is straightforward to existing chemical engineering simulation codes.

⎧ 0 ⎪ 1/2 ⎪ ⎛ d3 d 2 + α 2 ⎞ ⎟ ⎜ − · ⎪ 2 2 λ = ⎨ ⎝ αβ (α + γ ) ⎠ ⎪ 1/2 ⎪ ⎛ d3 d 2 + α 2 ⎞ ⎟ −⎜− · ⎪ ⎪ ⎝ αβ 2 (α + γ )2 ⎠ ⎩

therefore the spectral radius is ⎛ d d + α 2 ⎞1/2 ⎟ ρ(g′(u)) = ⎜ − 32 · 2 2 ⎝ αβ (α + γ ) ⎠

(28)

5. NUMERICAL EXPERIMENTATION In this part, the case for 1-pentene presented by Zhi and Lee,9 using the Peng−Robinson equation of state,3 was used to show the advantages of the iterative refinement suggested approach. As in the work of Zhi and Lee, it is not the accuracy of the

If α ≤ β ≤ γ, and given the definitions of d2 and d3, the spectral radius described by eq 28 is always less than unity. The smaller the spectral radius of the Jacobian matrix at the fixed point, the faster the nearby iterates will converge to it. 6975

dx.doi.org/10.1021/ie2023004 | Ind. Eng. Chem. Res. 2012, 51, 6972−6976

Industrial & Engineering Chemistry Research

Research Note

referrals to root or trigonometric functions when the cubic equation of state possesses three real roots. An iterative refinement approach to face these round-off errors was presented. 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.

model with which we are concerned, but rather the validity of the compressibility factor calculation. For this purpose, f(T) in eq 4 is used as f (T ) = ⎡⎣1 + m(1 −

2

Tr )⎤⎦

where m is given by Robinson et al.,16 as m = 0.3796 + 1.485ω − 0.1644ω2 + 0.01667ω3 (ω denotes the acentric factor). The critical properties needed and the acentric factor were taken from the compilation of Reid et al.17 The roots of the cubic equation associated were calculated by Cardano’s Method without and with the iterative refinement procedure described, and the results are listed in Table 2. From Table 2, it is easy to notice that, from point 1 to point 10, the solution of the equation of state with Cardano’s Method gave unreasonable results (for example, the liquid volume of point 10 is larger than that in point 11). Moreover, the solutions obtained gave negative liquid volumes in points 3 and 6. After applying the proposed iterative refinement procedure, note that the results are reasonable. Each root obtained after application of the proposed iterative refinement procedure was substituted back into its original cubic equation and the resulting residuals, as expected, are better than those obtained without the application of the proposed iterative refinement procedure that is without round-off errors. Table 3 shows the spectral radius for the proposed iterative refinement procedure, at the initial and final steps, respectively,

■

*Tel.: +525 5 5581 4982. Fax: +525 5 5581 4982. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

Spectral Radius ρ (g′(u)) initial

final

number of iterative refinement steps

0.2974 0.3038 0.0070 0.3167 0.3230 0.1060 0.1792 0.1763 0.1800 0.1834 0.1868 0.1904 0.1939 0.2304 0.2697 0.3136

0.1521 0.1556 0.1591 0.1625 0.1660 0.1694 0.1729 0.1764 0.1799 0.1834 0.1869 0.1904 0.1939 0.2304 0.2697 0.3136

3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0

REFERENCES

(1) Mathias, P. M.; Klotz, H. C. Take a Closer Look at Thermodynamic Property Models. Chem. Eng. Progress 1994, 6, 67− 75. (2) Soave, G. Equilibrium Constants from a Modified Redlich− Kwong Equation of State. Chem. Eng. Sci. 1972, 27 (6), 1197−1203. (3) Peng, D.-Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59−64. (4) Orbey, H.; Sandler, S. I. Modeling Vapor−Liquid Equilibria: Cubic Equations of State and Their Mixing Rules; Cambridge University Press: Cambridge, U.K., 1998. (5) Kamath, R. S.; Biegler, L. T.; Grossmann, I. E. An EquationOriented Approach for Handling Thermodynamics Based on Cubic Equation of State in Process Optimization. Comput. Chem. Eng. 2010, 34 (12), 2085−2096. (6) Tester, J. W.; Modell, M. Thermodynamics and Its Applications, 3rd ed.; Prentice−Hall PTR: Upper Saddle River, NJ, 1996. (7) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice−Hall: Upper Saddle River, NJ, 1998. (8) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: Cambridge, U.K., 2007. (9) Zhi, Y.; Lee, H. Fallibility of Analytic Roots of Cubic Equations of State in Low Temperature Region. Fluid Phase Equilib. 2002, 201 (2), 287−294. (10) Gosset, R.; Heyen, G.; Kalitventzeff, B. An Efficient Algorithm to Solve Cubic Equations of State. Fluid Phase Equilib. 1986, 25 (1), 51−64. (11) Deiters, U. K. Calculation of Densities from Cubic Equations of State. AIChE J. 2002, 48 (4), 882−886. (12) Eggeman, T.; Chafin, S. Beware the Pitfalls of CO2 Freezing Prediction. Chem. Eng. Progress 2005, 101, 39−44. (13) Salim, P. H. Comment on the paper of Zhi and Lee entitled “Fallibility of Analytic Roots of Cubic Equations of State in Low Temperature Region”. Fluid Phase Equilib. 2006, 240 (2), 224−226. (14) Ortega, J. M.; Rheinboldt, W. C. Iterative Solution of Nonlinear Equations in Several Variables; Society for Industrial Mathematics: Philadelphia, PA, 1987. (15) Burden, R. L.; Faires, J. D. Numerical Analysis, 9th ed.; Brooks/ Cole: Boston, MA, 2010. (16) Robinson, D. B.; Peng, D.-Y.; Chung, S. Y.-K. The Development of the Peng−Robinson Equation and Its Application to Phase Equilibrium in a System Containing Methanol. Fluid Phase Equilib. 1985, 24 (1−2), 25−41. (17) Reid, R.; Poling, B. E.; Prausnitz, J. M. The Properties of Gases and Liquids, 4th ed.; McGraw−Hill: New York, 1987.

Table 3. Behavior of the Proposed Iterative Refinement Procedure

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

AUTHOR INFORMATION

Corresponding Author

when an accuracy of 15 significant digits is used. Table 3 also shows the number of iterative refinement steps to fulfill this requirement using the Steffensen’s method. It is noteworthy that the spectral radius in all cases is far from 1, that is, a value close to 0.1; consequently, a small number of iterations is required to reach the final solution.

6. CONCLUSIONS With many cubic equations of state, widely used in chemical engineering process simulation, the calculation of molar volumes or compressibility factors, for a fixed pressure, temperature, and composition, using Cardano’s Method, can produce physically unreasonable values, because of several 6976

dx.doi.org/10.1021/ie2023004 | Ind. Eng. Chem. Res. 2012, 51, 6972−6976