Ind. Eng. Chem. Res. 2007, 46, 5429-5436
5429
GENERAL RESEARCH Improved Two-Sided Tangent Plane Initialization and Two-Phase-Split Calculations Wilson A. Can˜ as-Marı´n,* Julia´ n D. Ortiz-Arango,† and Uriel E. Guerrero-Aconcha‡ Grupo de InVestigacio´ n en Termodina´ mica AVanzada y Feno´ menos en Medios Porosos (Gitefep), Beyoil Ltd., Postal code T2P 2A6, Calgary, Canada
Sandra Hernandez-Ba´ ez Instituto Colombiano del Petro´ leo, ECOPETROL S.A., Piedecuesta, Colombia
The initialization procedure proposed by Michelsen for locating stationary points in the tangent plane distance function was applied to several cases considered as challenging in the literature concerning phase equilibrium calculations. For most of the systems studied in the literature, such a simple procedure is enough to detect the phase instability; thus it is unnecessary to resort to more robust, more elaborated, and more time-consuming initialization schemes. For the conditions and systems where problems were found, a simple procedure was proposed to extend Michelsen’s original scheme with the aim of increasing the probability of obtaining a negative value of the tangent plane distance function. Two algorithms for phase-split calculations were used. Both algorithms use a relaxation parameter to increase their convergence regions. No difficulties were found to calculate the equilibrium distributions in contrast with several convergence problems reported, in which other approaches were applied to the same two-phase systems. 1. Introduction The phase behavior prediction of fluid mixtures is very important for the design, optimization, and operation of chemical plants and the exploitation of hydrocarbon reservoirs. For a given mixture at certain conditions of temperature (T) and pressure (P), it is necessary to know the number of phases and their respective compositional distributions. At constant T and P conditions, the equilibrium state corresponds to the global Gibbs energy minimum. Nevertheless, it is not a trivial problem to find such a global solution due to the multiplicity of false solutions that can be found. The most conventional procedure is to suppose a priori the number of phases in equilibrium and to estimate initial values for the equilibrium factors (K-values), i.e. from saturation calculations when a liquid-vapor is assumed. Another procedure is to formulate the phase equilibrium problem as one of Gibbs energy minimization.1 However, both approaches can fail if the initial estimations of the iteration variables are too inaccurate, requiring sometimes a substantial amount of computational time only to arrive at a trivial solution, or simply to diverge in the case where the assumed number of phases is too high.2 For a review on methods for phase equilibrium calculations and extensive references, see the work of Teh and Rangaiah.3 Michelsen,2 based on the phase stability criterion, initially formulated by Gibbs4 and fully described by Baker et al.,5 presented the first numerical formulation for * To whom correspondence should be addressed. Tel.: +1-403472-0324. E-mail address:
[email protected]. † E-mail address:
[email protected]. ‡ E-mail address:
[email protected].
establishing whether or not a phase is thermodynamically stable, i.e., when a given mixture can split into multiple phases. For an N-component mixture, the tangent plane criterion can be expressed as
g(y) )
1
N
∑ yj( µj(y) - µj(z)) > 0
RT j)1
(1)
where g(y) is the difference between the reduced Gibbs free energy at the overall composition, z, and a certain test composition, y; µj is the chemical potential of the component j. The stability of the mixture with overall composition z is proven when g(y) g 0 for all y. The problem is normally solved either by finding all stationary points of such a function or by finding its global minimum. As Michelsen pointed out, it is sufficient for stability that g(y) is positive at all of its stationary points (especially its minimum), since it will then be positive everywhere.2 Cairns and Furzer6 presented a clear explanation of Michelsen’s algorithm. The main problem is to find all stationary points with complete certainty or, at least, to find a negative value for g(y). In simple terms, the problem is that a phase whose stability is in question must be assured to be at the lowest possible Gibbs energy before instability can be ruled out.7 As it is pointed out by Hua et al.,8 finding all the stationary points is not always necessary nor desirable in a phase stability analysis. If at any point the value of the tangent plane distance function is negative, the instability of the mixture can immediately be assured and then the computational procedure can be followed by phase-split calculations. Most of the existing methods for solving the phase stability problem are local in nature, and produce, at best, only local solutions. For example,
10.1021/ie061361b CCC: $37.00 © 2007 American Chemical Society Published on Web 06/30/2007
5430
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
Figure 1. Scheme of initialization for tangent plane used in this work (E-scheme).
Michelsen’s method solves tangent plane distance function (D) to find its stationary points using multiple initial guesses; but, if negative values are not found, it is not guaranteed that the system is indeed stable.1 Other alternative methods have been proposed: For example, Sun and Seider9 applied a homotopycontinuation method, where two types of initial points for homotopy paths were used; the stationary points are obtained along the homotopy paths. Eubank et al.10 introduced the area method, which is based on an exhaustive search over a grid, initially on a coarse grid and afterward refining the grid for trying to find the global minimum. McDonald and Floudas11,12 advocated global techniques involving branch and bound using convex underestimating functions. For several activity coefficient models, this method guarantees theoretically that the solution corresponding with the global minimum Gibbs energy is found. Later, this work was extended by Harding and Floudas to cubic equation of state models.13 Hua et al.8 presented the method interval Newton/generalized bisection (IN/GB) to locate all stationary points and applied it to systems containing two and three components. As it is pointed out by Harding and Floudas, it is not feasible to exactly localize all finite stationary points. Actually, it is only feasible to enclose all stationary points. The approach of Hua et al. encloses all stationary points of the plane tangent distance function within certain tolerance. Stochastic methods such as genetic algorithms (GA) and simulated annealing (SA) have also been successfully applied. However, these methods tend to be extremely slow for obtaining convergence.3,14,15 Obviously, the global methods increase the possibility for detecting phase instability, but their computational costs are very high. The initialization method proposed here is an extension of the two-sided tangent plane initialization procedure, seeking to increase the probability of obtaining the
correct solution of the tangent plane criterion using the local procedure proposed by Michelsen.2 The K-values obtained from stability analysis are excellent initial guesses for phase-split calculations. Explicitly, Hua et al.8 and Sofyan et al.16 reported difficulties using a commercial software (S1), based on Michelsen’s algorithm, for modeling the binary system consisting of hydrogen sulfur (1) and methane (2) at 190 K, 40.5 bar, and feed composition z1 ) 0.0187. Whereas the system is really unstable, S1 predicts erroneously a stable system. As it is stressed by Hua et al.,8 such a code, based on Michelsen’s tangent plane, can generally be considered as highly reliable. Obviously, faults are possible as was expressed above. It is well-known that a lot of commercial and in-house software currently use Michelsen’s scheme of initialization coupled with local convergence phase-split algorithms. We tested three commercial software packages widely used for modeling hydrocarbon reservoir fluids and designing chemical plants, and faults were found to detect phase instability in the simple binary system mentioned above. Sofyan et al.16 also reported that Michelsen’s algorithm incorrectly predicts a stable single phase for the binary system methane (1)/propane (2) at 277.6 K, 100 bar, and feed compositions z1 ) 0.68 and 0.73. Also, Blanchard et al.17 found problems with the same commercial routine mentioned above and also with another amply used commercial software (S2) for modeling the binary systems trans-2-hexen-1-ol/CO2 and tert-butyl alcohol (1)/ CO2 (2). The same problems have been recently reported by Souza et al.18 In this paper, a procedure to increase the probability of detecting phase instability based on Michelsen’s multiple guesses method is presented to be used with local convergence phase-split routines. As it was previously referred, Hua et al.8
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5431
established that it is not always necessary nor desirable to find all the stationary points to check the phase stability due to the high computational costs. It is only enough to detect a negative value for the tangent plane distance function (D). Obviously, similar to all methods mentioned above, if a negative value is not found, it is not guaranteed that the system is absolutely stable. In this work, two local convergence routines for phasesplit calculations are used. The first one is based on the Leibovici and Neoschil algorithm (called LN hereafter),19 and the second one is based on the Michelsen’s Gibbs energy minimization algorithm;20 but here, K-values are used instead fugacity coefficients in the formulation of the objective function to be minimized (called MM hereafter). For MM, an automatic step control (relaxation parameter) for the Newton-Raphson algorithm is also proposed to increase the convergence region (see Appendix A). The objective of the present work is 2-fold. First, to study the problems found by using the methodology of initial guesses to detect the stationary points in Michelsen’s tangent plane distance function. And second, to check the robustness of the two routines used here, MM and LN, for phase-split calculations in two-phase systems. This latter objective comes up because the problems have not only been found in detecting phase instability, but also in solving the phase-split problem.17,18
(2)
Applying eq 5, the methodology used here to generate initial guesses for searching stationary points in Michelsen’s tangent plane distance function is shown in Figure 1, and a more extensive explanation is given as follows: 1. Apply eqs 2 and/or 3. If instability is detected (D < 0), then go to step 3. Else, continue with step 2. 2. Take ψ equal to 0.1 (or other value preset) in eq 5. If one stationary point (nontrivial) is located, store it; but if it does not have negative tangent plane distance, take ψ ) ψ + 0.1 and use as Ki-values the ones obtained in such a stationary point. In contrast, if a trivial solution was found, then use the Ki-values from the Wilson’s correlation. If the stationary point located has a negative value for the tangent plane distance and if more exploration is desired, increase ψ again by 0.1 and return to step 2 until ψ takes the superior limit of 0.9. Then, select the stationary point with the more negative value, use its Ki-values as initial guesses, and go to step 3. Else, go to step 3. If no stationary points exhibiting a negative value for tangent plane distance can be located, then take Ki-values from the stationary point (nontrivial) with the least positive value as initial guesses and go to step 3. 3. Perform the phase-split calculation (if it is required). Obviously if instability is not found, the procedure could be repeated with ψ values of 0.05, 0.15, ..., 0.95. According to Hoteit and Firoozabadi,22 the K-values obtained in step 2 should be calculated from the mole numbers, Yi, instead from the normalized mole fractions, yi. The explanation is that by using the mole fractions, yi, instead of Yi, the estimation of the K-values is poorer for low pressures. After each iteration, Michelsen2 proposed to calculate the value of a modified distance function, g*(Y), and a convergence variable, r(Y),
(3)
g*(Y) ) 1 +
2. Scheme of Initialization for the Tangent Plane Used in This Work (E-Scheme) For liquid-vapor equilibrium, Michelsen2 proposed to initialize the tangent plane analysis using two types of initial guesses,
YVi ) Ki zi (vapor-type trial phase) Y Li )
zi (liquid-type trial phase) Ki
N
where Yi is defined as the mole number. In cases where the feed of composition, zi, can be unambiguously identified as a liquid or as a vapor, the search is for vapor-type incipient phase (eq 2) or liquid-type (eq 3), respectively. Only one of these guesses is then used. According to Michelsen, normally one of these initial estimations will converge to the trivial solution and the other one to the desired minimum for a specification in the two-phase region. Obviously, in some cases, both may lead to the trivial solution, and in some cases, two negative solutions will be obtained. In the latter case, the more negative solution will provide a much better guess of the K-values for the solution of the phase-split algorithm. However, whenever it is not feasible to define unambiguously a mixture as a vapor- or liquid-type, both guesses are generally used, e.g., near to the critical point or close to the three-phase region. For this reason, we prefer to define a mathematical function to take these possibilities into account. Our function is
Y hi )
ψY Li
+ (1 -
ψ)YVi
Y hi )
zi - 1) + 1] Ki
ln φi(Y) - ln Zi - ln φi(Z) - 1) (6) r(Y) )
(5)
Here, Ki are the equilibrium constants for each component (K-values), initially evaluated with the Wilson’s correlation.21
2g* β
(7)
with N
β)
∂g*
(Yi - Zi) ∑ ∂Y i)1
(8)
i
Here, r is calculated at each iteration, and if its value exceeds 0.8, the solution is considered to be convergent to the trivial solution; then, the search is abandoned for that guess and another guess is taken. Another simple equation for checking if a trivial solution is being approached is given by23 N
( 1 - Ki)2 < 10-4 ∑ i)1
(4)
where ψ is an interpolating parameter taking values within the interval [0,1]. If eqs 2 and 3 are introduced into eq 4 and then rearranged, eq 5 is directly obtained as
[ψ(K 2i
Yi(ln Yi + ∑ i)1
(9)
A convergence criterion for stability analysis calculation was taken as23 N
∑ i)1
( ) 1-
fi(z) N
fi(y)
2
< 10-12
(10)
Yi ∑ i)1
These criteria were enough for the systems modeled here. Equation 9 is evaluated at each iteration, and if the criterion is
5432
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 Table 1. Phase-Stability Analysis Using E-Scheme for the Unstable Feed Composition (z1 ) 0.0187) of the Binary System Hydrogen Sulfide (1) and Methane (2) at 190 K and 40.53 Bar Using the SRK Equation of State
Figure 2. Variation of tangent plane distance function (D) with the number of iterations for several ψ values: binary system hydrogen sulfide (1) and methane (2) for z1) 0.0187 at 190 K and 40.53 bar using the SRK equation of state. The legend for ψ values is as follows: (]) ) 0, (9) ) 0.1, (4) ) 0.2, ([×]) ) 0.3, (/) ) 0.4, (O) ) 0.5, ([+]) ) 0.6, (-) ) 0.7, (b) ) 0.8, ([) ) 0.9, (0) ) 1.
accomplished, then the solution is considered to be convergent to the trivial solution, the search is abandoned for that guess, and another guess is taken. This procedure avoids useless iterations and reduces computer time. 3. Results and Discussion The methodology proposed here was applied to several test problems involving vapor-liquid and liquid-liquid equilibrium in binary, ternary, and quaternary mixtures at several operation conditions. Some of them were also used by Hua et al.,8 Sofyan et al.,16 and Harding and Floudas13 to show the strengths of their respective methods. Both Soave-Redlich-Kwong (SRK) and Peng-Robinson (PR) equations of state are used. Critical properties, acentric factors, and binary interaction coefficient data used here are reported by them. 3.1. System H2S (1)/CH4 (2). Definitely, Michelsen’s scheme (ψ ) 0 and ψ ) 1) fails to localize the instability for the binary system H2S/CH4 at 190 K, 40.5 bar, and feed composition z1 ) 0.0187. Considering the feed as completely vapor and searching a liquid trial phase (applying eq 3) and then as completely liquid and searching a vapor trial phase (applying eq 2), the Michelsen’stangentplaneconvergesto D)0.0109and0.000 001 047 61 (here, the incipient phase is equal to the feed phase), respectively, thus predicting erroneously a stable mixture under such conditions. The results applying the methodology proposed are shown in Table 1. Results in bold type mean nontrivial solutions. It is noted that 55% of the calculations performed, varying the interpolating parameter ψ in eq 5, converged to the global minimum value reported by Hua et al.8 and Sofyan et al.16 as -0.004. Yi values coincided exactly with those shown by them. For this binary mixture, Figure 2 shows the variation of the tangent plane distance function (D) with the progress in the iterative process for several values of the interpolating parameter ψ. For certain ψ values, D is negative from the first iteration. It is important to point out the high sensitivity of the numerical scheme to ψ-values for this binary system. This high sensitivity increases the probability for detecting phase instability. Liquidliquid equilibrium was predicted at these conditions. 3.2. System CH4 (1)/C3H8 (2). Although Sofyan et al.16 reported difficulties with Michelsen’ scheme for the binary mixture methane (1) and propane (2) at 277.6 K, 100 bar, and feed compositions z1 ) 0.68 and 0.73 in predicting an erroneous
interpolating parameter value
tangent plane distance (D)
Y (H2S)
Y (CH4)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.010920981 -0.003940130 0.000000279 0.010921066 -0.003939701 -0.003939701 -0.003939701 -0.003939702 0.000000278 -0.003939741 0.000001048
0.884774766 0.076674045 0.018673397 0.884776349 0.076674138 0.076675199 0.076675199 0.076675199 0.018673407 0.076675094 0.018648518
0.115225234 0.923325955 0.981326603 0.115223651 0.923325862 0.923325862 0.923325862 0.923325862 0.981326593 0.923324906 0.981351482
Table 2. Phase-Stability Analysis Using E-Scheme for the Unstable Feed Composition (z1 ) 0.68) of the Binary System Methane (1) and Propane (2) at 277.6 K and 100 Bar Using the SRK Equation of State interpolating parameter value
tangent plane distance (D)
Y (CH4)
Y (C3H8)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-4.91998E-05 -4.88196E-05 -4.83298E-05 -4.76756E-05 -4.67588E-05 -4.53841E-05 -4.31024E-05 -3.86117E-05 -2.61686E-05 -0.000334574 -0.000334445
0.677479936 0.677496702 0.677518342 0.677547321 0.677588089 0.677649559 0.677752509 0.677958668 0.678557698 0.772445454 0.772465025
0.322520064 0.322503298 0.322481658 0.322452679 0.322411911 0.322350441 0.322247491 0.322041332 0.321442302 0.227554546 0.227534975
Table 3. Phase-Stability Analysis Using E-Scheme for the Unstable Feed Composition (z1 ) 0.73) of the Binary System Methane (1) and Propane (2) at 277.6 K and 100 Bar Using the SRK Equation of State interpolating parameter value
tangent plane distance (D)
Y (CH4)
Y (C3H8)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-0.000294813 -0.000294813 -2.07991E-05 -0.000297252 -2.34044E-05 -0.00199265 -0.000294806 -2.11022E-05 -0.00029481 -2.27886E-05 -2.34569E-05
0.650282923 0.650282926 0.757019969 0.650259539 0.757073055 0.729031364 0.650282994 0.757026041 0.650282952 0.757060373 0.75707414
0.349717077 0.349717074 0.242980031 0.349740461 0.242926945 0.270968636 0.349717006 0.242973959 0.349717048 0.242939627 0.24292586
stable thermodynamically system; we did not find such results. Tables 2 and 3 show the calculations for both compositions. In Table 2, if ψ ) 1 or 0.9, the methodology detects immediately the global minimum reported by Hua et al. and Sofyan et al. as -0.000 33. The trivial solution is the most obtained. This system is, under such conditions, near to its critical point, and it is considered to be one very difficult to model.3 For z1 ) 0.73, if ψ ) 0, the global minimum of -0.000 29 is easily detected (see Table 3). 3.3. System CH4 (1)/CO2 (2)/H2S (3). This system was modeled by Sofyan et al.16 Here, only results for two compositions presenting multiple local minima were chosen to be reported; no problems were found for the other cases. Such compositions and conditions are the following: [0.4989 (1), 0.0988 (2), 0.4023 (3)] at 208.5 K and 55.1 bar and [0.48 (1), 0.12 (2), 0.40 (3)] at 210 K and 57.5 bar. Table 4 shows the
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5433 Table 4. Phase-Stability Analysis Using E-Scheme for the Unstable Feed Compositions [0.4989 (1), 0.0988 (2), 0.4023 (3)] of the Ternary System Methane (1), Carbon Dioxide (2), and Hydrogen Sulfide (3) at 208.5 K and 55.1 Bar Using the PR Equation of State interpolating parameter
tangent plane distance (D)
Y (CH4)
Y (CO2)
Y (H2S)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-0.006823587 -0.006823559 -0.014704852 -0.00682343 -0.009740324 -0.006823464 -0.009740324 -0.006823427 -0.014704852 -0.006823467 -0.014704845
0.244858084 0.244858186 0.919378017 0.244858632 0.817478814 0.24485852 0.817478814 0.244858639 0.919378017 0.24485851 0.919377908
0.091253686 0.091253697 0.034319845 0.091253747 0.059930009 0.091253734 0.059930009 0.091253748 0.034319845 0.091253733 0.034319879
0.66388823 0.663888117 0.046302138 0.663887622 0.122591177 0.663887745 0.122591177 0.663887613 0.046302138 0.663887757 0.046302214
Table 5. Phase-Stability Analysis Using E-Scheme for the Unstable Feed Compositions [0.48 (1), 0.12 (2), 0.40 (3)] of the Ternary System Methane (1), Carbon Dioxide (2), and Hydrogen Sulfide (3) at 210 0.5 K and 57.5 Bar Using the PR Equation of State interpolating parameter
tangent plane distance (D)
Y (CH4)
Y (CO2)
Y (H2S)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-0.00223853 -0.002239119 -0.005478022 -0.002241149 -0.005478056 -0.002238578 -0.005478058 -0.002241149 -0.005478022 -0.002238525 -0.005478022
0.289843292 0.289841106 0.902920767 0.289832868 0.90292051 0.289843157 0.902920553 0.289832868 0.902920912 0.28981841 0.902920613
0.116635493 0.1166353 0.044562552 0.116634585 0.044562639 0.11663548 0.044562624 0.116634585 0.044562503 0.116633106 0.044562604
0.593521215 0.593523594 0.052516681 0.593532547 0.052516851 0.593521363 0.052516823 0.593532547 0.052516585 0.593548484 0.052516784
results for the first composition. As it can be observed, the Michelsen’s scheme (ψ ) 0 and ψ ) 1) detects the phase instability and localizes the global minimum (-0.015). The methodology proposed additionally locates the stationary points -0.0097 and -0.0068, also found by Sofyan et al. Table 5 shows the results for the second composition. As it can be observed, the Michelsen’s initialization (ψ ) 0 and ψ ) 1) detects phase instability and also locates the global minimum reported by Sofyan et al. as -0.0055. No trivial solutions were obtained with the methodology used here. All the liquid-vapor systems reported by Hua et al. and Sofyan et al. were also studied, and Michelsen’s liquid-vapor initialization procedure did not have problems detecting the instability of phases. These include the following: CO2/CH4, N2/CH4/C2H6, and CH4/CO2/H2S/H2O, at all conditions of temperature and pressure reported by them. 3.4. System trans-2-Hexen-1-ol (1)/CO2 (2). Blanchard et al.17 and Souza et al.18 found problems for modeling this system
using S1 and two modules from S2, so-called FLASH3 and RGIBBS (see Table 2 in the work of Blanchard et al. or Table 6 in the work of Souza et al.). In short, for a mole fraction of CO2 equal to 0.985 at 303.15 K and 63.819 bar, the three routines predict phase instability, but the phase-split calculations give erroneous solutions. If the pressure is 67.871 bar, FLASH3 and RGIBBS give correct results, but S1 does not converge after 1000 iterations. For a mole fraction of CO2 equal to 0.80 at 303.15 K and 71.007 25 bar, the three routines predict correctly an instable system, but the three routines give false phase-splits. If the pressure is 71.00826 bar, the three routines again predict phase instability, but FLASH3 and RGIBBS give erroneous solutions; only S1 gives a correct result. For a mole fraction equal to 0.7 at 303.15 K and 70.09 bar, the three routines predict erroneously a stable system, and thus, no phase-split calculation is feasible. For a mole fraction equal to 0.97 at 323.15 K and 97.15 bar, FLASH3 and RGIBBS give correct results, but S1 produces a numerical computation error. For a mole fraction equal to 0.742 at 323.15 K and 130 bar, FLASH3 and S1 give correct results, but RGIBBS predicts erroneously a stable system. If the pressure is 135 bar, S1 and RGIBBS give correct solutions, but now, FLASH3 predicts erroneously a stable system. Therefore, some problems have not only been found in detecting phase instability but also in solving the phase-split problem.17,18 The last case is mainly related to the use of deficient or bad conditioned numerical schemes. Generally, the K-values obtained from the stability analysis are excellent initial guesses for solving two-phase-split problems. We also modeled this system using the pure component properties reported by Stradi et al.24 and the two phase-split routines, LN and MM, mentioned above. Our results are shown in Table 6. In short, no problems were found to obtain the compositions in equilibrium. However, for a mole fraction equal to 0.7 at 303.15 K and 70.09 bar, it was not feasible to localize a stationary point with negative distance. But using the K-values, K(1) ) 0.003 244 and K(2) ) 1.424 602, obtained in the stationary point with D ) 0.001805, both routines converged after 83 iterations to accomplish the convergence criterion given N by 1/(N)∑i)1 [ln(fVi /fLi )]2 < 10-18, irrespective of the vaporphase fraction within the interval [0,1] used as initial guess. The equilibrium compositions reported by Blanchard et al.17 were always obtained. A vapor-liquid equilibrium was predicted with a vapor-phase fraction equal to 0.0576. For this system, it was intended to obtain convergence using directly the K-values from Wilson’s correlation,21 i.e., K(1) ) 0.000 057 4, K(2) ) 1.029 845, and then varying the initial guess for the vapor-phase fraction; but, all the attempts always produced numerical errors. Therefore, an algorithm for phase-split calculations based on asuming, a priori, a two-phase system
Table 6. Results of Phase-Stability and Phase-Split Calculations for the Binary Systems trans-2-Hexen-ol (1)/CO2 (2) and tert-Butyl Alcohol (1)/CO2 (2) Using E-Scheme and the LN or MM Algorithmsa binary mixture
this work
trans-2-hexen-1-ol (1)/CO2 (2)
S1
phase-stability
phase-split
phase-split
case
feed Z(2)
temp (K)
pressure (bar)
Y(1)
Y(2)
X(2)
Y(2)
X(2)
Y(2)
1 2 3 4 5 6 7 8
0.985b
303.150
0.8c
303.150
0.7b 0.97b 0.742b
303.150 323.150 323.150
63.819 67.871 71.00725 71.00826 70.090 97.750 130.000 135.000
0.0004216 0.0005735 0.0013013 0.0013021 0.0009750 0.0091621 0.0421317 0.0499334
0.999578 0.999426 0.998699 0.998698 0.999025 0.990838 0.957868 0.950067
0.6132015 0.6138379 0.6852428 0.6852447 0.6834999 0.6282918 0.7245074 0.7356604
0.9995342 0.9995324 0.9689110 0.9689097 0.9700945 0.9947671 0.9553362 0.9488434
0.6131 0.6683 0.6846 0.6846 0.6828 0.6281 0.7240 0.7352
0.9995 0.9993 0.969 0.969 0.9702 0.9947 0.9554 0.9489S1
74.660
0.003710
0.996291
0.993527
0.996234
0.9937
0.9963
tert-butyl alcohol (1)/CO2 (2) 9 a
0.995b
305.950
Comparisons with the results from S1 software.8 PR equation of state. b Vapor-liquid equilibrium. c Liquid-liquid equilibrium.
5434
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
and, moreover, using directly the equilibrium constants from Wilson’s correlation is prone to fail. The K-values in the equilibrium are K(1) ) 0.094 488 and K(2) ) 1.419 305. The deviation of the K-values for CO2 in the stationary point with positive distance is only 0.37% from the equilibrium value; whereas the equilibrium constants from the Wilson’s correlation have a very high deviation of 27%. 3.5. System tert-Butyl Alcohol (1)/CO2 (2). This binary system was also modeled by Blanchard et al.17 S1 did not converge after 1000 iterations, and RGIBBS did not detect phaseinstability; only FLASH3 converged to the correct solution. Our results are presented as case 9 in Table 6. We did not find problems predicting the phase instability and the compositions in equilibrium; both LN and MM converged easily and produced the same results. 4. Conclusions Any software with stability analysis based on Michelsen’ algorithm with his scheme of multiple initial guesses is expected to be very reliable for modeling two-phase systems. However, in some cases, it can fail in detecting phase instability, as it was shown in this work. A very simple and easy to implement extended scheme (E-scheme) for initializing the stability analysis with the Michelsen’s tangent plane was presented. This scheme increases the probability to find phase instabilities and still keeps the simplicity of Michelsen’s procedure. However, for the binary system trans-2-hexen-ol/CO2 with zCO2 ) 0.7, it was not possible to obtain a negative value for the tangent plane distance function. In cases like this, the advantage of other more elaborated algorithms such as Hua et al.’s is obvious. In addition, it was very interesting to find that by using the K-values from a stationary point with positive distance as initial guesses for both phase-split routines, it was easy to calculate the compositions in the equilibrium. In contrast, it was not possible to find such compositions using the K-values obtained directly from Wilson’s correlation as initial guesses. For the systems modeled here, both local convergence routines for phase-split calculations, LN and MM, did not show problems in finding the respective twophase equilibrium solutions. Appendix A Michelsen20
Mathematical Relations for the MM Routine. proposed the following objective function to be minimized NP
F(β) )
NC
βj - ∑ zi ln Ei ∑ j)1 i)1
(A1)
The mole fraction of the component i in the phase j is given by
yij )
zi φij Ei
where φij is the fugacity coefficient of component i in phase j. Equation A4 accomplishes the mass balance requirement and equal fugacity values in all phases in equilibrium.20 Using K-Values. In phase behavior calculations, the K-values are more frequently used than the fugacity coefficients. The K-values are defined as
Kij )
φir φij
Here, φir is the fugacity coefficient of the component i in a phase taken as reference, r. According to Leibovici,25 instead of using eqs A2 and A4, it is possible to define the following alternative relations in terms of the K-values NP
Ei )
Kij βj ∑ i)1
(A6)
zi K Ei ij
(A7)
and
yij )
Material Balance Restriction. The material balance for a multiphase multicomponent system is given by NP
βj yij ) zi, ∑ j)1
i ) 1,2, ..., NC
βj
∑ j)1 φ
(A8)
Replacing eq A7 into A8 and using eq A6, we obtain NP
zi
zi
i
i
NP
βj Kij ) ∑ Kij βj ) zi ∑ E E j)1 j)1
(A9)
This shows that the material balance is fulfilled. Isofugacity Restriction. For a multicomponent multiphase system in equilibrium, the fugacity of the component i is equal in all phases. In short,
fi1 ) fi2 ) ... ) fiNP, i ) 1,2, ..., NC
yi1φi1 ) yi2φi2 ) ... ) yiNPφiNP, i ) 1,2, ..., NC NP
(A5)
(A10a)
or using fugacity coefficients,
where,
Ei ≡
(A4)
(A10b)
(A2)
Taking the phase 1 as reference in eq A5 (using also eq A7) and replacing into eq A10b, the result is
(A3)
yiNP yi2 zi yi1 ) ) ... ) , i ) 1,2, ..., NC Ki1 Ki2 KiNP Ei
ij
β ) (β1, β2, ..., βNP)
β represents the mole amount of each phase, and zi is the overall composition of component i. F(β) is a convex function defined on the convex set β. When the number of components is equal or exceeds the number of phases, the Hessian matrix is positive defined, and F(β) is a strictly convex function. The minimum of a strictly convex function is unique. Therefore, the β vector can be obtained.
(A11)
This shows that the material balance is also fulfilled. Here, eq A1 in terms of the K-values was solved using Newton optimization. Where,
∆β ) -H-1 grad F(β)
(A12)
∆β is the search direction vector, H is the Hessian matrix, and
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5435
grad F(β) is the first derivative vector of the objective function. The respective elements are given by
∂F(β)
NC
)1-
∂βj
∑ i)1
zi Kij , j ) 1,2, ..., NP Ei
(A13)
and
∂2F( β)
NC
)
∂βS∂ βm
∑ i)1
Nomenclature
zi Ksi Kmi Ei
(A14)
Relaxation Parameter. We decided to use a relaxation parameter seeking to control the Newton-Raphson procedure. Defining λh as the relaxation parameter, eq A12 can be written as
βk(m+1) ) βk(m) + λh(m)∆βk(m)
(A15)
For that, positive compositions, yij > 0, are obtained. According to eq A7, the following must be fulfilled:
zi K > 0, i ) 1,2, ..., NC Ei ij
(A16)
)
Kij βj ∑ j)1
(m+1)
) NP
Kij(βj(m) + λ(m)∆βj(m)) > 0 ∑ j)1
[ ]
(A17)
NP
λi(m) > -
Kij βj(m) ∑ j)1 NP
, i ) 1,2, ..., NC
(A18)
Kij∆βj(m) ∑ j)1
There is a different λ for each component i. By definition, λ is positive. Therefore, only λi’s with negative denominators are taken into account. The maximum, but not acceptable, value is given by
λ(MAX)(m) ) min (λi(m)), i ) 1,2, ..., NC
(A18)
In eq A15, λh(m) ) pλ(MAX)(m), where p is a fraction, e.g., 0.9. If the positive compositional restriction is not violated, the iterative procedure is carried out using a λ value equal to 1 in eq A15. Obviously, a similar expression to eq A17 but in terms of fugacity coefficients is given simply by
[ ] ( ) ∑( ) NP
λ
Superscripts L, V ) liquid and vapor phases, respectively.
>-
∑ j)1
βj(m) φij
NP
∆βj(m)
j)1
φij
, i ) 1,2, ..., NC
β ) used for calculating r (eq 8) µ ) chemical potential φ ) fugacity coefficient ψ ) interpolating parameter Literature Cited
We finally obtain
(m)
i, j ) index for component
Greek Symbols
NP
Ei
D ) tangent plane distance f) fugacity g ) reduced Gibbs free energy K ) equilibrium constant N ) number of components r ) convergence criterion for trivial solution (eq 7) R ) ideal gas constant T ) temperature y ) trial phase mole fraction Y ) mole number z ) feed composition Subscripts
But, by definition, zi > 0 and Kij > 0. Therefore, Ei should also be larger than zero. Replacing eq A15 into eq A6, (m+1)
Multiplying and dividing eq A19 by φir, eq A17 is immediately obtained. Both LN and MM methods are designed to calculate phase distributions at fixed K-values. In an outer loop, equilibrium ratios are updated by successive substitution (first-order), while in an inner loop β is updated by Newton-Raphson (secondorder).
(A19)
(1) Gautam, R.; Seider, W. D. Computation of phase and chemical equilibrium. Part I. Local and constrained minima in Gibbs free energy. AIChE J. 1979, 25, 991. (2) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (3) Teh, Y. S.; Rangaiah, G. P. A Study of Equation-Solving and Gibbs Free Energy Minimization Methods for Phase Equilibrium Calculations. Chem. Eng. Res. Des. 2002, 80, 743. (4) Gibbs, J. W. A Method of Geometrical Representation of the Thermodynamic Properties of Substances by Means of Surfaces. Trans. Conn. Acad. Arts. Sci. 1873, 11, 382. (5) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. SPE J. 1982, 22, 731. (6) Cairns, B. P.; Furzer, I. A. Multicomponent Three-Phase Azeotropic Distillation. 2. Phase-Stability and Phase-Splitting Algorithms. Ind. Eng. Chem. Res. 1990, 29, 1364. (7) Heidemann, R. A. Computation of High Pressure Phase Equilibria. Fluid Phase Equilib. 1983, 14, 55. (8) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Enhanced interval analysis for phase stability: Cubic Equation of State models. Ind. Eng. Chem. Res. 1998, 37, 1519. (9) Sun, A. C.; Seider, W. D. Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilib. 1995, 103, 21. (10) Eubank, P. T.; Elhassan, A. E.; Barrufet, M. A.; Whiting, W. B. Area Method for the Prediction of Fluid Phase Equilibria. Ind. Eng. Chem. Res. 1992, 31, 942. (11) McDonald, C. M.; Floudas, C. A. Global optimization for the phase and chemical equilibrium problem: application to the NRTL equation. Comput. Chem. Eng.1995, 19, 1111. (12) McDonald, C. M.; Floudas, C. A. GLOPEQ: A new computational tool for the phase and chemical equilibrium problem. Comput. Chem. Eng. 1997, 21, 1. (13) Harding, S. T.; Floudas, C. A. Phase Stability with Cubic Equations of State: Global Optimization Approach. AIChE J. 2000, 46, 1422.
5436
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
(14) Zhu, Y. High pressure phase equilibrium through the simulated annealing algorithm: application to SRK and PR equation of state. Presented to the 2000 AIChE meeting, 2000. (15) Gomez, S.; Gosselin, O.; Barker, J. W. Gradient-based historymatching with a global optimization method. SPE J. 2001, 200. (16) Sofyan, Y.; Ghajar, A. J.; Gasem, K. A. M. Multiphase Equilibrium Calculations Using Gibbs Minimization Techniques. Ind. Eng. Chem. Res. 2003, 42, 3786. (17) Blanchard, L. A.; Xu, G.; Stadtherr, M. A.; Brennecke, J. F. Phase BehaVior and Its Effects on Reactions in Liquid and Supercritical CO2; Department of Chemical Engineering, University of Notre Dame: Notre Dame, IN. (18) Souza, A. T.; Cardozo-Filho, L.; Wolf, F.; Guirardello, R. Application of Interval Analysis for Gibbs and Helmholtz Free Energy Global Minimization in Phase Stability Analysis. Braz. J. Chem. Eng. 2006, 23, 117. (19) Leibovici, C. F.; Neoschil, J. A solution of Rachford-Rice equations for multiphase systems. Fluid Phase Equilib. 1995, 112, 217. (20) Michelsen, M. L. Calculation of multiphase equilibrium. Comput. Chem. Eng. 1994, 18, 545.
(21) Wilson, G. M. A Modified Redlich-Kwong Equation of State, Application to General Physical Data Calculation. Paper presented at the AIChE 69th National Meeting, Cleveland, Ohio, May 4-7, 1969. (22) Hoteit, H.; Firoozabadi, A. Simple phase stability-testing algorithm in the reduction method. AIChE J. 2006, 52, 2909. (23) Soreide, I. Improved Phase Behavior Predictions of Petroleum Reservoirs Fluids from a Cubic Equation of State. Doctoral Thesis, The Norwegian Institute of Technology, Norway, April 1989. (24) Stradi, B. A.; Stadtherr, M. A.; Brennecke, J. F. Multicomponent phase equilibrium measurements and modeling for the allylic epoxidation of trans-2-hexen-1-ol to (2R, 3R)-(+)-3-propyloxiramethanol in highpressure carbon dioxide. J. Supercrit. Fluids 2001, 20, 1. (25) Leibovici, C. F. Private Communication, 2003.
ReceiVed for reView October 23, 2006 ReVised manuscript receiVed April 9, 2007 Accepted May 22, 2007 IE061361B
