Ind. Eng. Chem. Res. 1998, 37, 1483-1489
1483
A Rigorous Mathematical Proof of the Area Method for Phase Stability Ahmed E. Elhassan,† Stefan G. Tsvetkov,‡ Robert J. B. Craven,† Roumiana P. Stateva,‡ and William A. Wakeham*,† IUPAC Thermodynamic Tables Project Centre, Department of Chemical Engineering and Chemical Technology, Imperial College, London SW7 2BY, U.K., and Institute of Chemical Engineering, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria
The paper introduces new developments of the original AREA method. A rigorous mathematical proof that the equilibrium points are the only ones which satisfy the maximum AREA criterion in the case of a two-component, two-phase system is given for the first time. A rigorous proof that the maximum AREA criterion is a necessary but not a sufficient condition for equilibrium in the case of an N-component, two-phase system is given also for the first time in the paper. Two test examples which reinforce the validity of the theoretical results obtained are presented and discussed. Introduction Although the criteria for phase stability have been known since the time of Gibbs (1873), the problem of locating the equilibrium compositions of a multiphase, multicomponent system at specified T and P has continued to exercise thermodynamicists to the present day. The main academic challenge comes from finding the true solution, which for mixtures corresponds with the global Gibbs energy minimum, using complex and nonlinear methods that can predict an unknown number of local minima in addition to the desired global one. The problem is not of academic interest only. Phase stability analysis is of crucial importance in many chemical engineering applications such as three-phase distillation, supercritical fluid extraction, petroleum and reservoir engineering, plant design, and many others. In all these areas large savings in terms of cost can be achieved if the stable phase configuration of a system can be predicted a priori with certainty. As dwindling resources continue to make increasing demands on plant efficiency, the necessity for accurate data, models, and phase equilibria algorithms is increasing. This explains the renewed interest in this problem in recent times; see, for example, Michelsen, 1992; McDonald and Floudas, 1995a, 1997; Sun and Seider, 1995, to name a few. Prior Work There are essentially two approaches that can be taken in order to demonstrate whether an equilibrium solution corresponds with a global minimum of the Gibbs free energy: direct minimization of the Gibbs free energy and minimization of the tangent-plane distance function (the tangent plane criterion). The second approach can be very effective in improving the chances * Corresponding author. Present address: Department of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, London SW7 2BY, U.K. Telephone: +44 171 594 5005. Fax: +44 171 594 8802. E-mail:
[email protected]. † Imperial College. ‡ Bulgarian Academy of Sciences.
of finding the true equilibrium solution, as has been discussed by McDonald and Floudas (1995b). For an N-component system the tangent plane criterion can be expressed as
g(y) )
1
N
∑yj(µj(y) - µj(z)) > 0 RTj)1
(1)
where g(y) is the difference between the reduced Gibbs energy of mixing at the overall composition z (the equation of the tangent hyperplane at z in the Ncomponent thermodynamic phase space) and µj is the chemical potential of species j. The tangent-plane distance function was introduced independently by Baker et al. (1982) and Michelsen (1982) and since then has been used extensively for phase stability analysis. Stability of the existing phase configuration is proven when g(y) g 0 for all y. In fact, as Michelsen pointed out, it is sufficient for stability that g(y) is positive at its stationary points (ysp), in particular its minima, since it will then be positive everywhere (Michelsen, 1982). If, however, the computational methods used cannot find with complete certainty all the stationary points, there is no theoretically based guarantee of stability. In other words, even if no negative solutions are obtained for eq 1, the postulated configuration may still be unstable. There have been several attempts to overcome this problem by introducing methods and algorithms based on different numerical techniques. For example, Sun and Seider (1995) applied a homotopy-continuation method in an attempt to locate all stationary points of the tangent-plane distance function. However, again no theoretical guarantee of identifying all of them can be given, as the formulation of the phase stability problem contains logarithmic terms. Such guarantees can only be made for polynomial systems of equations. Hua et al. (1996) used an interval Newton method, combined with a generalized bisection approach, which also lacks any theoretical guarantee, since interval methods can just enclose and not rigorously identify all stationary points. McDonald and Floudas advocated global optimization techniques which theoretically guar-
S0888-5885(97)00265-0 CCC: $15.00 © 1998 American Chemical Society Published on Web 03/19/1998
1484 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
anteed that the solution corresponding with the minimum Gibbs energy was found for an important but limited number of phase equilibria models (McDonald and Floudas, 1995a-c, 1997). However, in general there appears still to be a need for an efficient generalpurpose method that can perform phase stability calculations which guarantees mathematically the stable solution for any arbitrary equation of state or activity coefficient model. Theoretical Developments In 1992 Eubank et al. introduced a radically different approach to the solution of this problem which involved integration, rather than differentiation of the Gibbs energy surface. For a two-component system the objective function takes the form
(
)
φb - φa b F ) xa φ dx (x1 - xa1) 2 1
∫
xb1
|
G(a1 + t(b1 - a1), a2 + t(b2 - a2))l dt
is the area of the “curvilinear” trapezium with the same vertices and Gibbs energy surface, bounding it between the vertices (G(a1,a2), G(b1,b2)). Let the maximum of the above objective function be found by the integral technique as discussed originally by Eubank et al. (1992). The following will be proved: The points with coordinates a(a1,a2) and b(b1,b2) which maximize the objective function, are equilibrium points. In other words, they satisfy the equipotentiality criterion given by
(
)
(4)
(
)
(4a)
∂G(a1,a2) ∂G(b1,b2) ) ∂a1 ∂b1
(2)
(3)
where
S1 )
|
1
0
and
where geometrically F is the area between the straight line, connecting points xa1 and xb1, and the Gibbs energy curve. When F attains its maximum (positive) value, then xa1 and xb1 represent the equilibrium compositions of component 1 in phases a and b, respectively. This is known as the maximum AREA criterion. In their paper Eubank et al. (1992) showed briefly how the AREA method could be generalized to the N-component, twophase case. This concept was extended further by Elhassan et al. (1996) to the one-component (pure fluid) case and to the determination of dew- and bubble-point temperatures and pressures in addition to the isothermal flash problem. This paper analyzes the results obtained with the AREA method in the case of a two-component, twophase system and in the case of an N-component, twophase system (Eubank et al., 1992; Elhassan et al., 1996). The paper provides for the first time a rigorous mathematical proof that in the latter case equilibrium solutions satisfy the maximum AREA criterion. Furthermore, it provides for the first time a rigorous proof that, in the case of a two-component, two-phase system, the equilibrium points are the only ones which satisfy the maximum AREA criterion. Two-Component, Two-Phase Problem. To begin with we cast the objective function in terms of mole numbers rather than mole fractions for the twocomponent, two-phase system. Let G(θ1,θ2) be the function describing the Gibbs energy surface for a two-component system at constant T and P. Let θ1 and θ2 be the two independent variables. Let a(a1,a2) and b(b1,b2) be two points on the Gibbs energy surface with compositions, given in mole numbers, which maximize the following objective function:
F (a1,a2,b1,b2) ) S2 - S1
|∫
S2 )
|
G(a1,a2) + G(b1,b2) l 2
is the area of the trapezium with vertices (a, b, G(a1,a2), G(b1,b2)) and
∂G(a1,a2) ∂G(b1,b2) ) ∂a2 ∂b2
To prove that, we derive expressions for the difference of the chemical potentials of each component in the two phases namely, ∂G(a1,a2)/∂a1 - ∂G(b1,b2)/∂b1 and ∂G(a1,a2)/ ∂a2 - ∂G(b1,b2)/∂b2. We obtain the above expressions from the corresponding expression of the objective function, taking into consideration that at the maximum the stationarity conditions hold. After some trivial manipulations, the details for which are given as Supporting Information, we end up with
(
)
∂G(a1,a2) ∂G(b1,b2) 1 + l(b1 - a1) 2 ∂a1 ∂b1 ∂G(a1,a2) ∂G(b1,b2) 1 l(b - a2) )0 2 2 ∂a2 ∂b2
(
)
or
(
(b1 - a1)
)
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1 ∂G(a1,a2) ∂G(b1,b2) (b2 - a2) ) 0 (5) ∂a2 ∂b2
(
)
Equation 5 is a corollary of the maximum of the objective function at points a(a1,a2) and b(b1,b2). Hence, a(a1,a2) and b(b1,b2) must satisfy eq 5. It is obvious that a(a1,a2) and b(b1,b2) satisfy eq 5 if they are equilibrium points, for which eqs 4 and 4a hold. This is the first rigorous proof that the maximum AREA criterion corresponds with the equilibrium points, though this result has been used heuristically previously by Eubank et al. (1992). Since eq 5 admits the possibility of other solutions, it is important to demonstrate whether there are other points which maximize the objective function (and hence satisfy eq 5) but do not satisfy eqs 4 and 4a and which ipso facto are not equilibrium points. The following will be proved: For the particular case of a two-component, two-phase system, only the equilibrium points maximize the objective function (satisfy the AREA criterion). We seek a mathematical expression for the characteristics of the points which maximize the objective function and which can be used in addition to eq 5. To
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1485
do that, we analyze the case when the overall composition Z falls between a(a1,a2) and b(b1,b2). We derive expressions, at a constant slope D of the tie-line which passes through the point Z, for the following partial derivatives of the objective function at its maximum:
(
)
∂F (a1,a2,b1,b2) ∂b1
and
a1,a2,D
(
)
∂F (a1,a2,b1,b2) ∂a1
b1,b2,D
(
a1
(
a2
(
b1
( )
(
b2
)
(
(
)
)
∂F (a1,a2,b1,b2) db2 ) 0 (6) ∂b2
)
∂F (a1,a2,b1,b2) ∂b1
)
a1,a2,D
(
(
)
∂F (a1,a2,b1,b2) ∂b1
a1,a2,b1
(
)
∂F (a1,a2,b1,b2) ∂a1
)
b1,b2,D
(
(
a1,a2,b2
∂b2 ∂b1
)
∂F (a1,a2,b1,b2) ∂a1
a2,b1,b2
) 0 (7)
(
)
( )( )
+
∂a2 ∂a1
where f in this case is
) 0 (8)
f)
D
(
(
)
∂F (a1,a2,b1,b2) ∂a1
) )]
∂G(a1,a2) ∂G(b1,b2) ∂a2 ∂b2
-
[( ( 1
) )]
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
∂G(a1,a2) ∂G(b1,b2) b2 ∂a2 ∂b2 Further eqs 9 and 10 can be simplified to
(
)
(13)
(
(14)
)
2b2 - Z2 ∂2G ∂f ∂2G ) + ∂a2 ∂a1∂a2 2b1 - Z1 ∂a 2 2
Using the Gibbs-Duhem equation, which in the present case can be written as
)
x(1 + D2) b 2
) 0 (9)
)
2b2 - Z2 ∂2G ∂f ∂2G ) + ∂a1 ∂a 2 2b1 - Z1 ∂a2∂a1 1
[( (
b1,b2,D
(
Then
x(1 + D ) a ∂G(a1,a2) - ∂G(b1,b2) + 1 2 ∂a1 ∂b1 a2
) )(
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
2b2 - Z2 ∂G(a1,a2) ∂G(b1,b2) (12) 2b1 - Z1 ∂a2 ∂b2
) 2
)
Now:
a2,b1,b2
a1,a2,D
)(
da2 ∂f ∂f )/ da1 ∂a1 ∂a2
The corresponding expressions and relations obtained previously for the different derivatives are inserted into eqs 7 and 8, and after some trivial manipulations, the following results are obtained:
∂F (a1,a2,b1,b2) ∂b1
)
2b2 - Z2 ∂G(a1,a2) ∂G(b1,b2) (11) 2b1 - Z1 ∂a2 ∂b2
D
) ( )
∂F (a1,a2,b1,b2) ∂a2
)
(
-
and
(
) (
and eq 5 becomes
+
) ( )
∂F (a1,a2,b1,b2) ∂b2
) (
∂G(a1,a2) ∂G(b1,b2) ) ∂a1 ∂b1
From eq 6 it follows that
(
)
∂G(a1,a2) ∂G(b1,b2) ) 0 (10a) ∂a2 ∂b2
b2 - a2 2b2 - Z2 Z2 - 2a2 ) ) b1 - a1 2b1 - Z1 Z1 - 2a1
∂F (a1,a2,b1,b2) ∂F (a1,a2,b1,b2) da2 + db1 + ∂a2 ∂b1
(
)
For a fixed overall composition (Z1, Z2), the following holds:
∂F (a1,a2,b1,b2) da1 + ∂a1
(
)
∂G(a1,a2) ∂G(b1,b2) ) 0 (9a) ∂a2 ∂b2
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
where D ) (b2 - a2)/(b1 - a1). The total derivative of F (a1,a2,b1,b2) is
dF (a1,a2,b1,b2) )
)
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
a1
∂2G ∂2G + a )0 2 ∂a1 ∂a2 ∂a12
(15)
a1
∂2G ∂2G )0 + a2 ∂a2 ∂a1 ∂a 2
(16)
and
) 0 (10)
2
then
1486 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
a2 ∂2G ∂2G )2 a1 ∂a1 ∂a2 ∂a1
(15a)
a1 ∂2G ∂2G ) a2 ∂a2 ∂a1 ∂a22
(16a)
a2 a1 ) b2 b1
Substitution of these results into eqs 13 and 14 yields the following:
(
)
(
(
)
∂2G a2 2b2 - Z2 (17) ∂a2 ∂a1 a1 2b1 - Z1
)
2b2 - Z2 a1 ∂2G ∂f ∂2G ) ) ∂a2 ∂a1 ∂a2 2b1 - Z1 a2 ∂a2 ∂a1
(
)
Then:
[( ) ( )]
)
a2 a1
and
a2 ) a1C1
(19)
Analogous manipulation but using
f)
(
)
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
(
)(
)
b2 ) b1C2
(20)
Z2 - 2a2 ∂G(a1,a2) ∂G(b1,b2) Z1 - 2a1 ∂a2 ∂b2
yields the following relation:
Then eqs 9a and 10a can be rewritten as follows:
(
) ( ) (
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1
)
∂G(a1,a2) ∂G(b1,b2) + ∂a1 ∂b1 C2
)
∂G(a1,a2) ∂G(b1,b2) ) 0 (22) ∂a2 ∂b2
Subtracting eq 22 from eq 21 leads to
(C1 - C2)
(
)
∂G(a1,a2) ∂G(b1,b2) )0 ∂a2 ∂b2
(23)
To satisfy eq 23, either C1 - C2 should be zero or
(
)
∂G(a1,a2) ∂G(b1,b2) )0 ∂a2 ∂b2
)
(24)
Applying the result of eq 24 to eq 11, respectively leads to
(
)
(25)
The results of this analysis for the case of a twocomponent, two-phase system demonstrate unequivocally that with the exception of a single point, which is the trivial solution, the only points which maximize the objective function (satisfy the AREA criterion) are the equilibrium points. N-Component, Two-Phase Problem. In the following, the ideas outlined for the case of the twocomponent, two-phase system will be extended to the N-component, two-phase case. Since the methodology is very similar, most of the details concerning the algebra will be omitted, and only the results of the main derivations are given as Supporting Information. Let G(θ1,θ2,...,θN) be the function describing the Gibbs energy surface for an N-component system. Let θ1, θ2, ..., θN be the N independent variables. Let a(a1,a2,...,aN) and b(b,,b2,...,bN) be two points on the Gibbs energy surface with compositions in mole numbers, which maximize the objective function F in a tie-line plane, which passes through the point of the overall composition Z. The objective function is given as follows:
F (a1,a2,...,aN,b1,b2,...,bN) ) S2 - S1
(26)
where
|
S1 )
∂G(a1,a2) ∂G(b1,b2) ) 0 (21) C1 ∂a2 ∂b2
(
(
∂G(a1,a2) ∂G(b1,b2) )0 ∂a1 ∂b1
a1 ∂2G a2 2b2 - Z2 (18) a2 ∂a2 ∂a1 a1 2b1 - Z1
da2 ∂f ∂f )/ da1 ∂a1 ∂a2
which is the trivial solution. Taking this into consideration, then eq 23 can be satisfied if and only if
∂G(a1,a2) ∂G(b1,b2) )0 ∂a2 ∂b2
a2 ∂2G 2b2 - Z2 ∂2G ∂f )+ ) ∂a1 a1 ∂a1 ∂a2 2b1 - Z1 ∂a2 ∂a1 -
Let us suppose that C1 - C2 ) 0. Then it follows that
S2 )
|
|
G(a1,a2,...,aN) + G(b1,b2,...,bN) l 2
(26a)
|
∫01G[a1 + t(b1 - a1), ..., aN + t(bN - aN)]l dt
(26b)
The following will be proved: The points with coordinates a(a1,a2,...,aN) and b(b1,b2,...,bN), which maximize the objective function, are equilibrium points. Extending the derivations, made in the previous section, to the case of an N-component system leads to the following expression: N
(
(bi - ai) ∑ i)1
∂G(a1,a2,...,aN) ∂ai
-
)
∂G(b1,b2,...,bN) ∂bi
)0 (27)
Equation 27 is an analogue of eq 5 and is a corollary of the maximum of the objective function F at points a(a1,a2,...,aN) and b(b1,b2,..,bN) in a tie-line plane, which
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1487
passes through the point of the overall composition. Hence, a(a1,a2,...,aN) and b(b1,b2,...,bN) must satisfy eq 27. It is obvious that a(a1,a2,...,aN) and b(b1,b2,...,bN) satisfy eq 27 if they are equilibrium points for which the following holds:
∂G(a1,a2,...,aN) ∂G(b1,b2,...,bN) ) ∂ai ∂bi
i ) 1, 2, ..., N (28)
As in the case of a two-component, two-phase system, this result is a rigorous proof, given here for the first time. Like eq 5, eq 27 admits the possibility of other solutions. To study this possibility, we have continued further our analysis following a pattern analogous with that for the two-component, two-phase systems. The results obtained demonstrate that, unlike the latter case, there are other solutions (points), in addition to the trivial ones, which maximize the objective function. They satisfy eq 27, but they do not satisfy eq 28, and ipso facto are not equilibrium points. Hence, in the case of an N-component, two-phase system, the maximum AREA criterion gives only a necessary but not a sufficient condition for equilibrium. The insufficiency of the maximum AREA criterion in the case of threecomponent, two-phase systems is also briefly discussed by Elhassan et al. (1996), however, without any mathematical justification. We have not been able to show mathematically what characteristics these points possess. However, the benefit from obtaining (if it is possible at all) a mathematical interpretation of these characteristics is doubtful since they will not provide any useful information about the equilibrium solutions themselves. Test Examples For the case of a two-component, two-phase system, it has been demonstrated analytically that the equilibrium points are the only points satisfying the maximum AREA criterion. Examples of the robustness of the AREA method in locating the equilibrium solution have been given by Eubank et al. (1992) and by Elhassan et al. (1996). Therefore, this section of the paper gives two test examples for the more interesting N-component, two-phase case only. They will serve to reinforce the validity of the theoretical results obtained and to demonstrate explicitly, both numerically and graphically, the dangers of obtaining solutions which satisfy the maximum AREA criterion but are nonequilibrium. Two three-component systems are presented; the first (system A) is modeled with a simple excess Gibbs energy model, and the second (system B) is modeled with the Soave-Redlich-Kwong (SRK) cubic EOS. (1) System A. The Gibbs energy model used is
g ) gid + gE
(29)
where 3
gid )
xi ln xi ∑ i)1
(29a)
and
gE ) a12x1x2 + a13x1x3 + a23x2x3
(29b)
Table 1. Equilibrium and “False” Equilibrium Compositions (Mole Fractions), Phase a z
equilibrium compositions
false equilibrium
z1
z2
xa1
xa2
xa1
xa2
0.4 0.5 0.5 0.5 0.5
0.15 0.05 0.11 0.15 0.5
0.312 62 0.204 52 0.267 09 0.341 30 0.170 72
0.161 03 0.058 81 0.130 17 0.173 65 0.829 28
0.318 66 0.204 72 0.265 76 0.336 38 0.170 72
0.162 75 0.057 38 0.129 00 0.171 57 0.829 28
The following values for the cross coefficients are used: a12 ) 2.3, a13 ) 2.4, and a23 ) 1.8 (see Shyu et al., 1995). Since the derivations leading to eq 27 are given in mole numbers rather than in mole fractions, it is necessary to express eq 27 in terms of mole fractions, in order to apply directly the Gibbs energy model given by eqs 29. The corresponding transformations are given in Derivations III and IV of the Supporting Information, respectively. In order to broaden the perspective and give a numerical interpretation of the rigorous results provided in the previous sections of the paper, calculations are performed for system A, applying the following representative scheme: 1. Choose an overall composition z. 2. Fix the slope D of the tie line and determine the compositions which maximize the objective function. 3. Check whether these compositions satisfy eqs 28 as well. If the compositions satisfy eq 28, they are equilibrium. If the compositions do not satisfy eq 28, go to step 4. 4. Change the slope of the tie line and go back to step 2. As pointed out in the previous paragraph, the maximization of the objective function is carried out in a tieline plane, which passes through the point of the overall composition. For an N-component system there is an infinite number of planes, all of which pass through z, but only one of these coincides with the equilibrium tie line. Taking into consideration that the maximum AREA criterion (eq 27) will be satisfied in each of the above planes and that the orientation of the equilibrium tie-line plane is not known a priori, it is necessary to check whether the compositions which maximize the area in a given plane (at a fixed D) are equilibrium, i.e., whether they satisfy eq 28. The representative scheme follows this pattern and determines the slope of the equilibrium tie line and the equilibrium compositions in an iterative manner. Furthermore, the compositions which correspond with the maximum area, among all maxima (obtained at different D), are determined as well. The former “global maximum area” is determined at some DM * DE. These compositions satisfy eq 27 but do not satisfy eq 28 and hence are not equilibrium points. That is why we call them “false” equilibrium points. Some of the numerical results obtained are presented in Tables 1 and 2. For a given overall composition z, the tables give values of the equilibrium compositions and values of the compositions which correspond with the global maximum area. It is demonstrated numerically that when the third component is not present (i.e., at the basis of the binodal where x3 ) 0), the compositions which correspond with this maximum among all maxima and the compositions which correspond with the equilibrium are identical. When a graph of the true and false equilibrium points is prepared, it shows that the false equilibrium points
1488 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
Figure 1. Diagram of the maximum value of F, calculated at constant slope D of the tie line. The overall compositions are z ) 0.48, 0.14, and 0.38.
Figure 2. Variation of the chemical potential of component 2 in phases a and b as a function of D, calculated at the maximum of the objective function F. The overall compositions are z ) 0.48, 0.14, and 0.38. Table 2. Equilibrium and “False” Equilibrium Compositions (Mole Fractions), Phase b z
equilibrium compositions
z1
z2
xb1
0.4 0.5 0.5 0.5 0.5
0.15 0.05 0.11 0.15 0.5
0.626 53 0.764 61 0.676 22 0.598 45 0.829 28
false equilibrium
xb2
xb1
xb2
0.121 40 0.042 11 0.094 74 0.135 33 0.170 72
0.641 93 0.765 46 0.674 12 0.572 77 0.829 28
0.112 08 0.041 57 0.095 87 0.137 77 0.170 72
lie totally inside the binodal, thus forming a “binodal of false equilibrium points”. The binodal of the true equilibrium points and the binodal of the false equilibrium points coincide at x3 ) 0. For example, Figure 1 shows the values for the maximum of the objective function F obtained at different values of the slope D of the tie line at a specified z. Figure 2 shows the variation of the chemical potential of component 2 in phase a and phase b, at the specified z, as a function of D. As discussed above, the maximum of all maxima area does not correspond with the equality of the chemical potentials of the components in the system, or in other words the equilibrium occurs at DE * DM.
Figure 3. Diagram of the maximum value of F, as a function of the slope D of the tie line, for the methane + carbon dioxide + hydrogen sulfide system, using the SRK cubic EOS. T ) 212.4 K, P ) 6.06 MPa. The overall compositions are z ) 0.73, 0.07, and 0.20.
Figure 4. Variation of the chemical potential of carbon dioxide in the liquid and vapor phases for the methane + carbon dioxide + hydrogen sulfide system, at T ) 212.4 K and P ) 6.06 MPa. The overall compositions are z ) 0.73, 0.07, and 0.20. The chemical potentials are evaluated at the maximum value of the objective function F for a given slope D of the tie line. (s) liquid phase. (-‚-) vapor phase.
(2) System B. System B, namely, methane + carbon dioxide + hydrogen sulfide, is modeled with the SRK cubic EOS. Calculations, analogous to those performed for system A, are run for system B at T ) 212.4 K and P ) 6.06 MPa. The system exhibits a very narrow three-phase region in the temperature plane, but for the purposes of illustrating the above theoretical developments, the overall compositions studied are outside the three-phase region. The binary interaction parameters used by the SRK CEOS are k12 ) 0.0933, k13 ) 0.0823, and k23 ) 0.0989. Figure 3 gives the maximum of F as a function of the slope D and at z ) 0.73, 0.07, and 0.20. Figure 4 shows the variation of the chemical potential of carbon dioxide in the liquid and vapor phases as a function of D and at the same overall composition. As in the case of system A, the maximum of all maxima of the objective function at a specified overall composition does not correspond with the equilibrium compositions.
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1489
Conclusions
Greek Letters
A rigorous mathematical proof that the maximum AREA criterion gives a necessary and sufficient condition for equilibrium in the case of a two-component, twophase system is presented for the first time in this paper. Thus, the application of the AREA criterion identifies the system as stable or unstable and determines the equilibrium solution. A rigorous mathematical proof that the maximum AREA criterion is a necessary but not a sufficient condition for equilibrium in the case of an N-component, two-phase system is also given for the first time in this paper. The mathematical derivations, and the conclusions arising from them, lead to a wide variety of possible numerical implementations of the AREA criterion to determine an equilibrium solution. However, in order to find the global solution among all possible equilibria, a combined algorithm, which incorporates the integral scheme, will be necessary. Moreover, before such an algorithm can be successfully implemented, a theoretical proof that the procedure provides the global minimum of the Gibbs energy will be necessary. This task is beyond the scope of the present paper and will be the subject of further study.
φ ) reduced Gibbs energy of mixing µi ) chemical potential of component i, J θ ) independent variable Superscripts a, b ) phase identification id ) ideal E ) excess Subscripts a, b ) phase identification E ) corresponding to the equilibrium M ) corresponding to the global maximum area s.p. ) stationary point T ) total
Supporting Information Available: Details of derivations for the cases of two-component, two-phase and N-component, two-phase systems and transformations of mole numbers to mole fractions (12 pages). Ordering information is given on any current masthead page. Literature Cited
Acknowledgment The authors thank the United Kingdom Department of Trade and Industry for the major part of the funding of the Project Centre and the British Council in Sofia, Bulgaria, for funding the travelling expenses of R.P.S. and S.G.T.’s visit to the U.K. In addition, R.P.S. and S.G.T. acknowledge the financial support of the National Science Fund of Bulgaria, under Grant 517/95. Nomenclature ai ) number of moles of component i in phase a aij ) cross coefficients in the Gibbs energy model bi ) number of moles of component i in phase b C1, C2 ) constants D ) slope of the tie line f ) function, eq 12 F ) objective function, eq 2 F ) objective function, eqs 3 and 26 g ) molar Gibbs energy, J mol-1 g(y) ) tangent-plane distance function, eq 1 G ) total Gibbs energy, J kij ) binary interaction coefficient, SRK CEOS l ) length of segment N ) number of components in the system n ) number of moles P ) pressure, MPa R ) universal gas constant, J K-1 mol-1 S1 ) area of the trapezium S2 ) area of the “curvilinear” trapezium T ) temperature, K t ) parameter xi ) mole fraction of component i in the liquid phase yi ) mole fraction of component i in the vapor phase zi ) mole fraction of component i in the feed Zi ) number of moles of component i in the feed
Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982, 22, 731. Elhassan, A. E.; Lopez, A. A.; Craven, R. J. B. Solution of the Multiphase Problem for Pure Component, Binary and Ternary Systems Using the Area Method. J. Chem. Soc., Faraday Trans. 1996, 92, 4419. 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. 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. Hua, J. Z.; Brenneke, J. F.; Stadther, M. A. Reliable Prediction of Phase Stability Using an Interval Newton Method. Fluid Phase Equilib. 1996, 116, 52. McDonald, C. M.; Floudas, C. A. Global Optimisation and Analysis for the Gibbs Free Energy Function Using the UNIFAC, Wilson and ASOG Equation. Ind. Eng. Chem. Res. 1995a, 34, 1674. McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase Stability Problem. AIChE J. 1995b, 41, 1798. McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase and Chemical Equilibrium Problem. Application to the NRTL Equation. Comput. Chem. Eng. 1995c, 19, 111. McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1997, 21, 1. Michelsen, M. L. The Isothermal Flash Problem. Part I: Stability. Fluid Phase Equilib. 1982, 9, 1. Shyu, G.-S.; Hanif, N. S. M.; Alvarado, J. F. J.; Eubank, P. T. Equal Area Rule Methods for Ternary Systems. Ind. Eng. Chem. Res. 1995, 34, 4562. Sun, A. C.; Seider, W. D. Homotopy Continuation Method for Stability Analysis in the Global Minimisation of the Gibbs Free Energy. Fluid Phase Equilib. 1995, 103, 213.
Received for review April 4, 1997 Revised manuscript received November 18, 1997 Accepted November 26, 1997 IE970265V
