4348
Ind. Eng. Chem. Res. 1996, 35, 4348-4353
Maximum Partial Area Rule for Phase Equilibrium Calculations Guor-Shiarn Shyu, Nishawn S. M. Hanif, Kenneth R. Hall, and Philip T. Eubank* Chemical Engineering Department, Texas A&M University, College Station, Texas 77843-3122
Eubank and Hall have shown recently that the tangent-line criterion reduces to an equal area construction for the derivative of the total Gibbs energy plotted against composition. This paper presents another criterion along with a thermodynamic proof, the maximum partial area rule (MPAR), which establishes an area on the plot of the derivative of the Gibbs energy with respect to composition which is always a maximum at equilibrium. The method based upon this criterion is especially powerful for complex phase equilibria involving multiple components and phases. Introduction The phase equilibrium behavior of fluid mixtures is an important design consideration for both chemical processes and oil production. Phase equilibrium calculations commonly employ one of two techniques, K value or Gibbs energy minimization (GEM), as noted by Ammar and Renon (1987). The first procedure solves simultaneously a set of material balances and thermodynamic equations. The second procedure minimizes the Gibbs energy with respect to the composition of the different components in the mixture. The importance of using GEM in phase equilibria calculations has received considerable attention from authors including Baker et al. (1982), Michelsen (1982a,b), and Eubank et al. (1992). Techniques for using GEM to solve phase equilibrium problems from basic thermodynamic equations are also common, such as those of Heidemann (1974), Gautam and Seider (1979), Trangenstein (1985), and Cairns and Furzer (1990). The two-phase, n-component equilibrium problem is geometrically a determination of where a tangent (hyper)plane touches a (hyper)surface at two points (compositions) without intersection of the surface, the latter condition ensuring a global minimum solution to the equilibrium problem. This geometry is equivalent to having fugacity balances for all the components. The number of degrees of freedom on where to put the plane at a given P and T is n - 2. The problem is closed when, at P and T, one specifies a feed composition, i.e., the degrees of freedom are then zero. Geometrically, this is equivalent to requiring that the line connecting the two compositions mentioned above also passes through the specified feed composition. It must be understood that we are only taking a new view angle with respect to a well-known phase equilibrium geometry; we are not changing it. In essence, we are suggesting a new algorithmic attack on a wellknown phase equilibrium problem. Eubank and Hall (1995) have shown that the geometry of the tangent-line criterion proposed by Baker et al. (1982) can be recast so that loops are emphasized which enable the execution of algorithms with easy implementation and fast convergence. Their method is an equal area construction for the derivative of the Gibbs energy with respect to composition. Shyu et al. (1995) generated a more general procedure and applied * To whom correspondence should be addressed. Phone: (409) 845-3339. FAX: (409) 845-6446. E-mail: p-eubank@ tamu.edu.
S0888-5885(96)00203-5 CCC: $12.00
Figure 1. Illustrative plot of the gi′ function for a gE model.
the method to ternary systems, including two- and three-phase systems. They showed that the equal area rule (EAR) has decided advantages compared to the tangent plane method for calculations near critical points, for critical end points, and in the retrograde region because of its additional sensitivity. Shyu et al. (1995) use a tie-line vector procedure to search the equilibrium slope for the equilibrium compositions. Their method seeks the roots of an orthogonal function. In this paper, we introduce a different criterion to establish the equilibrium slope. This criterion is the maximum partial area rule (MPAR), which we have based upon the observation that a certain area on the plot of the derivative of Gibbs energy with respect to composition against composition is always a maximum at equilibrium. We also present a thermodynamic proof for the new criterion. Equal Area Rule (EAR) Figure 1 presents a typical plot of the derivative of the Gibbs energy with respect to composition against x1 in which the Gibbs energy has been calculated from a gE model. The equilibrium compositions are xR1 and xβ1. Equilibrium compositions result from satisfying the tangent-line criterion:
( ) ( ) ∂g
R
∂x1
T,P
)
∂g
β
∂x1
T,P
gβ - gR )
xβ1 - xR1
≡ gj 1′
(1)
where g is the total Gibbs energy or the Gibbs energy of mixing reduced by RT, R is the gas constant, T is © 1996 American Chemical Society
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4349
Figure 3. Illustrative phase diagram for a ternary system.
upon only two variables, x1 and x2, in a ternary system:
g ) g(x1, x2) Figure 2. Flow diagram applying the EAR in a ternary system.
temperature, P is pressure, and superscripts R and β denote separate equilibrium phases. Eubank and Hall (1995) derive the equal area rule which satisfies eq 1:
∫xx
( )
β 1 R 1
∂g
∂x1
dx1 ) gβ - gR ≡ gj 1′(xβ1 - xR1 )
(2)
T,P
where gj 1′ is the value of (∂g/∂x1)T,P at equilibrium. In Figure 1, equilibrium occurs when a specific value of g1′ equalizes the upper area, U, and the lower area, L. Eubank and Hall (1995) also have proposed an efficient algorithm which normally converges in one to three iterations for equal area problems assuming reasonable, if not accurate, initial guesses. Examples using this method with binary systems can be found in Eubank and Hall (1995) and Hanif et al. (1996). Equal Area Rule for Ternary Systems Application of the EAR to ternary systems requires the introduction of a slope variable defined by
D12 ≡
( ) ∂x2 ∂x1
(3)
T,P
The value of this variable comes from the tie-line vector procedure by which we mean the external loop of the procedure used in binary or essential-binary systems to search for gj 1′ as the internal loop. The flow chart for using the EAR method with a ternary system appears as Figure 2. The reason for using the tie-line vector procedure is to solve problems involving a multiple-component system using a procedure similar to that used for a binary system. In fact, the current procedure reduces treatment of multicomponent systems effectively to that of a binary system, thus the nomenclature essentialbinary. For example, at constant temperature and pressure, the restriction that the mole fractions sum to unity implies that the Gibbs energy of mixing depends
(4)
Figure 3 illustrates that the compositions of the two phases in equilibrium lie on a straight line which passes through the overall composition of the mixture (z1, z2 ) and satisfies the material balances as in a standard flash calculation. The constant slope of the tie line represents the relative change between two components with respect to their overall composition (z1, z2). For multicomponents, Shyu et al. (1995) define the relative change function, Dij, for components (i, j) at an overall composition (zi, zj) along a tie line as
Dij ≡
()
xj - zj
∂xj
∂xi
)
T,P
xi - zi
(i * j)
(5)
For binary systems, D12 ) -1, whereas for general multicomponent systems, N-1
N
∑ ∑ Dij ) -1 i)1 j)i+1
(6)
where N is the number of components. For ternary systems, Dij is the slope of a tie line. If Dij and the overall composition, [zi], are available, x2 and x3 can be written in terms of x1 using
x2 ) D12(x1 - z1) + z2
(7)
x3 ) D13(x1 - z1) + z3
(8)
D12 + D13 ) -1
(9)
Substituting eqs 7, 8 and 9 into eq 4 produces
g ) g(x1, D12)
(10)
and its derivative becomes
g 1′ )
( ) ∂g
∂x1
(11)
T,P,D12
Using D12 as the search variable, the g function has the same complexity as a binary system and the basic EAR method of eqs 1-3 applies.
4350 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
Shyu (1995) presents the expression for calculating g′ of multicomponent systems:
( )
gi′ ) ln
fˆi fN
N-2 N-1
+
fi fˆN
∑ ∑ i)1 j)i+1
Dij ln
( ) fˆj fN
fj fˆN
(12)
where f is a pure-component fugacity and fˆi is the fugacity of component i in the mixture. In ternary systems using an activity coefficient model, eq 12 becomes:
( )
g1′ ) ln
( )
x1γ1 x2γ2 + D12 ln x3γ3 x3γ3
(13)
where γi is the activity coefficient of component i in the mixture. Using an EOS/MCR (equation of state/mixture combining rule) model, the equivalent expression is
( )
g1′ ) ln
x1φˆ 1φ3
x3φˆ 3φ1
+ D12 ln
( ) x2φˆ 2φ3
(14)
x3φˆ 3φ2
where φi is the fugacity coefficient of pure component i and φˆ i is the fugacity coefficient of component i in the mixture. Using the tie-line vector procedure for the ternary system enables us to observe van der Waals loops as in binary systems. Shyu et al. (1995) propose using -1 as the initial guess for D12 in the dilute third-component region. The initial guess for additional D12 is an extrapolation from the converged D12 of the previous equilibrium determination. The procedure is based partially upon the law of rectilinear diameters and is effective when approaching critical points (plait points). The equilibrium D12 occurs when the function N-1
Ξ≡
∑ i)1
[( ) ( ) ] ∂g
∂g
R
β
-
∂xi
T,P,xj
∂xi
Figure 4. Calculation of Ξ as a function of D12 in the two-phase region. Equilibrium occurs when Ξ equals zero. This diagram uses a LLE model, GE/RT ) 3.2x1x2 + 1.5x1x3, with an overall composition of z1 ) 0.35 and z2 ) 0.35.
(j * i, j < 3) (15)
T,P,xj
approaches zero. Ξ represents the sum of the differences of gradients between two equilibrium points. The equilibrium compositions result from a search along the tie-line slope D12, which makes Ξ less than a tolerance value (e.g., 1 × 10-7). Because Ξ depends only upon the slope D12 in a ternary system, it is easy to find using any root searching procedure (e.g., Newton-Raphson or bisection). Figure 4 illustrates Ξ calculated using a Margules model for the excess Gibbs energy plotted against the assumed slope D12. Phase loops do not appear if the assumed D12 is either too small or too large.
Figure 5. Illustrative plot demonstrating the areas U, L, and AOC.
the same composition. On Figure 5, EAR requires that AOC + U2 equal |L|. AOC in the Figure 5 can be shown to equal Lv in Figure 6 from the derivation
AOC )
ML 1 R
1
∂g
∂x1
dx1 - gj 1′(z1 - xR1 ) )
T,P
gML - gR -
Maximum Partial Area Rule When dealing with multicomponent systems, we introduce here another criterion to establish the equilibrium value of Dij. It is called the maximum partial area rule or MPAR (Shyu, 1995). Figure 5 shows an area (AOC) related to the overall composition (z1). This area is calculated from one equilibrium composition to the input overall composition. The numerical value of AOC equals that of the distance between the Gibbs curve and its tangent line at the fixed input overall composition as illustrated by Figure 6 on which ML denotes a value calculated from the model at the overall composition and TL denotes the tangent-line value for
()
∫xx
ML
g
xβ1 TL
- g - (g R
[ ] gβ - gR -
xR1 R
(z1 - xR1 ) )
- g ) ) gML - gTL ≡ Lv (16)
In the vector procedure, we introduce the known composition [zi] and use an assumed set [D1j] to calculate the g1′ curve. All compositions, except for the overall composition, depend upon the selected [D1j]. Thus, the Gibbs energy varies depending upon the values of [xi], except at the overall composition, which is fixed. All [D1j] must pass through this overall composition. At equilibrium, the tangent-line traces the lowest value of the Gibbs energy (under the constraint that it
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4351
Figure 6. Illustration of MPAR for an activity coefficient model.
not cut the g curve, which implies that it is the global tangent line). In the vector procedure, we compare the Gibbs energy at the overall composition for curves with different [D1j]. The equilibrium [D1j] is that which has the lowest Gibbs energy at the overall composition. When this occurs, [gML - gTL] must be a maximum, and thus, AOC also must be a maximum from eq 16. In addition, we could have based eq 16 upon the β phase rather than upon the R phase without changing the results. The EAR as proposed by Eubank and Hall (1995) is useful for binary systems because Gibbs energies calculated from a model at a fixed x1 are invariant. Thus, we can choose any composition about which to enforce the EAR, the most convenient of which is the middle root. However, for multicomponent systems, only the point of the overall composition has an invariant Gibbs energy for all possible tie lines. Thus, AOC is a more useful construct for the EAR in multicomponent systems. Procedure for MPAR As Applied to a Ternary System. 1. Input P, T, and feed composition, (z1, z2). 2. Choose initial D12 ) -1. 3. Use eq 12 to determine g′. 4. Plot g′ vs x1, and look for equal area (see Figure 1). Get an approximation for the initial guess of gj 1′. 5. If no equal area exists, then pick another D12 and repeat steps 3 and 4. 6. If equal area exists, then find AOC by eq 16.
Figure 7. Comparison for the positive area and partial equal area calculated at different D12. The model is gE ) 2.8x1x2 + 0.9x1x3 + 1.2x2x3. The overall composition is z1 ) 0.4, z2 ) 0.37.
7. Use additional D12 to repeat steps 3 and 6. 8. Determine trend for AOC vs D12 (see Table 1). 9. Find optimum D12 which maximizes AOC. Equilibrium compositions are determined by step 6. Applications of MPAR At equilibrium in multicomponent systems, we have observed that AOC is always a maximum. This observation is especially useful when coupled with interactive graphics and the effective-binary g1′ vs x1 diagrams. For example, a four-component system contains the vectors D12 and D13. A search for the equilibrium value of these two unbounded elements employs a rough grid in (D12, D13) space which presents on the computer screen a number of g1′ vs x1 curves. Many of the curves do not exhibit phase loops, but among those that do, the one which satisfies the constraint of the EAR with a maximum AOC is easy to select. AOC has physical meaning and is easier to use in phase equilibria calculations than areas U or L.
Table 1. Application of MPEAR Using the EOS/MCR Model for the CO2 (1)/Methanol (2)/Propane (3) System at P ) 12.1 bara D12 -2.10 -2.12 -2.14 -2.16 -2.17 -2.18 -2.20 -2.22 -2.24 -2.26 -2.28 -2.30 -2.32 -2.34 -2.36 -2.38 -2.40
xR1 0.0278 0.0288 0.0299 0.0311 0.0317 0.0324 0.0339 0.0355 0.0373 0.0392 0.0412 0.0433 0.0455 0.0477 0.0501 0.0525 0.0549
xβ1 0.4414 0.4399 0.4384 0.4369 0.4362 0.4355 0.4341 0.4327 0.4313 0.4300 0.4286 0.4273 0.4260 0.4247 0.4235 0.4223 0.4210
gR -0.1353 -0.1343 -0.1334 -0.1325 -0.1320 -0.1315 -0.1306 -0.1298 -0.1289 -0.1281 -0.1273 -0.1265 -0.1258 -0.1250 -0.1243 -0.1236 -0.1228
gβ -0.6774 -0.6771 -0.6769 -0.6766 -0.6765 -0.6763 -0.6760 -0.6757 -0.6754 -0.6751 -0.6748 -0.6745 -0.6742 -0.6738 -0.6735 -0.6732 -0.6728
gj 1′ -1.311 -1.320 -1.330 -1.341 -1.346 -1.352 -1.363 -1.375 -1.387 -1.400 -1.413 -1.427 -1.441 -1.456 -1.471 -1.486 -1.502
Ξ
PAb
AOCc
0.259 0.179 0.094 0.005 -0.040 -0.086 -0.181 -0.277 -0.374 -0.471 -0.567 -0.662 -0.754 -0.844 -0.931 -1.014 -1.094
4.514 × 4.513 × 10-1 4.511 × 10-1 4.508 × 10-1 4.507 × 10-1 4.505 × 10-1 4.502 × 10-1 4.497 × 10-1 4.492 × 10-1 4.487 × 10-1 4.481 × 10-1 4.474 × 10-1 4.467 × 10-1 4.460 × 10-1 4.452 × 10-1 4.444 × 10-1 4.435 × 10-1 10-1
3.887 × 10-1 3.889 × 10-1 3.890 × 10-1 3.891 × 10-1 3.890 × 10-1 3.890 × 10-1 3.889 × 10-1 3.887 × 10-1 3.884 × 10-1 3.881 × 10-1 3.877 × 10-1 3.872 × 10-1 3.866 × 10-1 3.860 × 10-1 3.853 × 10-1 3.846 × 10-1 3.838 × 10-1
a The overall composition for this system is z ) 0.28 and z ) 0.37 and gz ) -0.0771. gz is the value of Gibbs energy at the overall 1 2 composition. b PA means positive area. c AOC means area to overall composition. Note: Equilibrium occurs between D12 ) -2.16 and -2.17 when Ξ ) 0.
4352 Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996
CO2 + methanol + propane. The required binary parameters come from Wong and Sandler (1992). Figure 8 is analogous to Figure 7. Figure 9 presents the phase diagram for this system using MPAR with the procedure suggested by Shyu et al. (1995) for locus calculations. Conclusion
Figure 8. Illustration of MPAR using an EOS/MCR model. This diagram represents CO2 (1) + methanol (2) + propane (3) at P ) 12.1 bar. The overall composition is z1 ) 0.28, z2 ) 0.37.
Earlier papers by Eubank and Hall (1995), Shyu et al. (1995), and Hanif et al. (1996) have proposed EAR and demonstrated its efficacy for phase equilibria calculations in binary systems, multicomponent systems, and multiphase systems. This paper contains a major enhancement of EAR with the observation and proof that the partial area from the equilibrium composition up to the overall composition maximizes at equilibrium. We have developed a method which utilizes this fact as a powerful tool for performing equilibrium calculations. MPAR is a rigorous Gibbs minimization procedure which, in combination with the EAR concept, provides a necessary and sufficient condition for equilibrium involving multicomponent systems. The application technique incorporates stability analyses. In addition, the method lends itself well to use of interactive graphics, which greatly simplifies successful implementation of the procedures for efficient calculations. MPAR appears to add significantly to the already existing advantages EAR possesses compared to other methods for calculating phase equilibria. Acknowledgment The National Science Foundation (Grant CTS9317812), the U.S. Department of Energy (Grant DEFG03-93AR14357), Texaco (E & P Technology), Marathon Oil Company, Chevron Oil Company, the Center for Energy and Mineral Resources of Texas A&M University, and the Texas Engineering Experiment Station have provided financial assistance for this effort. Nomenclature
Figure 9. Ternary diagram for the mixture carbon dioxide (1) + methanol (2) + propane (3) at 313 K and 12.1 bar. The circles are the experimental data reported by Galivel-Solastiouk et al. (1986). The continuous line is the predicted locus of the phases in equilibrium.
For simplicity, we use ternary systems to illustrate applications of MPAR. Figure 7 contains two-phase results using a Margules activity coefficient model with overall composition of z1 ) 0.4 and z2 ) 0.37. When approaching equilibrium, Ξ approaches zero and the AOC approaches a maximum. In this example, the equilibrium D12 lies between -1.10 and -1.11. The positive EAR area is not a maximum. The maximum of AOC occurs at the same D12 for which Ξ equals zero. However, the D12 for the maximum of the positive EAR area is different from the equilibrium D12. Table 1 and Figure 8 contain the results of using an EOS/MCR model to demonstrate the use of MPAR. We use the EOS of Peng and Robinson (1976) with the MCR of Wong and Sandler (1992) as applied to the system
AOC ) area related to the overall composition Dij ) composition gradient EAR ) equal area rule EOS ) equation of state f ) fugacity G ) Gibbs energy GEM ) Gibbs energy minimization g ) Gibbs energy or Gibbs energy of mixing (dimensionless) (G/RT) gi′ ) derivative of g with respect to composition Lv ) graphical representation of AOC (see Figure 6) MCR ) mixture combining rule MPAR ) maximum partial area rule N ) total number of components P ) pressure R ) gas constant T ) temperature x ) composition (mole fraction) zi ) overall composition of component i in a system Greek Letters Ξ ) function defined in eq 15 φi ) fugacity coefficient of pure component i φˆ i ) fugacity coefficient of the component i in the mixture γi ) activity coefficient of component i in the mixture
Ind. Eng. Chem. Res., Vol. 35, No. 11, 1996 4353 Superscripts R, β ) denote phases E ) excess function ML ) g at the overall composition calculated from a model TL ) g at the overall composition on the tangent line Subscripts 1,2,3 ) component 1, 2, 3 i,j ) molecular species min, max ) minimum or maximum, respectively
Literature Cited Ammar, M. N.; Renon, H. The Isothermal Flash Problem: New Methods for Phase Split Calculations. AIChE J. 1987, 33, 926939 Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982 (Oct), 731-741. Cairns, B. P.; Furzer, I. A. Multicomponent Three-Phase Azeotropic Distillation. Part II. Phase- Stability and Phase-Splitting Algorithms. Ind. Eng. Chem. Res. 1990, 29, 1364-1382. Eubank, P. T.; Elhassan, A. E.; Barrufet, M. A.; Whiting, W. B. Area Method for Prediction of Fluid Phase Equilibria. Ind. Eng. Chem. Res. 1992, 31, 942-949. Eubank, P. T.; Hall, K. R. An Equal Area Rule and Algorithm For Determining Phase Compositions. AIChE J. 1995, 41, 924-927. Galivel-Solastiouk, F.; Laugier, S.; Richon, D. Vapor-Liquid Equilibrium Data for the Propane-Methanol and PropaneMethanol-Carbon Dioxide System. Fluid Phase Equilib. 1986, 18, 73. 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-999.
Hanif, N. S.; Shyu, G.-S.; Hall, K. R.; Eubank, P. T. Calculation of Multi-Phase Equilibria Using the Equal Area Rule with Application to Hydrocarbon/Water Mixtures. Fluid Phase Equilib. 1996, in press. Heidemann, R. A. Three-Phase Equilibria using Equations of State. AIChE J. 1974, 20, 847-855. Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982a, 9, 1-19. Michelsen, M. L. The Isothermal Flash Problem. Part II. PhaseSplit Calculation. Fluid Phase Equilib. 1982b, 9, 21-40. Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59-64. Shyu, G.-S. Studies of Methods Used in Phase Equilibrium Calculations. Ph.D. Dissertation, Texas A&M University, College Station, TX, 1995. Shyu, G.-S.; Hanif, N. S. M; Hall, K. R.; Eubank, P. T. Equal Area Rule Methods For Ternary Systems. Ind. Eng. Chem. Res. 1995, 34, 4562-4570. Trangenstein, J. A. Minimization of Gibbs Free Energy in Compositional Reservoir Simulation. Soc. Pet. Eng. J. 1985 (Sept), 233-246. Wong, D. S.; Sandler, S. I. A Theoretically Correct New Mixing Rule for Cubic Equations of State for Both Highly and Slightly Non-ideal Mixtures. AIChE J. 1992, 38, 671-680.
Received for review April 10, 1996 Revised manuscript received July 31, 1996 Accepted July 31, 1996X IE960203P
X Abstract published in Advance ACS Abstracts, October 1, 1996.
