Ind. Eng. Chem. Res. 2006, 45, 9161-9170
9161
An Algorithm to Calculate K- and L-Points Samir H. Mushrif† and Aaron V. Phoenix* Department of Chemical Engineering, UniVersity of Saskatchewan, 57 Campus DriVe, Saskatoon, Saskatchewan, S7N 5A9, Canada
An algorithm to directly locate K- or L-points (which are often called critical endpoints) is proposed that utilizes a robust critical point solver and a standard phase stability test within a nested-loop structure. Successive substitution (SS), accelerated successive substitution, and Newton-Raphson (NR) procedures were used to perform the phase stability test, and it was determined that a hybrid between the SS and the NR methods balances the robustness of the SS procedure against the convergence rate of the NR procedure. The algorithm was tested using the Peng-Robinson equation of state for a variety of binary systems, including an ethane + ethanol system, three methane + n-alkane systems, and three propane + polyaromatic systems. An analysis of convergence difficulties is discussed, relative to the Gibbs free-energy surface topology. It was determined that an accurate initial critical composition was necessary to avoid convergence to a trivial solution. Convergence was not sensitive to the initial guess used to locate the noncritical equilibrium phase. Introduction
N
When a liquid phase and a vapor phase become critical in the presence of a noncritical, heavier equilibrium liquid phase, it is called a K-point; when two liquid phases become critical in the presence of a noncritical equilibrium vapor phase, it is called an L-point.1 Conventionally, the K- and L-points are also called upper critical endpoints (UCEPs) and lower critical end points (LCEPs). However, there is a minor difference in the terminology: the criterion for distinguishing an LCEP from a UCEP is not the type of critical point (liquid-liquid or liquidvapor) but its location in P-T space. An UCEP is always located at a higher temperature and pressure than an LCEP, and when there is only one critical endpoint present, it is referred to as an UCEP. Hence, a K-point will always be a UCEP but an L-point can be a UCEP or an LCEP, depending on the global phase behavior of the system. Within the framework of the Scott and van Konynenburg2,3 nomenclature for binary systems, K- and L-points exist in Type II and higher systems. These K- and L-points can be thought of as landmarks within the system phase space, and experimental determination of the points sometimes requires advanced facilities and equipment.4 To mitigate experimental costs or to examine complex phase behavior of model systems pertinent to industry, K- and L-points are often calculated from thermodynamic models. Conversely, these landmarks may be used to tune the equation-of-state parameters, so that thermodynamic models perform better in regions of complex phase behavior. The criteria for critical points have been well-documented. Reid and Beegle5 and Peng and Robinson6 provide a review of the critical point criteria used in calculations prior to 1977. These techniques involved expensive computations of many determinants and their derivatives. Heidemann and Khalil7 proposed a method of calculating critical points of mixtures using much less costly critical criteria. Their development was based upon a Taylor series expansion of the Helmholtz free energy about the critical point at constant temperature and volume: * To whom correspondence should be addressed. Tel: +1-306-9664190. Fax: +1-306-966-4777. E-mail:
[email protected]. † Current address: Department of Chemical Engineering, McGill University, Montreal, PQ, Canada.
[A - A0 1
∑ i)1
µi0∆ni]T0,V0 )
N
N
N
∑∑ ∑
(
1
N
N
∑∑
( ) ∂2A
2! i)1 j)1 ∂ni∂nj
∂3A
3! i)1 j)1 k)1 ∂ni∂nj∂nk
)
∆ni∆nj +
T0,V0
∆ni∆nj∆nk + O(∆n4) (1)
T0,V0
Briefly, the two equations that defined their critical point method were
det(Q) ) 0 and
C)
(2)
( ) ∑ ∑∑ ( ) V-b 2b
2 N
N
N
2
nT
k)1 j)1 i)1
∂2 ln(fi) ∂nk∂nj
∆ni∆nj∆nk ) 0 (3)
where the elements of Q, qij, were defined as
qij )
( )( ) ∂ ln(fi) T nT 100 ∂nj
(4)
The elements qij defined in eq 4 were scaled coefficients of the quadratic form of eq 1, whereas eq 3 represents a scaled version of the cubic form of eq 1. The ∆ni values within eq 3 were determined by solving the linear system
Q∆n ) 0
(5)
where ∆n ) (∆n1, ..., ∆nN)T and ∆nN was arbitrarily set equal to 1. Equations 2 and 3 were solved using two nested one-dimensional searches and did not require user supplied initiation points. The inner loop solved eq 2 for T at a fixed volume, whereas the outer loop used eq 3 to update the volume. It was shown to be reasonably robust for typical critical-point calculations. Michelsen8 modified Heideman and Khalil’s method of solving eqs 2 and 3 by recognizing that the cubic form of the Taylor series expansion that gave eq 3 could be expressed in terms of a double summation, using the numerical derivatives of Q in the direction of ∆n. Michesen’s alternate form of eq 3 was
10.1021/ie060359h CCC: $33.50 © 2006 American Chemical Society Published on Web 11/14/2006
9162
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006
C)
∑i ∑j
T ∆ni∆njQ* ij ) ∆n Q*∆n ) 0
(6)
Q* )
∂ (Q(n + s∆n,T,V))s)0 ∂s
(7)
where
The variable s in this development represents a “distance” from the critical point. Michelsen then used his phase stability criteria9 as the basis for another alternative formulation of the critical-point criteria.10 His criteria were alternative versions of eqs 2 and 6 that were expressed in terms of “deviation” variables from the critical point. The direction of the deviation was defined as the eigenvector of the matrix used in lieu of Q in eq 2. Michelsen’s development used temperature and pressure as independent variables and the values calculated during the solution procedure were shown to give important information about the free-energy surface topography around the critical point. Michelsen discussed the topography of a three-phase boundary around the critical point,10 or, in other words, around a K- or L-point. A more-detailed development and analysis of Michelsen’s method for calculating critical points followed in a later paper.11 This paper redeveloped the criteria expressing the Taylor series expansion of the distance from the Gibbs free energy to a plane tangent to the free-energy surface at the critical point, in terms of a scalar distance parameter s and a direction vector u. The paper emphasized that a good initial guess had to be provided for the method to work. Michelsen and Heidemann12 extended Michelsen’s earlier work11 to calculate tricritical points by expanding the Taylor series expansion to more terms. Two additional criteria more than those for a critical point were needed to define a tricritical point. For the purposes of this paper, the specific criteria are not needed, but Michelsen and Heidemann’s nested-loop solution procedure involved solving for a regular critical-point temperature and pressure (or volume) at a fixed composition using Michelsen’s method11 (Newton-Raphson technique), whereas an outer loop updated two compositions using an additional Newton-Raphson technique with finite difference derivatives. Alternatively, Michelsen and Heidemann also proposed solving for the tricritical temperature, pressure, and composition using a full, four-variable Newton-Raphson technique. Continuing the work in locating tricritical points, Kohse and Heidemann13 modified Michelsen and Heidemann’s iterative method. They utilized Heidemann and Khalil’s critical-point algorithm to locate the critical-point temperature and volume in two inner loops and updated the tricritical compositions using a two-dimensional Newton-Raphson method in an outer loop. An alternative three-dimensional search technique in the outer loop was attempted but is not relevant to the present work. Kohse and Heidemann found that their two-dimensional search algorithm was better able to address poor initial composition guesses than their three-dimensional search algorithm. More recently, Gauter et al.4 calculated critical endpoints for ternary systems using the approach taken by Heidemann and Khalil7 to determine a critical point and then tracking the critical locus in search of an equilibrium phase. Unlike Gauter et al.’s work,4 the purpose of this work was to develop and analyze an algorithm to directly calculate Kand L-points without following the critical loci. The algorithm was developed for multicomponent systems, but the analysis and examples in this paper will be limited to binary and ternary
Figure 1. Flowsheet of the K- and L-point algorithm.
systems. The algorithm can use any equation of state that relates the pressure, temperature, and volume of fluid mixtures; however, for the purposes of this study, only the Peng-Robinson equation of state was used with standard van der Waals mixing rules.14 The structure of the algorithm is based on the nestedloop structure used by Kohse and Heidemann to find tricritical points:13 two inner loops use the methods of Heidemann and Khalil7 and Michelsen11 to find a critical-point temperature and volume, respectively, at a fixed composition; an outer loop modifies the critical-point composition using a one-dimensional search for an equilibrium phase. The Algorithm The algorithm used to calculate K- and L-points consists of two nested inner loops to calculate a critical-point temperature and volume at a fixed composition z. An outer loop uses the critical point as a test phase, searches for an incipient phase at a trial composition nˆ , and updates the critical composition to iteratively decrease the tangent plane distance of the incipient phase to zero. Both the inner and outer loops use a NewtonRaphson method with numerical derivatives. A simplified flow sheet is shown in Figure 1. Following the developments of Michelsen and Heidemann,12 the critical-point criteria solved within the inner loops (Figures 1a and 1b) of the code are based on the tangent-plane criteria. The condition for the stability of the mixture, with respect to a trial phase, nˆ , at a specified temperature and pressure (T0, P0), is
F(nˆ ) ) 1+
∑i nˆ i(ln nˆ i + ln φi(nˆ ) - ln zi - ln φi(z) - 1) g 0
(8)
where F is a modified “tangent plane distance” from the Gibbs
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9163
free-energy surface to the hyperplane tangent to the Gibbs freeenergy surface at the critical-point composition, z. The parameter φi represents the fugacity coefficient of component i, calculated with an equation of state. Introducing a scaled, translated mole number, Xi ) (nˆ i - zi)/ zi1/2, the first and second partial derivatives of F, with respect to Xi, are
( )
∂F gˆ i ) ) zi1/2(ln nˆ i + ln φi(nˆ ) - ln zi - ln φi(z)) (9) ∂Xi and
Bij )
( )
( )
∂(ln φi) ∂2F ) σij + (zizj)1/2 ∂Xi∂Xj ∂nˆ j
σij ) 1
(for i ) j)
(11a)
)0
(for i * j)
(11b)
gˆi is an element of the vector gˆ , and Bij is an element of the matrix B. The matrix B is a modified form of the matrix Q given in eq 2. A distance parameter, s, is introduced such that
(12)
where u is a vector with property ) 1. Michelsen and Heidemann12 chose u to be an eigenvector corresponding to the minimum eigenvalue of the singular matrix B. As in their work, in this study, the vector u is calculated from the normalized solution of uTu
Bu ) 0
(13)
after setting uN ) 1. Considering the dependence of F on s, the Taylor series expansion of F about the critical phase can be written as
F ) b′s2 + c′s3 + d′s4 + ...
(14)
where the coefficients b′, c′, etc. are related to the derivatives of F, with respect to s, and evaluated at s ) 0. At a critical point, b′ ) c′ ) 0. For the critical point to be stable, d′ g 0. Details of the derivation can be found in Michelsen and Heidemann’s work.12 The criterion of b′ ) 0 is met when the matrix B is singular and is used to find the critical temperature at a fixed composition and volume. Using the scaling factor proposed by Heidemann and Khalil,7 the equation used to find the critical temperature is given by
Q(T;V,z) )
(100T ) det(B) ) 0
X ) su + s2w + ...
(15)
where the determinant of B can be calculated by inverse iteration.11 Equation 15 is solved using a Newton-Raphson procedure with a numerical derivative (see loop (a) in Figure 1). Numerical differentiation was achieved with a central differencing scheme, using a perturbation in temperature of dT ) 10-5 T. Convergence was defined as a relative change in absolute temperature of 0, the trial phase is an incipient phase, and if θ < 0, the reference phase is unstable. Equations 24 and 25 can be combined to give N working equations used to solve for a stationary point:
gi ) Yi - exp{ln(fcrit i ) - ln(φi(y)) - ln P} ) 0 (for 1 e i e N) (26) Graphically, a stationary point that is given by a solution to eq 24 or eq 26 would be identified by a tangent to the free-energy surface at the trial phase parallel to the tangent at the reference phase. The search for an equilibrium phase begins with an initial guess for the equilibrium phase composition, and eq 26 is solved iteratively using an accelerated successive substitution method15 or a Newton-Raphson method with analytical derivatives (loop (c) in Figure 1). In the case of a successive substitution method, eq 26 is rewritten as
Yi )exp{ln(f crit i )-ln(φi(y))- ln P}
(for 1 e i e N) (27)
A K- or L-point has been found if
θ)0
(28)
If θ, as calculated by eq 25, does not meet the convergence criterion (θ2 < 10-12), the critical-phase composition is updated using a Newton-Raphson method:
) z(k) z(k+1) 1 1 -
θ(k) θ′(k)
(29)
where θ′ is approximated by perturbing the critical composition, recalculating θ, and calculating the finite difference analogue of the derivative. Discussion There is no a priori guarantee that a K- or L-point exists in a multicomponent system; however, if a K- or L-point exists, it will have (N - 2) degrees of freedom16 in a nonreacting system. Therefore, K- and L-points have zero degrees of freedom in a binary system and no system variables can be used as parameters in the problem. For nonreacting systems with more than two components, (N - 2) system variables or relationships must be specified to close the set of equations. This study did not examine which variables were best to specify for multicomponent systems.
Figure 2. Gibbs energy plot for the ethane + ethanol system at 315 K and 5.34 MPa: (9) critical phase and (2) equilibrium phase.
Initial compositions, temperatures, and volumes for the critical phase and equilibrium phase were required within the algorithm, along with the pure component data and the equation-of-state binary interaction parameters. Initial guesses for critical temperatures and volumes were the same as those used by Heidemann and Khalil:7
T(0) ) 1.5 V(0) )
{
∑i ziTCi
3.85b (K-point) 1.75b (L-point)
(30a) (30b)
where b is the equation-of-state co-volume parameter. The success of the algorithm to locate a K- or L-point was largely dependent on two factors: the binary interaction parameters (δij) and the initial critical composition (z(0)). The existence of a mathematical K- or L-point was sensitive to the choice of binary interaction parameters, and when a K- or L-point did exist, a good choice of initial critical composition was paramount to avoiding trivial solutions. The importance of binary interaction parameters can be illustrated by considering the dimensionless Gibbs free energy of mixing surface for the ethane + ethanol system. Figure 2 shows two free-energy surfaces calculated at the experimental K-point temperature,17 314.66 K, and a pressure of 5.34 MPa. When δij ) 0.04, the free-energy surface is concave upward over the entire composition range, making the existence of a K- or L-point impossible. If δij ) 0.135, the free-energy surface is concave down over the midrange of compositions, allowing for a possible K- or L-point. Indeed, the calculated K-point is graphically shown by the tangent plane joining the ethane-rich critical phase (represented by a solid square (9) in Figure 2) and the ethanol-rich equilibrium phase (represented by a solid triangle (2) in Figure 2). Similarly, for ternary systems, Figures 3 and 4 show the freeenergy surface contours for a ternary system of CO2 (1) + 1-pentanol (2) + tridecane (3) at one of the experimental K-points of the system (317.43 K and 9.22 MPa).4 Figure 3 shows the contour plot with all binary interaction parameters set to zero. The free-energy surface is concave upward
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9165
Figure 5. Free-energy surfaces of n-pentane (1) + anthracene (2) binary mixtures: (a) zn-pentane ) 0.99 and (b) zn-pentane ) 0.975. Legend of symbols: (9) critical phase and (2) stationary point. Legend of lines: (s) critical-phase tangent and (- - -) stationary-point tangent. Figure 3. Gibbs energy contour plot for the ternary system of CO2 (1) + 1-pentanol (2) + tridecane (3) at 317.43 K and 9.22 MPa, which is an experimental K-point for the system.4 All δij ) 0.0.
Figure 4. Gibbs energy contour plot for the ternary system of CO2 (1) + 1-pentanol (2) + tridecane (3) at 317.43 K and 9.22 MPa, which is an experimental K-point for the system.4 δ12 ) 0.17, δ13 ) 0.15, and δ23 ) 0.10.
throughout the composition range, indicating that no K- or L-point is present. The surface contour in Figure 4 was created with δ12 ) 0.17, δ13 ) 0.15, and δ23 ) 0.10, and it clearly shows a possible phase split if the feed composition lies near x2 ) 0.7. Cleary, it is necessary for the Gibbs free-energy surface to be concave downward over a portion of the composition range if a K- or L-point can exist. In other words, there must be some composition for which a single phase (critical or not) is unstable. If this property does not exist for a given interaction parameter, the equilibrium phase search in the algorithm will converge to a trivial solution, i.e., the critical composition. An exploration
of the effects of δij on the overall P-T phase diagram and on the location of K- and L-points has been given in a related paper.18 The convergence issues associated with locating a critical point at a given composition have been discussed by Heidemann and Khalil7 and will not be discussed in this paper. However, the initial critical composition, z(0), was observed to affect the successful convergence of the proposed algorithm to a K- or L-point significantly. An evaluation of some experimental critical endpoints of typical hydrocarbon binary systems indicated that the K- and L-points were closer to the critical point of the more-volatile component along the critical-point loci in the P-T space. Although specific to the hydrocarbon systems that have been examined, and likely to be dependent on the relative volatilities of the mixtures studied, this methodology directed the initial critical compositions to be chosen with a volatile component mole fraction of 0.95 < z1 < 1.0. Even with this range of initial critical compositions as a guide, the topography of the free-energy surface at the critical point often showed characteristics resulting in convergence to a trivial solution where the equilibrium phase was identical to the critical phase. Analysis of the free-energy surface indicated that, although a mathematical K- or L-point existed, and critical points could be located, the domain of convergence, with respect to the critical-point composition, was limited. This is illustrated by the free-energy surface of the n-pentane (1) + anthracene (2) binary shown in Figure 5, which was generated with δij ) 0.1. With the critical composition set to z1 ) 0.99, as shown in Figure 5a, no nontrivial solution to eq 24 can be determined and the algorithm fails to find a K- or L-point. A critical composition with z1 ) 0.975 gives the free-energy surface shown in Figure 5b. A stationary point exists under these conditions, as indicated by the dashed line drawn parallel to the tangent drawn at the critical composition. Equation 24 will locate this point, giving a negative value for θ and providing a nontrivial solution with which the critical composition can be updated. Thus, it can be seen that failure of the algorithm can occur when a critical composition is updated to a value where no stationary points exist other than the trivial solution. This frequency of this occurring could be reduced, but not eliminated, by limiting the percentage change of the critical composition on each
9166
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006
Figure 6. Free-energy surfaces of a n-pentane (1) + anthracene (2) binary mixtures with zn-pentane ) 0.95. Legend of symbols: (9) critical phase and (2) stationary point. Legend of lines: (s) critical-phase tangent and (- - -) stationary-point tangent.
Figure 7. Iteration counts for four methods of solving for a stationary point for the n-pentane (1) + anthracene (2) system at a critical composition of z1 ) 0.95. Solid symbols indicate that the stationary point was foun; open symbols indicate that the trivial solution was found.
iteration of the outer loop. No standard method of avoiding this scenario during initiation of the algorithm was determined; however, caution is suggested if initial critical compositions near to the composition boundary will be used. Equation 24 required an initiation protocol to solve for a stationary point iteratively. Initiation was of particular importance during the first iteration of the outer loop, because, in subsequent outer-loop iterations, eq 24 quickly converged using the previous solution. Because the critical phase was predominantly rich in the more-volatile component, an initial trial composition rich in the less-volatile component worked for most systems. However, if the solution of eq 24 was heading toward a trivial solution, computations were terminated and new guess for the trial composition was chosen manually. Convergence toward a trivial solution was determined by examining the determinant of the matrix
Figure 7 shows the iteration counts required by each of the four methods to converge to a solution, using 0 < y (0) 1 < 1 as the initial guess. Convergence to the stationary point or the trivial solution is indicated by the solid symbols or open symbols, respectively. Convergence to a stationary point was defined by
[ ]
∂g1 ‚‚‚ ‚‚‚ ∂Y1 • l ‚• M) • l ‚• ∂gn ‚‚‚ ‚‚‚ ∂Y1
∂g1 ∂Yn l l ∂gn ∂Yn
(31)
M becomes singular at the trivial solution, the critical composition, and its determinant approaches zero. Using the n-pentane (1) + anthracene (2) system as a test case, the convergence behavior of the search for a stationary point was explored. Figure 6 shows the free-energy surface, the critical point, and the stationary point where the critical composition was fixed at z1 ) 0.95. The stationary point in this system occurs at a composition of y1 ) 0.361, corresponding to θ ) -0.08. Four methods to solve eq 24 for the stationary point were analyzed: a Newton-Raphson method (eq 26), a successive substitution method (eq 27), and two accelerated successive substitution methods that utilize the method 2 and method 3 techniques documented by Mehra et al.15
N
∑i gi2 N
< 10-12
(32)
All methods converged to the stationary point when the procedure was initiated with a composition of y1 < 0.58. As expected, the Newton-Raphson method showed the smallest domain of convergence to the stationary point while also showing the highest rate of convergence with iteration counts of 10-4 were attempted, convergence to the stationary point in the region y1 > 0.6 was not observed as the composition at the time of switching was outside the domain of convergence for the Newton-Raphson method shown in Figure 7. The remainder of the discussion is based on a combined successive substitution/Newton-Raphson method with a switching tolerance of 10-6. The successive substitution was used only for the first iteration of the outer loop. Because the noncritical equilibrium phase in a K-point is a liquid and in an L-point is a vapor, different initial compositions were used to start a search for each of these points. It was determined that using a value of n1 ) 0.001 for an L-point calculation worked well. To find a K-point, it was determined that using an initial stationary point at the opposite end of the composition spectrum was required (n1 ) 0.999). Although convergence to an L-point stationary point was relatively insensitive to the initial guess for n1, convergence to the nontrivial solution for a K-point was sensitive to the initial guess.
Figure 9. Iteration counts for the combination of successive substitution and Newton-Raphson methods of solving for a stationary point for the n-pentane (1) + anthracene (2) system at a critical composition of z1 ) 0.95. Solid symbols indicate that the stationary point was found; open symbols indicate that the trivial solution was found.
Figure 10. Gibbs energy plot for ethane (1) + ethanol (2) system at the calculated L-point: (9) critical phase and (2) equilibrium phase.
In the cases studied, the topography of the free-energy surface was such that the features of the surface corresponding to the vapor phase covered only a small portion of the composition space. Figure 10 illustrates this point for the ethane (1) + ethanol (2) system, where the L-point vapor phase has an ethane mole fraction of 0.989. Along the free-energy surface, the ethanerich vapor phase exists for ethane compositions of >0.96. Using an initial guess for the stationary point less than this value quickly results in the trivial solution. This sensitivity was more pronounced for the L-point calculations that involve mixtures of light alkanes with heavier alkanes or aromatics. After a stationary point is found, the critical-phase composition is updated in the outer loop (d) of Figure 1 until θ2 < 10-12.
9168
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006
Table 1. Calculated K- and L-Points Calculated point
z (0) 1
δij
iterations
temperature, T (K)
pressure, P (MPa)
Experimental z1
n1
K L
0.135 0.036
0.950 0.680
9 7
314.8 308.4
Ethane (1) + Ethanol (2) 5.341 0.9727 0.1894 4.277 0.6856 0.9893
K L
0 0
0.990 0.960
6 4
195.8 194.5
Methane (1) + n-Pentane (2) 5.212 0.9935 0.9149 4.994 0.9575 0.9983
temperature, T (K)
pressure, P (MPa)
314.7 309.0
5.556 4.98
z1
n1
reference 17 17
K L
-0.01 0.014
0.998 0.935
5 10
192.7 181.8
Methane (1) + n-Hexane (2) 4.862 0.9986 0.8769 3.415 0.9311 0.99998
195.9 182.5
5.206 3.415
K
-0.02
0.999
11
191.4
Methane (1) + n-Heptane (2) 4.719 0.9996 0.8647
191.7
4.785
20
374.0
Propane (1) + Phenanthrene (2) 4.501 0.9983 0.8758
377.3
4.67
21
375.6
Propane (1) + Fluorene (2) 4.583 0.9968 0.8059
385.5
5.11
21
378.8
4.76
21
K K K
-0.03
0.998
-0.07
0.997
0
0.997
4 3 8
Propane (1) + Triphenylmethane (2) 4.520 0.9983 0.8989
374.3
Table 2. Pure-Component Critical Properties and Acentric Factorsa
component
critical temperature, TC (K)
critical pressure, PC (MPa)
ω
methane ethane ethanol propane 1-pentanol n-pentane n-hexane n-heptane n-tridecane anthracene phenanthrene fluorene triphenylmethane carbon dioxide
190.4 305.4 513.9 369.8 588.2 469.7 507.5 540.3 676 873 869.3 870 865 304.1
4.6 4.88 6.14 4.25 3.91 3.37 3.01 2.74 1.72 2.9 2.9 4.7 2.2 7.38
0.0109 0.0979 0.643 0.1518 0.5784 0.2522 0.299 0.3494 0.6203 0.489 0.495 0.3493 0.5735 0.239
a
Taken from ref 22.
The updating procedure modifies one of the zi using a NewtonRaphson method with numerical derivatives. This work and its related publication were focused on binary systems; therefore, the choice of which critical composition variable to change was not investigated. Instead, the mole fraction of the most-volatile component, z1, was treated as the independent variable in the outer loop. Algorithm failure due to the outer loop was most often due to the critical-phase composition changing to a value where no stationary point existed, as discussed previously and shown in Figure 5a. In such a case, the algorithm results in a trivial solution. To help mitigate this possibility, changes in the independent zi were limited to 10% or less of the distance to the composition boundary, zi ) 1. Although the aforementioned discussion of the obstacles associated with the proposed algorithm has focused on the n-pentane + anthracene system, the issues are representative of those observed in all the hydrocarbon systems tested. Table 1 summarizes the K- and L-points found for various hydrocarbon
0.998 0.929
0.7177 0.99995
19 19
binary systems using a binary interaction parameter that best reproduced the experimental K- or L-point. The existence of a K- or L-point was highly dependent on the binary interaction parameter used, and the listed parameter was that which gave the nearest match to the experimental K- or L-point. The effects of the binary interaction parameter on calculated K- and L-points and on the global P-T phase diagram are discussed by Mushrif and Phoenix.18 Table 1 also lists the initial critical compositions used and the outer-loop iteration counts. The pure component properties used for the generalized Peng-Robinson equation of state can be found in Table 2. The generic initiation for the stationary-point loop was effective for locating all the stable K- and L-points. However, no generic means of determining an initial composition for the critical-point calculation was observed. In the methane + n-alkane series, it was observed that, as the disparity in molecular size increased, the initial critical-point composition required for convergence became richer in the volatile component. Convergence to the K-point required a precision of three significant figures in the initial composition. A similar sensitivity to the initial composition was observed for the propane + polyaromatic systems. In all cases, this sensitivity was due to the small range of critical compositions for which a nontrivial stationary point existed, as illustrated in Figure 5. The outerloop iteration counts for all test calculations were reasonably low, because of the precision required for the initial critical composition guess. Convergence failure also occurred in situations when the critical composition lead to a negative pressure, as evaluated from the Peng-Robinson equation of state. Because the equations used to locate the stationary point were written in terms of the logarithm of fugacities, these negative pressures were not tolerated by the computational procedure. Although a critical point at negative pressures is unrealistic, it should be possible to recast the stability equations using temperature and
Table 3. Unstable K- and L-Points system
point
δij
z1(0)
n1(0)
iterations
T (K)
P (MPa)
z1
n1
ethane (1) + ethanol (2) methane (1) + n-pentane (2) methane (1) + n-hexane (2) propane (1) + phenanthrene (2) propane (1) + triphenylmethane (2)
L K L K K
0.036 0 0 -0.01 -0.03
0.800 0.980 0.970 0.980 0.980
0.100 0.001 0.990 0.001 0.001
90 15 36 5 8
323.5 192.8 196.2 177.9 350.0
4.699 4.604 5.266 1.221 1.34
0.7936 0.9724 0.9883 0.9778 0.977
0.7352 0.9549 0.9939 0.9648 0.964
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9169
volume as independent variables and improve the robustness of the algorithm. This option is currently being investigated. By changing the initial guesses for the critical phase and for the stationary point calculation, numerous unstable K- and L-points were located, as determined by d < 0. These points are listed in Table 3. Note that the initial critical composition used to find these points did not differ significantly from the initial compositions used to find the stable points. In the case of the ethane + ethanol and the methane + n-pentane systems, a different initial stationary-point composition was required. Convergence to these unstable points with initial guesses similar to those required to find the stable points illustrates the importance of checking the stability of a converged solution. The unstable K-points reported for the propane + polyaromatic systems could only be determined for compositions within three significant digits. These K-points were confirmed by varying the critical- and stationary-point initial compositions; however, numerical uncertainties along with the “flat” nature of the freeenergy surface around the solution limited the final precision of the results. Conclusions An algorithm was proposed to calculate the K- and L-points of binary systems, which are also termed critical endpoints. The algorithm used the robust critical-point algorithm of Heidemann and Khalil7 to determine a critical point and Michelsen’s stability test9 to locate a stationary point representative of a free-energy tangent plane parallel to the critical-point tangent plane. Test runs on several systems representative of hydrocarbon binaries showed that the algorithm was able to locate both stable and unstable K- and L-points; however, the existence of a critical endpoint was dependent on the binary interaction parameter used in the Peng-Robinson equation of state. Furthermore, convergence to a nontrivial solution required initial critical-point compositions that were similar to the final solution. Initiation of the stationary-point calculation similar to the final solution was not determined to be as important, and initiation near to the borders of composition space is recommended. Not surprisingly, a hybrid successive substitution/Newton-Raphson method gives a balance between robustness and speed when solving the stability equations for a stationary point. The topography of the free-energy surface at a critical endpoint automatically limits the radius of convergence of the problem as formulated. The difficulty of initiating this numerical problem for a binary system is further compounded by its zero degrees of freedom. However, locating these landmark points in P-T-x space may be an important step toward determining the suitability of equations of state in predicting complex phase behavior of petroleum systems and in tuning equation-of-state binary interaction parameters for proper phase behavior predictions in complex multiphase regions. Nomenclature A ) Helmholtz free energy (J/mol) a′ ) expansion coefficient for tangent-plane distance B ) Hessian matrix in critical-point calculations Bij ) element of B b ) Peng-Robinson equation of state co-volume parameter b′ ) expansion coefficient for tangent-plane distance C ) cubic form of the Taylor series expansion of Helmholtz energy c′ ) expansion coefficient for tangent-plane distance d′ ) expansion coefficient for tangent-plane distance
e′ ) expansion coefficient for tangent-plane distance F ) tangent-plane distance Fs ) gradient of F, with respect to parameter s fi ) fugacity of component i (Pa) f ′ ) expansion coefficient for tangent-plane distance G ) Gibbs free energy (J/mol) g ) function to be solved in equilibrium-phase calculations gi ) element of g gˆ ) gradient of F in X gˆ i ) element of gˆ M ) matrix defined by eq 29 N ) number of components n ) mole-number vector ni ) number of moles of the ith component nˆ ) vector representing trial-phase mole numbers in criticalpoint calculations nˆ i ) element of nˆ P ) pressure (Pa) Q ) matrix in the quadratic form of expansion of A Q* ij ) element of Q* Q* ) gradient of Q in the direction of ∆n as defined in eq 7 q ) coefficient in expansion of X qij ) element of Q R ) gas constant (J/(mol K)) r ) mole number vector ri ) element of r s ) distance parameter T ) temperature (K) u ) coefficient in expansion of X ui ) element of u V ) volume (m3) V ) molar volume (m3/mol) w ) coefficient in expansion of X xi ) mole fraction of component i X ) translated mole-number vector in critical-point calculation Xi ) element of X y ) normalized mole-fraction vector y1, y2, ... ) elements of y Y ) un-normalized mole-fraction vector Y1, Y2, ... ) elements of Y z ) critical-phase composition vector zi ) element of z Greek Symbols γ ) Lagrangian multiplier in critical-point calculations δij ) binary interaction parameter for components i and j φi ) fugacity coefficient of component i λ ) acceleration factor for successive substitution method µi ) chemical potential of component i (J/mol) θ ) dimensionless tangent plane distance σij ) Kronecker delta Subscripts and Superscripts crit, C ) critical i, j, k, ... ) component indices 0 ) test phase
9170
Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006
T ) total (0) ) initial guess Acknowledgment Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Department of Chemical Engineering, University of Saskatchewan, is greatly appreciated. Literature Cited (1) Minicucci, D.; Shaw, J. M.; Zou, X.-Y. The Impact of LiquidLiquid-Vapour Phase Behaviour on Coke Formation from Model Coke Precursors. Fluid Phase Equilib. 2002, 194-197, 353-360. (2) Scott, R. L.; van Konynenburg, P. H. van der Waals and Related Models for Hydrocarbon Mixtures. Discuss. Faraday Soc. 1970, (49), 8797. (3) van Konynenburg, P. H.; Scott, R. L. Critical Lines and Phase Equilibria in Binary van der Waals Mixtures. Philos. Trans. R. Soc. London, A 1980, 298, 495-540. (4) Gauter, K.; Heidemann, R. A.; Peters, C. J. Modeling of Fluid Multiphase Equilibria in Ternary Systems of Carbon Dioxide as a NearCritical Solvent and Two Low Volatile Solutes. Fluid Phase Equilib. 1999, 158-160, 133-141. (5) Reid, R. C.; Beegle, B. L. Critical Point Criterion in Legendre Transform Notation. AIChE J. 1977, 23, 726-732. (6) Peng, D.-Y.; Robinson, D. B. A rigorous method for predicting the critical properties of multicomponent systems from an equation of state. AIChE J. 1977, 23 (2), 137-144. (7) Heidemann, R. A.; Khalil, A. M. The Calculation of Critical Points. AIChE J. 1980, 26, 769-779. (8) Michelsen, M. L. Calculation of phase envelopes and critical points for multicomponent mixtures. Fluid Phase Equilib. 1980, 4 (1-2), 1-10. (9) Michelsen, M. L. The Isothermal Flash Problem. Part I: Stability. Fluid Phase Equilib. 1982, 9, 1-19. (10) Michelsen, M. L. The Isothermal Flash Problem. Part II: Phase split calculation. Fluid Phase Equilib. 1982, 9, 21-40.
(11) Michelsen, M. L. Calculation of Critical Points and of Phase Boundaries in the Critical Region. Fluid Phase Equilib. 1984, 16, 57-76. (12) Michelsen, M. L.; Heidemann, R. A. Calculation of Tri-critical Points. Fluid Phase Equilib. 1988, 39, 53-74. (13) Kohse, B. F.; Heidemann, R. A. Computation of Tricritical Points in Ternary Systems. AIChE J. 1993, 39, 1242-1256. (14) Peng, D.-Y.; Robinson, D. B. A New Two Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59-64. (15) Mehra, R. K.; Heidemann, R. A.; Aziz, K. An Accelerated Successive Substitution Algorithm. Can. J. Chem. Eng. 1983, 61, 590596. (16) Knobler, C. M.; Scott, R. L. Multicritical Points in Fluid Mixtures: Experimental Studies. In Phase Transitions and Critical Phenomena; Domb, C., Lebowitz, J. L., Eds.; Academic Press: London, 1984; Vol. 9, pp 163-231. (17) Lam, D. H.; Jangkamolkulchai, A.; Luks, K. D. Liquid-LiquidVapor Phase Equilibrium Behavior of Certain Binary Ethane + n-Alkanol Mixtures. Fluid Phase Equilib. 1990, 59, 263-277. (18) Mushrif, S. H.; Phoenix, A. V. The Effect of Peng-Robinson Binary Interaction Parameters on the Predicted Multiphase Behavior of Selected Binary Systems. To be submitted to Ind. Eng. Chem. Res. (19) Lin, Y.-N.; Chen, R. J. J.; Chappelear, P. S.; Kobayashi, R. VaporLiquid Equilibrium of the Methane-n-Hexane System at Low Temperature. J. Chem. Eng. Data 1977, 22, 402-408. (20) Chang, H. L.; Hurt, L. J.; Kobayashi, R. Vapour-Liquid Equilibria of Light Hydrocarbons at Low Temperatures and High Pressures: The Methane-n-Heptane System. AIChE J. 1966, 12, 1212-1216. (21) Peters, C. J.; Rijkers, M. P. W. M.; Roo, D. J. L.; Arons Swaan, J. D. Phase Equilibria in Binary Mixtures of Near Critical Propane and PolyAromatic Hydrocarbons. Fluid Phase Equilib. 1989, 52, 373-387. (22) Yaws, C. L. Chemical Properties Handbook; McGraw-Hill: New York, 1999.
ReceiVed for reView March 23, 2006 ReVised manuscript receiVed August 29, 2006 Accepted September 18, 2006 IE060359H
