The Area Method for Phase Stability Analysis Revisited: Further

Ind. Eng. Chem. Res. 2007, 46, 1611-1631

1611

The Area Method for Phase Stability Analysis Revisited: Further Developments. Formulation in Terms of the Convex Envelope of Thermodynamic Surfaces Janos Balogh,*,† Robert J. B. Craven,‡ and Roumiana P. Stateva§ Department of Computer Science, Juhasz Gyula Teachers’ Training College, UniVersity of Szeged, Szeged H-6701, Hungary, ESDU International plc, London N1 6UA, UK, and Institute of Chemical Engineering, Bulgarian Academy of Sciences, Sofia, Bulgaria

The Area method for phase stability analysis is re-examined from the viewpoint of the topology of thermodynamic surfaces. A new generalization of the Area method objective function and the maximum area criterion is proposed in terms of the convex envelope of thermodynamic surfaces. Several special cases for two-component, single- and two-phase systems are analyzed and the results generalized to hypersurfaces, characteristic of the Gibbs energy of multicomponent, multiphase systems. It is shown conclusively that the area (or more generally hypervolume) between the free energy hypersurface and a hyperplane that intersects the hypersurface is maximal when the hyperplane coincides with the convex envelope of the hypersurface. Extension of the maximum area criterion to four basic thermodynamic surfaces for pure fluids and mixtures is discussed. The authors believe the approach adopted here sets the Area method on secure theoretical foundations, answers criticisms concerning its validity, and will stimulate fresh research. Introduction Thermodynamic characterization of multiphase, multicomponent systems continues to be of great importance in many branches of the chemical processing industry. While the thermodynamic principles involved have been well understood since the time of Gibbs, the mathematical computations involved are time-consuming and nontrivial. Before the development of powerful and fast computers, the necessary mathematical computations were only possible for the simplest thermodynamic models and systems containing few components. However, in the last 25 years, the development of fast and powerful computers has led to renewed interest in the subject. It is now possible to use sophisticated numerical procedures for determining the stable phase configuration at equilibrium of complex systems and better models for estimating the thermodynamic properties of the phases. Thermodynamic characterization can be split into three stages: Phase stability analysis determines whether or not the system is stable as a single, homogeneous phase. If the system is unstable, the second step is to determine the stable phase configuration and the component distribution among the equilibrium phases. If desired, the final step is to determine the thermodynamic properties for all the equilibrium phases, usually from some prediction scheme, which may be different for each phase. This paper is primarily concerned with the first of these. Many effective procedures for determining the stable phase configuration of a thermodynamic system rely on the geometrical concept of the tangent plane. As its name implies, this is a hyperplane (or line for the simplest systems) that touches a thermodynamic hypersurface at one or more points. Clearly for a single-phase system, there is infinity of such planes, but the occurrence of multiple phases restricts their number for thermodynamically stable systems. These restrictions were investigated in detail by Baker et al.1 and by Michelsen,2 who introduced the concept of the tangent plane distance function. This function and variants of it3,4 have found widespread use in phase stability analysis. * To whom correspondence should be addressed. Tel.: +36-62546080. Fax: +36-62-420953. E-mail: [email protected]. † University of Szeged. ‡ ESDU International plc. § Bulgarian Academy of Sciences.

Perhaps the most important property of any tangent plane is that for stability, the plane should not cross the thermodynamic surface at any other point on the surface. This condition imposes severe restrictions on surface topology, namely, that a stable surface must either be entirely concave or entirely convex throughout. It is known that, at equilibrium, the stable phase configuration corresponds with the minimum Gibbs energy at a given temperature and pressure and the minimum Helmholtz energy at a given temperature and volume. The requirements of the minimum free energy and the “noncrossing rule” taken together imply that these free energy surfaces must be convex throughout their respective domains. Recently, several papers have been published that exploit the properties of convex functions in phase stability analysis.5-10 In 1992, Eubank et al.11 introduced a radically different approach to phase stability analysis. They investigated the properties of the area between the thermodynamic surface and the tangent plane and concluded that, for a two-component, twophase system, this area was maximal. In terms of the Gibbs energy expressed as a function of composition at a given pressure and temperature, the area is given by

(xb1 - xa1) (G(xb1) + G(xa1)) II ) a G(x) dx x1 2

∫

x1b

(1)

where, geometrically, II is the area between the straight line connecting points xa1, xb1 and the Gibbs energy surface. When II is maximal (and positive in value), then xa1 and xb1 are the equilibrium compositions of component 1 in phases a and b, respectively. This condition is known as the maximum area criterion, and in their paper, Eubank et al.11 showed briefly how the Area method could be generalized to the Nc-component, two-phase case. The Area method concept was extended further by Elhassan et al.12 to the one-component (pure fluid) case and to the determination of dew- and bubble-point temperatures and pressures for ethane-n-octane mixtures. Formal proofs that the maximum area criterion corresponds with global phase stability are however lacking. It is quite straightforward to show that the necessary conditions for phase stability are satisfied when eq 1 is maximal,12 but it is more difficult to prove that the sufficient criteria are met. In a further contribution, Elhassan et al.13 and later Elhassan et al.14 sought to provide such a proof. However, their methods relied very

10.1021/ie060949f CCC: $37.00 © 2007 American Chemical Society Published on Web 02/01/2007

1612

Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007

Figure 1. Simple example of a convex set.

Figure 3. Convex hull (Clim) and convex envelope (conv(S)) of a nonconvex set S.

complex, multicomponent, multiphase cases where the Gibbs energy is in general a hypersurface. First, properties of convex sets and functions will be examined and the concept of the convex envelope will be related to thermodynamic stability. Second, the connection between the maximum area criterion and the convex envelope will be established for cases characteristic of thermodynamic systems in general. Particular situations will be examined and for each one of these, a rigorous proof that a tangent line to a curve maximizes the area between the tangent and the curve will be given. All essential elements of the discussion will be based on considerations of surface topology. Not only does this lead to a considerable simplification of the mathematics but also to conclusions that are quite general and do not depend on any thermodynamic assumptions. Finally, we shall apply these conclusions to a variety of thermodynamic surfaces for pure fluids and multicomponent mixtures. Figure 2. Simple example of a nonconvex set with elements P1 and P2 which do not satisfy eq 2.

heavily on the differentiation of functions of the form of eq 1, and as derivatives will, in general, only identify local stationary or saddle points and not the global minimum of the free energy, these procedures were open to criticism. Doubts were expressed by Barton and Hwang15 about their validity and, more importantly, about the existence of any explicit relationship between phase stability and the maximum area criterion. For this reason perhaps, the Area method has not been applied extensively to phase stability problems, apart from the continuing research by those at Texas A&M University (refs 16-19 to mention just a few). From the above, it is clear that there are certain points that require further clarification. While it is our contention that previous work has demonstrated unequivocally a clear connection between phase stability and the maximum area criterion, it is unclear whether the maximum area criterion is synonymous with global free energy minimization. To this extent, at least, the criticisms of Barton and Hwang15 are valid. It is possible, for example that the maximum area criterion is only capable of identifying local and not global free energy minima. Also, there may be situations where application of the maximum area criterion does not produce correct answers without the imposition of other constraints. These considerations have led us to examine the Area method again. For clarity, the topological arguments presented here are for the two-component case, where the Gibbs energy surface is a curve. Subsequently, the arguments are extended to more

Convex Hull and Envelope of a Function Let us consider a set S, defined as follows: S is convex if and only if for any pair of elements P1 and P2 belonging to S and for any λ ∈ [0,1], the following holds:

λP1 + (1 - λ)P2∈ S

(2)

Figure 1 is a simple example of a convex set. However, in general, most sets are not convex as exemplified by Figure 2. In Figure 2, a pair of elements P1 and P2 can be selected as shown such that eq 2 is not satisfied and by definition the set is not convex. Now consider Figure 3. The set C is convex and it contains a set S, which is not convex, as S contains pairs of points that do not satisfy eq 2. It is clear from Figure 3 that a convex set C1, a subset of the convex set C, but which contains the nonconvex set S, can be selected, such that:

S ⊆ C1 ⊆ C

(3)

In a similar way, we can select a third convex set C2, which is a member of C1 and contains S. In the limit, a set Clim can be selected, which is the smallest possible convex set that contains the nonconvex set S, so that finally

S ⊆ Clim ⊆ ... ⊆ C2 ⊆ C1 ⊆ C

(4)

Clim is called the convex hull of the set S. It is the smallest convex set that can be constructed which contains set S. In general, S may be either convex or not. However, in the special

Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1613

Figure 4. Relationship between the convex envelope (conv(G(x))) and the Gibbs energy curve for a two-component, two-phase system. The hatched area denotes the epigraph of the curve.

case of a convex set, the convex hull of the set is identical to the set itself. The conVex enVelope of the set S is the boundary of its convex hull. In this paper, we denote it by conv(S), and it is shown in Figure 3 by the bold line. The area bounded by and contained within the bold line represents the convex hull of S. A convex hull of a function in an interval, where the function is defined, is the convex hull of its epigraph. The epigraph of a function F(x), denoted epi(F(x)), represents the set of points on the same interval, lying on or above the function in the interval. By analogy, the conVex enVelope of a function is the boundary function of its convex hull. In other words, the convex envelope of a function is the greatest convex function that underestimates the function in the same interval. Further discussion follows in the next paragraph.

hyperplanes. For compositions expressed in terms of mole fractions, a composition hyperspace of dimension Nc - 1 must be considered where Nc is the number of system components. Points on or above the multidimensional hypersurface, the epigraph of the hypersurface, may in general be regarded as a nonconvex set S. The convex envelope of S consists of two parts: tangent hyperplanes to the hypersurface and the hypersurface itself. Any tangent hyperplanes must lie everywhere on or below the surface. If this were not the case, the hyperplane would cut the hypersurface at some point, thereby connecting pairs of elements of S for which the convexity condition, eq 2 is invalid. However, the condition that any tangent hyperplane must be everywhere touching or below the surface is the condition for global stability as given by Baker et al.1 Therefore, the convex envelope of a multidimensional Gibbs hypersurface is identical to the surface of dissipated energy as defined by Gibbs.20 By analogy with the two-component case, regions where the surface coincides with the convex envelope are in the single phase. Multiphase regions lie on a tangent hyperplane G′(x), which gives the overall Gibbs energy of the system of composition x at constant temperature and pressure. On the tangent hyperplane, the system splits or dissipates into P phases in such a way that the equilibrium compositions are at points where G′(x) is tangential to G(x) and the chemical potentials of any component in the phases at equilibrium are equal. The above discussion leads to the conclusion that for a multidimensional Gibbs surface G(x), the curve of global stability is the convex envelope of G(x). This conclusion can be summarized with the following notation:

Gstab(x) ) conv(G(x))

(5)

Since conv(G(x)) lies everywhere on or below G(x), it follows that

G(x) - conv(G(x)) g 0

(6)

Relationship between the Gibbs Energy Function of the Surface and Its Convex Envelope

Relationship between the Convex Envelope and the Maximum Area Criterion

Figure 4 is a typical curve for the Gibbs energy of a twophase two-component system at constant temperature and pressure. It is in general a continuous function in the interval [0,1] and has one maximal value. Points along and above the curve, the epigraph of the curve, can be viewed as a nonconvex set S. The convex envelope of S is indicated by the bold line. However, the bold line also represents the “curve” of dissipated energy, a term introduced originally by Gibbs as the “surface of dissipated energy” for multicomponent systems or the curve of dissipated energy for two-component systems.20 In other words, for a two-phase, two-component system, the surface of dissipated energy, which corresponds with the global minimum value of the total Gibbs energy, is the convex envelope of the Gibbs energy curve. Any point on conv(G(x))sthe convex envelope of G(x)scorresponds with the Gibbs energy of the globally stable phase configuration for a system of overall composition x. Where G(x) and conv(G(x)) coincide, the system is stable in the single phase. Elsewhere, the system splits or dissipates into two phases in such a way that the equilibrium compositions are given by two points where a line G′(x) is tangential to G(x) at both points and the chemical potentials of any component in the two phases are equal. Figure 4 makes clear that the common tangent line completes the convex envelope of G(x). The extension of the above to multicomponent, multiphase systems is straightforward. In the general case, the Gibbs energy curve becomes a hypersurface and tangent lines become

The connection between the maximum area criterion, eq 1, and global phase stability can be established through the concept of the convex envelope of G(x) as follows. To proceed further, it is helpful to rewrite eq 1 as

II )

∫uV(G(x) - L′′(x)) dx

(7)

or, more succinctly, as

II )

G(x) dG dx ∫uV ∫L′′(x)

(8)

where L′′(x) is a straight line bounded by the elements u and V which lie on G(x) (eq 9)

L′′(x) )

(

)

G(V) - G(u) VG(u) - uG(V) x+ V-u V-u

(9)

A straightforward integration of eq 7 reproduces eq 1, where xa1, xb1 replace u, V, respectively. In eq 7, 0 e u e V e 1 and x satisfies 0 e x e 1. Therefore, the algebraic sign of II is determined by that of the integrand. If L′′(x) is set to conv((G(x)), it follows immediately from eq 6 that II must be positive or zero. If it can be shown that this setting also generates the global maximum of II, the connection between the convex envelope of the Gibbs energy curve and the maximum area criterion is proved.

1614


follows directly from the definition of convexity that

L′′(x) g G(x)

(10)

for any arbitrary x ∈ [u,V]. Hence, as u e V

II )

Figure 5. Area between the convex envelope of a curve and the curve itself.

∫uV(G(x) - L′′(x)) dx e ∫uV0 dx ) 0

(11)

As the integrand is e 0, II is in general negative. The equality in eq 11 holds if and only if u ) V. In this case, the lower and upper limits of integration coincide and the line L′′(x) becomes tangential to G(x). As there is no restriction on the position of u and V in the interval [0,1] other than 0 e u e V e 1, it follows that the objective function, eq 7, is maximal for all tangents to G(x) in the interval [0,1], and therefore, the envelope of all tangents to G(x) in the same interval represents a curve of maxima for II. However, the envelope of all tangents to the curve G(x) is also the convex envelope of G(x), and so, in this case at least, the maximal condition for eq 7 can be written as follows:

IImax )

∫01(G(x) - conv(G(x))) dx ) 0

(12)

In other words, the objective function, eq 7, is maximal over the whole composition range, when the line L′′(x) is replaced by the convex envelope of G(x). The maximal value of II is zero, and for any other line within the epigraph of G(x), the value of II is always less than zero. It is now straightforward to generalize this result to multidimensional hypersurfaces with only one local minimizer. Such surfaces are typical of the Gibbs energy hypersurface for multicomponent single-phase systems. In this case, the curve G(x) is a convex function of the vector space composition coordinates x. The line L′′(x) becomes a hyperplane, which in general lies within the epigraph of G(x). Equation 7 becomes Figure 6. Typical Gibbs energy curve for a two-component, one-phase system.

Let us suppose that G(x) is continuous in the interval [0,1], it has m local minimizers (denoted by min1, ..., minm) and is convex in the subintervals [0,min1] and [minm,1]. (See Figure 5). In order to prove rigorously that the area enclosed between the tangent lines to the curve G(x) and the curve maximizes the area defined by eq 7, the following typical cases for the Gibbs energy of two-component systems will be considered in detail and the results will be generalized. Case 1: Function G(x) Has One Local Minimizer in the Interval [0,1]. The one local minimizer case (Figure 6) is very simple and almost trivial, but nevertheless, it is instructive to consider it in some detail. Figure 6 is a typical curve for the Gibbs energy of a single-phase, two-component system. Furthermore, it is possible to generalize the results from this example to multicomponent, single-phase systems and to other, more complex cases. Statement. If G(x) is a convex function with just one local minimizer (min1) in the interval [0,1] and L′′(x) is a straight line within the epigraph of G(x) and touches G(x) at two points u and V, such that 0 e u e V e 1, then the value of eq 7 is maximal if and only if L′′(x) is a tangent line, L′(x), of G(x) at any point x ∈ [0,1]. Proof. Since the function G(x) is convex in the interval [0,1], then, for any line within the epigraph of the curve G(x) that touches G(x) at two points, denoted by u and V (u e V), it

II )

∫∫∫(G(x) - L′′(x)) dx

(13)

x

and since L′′(x) g G(x) within the epigraph of G(x), the hypervolume defined through eq 13 is in general negative. When L′′(x) ) G(x), the lower and upper limits of integration coincide and II vanishes. In this case, the hyperplane becomes a tangent hyperplane to the hypersurface and by similar reasoning to the two-component case, the envelope of all tangents to G(x) in the same interval represents a surface of maxima for II. However, the envelope of all tangents to the surface is simply the convex envelope of G(x) so that finally

IImax )

∫∫∫(G(x) - conv(G(x))) dx ) 0

(14)

x

represents the surface of maxima of II for multicomponent single-phase systems where the integration is taken over the whole of the vector composition space. Case 2: Function G(x) Has One Local Minimizer and One Local Maximizer in the Interval [0,1]. Case 2 is illustrated in Figure 7. Figure 7 is a typical curve for the Gibbs energy of a two-phase, two-component system when one of the phases is a solid. The Gibbs energy of the solid component, component 1, is a point at (1,G(1)). In this example, there is a discontinuity in the Gibbs energy curve between the fluid phase and the pure solid at x ) 1. As the Gibbs energy of the pure solid lies below that of the (hypothetical) solution at G(1), the Gibbs energy curve is not convex. Therefore, the system as exemplified by


Further insight into the nature of the objective function II in this case can be found by splitting the integral of eq 7 into two parts so that

∫u1(G(x) - L′′(x)) dx ) ∫uV(G(x) - L′′(x)) dx + ∫V1(G(x) - L′′(x)) dx

Figure 7. Typical Gibbs energy curve for a two-component, two-phase system, where one of the phases is a pure solid.

this curve is not stable in the single phase for all x, and unlike case 1, the convex envelope of G(x) will not coincide everywhere with the function itself. This is stated more succinctly as follows: Let the function G(x) be continuous in the interval [0,1], discontinuous at the point x ) 1, G(x ) 1) with one local minimizer, and convex in the subintervals [0,min1] and [min1,1]. Also let

G(1) < lim G(x) xv1

(15)

where, x v 1 denotes that x tends to 1 from the direction 0 f 1. Statement. If G(x) is a function with just one local minimizer (min1) and one local maximizer (max1) in the interval [0,1] and L′′(x) is a straight line within the epigraph of conv(G(x)) and crosses (or touches) G(x) at two points u and V such that 0 e u e V e 1, then the value of eq 7 is maximal if L′′(x) is a line, L′(x), which is tangential to G(x) at some point x ∈ [0,1] and passes through the point (1,G(1)); i.e., L′(x) is a supporting line of the curve G(x), x ∈ [0,1], from the point G(1). Proof. A proof of the above can be easily derived from consideration of the following three integrals:

∫u1(G(x) - L′′(x)) dx

(16)

∫u1(G(x) - conv(G(x))) dx

(17)

∫u1(conv(G(x)) - L′′(x)) dx

(18)

II1 ) II2 ) II3 )

where V is such that t e V e 1 and [t,G(t)] is the point at which the supporting line is tangential to G(x). The first integral on the right-hand side is analogous to the single-phase example, case 1. In general, it is negative but vanishes when u tends to V, where L′′(x) becomes a tangent line, L′(x), to G(x). All tangents in the interval [u,V] are acceptable solutions, but the solutions of interest for this example are those that lie in the interval [0,t] where the tangent line through t also passes through the point [1,G(1)]. The second integral on the right-hand side is in general positive since the crossing point at V reverses the sign of G(x) - L′′(x). This integral is maximal when G(x) - L′′(x) is maximal, and the interval [V,1] is maximal simultaneously. The first condition is met for the supporting line L′(x) connecting [t,G(t)] and [1,G(1)], the convex envelope of G(x) in the interval [t,1]. The second condition is also met since the convexity condition, ∂2G(x)/∂x2 g 0, implies that either side of [t,G(t)], G(x) will lie everywhere above the tangent line at t. Lines drawn from the point [1,G(1)] through points either side of [t,G(t)] will therefore lie above the supporting line leading to a reduction in the value of the integral. Finally, therefore as for case 1, the objective function II, eq 7 for case 2, is maximal along the convex envelope of G(x) so that eq 12 becomes

IImax )

∫01(G(x) - conv(G(x))) dx ) ∫t 1(G(x) - conv(G(x))) dx > 0

(21)

However, unlike case 1 where the line L′′(x) is everywhere confined within the epigraph of G(x), the integration limits for IImax do not coincide. Consequently, the maximal value of II, IImax, is always greater than zero (see Figure 7). Generalization to hypersurfaces follows the same line of reasoning as that adopted for case 1. It is only necessary to replace x with a composition vector x and set the limits of integration to the whole of the composition space. This leads immediately to eq 22:

IImax )

∫∫∫(G(x) - conv(G(x))) dx > 0

(22)

x

From the above, it follows immediately that

II3 ) II1 - II2

(20)

(19)

Since the line L′′(x) lies within the epigraph of the convex function conv(G(x)), then from the results of case 1, II3 is in general negative, with a maximal value of zero for the special case L′′(x) ) conv(G(x)), a result in agreement with eq 19. However, since L′′(x) is by definition g conv(G(x)) over the interval [u,1], then the integrand of eq 16 is e that of eq 17 in the same interval and II1 e II2. Therefore, the integral II2 is the maximal value of II1, and at this maximum, the function L′′(x) coincides with the convex envelope of G(x). From case 1, it follows immediately that this result holds over the whole range of x [0,1].

where the integration is taken over a region of composition space for which the integrand is positive. The integral over regions for which the integrand is negative vanishes when II attains its maximal value. Case 3: Function G(x) Has Two Local Minimizers and One Local Maximizer in the Interval [0,1]. Case 3 is illustrated in Figure 8. Figure 8 is a typical curve for the Gibbs energy of a two-phase, two-component system when both phases contain appreciable amounts of both components. In this example, the function G(x) is not convex as ∂2G(x)/∂x2 < 0 at the maximum and the system is not stable in the single phase for all x. Therefore, the convex envelope of G(x) will not coincide everywhere with the function itself. Statement. If G(x) is a function with two local minimizers (min1, min2) and one local maximizer (max1) in the interval [0,1], and L′′(x) is a straight line within the epigraph of conv(G(x)) and crosses (or touches) G(x) at four points u1, V1,

1616


Figure 8. Typical Gibbs energy curve for a two-component, two-phase system, where both phases contain significant amounts of both components.

u2, and V2 such that 0 e u1 e V1 e u2 e V2 e 1, then the value of the objective function integral, eq 7, is maximal when L′′(x) is a common tangent, L′(x), to G(x) at two points within the interval [0,1]. Proof. A proof of the above follows exactly that for case 2. Consideration of the three integrals

∫uV (G(x) - L′′(x)) dx

(23)

II2 )

∫uV (G(x) - conv(G(x))) dx

(24)

II3 )

∫uV (conv(G(x)) - L′′(x)) dx

(25)

II1 )

2

1

2

1

2

1

where L′′(x) lies within the epigraph of conv(G(x)) leads to the same conclusions as those found for case 2, namely, that II3 has a maximal value of zero when L′′(x) ) conv(G(x)) and II2 is the maximal value of II1 for the same condition. By reasoning similar to that for cases 1 and 2, it follows immediately that this result holds when the integration limits for eqs 23-25 extend over the whole of the composition range [0,1]. Further insight into the nature of the objective function II in this case can be found in a manner analogous to that used for case 2. The objective function integral, eq 7, is split into three parts so that

∫uV (G(x) - L′′(x)) dx ) ∫uV (G(x) - L′′(x)) dx + ∫Vu (G(x) - L′′(x)) dx + ∫uV (G(x) - L′′(x)) dx 2

1

1

1

2

and the interval [V1,u2] is maximal simultaneously. The first condition is met for the common tangent L′(x) connecting [t1,G(t1)] and [t2,G(t2)], the convex envelope of G(x) in the interval [t1,t2]. The second condition is also met since the convexity condition, ∂2G(x)/∂x2 g 0, implies that either side of [t1,G(t1)] and [t2,G(t2)], G(x) will lie everywhere above the common tangent. Therefore, lines connecting such points will lie above L′(x), leading to a reduction in the value of the integral. Therefore as for cases 1 and 2, the objective function II, eq 7 for case 3, is maximal along the convex envelope of G(x) so that eq 12 becomes the following:

IImax )

∫01(G(x) - conv(G(x))) dx ) ∫t t (G(x) - conv(G(x))) dx > 0 2

(26)

(27)

1

However, unlike case 1 where the line L′′(x) is everywhere confined within the epigraph of G(x), the integration limits for IImax do not coincide. Consequently, the maximal value of II, IImax, is always greater than zero (see Figure 9). Generalization to hypersurfaces follows the same line of reasoning as that adopted for case 1. It is only necessary to replace x with a composition vector x and set the limits of integration to the whole of the composition space. This leads immediately to eq 28 identical in form to eqs 22 and 14:

IImax )

1

2

Figure 9. Alternating signs of the areas of segments between a nonconvex curve and a line that cuts the curve in four places: (+) regions with positive area; (-) regions with negative area.

∫∫∫(G(x) - conv(G(x))) dx > 0

(28)

x

2

The terms on the right-hand side of eq 26 alternate in sign as the integrand changes sign at points where L′′(x) crosses the Gibbs energy curve G(x). The first and the third integrals on the right-hand side of eq 26 are analogous to the single-phase example, case 1. In general, they are negative but vanish when u1 tends to V1 and u2 tends to V2. At these limits, L′′(x) are tangent lines, L′(x), to G(x). All tangents in the intervals [u1,V1] and [u2,V2] are acceptable solutions, but the solutions of interest for this example are those that lie in the intervals [0,t1] and [t2,1], where t1 and t2 are chosen such that 0 < u1 e t1 e V1 < u2 e t2 e V2 < 1 and the line L′(x) is the common tangent to G(x) through the two points [t1,G(t1)] and [t2,G(t2)] (see Figure 9). The second integral of the three is in general positive. As for case 2, this integral is maximal when G(x) - L′′(x) is maximal

where the integration is taken over a region of composition space for which the integrand is positive. The regions (two in this case) for which the integrand is negative vanish when II is maximal. Case 4. Function G(x) Has Three Local Minimizers and Two Local Maximizers in the Interval [0,1]. Case 4 is illustrated in Figures 10 and 11. These figures are typical curves for the Gibbs energy of three-phase, two-component systems when a third phase is either incipient (Figure 10) or stable (Figure 11) over part of the composition range at a given temperature and pressure. As for cases 2 and 3, the convex envelope of G(x) does not coincide everywhere with the function itself but lies along the common tangent(s), L′(x), over part of the composition range. Statement. If G(x) is a function with three local minimizers (min1, min2, min3) and two local maximizers (max1, max2) in


Figure 10. Typical Gibbs energy curve for a two-component, three-phase system, where one of the phases is incipient over part of the composition range: (+) regions with positive area; (-) regions with negative area.

the interval [0,1] and L′′(x) is a straight line within the epigraph of conv(G(x)) and crosses (or touches) G(x) at six points u1, V1, u2, V2, u3, and V3 such that 0 e u1 e V1 e u2 e V2 e u3 e V3 e 1, then the value of eq 7 is maximal when L′′(x) are common tangents, L′(x), to G(x) at two (Figure 10) or four (Figure 11) points within the interval [0,1]. Proof. A proof of the above follows exactly that for case 2. Consideration of three integrals analogous to eqs (23-25), but with appropriate limits of integration leads to the same conclusions as those found for cases 2 and 3, namely, that II3 has a maximal value of zero when L′′(x) ) conv(G(x)) and II2 is the maximal value of II1 for the same condition. From cases 1-3, it follows immediately that this result holds when the integration limits extend over the whole of the composition range [0,1]. As for cases 2 and 3, further insight into the nature of the objective function II can be found by splitting the objective function integral into component parts such that the integrand alternates in sign at points where L′′(x) crosses the Gibbs energy curve G(x) so that

∫uV (G(x) - L′′(x)) dx ) ∫uV (G(x) - L′′(x)) dx + ∫Vu (G(x) - L′′(x)) dx + ∫uV (G(x) - L′′(x)) dx + ∫Vu (G(x) - L′′(x)) dx + ∫uV (G(x) - L′′(x)) dx 3

1

1

1

2

2

2

1

3

2

3

(29)

3

The first, third, and fifth integrals on the right-hand side of eq 29 are analogous to the single-phase example, case 1. In general, they are negative but vanish when the intervals between the limits of integration vanish. Under these conditions, L′′(x) are tangent lines, L′(x), to G(x). All tangents in the intervals [u1,V1], [u2,V2], and [u3,V3] are acceptable solutions, but like case 3, the solutions of interest for this example are those for which there exist common tangents to G(x) that lie in an interval between [u1,V3]. Three cases are of interest. Case 4a: Convex Envelope of G(x) Has Only One Common Tangent. This case is illustrated in Figure 10. The second and fourth integrals on the right-hand side of eq 29 are in general positive. By analogy with case 3, these integrals are maximal when G(x) - L′′(x) is maximal and the intervals [V1,u2] and [V2,u3] are maximal. However, for this example, the interval [u2,V2] vanishes at some point above the curve conv(G(x)). Below this point, the two separate positive terms in eq 29 coalesce and this example becomes identical to case 3 (vide supra).

Figure 11. Typical Gibbs energy curve for a two-component, three-phase system, where all the phases are stable over part of the composition range: (+) regions with positive area; (-) regions with negative area.

Case 4b: The convex envelope of G(x) has two common tangents. This case is illustrated in Figure 11. In this example, the interval [u2,V2] vanishes at some point below the curve conv(G(x)). The two separate positive terms in eq 29 retain their separate identity, and the arguments deployed for case 3 apply to each of the two terms separately. These lead to the same conclusion as before, namely, that the objective function, eq 7, is maximal when L′′(x) is equated to conv(G(x)), but in this case, the convex envelope consists in part of two common tangents L′(x) to the Gibbs surface, Figure 11. Equation 12 becomes

∫01(G(x) - conv(G(x))) dx ) ∫t t (G(x) t conv(G(x))) dx + ∫t (G(x) - conv(G(x))) dx > 0

IImax )

2

1

4

(30)

3

where t1, t2, t3, and t4 are points where the two common tangents touch the Gibbs surface. As for case 3, the limits of integration are not identical in eq 30 since the nonzero terms of IImax lie outside epi(G(x)). Generalization to hypersurfaces follows along the same lines adopted for case 3 and leads immediately to eq 28. While the integration extends in principle over the whole of the composition space, the nonzero terms are confined to regions outside epi(G(x)) and above the common tangent hyperplanes (two in this case). Case 4c: Convex Envelope of G(x) Has One Tangent Line That Touches G(x) at Three Points. This example is a limiting case of 4b where t2 tends to t3. The two separate common tangent lines become one that touches the Gibbs energy curve at three points t1, t2-3, and t4. Physically, this corresponds with the Gibbs energy curve for a two-component, three-phase system. Case 5: Function G(x) Has m Local Minimizers and m - 1 Local Maximizers in the Interval [0,1]. Statement. If G(x) is a function with m local minimizers (min1, min2, ...., minm) and m - 1 local maximizers (max1, max2, ...., maxm-1) in the interval [0,1] and L′′(x) is a straight line within the epigraph of conv(G(x)) and crosses (or touches) G(x) at 2m points u1, V1, u2, V2, u3, V3, ...., um, Vm, such that 0 e u1 e V1 e u2 e V2 e .... e um e Vm e 1, then the value of the objective function II (eq 7) within the interval [0,1] is maximal when L′′(x) comprises all possible common tangents, L′(x), to G(x) and the common tangents form part of the convex envelope of G(x).

1618


Proof. The proof proceeds along lines identical to that for case 3. Consideration of three integrals analogous to eqs 2325, but with appropriate limits of integration (eqs 31-33), leads to the same conclusions as those found for case 3, namely, that the maximal value of II3 is zero when L′′(x) ) conv(G(x)) and II2 is the maximal value of II1 for the same condition. From cases 1-3, it follows immediately that this result holds when the integration limits extend over the whole of the composition domain [0,1]

II1 )

∫u

Vm

(G(x) - L′′(x)) dx

(31)

1

∫uV (G(x) - conv(G(x))) dx

(32)

II3 )

∫uV (convG(x) - L′′(x)) dx

(33)

1

m

1

As for case 3, further insight into the nature of the objective function II can be found by splitting the objective function integral into component parts such that the integrand alternates in sign at points where L′′(x) crosses the Gibbs energy curve G(x). The alternating positive and negative terms can be grouped so that finally m

II )

∫uV (G(x) - L′′(x)) dx ) ∑∫uV (G(x) - L′′(x)) dx + m

i

1

i)1 m-1

i

∫V ∑ i)1

ui+1

(G(x) - L′′(x)) dx (34)

i

All the terms in the first summation satisfy the condition L′′(x) ⊆ epi(G(x)). They are negative with maximal values of zero when G(x) ) L′′(x) ) conv(G(x)) where ui ) Vi. They vanish when II is maximal, and make no contribution to the maximal value of II, IImax. The terms in the second summation satisfy the condition L′′(x) is not contained in epi(G(x)). Their integrands are all positive and are maximal when L′′(x) corresponds with conv(G(x)) over their respective ranges of integration (See cases 3 and 4). Finally therefore, the maximal value of the objective function II, eq 7, is given by

IImax )

∫01(G(x) - conv(G(x))) dx ) m-1

∫t ∑ i)1

t2i

IImax )

∫∫∫(G(x) - conv(G(x))) dx ) x

∑∫∫∫(G(x) - conv(G(x))) dx

(36)

x

II2 )

m

composition vector over the whole of the composition space and the L′′(x) becomes a hyperplane that cuts the hypersurface G(x) at 2m points, where m is the number of local minimizers. The integrals of eqs 31-35 become multiple integrals with limits of integration corresponding with their respective composition domains. Finally, the maximal value of the objective function IImax becomes

(G(x) - conv(G(x))) dx > 0 (35)

2i-1

where the [ti,G(ti)] are points where the common tangents L′(x) touch the curve G(x). Not all the terms in eq 35 may be active for the condition L′′(x) ) conv(G(x)). Some may vanish within the epigraph of conv(G(x)) before this condition is satisfied (see case 4a). The number of active terms will be identical to the number of common tangent lines that lie along conv(G(x)) within the interval [0,1]. As for case 4c, where the lower and upper integration limits for adjacent terms of the summation in eq 35 are the same, two separate common tangent lines become one that touches the Gibbs energy curve at three points. Physically, this situation corresponds with the Gibbs energy curve for a stable three-phase system. A similar result can be derived for a four-phase system, but the Gibbs phase rule limits the number of points where a common tangent can touch the Gibbs energy curve in this way. These results can be generalized to multiphase, multicomponent systems using the arguments developed for cases 3 and 4. The composition scalar in eqs 31-35 is replaced by a

where the first integration is over the whole of the vector composition space and the summation applies to integration over regions for which G(x) - conv(G(x)) > 0. Examples 1-5 have demonstrated that the area enclosed between a curve G(x), (or more generally a hypersurface) and a line (or hyperplane) through the epigraph of conv(G(x)) is maximal when the line (or hyperplane) coincides with the convex envelope of G(x). Although this result has been derived for Gibbs energy surfaces of two- and multicomponent mixtures at constant temperature and pressure, its application is more general. As the foregoing results have been derived from considerations of surface topology only, without, for example, any application of thermodynamic constraints such as the equality of chemical potentials criterion, they are applicable to all curves and surfaces that behave in a similar way. From a thermodynamic viewpoint, these results should not only be applicable to the Gibbs energy expressed as a function of composition at constant temperature and pressure, but also to the internal energy expressed as a function of composition at constant entropy and volume, Helmholtz energy expressed as a function of composition at constant temperature and volume, and the enthalpy expressed as a function of composition at constant entropy and pressure. Application to these surfaces is discussed in some detail below. Comparisons with Earlier Formulations of the Maximum Area Criterion In their pioneering paper, Eubank et al.11 focused attention exclusively on the Gibbs energy as a function of composition at constant pressure and temperature. They proved that when the area as defined by eq 7 was maximal, the equality of chemical potentials criterion was satisfied. In the Appendix, it is shown that, in general, the area between a common tangent line and a curve is a maximal value of eq 7. The proof is similar to that of Eubank et al. but makes no reference to equality of chemical potentials or any other thermodynamic criteria. This provides further evidence that the maximum area criterion, eq 7, should be applicable to surfaces other than those for the Gibbs energy of two-component systems at constant temperature and pressure. This point is illustrated below by applying the maximum area criterion to pure fluid phase equilibria and multicomponent, two-phase equilibria phenomena for the USV, HSp, AVT, and GpT surfaces. Eubank et al.11 discussed the curve depicted in Figure 12 at some length. This example corresponds with case 4a discussed above. Here, there are three local maximizers for eq 7, all of which satisfy the equality of chemical potentials criterion, areas II1, II2, and II3, all of which are positive. However, among these, area II3, where the common tangent line forms part of the convex envelope of G(x), is the largest since II1 ⊆ II3 and II2 ⊆ II3. In this example therefore, search procedures designed to locate the “maximum of the maxima” (see also ref 13) will correctly identify the stable equilibrium phase configuration.


Figure 12. Three maximal areas for the objective function, II, for the case illustrated in Figure 10: (- - - -) II1 or II2; (hatched region) II3.

Figure 13. Three maximal areas II4, II5, and II6 for the objective function, II, for the case illustrated in Figure 11: (- - - -) II1, II2, or II3; (hatched region) II5 or II6. II4 is the sum of II1, II2, and II3, which is less than the sum of II5 and II6. See text.

However, this last conclusion must be modified when conv(G(x)) has more than one common tangent line. An example of such a situation is shown in Figure 13, which corresponds with case 4b discussed above. In this example, search procedures designed to locate the local maximizers of eq 7 will identify areas II4, II5, or II6, all of which are positive and all of which satisfy the equality of chemical potentials criterion. There is no a priori reason why II4 (II4 ) II1 - |II2| + II3) should not be greater than either II5 or II6 (see Figure 13), so application of the maximum of the maxima principle may identify area II4 as the stable equilibrium solution. However, the tangent line corresponding with II4 intersects the curve G(x) at two points and hence does not correspond with the stable equilibrium phase configuration. Therefore the selection of L′(x)I4 (the tangent line corresponding with the area II4) as the stable tangent line is incorrect, even though II4 is the maximum of the maxima, when compared with either II5 or II6, separately. From this result, it is clear that blind application of the maximum area criterion in terms of the maximum of the maxima might lead to wrong results. In this case, common tangents corresponding with both II5 and II6 represent the stable phase configuration since as Figure 13 shows, these tangents do not intersect the curve G(x). However since II4 lies within the epigraph of conv(G(x)), while II5 and II6 lie within the regions bounded by the convex envelope of G(x) and G(x) itself, II5 + II6 > II4, a result epitomised by eq

Figure 14. Typical Gibbs energy curve for a two-component system with three separate two-phase regions: (hatched regions) II1, II2, or II3. The sum of II1, II2, and II3 is greater than any other possibility, e.g., the area bounded by the (- - - -) line and the curve G(x).

30. Therefore in principle, the maximum area criterion can still be retained, provided that it is modified to include the sum of all positive areas that lie between conv(G(x)) and G(x) over the whole composition domain, not just the areas associated with isolated regions. This last conclusion can be applied to more complex situations such as that depicted in Figure 14, for which there are three separate tangent lines that do not cross the curve G(x) at any point. Since all other common tangents for which eq 7 is maximal fall within the epigraph of conv(G(x)), the required tangent lines are those associated with II1, II2, and II3, each of which represents a stable equilibrium two-phase configuration. The sum of these is the required maximal value of eq 7 rather than any of the three areas separately. If only one (or the sum of any pair) of the three were selected, solutions corresponding with common tangents falling within the epigraph of conv(G(x)) may be associated with areas greater than any of these separately. These considerations imply that the curve G(x) has to be considered in its entirety when the maximum area criterion is used in phase stability analysis, a result consistent with eq 35. Finally, the modification of the maximum area criterion leads to the conclusion that the required global maximum of the m-1 objective function, eq 7, is Σi)1 IIi where m is the number of local minimizers in the function G(x) for ∀ xi ∈ [0,1] and not necessarily any of the separate terms of the summation. Within the summation, the individual terms correspond with stable twophase (or more for terms where the upper and lower bounds of the summation coincide) systems, but their individual values may well be less than those for metastable possibilities. Application of the Maximum Area Criterion to the Internal Energy, Enthalpy, Helmholtz Energy, and Gibbs Energy Surfaces for Pure Fluids and Multicomponent Mixtures USV Surface for Pure Fluids. The USV surface is the fundamental thermodynamic surface from which all the other surfaces may be derived. For stability, the internal energy should be minimal for fixed values for the canonical variables entropy and volume. Since entropy and volume are extensive properties and in general differ in coexisting phases, the objective function integral has to be taken along a line l in the plane U ) 0, that connects the state points (Sa,Va), (Sb,Vb) (Figure 15).

1620


( ( ( ( ( (

) ) ) ) ) )

∂la1

la1

(U(la2) - U(lb1)) (la2 - lb1) ∂U(lb1) 2 2 ∂lb1

la2

(U(la2) - U(lb1)) (la2 - lb1) ∂U(la2) ) 2 2 ∂la2 lb

lb1

(U(lb2) - U(la2)) (lb2 - la2) ∂U(la2) 2 2 ∂la2

lb2

(U(lb2) - U(la2)) (lb2 - la2) ∂U(lb2) ) 2 2 ∂lb2 la

la2

)

la1

∂II2(lb1,la2)

)

la2

∂II2(lb1,la2) ∂la2

1

∂II3(la2,lb2)

Figure 15. Typical internal energy curve and the projection of the tangent line on the plane U ) 0, for a one-component, two-phase system.

∂la2

)

lb2

∂II3(la2,lb2)

To proceed, assume for simplicity that the USV surface is similar to case 3 in the text. It follows immediately from eq A-14, that the objective function II in this example is given by

∂lb2

)} )} )} )} )} )}

(U(lb1) - U(la1)) (lb1 - la1) ∂U(lb1) 2 2 ∂lb1

∂II1(la1,lb1)

∂lb1

( ( ( ( ( (

lb1

)

lb1

∂lb1

{ { { { { {

(U(lb1) - U(la1)) (lb1 - la1) ∂U(la1) 2 2 ∂la1

∂II1(la1,lb1)

2

(40)

2

II )

II2i-1 + II2 ∑ i)1

(37)

It is now convenient to express l in terms of the entropy, S, and the volume, V. The procedure is similar to that given by Eubank et al.11 From Figure 15, it follows that

where the II2i-1 and II2 are given by eq A-18:

II1(la1,lb1) )

∫l

a 1

U(l) dl -

[

lb1U(la1) - la1U(lb1) +

II2(lb1,la2)

)

(lb1 + la1) (U(lb1) - U(la1)) 2

∫l U(l) dl -

]

b 1

[

+ 2

(la2

lb1)

(U(la2) - U(lb1))

]

and

∫l l U(l) dl -

[

(lb2

+ 2

dl ) dV(1 + DVS2)1/2

d(II2i-1) + dII2 ) ∑ ∑ i)1 i)1 2

i)1

∂II

∂lbi

∂II

∂lai

dlbi

lai

+

dlai

la

replaces u and

lb

(43)

Substitution of eqs 41-43 into eq 40 gives finally the following:

(

+

∂II

dlb1

la2

)

∂II1(Va1,Sa1,Vb1,Sb1)

lbi

∂lb1

Sb - Sa Vb - V a

(U(lb2) - U(la2)) (38)

() ∑( ) () () 2

(42)

where, for a given line l, DVS is dS/dV, the constant slope of the line connecting the points (Sa,Va), (Sb,Vb) in the U ) 0 plane. Therefore,

DVS ) dS/dV )

]

la2)

where, with a slight change of notation, replaces V. The total derivative of II is given by:

dII )

and

b 2

a 2

lb2U(la2) - la2U(lb2) +

2

)]

Sb - Sa 2 1/2 ) Vb - V a (Vb - Va)[1 + DVS2]1/2 (41)

la2

la2U(lb1) - lb1U(la2) +

II3(la2,lb2) )

[ (

lb - la ) (Vb - Va) 1 +

lb1

+

∂II

∂la2

dla2

(39)

lb1

For stationary values of II, the six partial derivatives of eq 39 should vanish. From eq 38, these derivatives are given as follows:

(

∂Va1

{

) (1 + DVS2)1/2

Vb1

(

)}

(

)}

(U(Vb1,Sb1) - U(Va1,Sa1)) (Vb1 - Va1) ∂U(Va1,Sa1) 2 2 ∂Va1

)

∂II1(Va1,Sa1,Vb1,Sb1) ∂Vb1

{

Va1

(U(Vb1,Sb1)

Vb1

) (1 + DVS2)1/2 - U(Va1,Sa1)) (Vb1 - Va1) ∂U(Vb1,Sb1) 2 2 ∂Vb1

Va1

(

)

∂II2(Vb1,Sb1,Va2,Sa2) ∂Vb1

{

)

∂II2(Vb1,Sb1,Va2,Sa2) ∂Va2

{

( (

)

∂II3(Va2,Sa2,Vb2,Sb2) ∂Va2

{

2

U(Vb1,Sb1))

-

(Va2

2

∂U(Vb1,Sb1) ∂Vb1 Va2

(

)}

(

)}

- U(Vb1,Sb1)) (Va2 - Vb1) ∂U(Va2,Sa2) 2 2 ∂Va2 ) (1 + DVS2)1/2

Vb2

(U(Vb2,Sb2) - U(Va2,Sa2)) (Vb2 - Va2) ∂U(Va2,Sa2) 2 2 ∂Va2

)

∂II3(Va2,Sa2,Vb2,Sb2) ∂Vb2 Va2

{

(U(Vb2,Sb2)

2

Vb1

Vb2

) (1 + DVS )

2 1/2

U(Va2,Sa2))

-

(Vb2

2

(

Va2)

)}

∂U(Vb2,Sb2) ∂Vb2 Va2

(44)

If the objective function II is to be stationary, the six partial derivatives of eqs 44 must go to zero. As shown in the Appendix, these derivatives vanish when (Ua,Sa,Va) and (Ub,Sb,Vb) are the same state point, in which case the area of the segment goes to zero or when (Ua,Sa,Va) and (Ub,Sb,Vb) are connected by a common tangent line to the USV surface. Only the first of these conditions is consistent with necessary thermodynamic criteria for equilibrium since at the same state point pressure, temperature, and chemical potential are the same. If a common tangent line is to be a line connecting points in coexisting phases, then the surface temperatures, pressures, and chemical potentials must be the same at both points. The six equations comprising eq 44 consist of three equation pairs, with each pair a function of Va, Vb, Sa and Sb. If one of the equation pairs is set to zero, this pair will constitute two equations with four unknowns Va, Vb, Sa and Sb. If one of these unknowns is fixed, the Gibbs phase rule states that it should be possible to determine the remaining three. However, in general, it is not possible to determine these unambiguously from a pair of equations containing four unknowns. Therefore, a condition or constraint has to be imposed, to provide essentially a third equation for the three unknowns to be found, e.g., from an equation of state giving U in terms of S and V. It is natural that the constraint should conform to the necessary criteria for phase equilibrium. With the addition of such a constraint, it can be shown that setting eqs 44 to zero satisfies the two remaining necessary criteria. A straightforward way of constraining eq 44 is to express the derivatives ∂U/∂V along the line l in terms of temperature and pressure. In general

dU ) T dS - p dV ∂S ) T( ) - p (∂U ∂V ) ∂V l

l

)

∂U Sb - Sa )T b -p ∂V l V - Va l

)}

(

Vb1)

) (1 + DVS2)1/2

Vb1

(U(Va2,Sa2)

( ) (

) (1 + DVS2)1/2

Va2

(U(Va2,Sa2)

(


) TD (∂U ∂V ) l

VS

-p

(45)

where the subscript l denotes differentiation in the direction of the line l and DVS has been substituted for (∂S/∂V)l from eq 43. In eq 45, both DVS and (∂U/∂V)l are constant along a given tie line connecting points (Ua,Sa,Va) and (Ub,Sb,Vb). Therefore, if a pair of terms in eq 44 is constrained by equating either the temperature or pressure at both points, the other intensive variable is similarly constrained through eq 45. Of course, no constraints are necessary when (Ua,Sa,Va), (Ub,Sb,Vb) are the same state point, since all thermodynamic properties will be identical for this condition. It now only remains to substitute the right-hand side of eq 45 for the ∂U/∂V in eq 44. For the first pair of equations, direct substitution gives the following:

(

)


(

{

∂Va1

) (1 + D 2VS)1/2

Vb1

(

1

(

1

)}

(U(Vb1,Sb1) - U(Va1,Sa1)) (Vb1 - Va1) Ta1(Sb1 - Sa1) -pa1 + 2 2 (Vb - Va )

)


{

∂Vb1

(U(Vb1,Sb1)

1

) (1 + D 2VS)1/2

Va1

)}

Tb1(Sb1 - Sa1) - U(Va1,Sa1)) (Vb1 - Va1) -pb1 + 2 2 (Vb - Va ) 1

(46)

As discussed above, the equality pa1 ) pb1 follows automatically if eq 46 is constrained by equating Ta1 and Tb1. Rearranging eqs 46 and equating them to zero to find the maximal condition gives

0 ) (1 + DVS2)1/2

{

([U(Vb1,Sb1) + pb1 Vb1] - [U(Va1,Sa1) + pa1 Va1]) 2 b b T1 S1 - Ta1 Sa1 2

(

)}

0 ) (1 + DVS2)1/2

{

([U(Vb1,Sb1) + pb1 Vb1] - [U(Va1,Sa1) + pa1 Va1]) 2 Tb1 Sb1 - Ta1 Sa1 2

(

)}

(47)

or

Ga1 ) Gb1

(48)

from both equations. For pure fluids, eq 48 is simply the equality of chemical potentials criterion for phase equilibrium. Therefore, the stationary values of eqs 44 satisfy the necessary criteria for phase equilibrium for the pure fluid USV surface, if the stationary values are constrained by equating the temperatures (or pressures) in both phases. It can be shown that the application of the equality of chemical potentials constraint, a priori, is also

1622


sufficient for providing satisfactory stationary values of the objective function, II. In conclusion therefore, stationary values of the objective function II for the pure fluid USV surface satisfy the necessary phase equilibrium criteria, if the stationary values are constrained by equating one of the three intensive variables, temperature, pressure, or chemical potential in the phases at equilibrium. The connection between the maximum area criterion and minimal internal energy follows in part from consideration of the second derivatives of the II2i-1 and II2. Comparisons with eqs A-19 and A-20 in the Appendix lead immediately to the following:

( ( ( ( ( (

∂

) ) ) ) ) )

2

II1(la1,lb1) 2

∂la1

lb1

∂2II1(la1,lb1) 2 ∂lb1

1

) -

la2

∂2II2(lb1,la2) 2

∂la2

2

∂la2

∂2II3(la2,lb2) ∂lb2

) -

lb1

∂2II3(la2,lb2)

2

( ( ( ( ( (

2

∂

)} )} )} )} )} )}

2

2

∂la1

lb1

la2

(la2 - lb1) ∂2U(la2) 2 2 ∂la2

lb1

(lb2 - la2) ∂2U(la2) ) 2 2 ∂la2 lb 2

(lb2 - la2) ∂2U(lb2) ) 2 2 ∂lb2 la 2

(

2

∂Va1

)

-(1 + DVS2)1/2

(

∂

2

)

II1(Va1,Sa1,Vb1,Sb1) 2 ∂Vb1

-(1 + DVS2)1/2

(

)

∂2II2(Vb1,Sb1,Va2,Sa2) 2

∂Vb1

)}

{ (

)}

-(1 + DVS2)1/2

)

∂2II3(Va2,Sa2,Vb2,Sb2) ∂Vb2

2

(Vb2 - Va2) ∂2U(Va2,Sb2) 2 2 ∂Va2

{ (

( )

)}

(Vb2 - Va2) ∂2U(Vb2,Sb2) 2 2 ∂Vb2

dp ) (49)

la2

() ()

(51)

l

(∂V∂p) dV + (∂T∂p) dT T

V

then

)} Vb1

)}

(Vb1 - Va1) ∂2U(Vb1,Sb1) 2 2 ∂Vb1

(∂V∂p) ) (∂V∂p) + (∂T∂p) (∂V∂T) l

T

V

l

Substitution in eq 51 yields the following:

( ) ( )(

( ))

()

∂2U ∂T ∂p ∂p ) DVS 2 ∂V ∂T V ∂V l ∂V l

dS )

(∂V∂S) dV + (∂T∂S) dT T

V

then

( ) ( ) ( )( ) ∂S ∂S ) ∂V l ∂V

+

T

cv ∂T ) DVS T ∂V l

Since

(∂T∂p) ) (∂V∂S)

Va1

V

T

direct substitution for (∂S/∂V)T gives

{

-(1 + DVS2)1/2

(

(50)

Va2

Since, in general

lb2

(Vb1 - Va1) ∂2U(Va1,Sa1) 2 2 ∂Va1

{ (

Vb2

)

Va2

Since, in general

{ (

Vb1

)

Vb2

∂ 2U ∂T ∂p ) DVS ∂V l ∂V ∂V2 l

)

Va2

)

2

∂Va2

)

Va1

{ (

(Va2 - Vb1) ∂2U(Va2,Sa2) 2 2 ∂Va2

For the objective function II to be maximal, the six second derivatives in eqs 50 must be negative. This implies that the six derivatives of the form (∂2U(V,S)/∂V 2) should be positive. The necessary relationships can be found from differentiation of eq 45 in the direction of l so that

)

Vb1

(

-(1 + DVS2)1/2

The same result may be obtained by differentiation of eqs 44. Substitution of V for l from eq 41 into eq 49 gives the following:

∂2II1(Va1,Sa1,Vb1,Sb1)

)

Vb1

∂2II3(Va2,Sa2,Vb2,Sb2)

la1

(la2 - lb1) ∂2U(lb1) 2 2 ∂lb1

)

2

∂Va2

-(1 + DVS2)1/2

(

U(la1)

(lb1 - la1) ∂2U(lb1) ) 2 2 ∂lb1 la

∂2II2(lb1,la2) 2 ∂lb1

{ { { { { {

) -

(lb1 - la1)

(

∂2II2(Vb1,Sb1,Va2,Sa2)

)}

(Va2 - Vb1) ∂2U(Vb1,Sb1) 2 2 ∂Vb1

Va2

(∂V∂T) ) cT (D - (∂T∂p) ) l

V

VS

and further substitution into eq 52 yields

V

T

(52)

( ) (

( ))

() ( ) (

∂2U ∂p 2 ∂p T ) DVS , 2 ∂T V ∂V T ∂V l cV T ∂ 2U ∂p 2 1 ) DVS + (53) 2 ∂T V VκT ∂V l cV

( ))

( )

where κT is the isothermal compressibility. For thermodynamically stable systems, the terms on the right-hand side of eq 53 are positive. It follows therefore that (∂2U/∂V2)l is always positive for stable thermodynamic systems and all the second derivatives on the left-hand side of eqs 49 and 50, i.e., of the objective function II, are negative. Therefore for the USV surface, the area enclosed between the tie line connecting points (Ua,Sa,Va), (Ub,Sb,Vb), and the internal energy must be maximal for stability in both the single phase and for coexisting phases. For coexisting phases, the maximal condition must be subject to an equality of intensive variables constraint, where the intensive variable is one of temperature, pressure, and chemical potential. It is possible to express eq 42 as

dl ) dS(1 + [dV/dS]2)1/2

(54)

where it is implicit that the dl are to be expressed in terms of dS rather than dV. The derivations follow through in exactly the same way and lead to identical conclusions, viz., that thermodynamic stability is associated with maximal values of the objective function II for the USV surface. Like the derivations in terms of volume, the maximal values are subject to an equality of intensive variables constraint. HSp and AVT Surfaces for Pure Fluids. In the case of pure fluids, the equilibrium constraints of equality of temperatures and equality of pressures in coexisting phases can be treated more naturally with the HSp and AVT surfaces, respectively. For stability, the enthalpy should be minimal for a given entropy and pressure and the Helmholtz energy should be minimal for a given volume and temperature. It can be shown that the maximum area criterion applies along isobars of enthalpy, H, plotted against entropy, S, and along isotherms of Helmholtz energy, A, plotted against volume, V. For case 3, equations analogous to eqs 38, 44, and 49 can be derived from the procedure outlined in the Appendix by replacing F(x) with H(S) for the HSp surface and F(x) with A(V) for the AVT surface. For these surfaces, DSp and DVT vanish as the path of integration coincides with the free independent variable, entropy, or volume. Further insight into the nature of the temperature and pressure constraints and their relation to the HSp and AVT surfaces may be obtained by applying the constraints directly to eq 46, which can be rearranged to give the following:

(

)

∂II1(Va1,Sa1,Vb1,Sb1) ∂la1

{

)

Vb1

([U(Vb1,Sb1)

+

pb1

Vb1]

2

[U(Va1,Sa1)

+

pa1

(

Va1])

-

)}

Tb1 Sb1 - Ta1 Sa1 2


(

)

∂II1(Va1,Sa1,Vb1,Sb1) ∂lb1

{

)

Va1

([U(Vb1,Sb1) + pb1 Vb1] - [U(Va1,Sa1) + pa1 Va1]) 2

(

)}

Tb1 Sb1 - Ta1 Sa1 2

(55)

or

(

)

∂II1(Va1,Sa1,Vb1,Sb1) ∂la1

)

Vb1

{

([U(Vb1,Sb1) - Tb1 Sb1] - [U(Va1,Sa1) - Ta1 Sa1]) + 2

(

)

∂II1(Va1,Sa1,Vb1,Sb1) ∂lb1

)}

(

pb1 Vb1 - pa1 Va1 2

)

Va1

{

([U(Vb1,Sb1) - Tb1 Sb1] - [U(Va1,Sa1) - Ta1 Sa1]) + 2

)}

(

pb1 Vb1 - pa1 Va1 2

(56)

The terms enclosed in small parentheses on the right-hand side of eqs 55 and 56 are simply the enthalpy and Helmholtz energy differences, respectively. The path of integration, l, is along the entropy (HSp) and volume (AVT) axes, respectively. Along isobars, the pressure is constant and the temperature T is (∂H/∂S). Along isotherms, the temperature is constant and the pressure is -(∂A/∂V). Equations 55 and 56 become the following:

( (

) { ) {

∂II1(Sa1,Sb1) ∂Sa1

)

Sb1

∂II1(Sa1,Sb1) ∂Sb1

)

Sa1

( )} ( )}

(H(Sb1) - H(Sa1)) (Sb1 - Sa1) ∂H(Sa1) 2 2 ∂Sa1

Sb1

(H(Sb1) - H(Sa1)) (Sb1 - Sa1) ∂H(Sb1) 2 2 ∂Sb1

Sa1

(57) and

(

)

∂II1(Va1,Vb1)

(

∂Va1

)

{

Vb1

)

∂II1(Va1,Vb1) ∂Vb1

( )}

(A(Vb1) - A(Va1)) (Vb1 - Va1) ∂A(Va1) 2 2 ∂Va1

Vb1

)

{

Va1

( )}

(A(Vb1) - A(Va1)) (Vb1 - Va1) ∂A(Vb1) 2 2 ∂Vb1

Va1

(58)

1624


Equations 57 and 58 are simply first derivatives of the Area method objective function along isobars for the HSp and isotherms for the AVT surfaces. Equating these to zero to obtain the stationary conditions and rearranging gives the following:

( ) { ∂H(Sa1)

)

∂Sa1

} ( ) {

(H(Sb1) - H(Sa1))

,

(Sb1 - Sa1)

Sb1

∂H(Sb1) ∂Sb1

( ) { ∂A(Va1)

}

(A(Vb1) - A(Va1))

)

∂Va1

Vb1

(Sb1 - Sa1)

Sa1

and

(Vb1 - Va1)

(59)

}

(60)

,

( ) { ∂A(Vb1) ∂Vb1

}

(H(Sb1) - H(Sa1))

)

)

(A(Vb1) - A(Va1)) (Vb1 - Va1)

Va1

Equations 59 and 60 are equations for the slopes of common tangents to points a and b along isobars for the HSp surface and isotherms along the AVT surface. These slopes are the system temperature and the negative of the pressure, respectively. Therefore, application of the maximum area criterion along HSp isobars and AVT isotherms automatically constrains both temperature and pressure to be equal in phases a and b without the need for “external” constraints. Furthermore, if the common values of temperature are substituted into eq 57 and those for pressure are substituted into eq 58, equations of the form Ga ) Gb are recovered. For the HSp and AVT surfaces therefore, stationary values of the objective function II, along isobars and isotherms, respectively, satisfy the basic necessary thermodynamic criteria for phase equilibrium, viz, the three intensive variables, pressure, temperature, and chemical potential should be equal in coexisting phases. The connection between the maximum area criterion, minimal enthalpy, or Helmholtz energy and stationary values for the objective function II follows, in part, from consideration of the second derivatives of eqs 57 and 58. Along an isobar (∂2H/ ∂S2)p is T/cp, where cp is the isobaric heat capacity, and along an isotherm, (∂2A/∂V2)T is [-(∂p/∂V)T] or 1/VκT, where κT is the isothermal compressibility. Therefore, differentiation of eqs 57 and 58 gives finally the following:

(

) { (

∂2II1(Sa1,Sb1) 2 ∂Sa1

) -

Sb1

( )} ) {

(Sb1 - Sa1) Ta1 , 2 cap1

∂2II1(Sa1,Sb1)

and

(

) { (

∂2II1(Va1,Vb1) 2

∂Va1

2

∂Sb1

lb - la )

( )}

(Sb1 - Sa1) Tb1 ) 2 cap1 Sa 1

( )}

(Vb1 - Va1) 1 ) 2 VκT b

V1

) {

∂2II1(Va1,Vb1) 2 ∂Vb1

a

(nbj (61)

-

naj )

[( ) ( ) ( ) ] Sb - Sa

nbj - naj

2

+

Vb - V a

2

+

nbj - naj

∑i

nbi - nai

2 1/2

)

nbj - naj

(nbj - naj )[DnjS2 + DnjV2 +

∑i Dn n 2]1/2 j i

(63)

and

,

1

) -

Va1

be negative and the objective function II must be maximal for stable systems. Therefore, only those stationary values of the objective function that correspond with maximal II need be considered. For example, in case 3, II1 and II3 are e 0, and II2 is g 0; so, the maximal condition is achieved when II1 and II3 vanish and II2 is a maximal, positive value. Although the USV surface is the fundamental thermodynamic surface, it is conceptually easier to apply the maximum area criterion to the HSp and AVT surfaces. For the latter surfaces, an intensive or field variable is one of the canonical variables and this can be fixed a priori, when applying the criterion without the necessity for other constraints. The canonical variables for the USV surface are however extensive and are not, in general, equal for coexisting phases. Therefore, the criterion must be applied along a general path in the plane U ) 0, not in the direction of one of the independent variables. GpT Surface for Pure Fluids and Thermodynamic Surfaces for Mixtures. To conclude the discussion on pure fluids, some remarks concerning the GpT surface are appropriate. In general, it is not possible to apply the maximum area constraint to phase equilibria for pure fluids with the Gibbs surface, since all three variables of the surface are equal here. The area between a tangent line to the surface and the surface itself is therefore zero for the whole surface. At a phase boundary, the first derivatives of the Gibbs energy with respect to pressure and temperature, which yield the volume and entropy, are discontinuous, since these are different in coexisting phases. Therefore for the GpT surface, slopes of the tangent lines are discontinuous at a phase boundary. These considerations do not, in general, apply to mixtures since, with the exception of azeotropes, mixture compositions are different in coexisting phases. For fixed temperature and pressure, plots of the Gibbs energy against composition yield curves similar to those discussed under cases 2-4. However, where compositions are the same in both phases, it may be advantageous to use a different thermodynamic surface. The relevant maximum area criteria may be derived in exactly the same way as for pure fluids and is given below. For case 3 in the text, eqs 37-40 apply to the fundamental USVn surface for multicomponent mixtures, where the path of integration l lies along the tangent hyperplane U ) 0 that connects the points (Sa,Va, na), (Sb, Vb, nb). Following Eubank et al.,11 eqs 41 and 42 become for a general mixture component j the following:

( )}

(Vb1 - Va1) 1 2 VκT

b

dl ) dnj[DnjS2 + DnjV2 + (62)

1

In eqs 61 and 62, the entropy and volume differences are positive by definition. Since temperature, isobaric heat capacity, and isothermal compressibility must be positive for thermodynamic stability, the second derivatives of eqs 61 and 62 must

∑i Dn n 2]1/2 j i

Since for a multicomponent mixture

dU ) T dS - p dV + then

∑i µi dni

(64)


( ) ( ) ( ) ∑( ) ∂U ∂nj

()

∂S

)T

-p

∂nj

l

∂V

+

∂nj

l

µi

i

l

S b - Sa

Vb - Va )T -p + ∂nj l nbj - naj nbj - naj

∂U

() ∂U ∂nj

) TDnjS - pDnjV +

l

∑i

components, then the equivalent of Nc - 2 chemical potentials must be specified. With fixed values for Nc - 2 chemical potentials, eqs 68 reduce to

∂ni

∂nj nbi

l

- nai

0 ) {naj (µbj - µaj ) + nak(µbk - µak)}

µi nbj - naj

0 ) {nbj (µbj - µaj ) + nbk (µbk - µak)}

∑i µiDn n

(65)

j i

and by analogy with eq 47, a typical equation pair for the maximal condition of the Area method objective function, II, becomes the following:

∑i Dn n ]

0 ) [DnjS + DnjV + 2

{ {

pb1

Vb1]

-

+

[U(Va1,Sa1)

2 Tb1 Sb1 - Ta1 Sa1

(

)

2

pa1

-

Va1])

-

∑i µai (nbi - nai )

∑i Dn n 2]1/2

}

j i

([U(Vb1,Sb1) + pb1 Vb1] - [U(Va1,Sa1) + pa1 Va1]) 2 b b T1 S1 - Ta1 Sa1

(

)

2

-

-

∑i µbi (nbi - nai )

}

(66)

For the determination of the stable phase configuration of mixtures, the T, p flash problem, where the temperature and pressure have fixed values, is of considerable importance. For fixed temperature and pressure, eqs 66 can be written more simply as

0)

0)

{ {

µbj ) µaj and

µbk ) µak

j i

+

0 ) [DnjS2 + DnjV2 +

where the chemical potentials for components j and k are not specified a priori. Dividing throughout by the total number of moles and some rearranging yields the conditions

2 1/2

2

([U(Vb1,Sb1)

(69)

(Gb - Ga)

-

2

(Gb - Ga)

-

2

} }

∑i

dGa b (ni - nai ) a dni

∑i

dGb b (ni - nai ) b dni

(67)

∑i µi ni

eqs 67 reduce to the following:

∑i nai (µbi - µai )}

0){ 0){

∑i nbi (µbi - µai )}

for acceptable solutions of eqs 69 with finite values of the mole fractions xj, xk. Equations 70 are the outstanding necessary conditions for equilibrium between two phases. This result shows, that for multicomponent, two-phase systems, the Area method is entirely consistent with the Gibbs phase rule. Apart from temperature and pressure, no degrees of freedom need be specified for binary systems since eqs 68 collapse to eqs 69 in this instance. This implies that, for twocomponent systems, it is not strictly necessary to know the overall composition to find the compositions of the phases at equilibrium. If the mole fractions in the two phases are the same, eqs 69 are no longer independent. In this instance, the Gibbs energy surface is unsuitable. However, eqs 66 can be expressed either in terms of the more suitable HSpn or AVTn surfaces, since neither entropy nor volume are the same in coexisting phases. Finally, if the fundamental USVn surface is used, the Area method should serve equally well for mixtures and pure fluids, but like the case for pure fluids, the solutions have to be constrained so that temperatures in coexisting phases, pressures, and possibly some chemical potentials should be the same. Concluding Remarks and Future Prospects

which is simply the Area method maximal condition in terms of the Gibbs energy. Since at constant temperature and pressure the Gibbs energy is given by

G)

(70)

(68)

For a multicomponent, two-phase system, there are Nc degrees of freedom, where Nc is the number of components. Since temperature and pressure are assumed to be fixed, Nc - 2 intensive variables must be fixed for the system to be completely specified. As the summations are over the total number of

A major objective of this paper was to remove any remaining doubts concerning the validity of the Area method maximum area criterion for phase stability analysis, and, in the opinion of the authors, this objective has been achieved. Subject to the modifications discussed above, not only is the method sound for phase stability analysis using the Gibbs and Helmholtz energy surfaces, but it is also valid for other surfaces with similar topology. Detailed analysis of typical curves for the Gibbs energy of two-component systems has shown that the original concept of the Area method for phase stability analysis, namely, that the stable equilibrium phase configuration corresponds with the maximum value for the area from the set of all possible areas must be modified to include the sum of all maximal areas of regions between the convex envelope of the Gibbs energy curve and the curve itself. Moreover, the summation should extend in principle over the whole of the composition domain. While the number of these may, at first sight, appear to be infinite, in practice, only regions outside the epigraph of the Gibbs energy curve (or, more generally, the Gibbs energy hypersurface) need be considered, since those within make no contribution to the maximal value of the Area method objective function. While the convex envelope formalism has made it possible to unite the phase behavior of pure fluids, one-, two-, and multicomponent mixtures within a single framework embodied

1626


by the maximum area criterion as encapsulated in eq 26, the capabilities (and limitations) of the Area method criterion, or more simply the Area method, have yet to be fully explored. The new formulation in terms of the convex envelope and the imposition of constraints where necessary, for example, for the USV surfaces of pure fluids and binary mixtures, removes some deficiencies in the old. The determination of the convex envelope of a function over a multidimensional hyperspace is a difficult and challenging task. There exist a number of publications in the literature for the approximate determination of the convex envelope, and this is a field of continuing research (see, for example, refs 21 and 22). If some characteristics of the function are known, such as Lipschitz continuity, it seems to be possible to use either some piecewise successive techniques or numerical integration techniques with a very fine grid. However, it should be remembered that what is required are parameters (mole fractions, volumes, etc.) that maximize the objective function, II. It may not be necessary in all instances, therefore, to use numerical procedures to locate the convex envelope. In some instances, it may be simpler to solve constrained sets of equations obtained by setting the first derivatives of II to zero and to examine the solutions obtained for maximal values (vide supra). This is essentially the derivative scheme, first proposed by Eubank et al.11 and discussed further by Elhassan et al.12 Within the general framework of the convex hull, it has been shown that the area between common tangent lines and the Gibbs energy function itself is equivalent to the area between the convex envelope of the function G(x) and G(x). This area is always maximal, not just for the Gibbs energy but for other thermodynamic functions such as the internal energy, enthalpy, and Helmholtz energy expressed in terms of their canonical variables. However, it would be interesting to see if the scheme can be extended to other mixture surfaces where one of the canonical variables is the dependent variable and to reacting systems. Of particular interest would be the S(UV) surface since the entropy, S, should be maximal for given internal energy and volume. For this surface therefore, the maximum area criterion should become the minimum area criterion and similar conclusions should be evident for other surfaces where the entropy is expressed as a function of other thermodynamic quantities. To the authors’ knowledge there has been no implementation of the Area method procedure in terms of hyperplanes and hypervolumes, a potentially interesting field of research as theoretical connections between the two-dimensional area and n-dimensional hypervolumes may exist for multiphase systems. The extension of the Area method criterion to the various thermodynamic surfaces (eqs 37-70), that followed from the procedure given in the Appendix, has produced some noteworthy results. A significant advantage of the scheme is its conceptual simplicity. Following the procedures given in the Appendix, Area method objective functions, and their first and second derivatives can be written down almost by inspection for all four basic thermodynamic surfaces. Clear connections have been established between the objective functions, for the several surfaces and the necessary (but not sufficient) criteria for thermodynamic stability. In general, stationary values of II, obtained by setting the first derivatives to zero, sometimes with the imposition of a constraint, satisfy the equality of intensive (or field) variables criteria but, not surprisingly, it is not possible to distinguish between the stationary values with these insufficient criteria alone. If the set of stationary values is divided into three subsets containing maximal, minimal, and saddle point values, it has been shown conclusively that only the subset containing the maximal values corresponds with realistic values

for thermodynamic properties such as the isobaric heat capacity and isothermal compressibility. If the maximal values are split into groups containing values associated with incipient, metastable, and stable phase regions, topological arguments, based on the convex envelope formalism, have shown that the sum of the areas associated with the stable phases is greater than any other single possibility. Until now, a major question to be resolved in the application of the Area method criterion is the problem of thermodynamic constraints. When are they necessary and how should they be applied for correct results? As originally formulated by Eubank et al.,11 constraints were unnecessary, as a binary system at constant temperature and pressure has no degrees of freedom. However, when the procedure was extended to ternary systems,12 the extra degree of freedom necessitated a thermodynamic constraint for satisfactory solutions to be obtained. It seems clear that an intimate connection between the Gibbs phase rule and the number of constraints exists, but this requires further elaboration for the multicomponent, multiphase case. The Gibbs phase rule certainly places restrictions on the nature of thermodynamic surfaces, in that tangent lines can touch the Gibbs curve at a limited number of points only, but to the authors’ knowledge, these restrictions and the consequences that flow therefrom have not been investigated in detail. The derivations for the USVx surface for pure fluids and binary mixtures for a two-phase system and their extension to the corresponding HSpx, AVTx, and GpTx surfaces have, however, shed some light on this problem. For the USVx surface, the fundamental thermodynamic surface, the area between the surface and the tangent line, has to be maximized subject to equality of temperature or pressure (pure fluids) or temperature and pressure (binary mixtures) criteria. It is not strictly necessary to specify what the values of temperature and pressure should be, only that they should be the same for coexisting phases. If, however, temperature (pure fluids) or temperature and pressure (binary mixtures) are known a priori, it is conceptually simpler to work in terms of the AVT (pure fluids) or the GpTx surface, since the necessary coordinate transformations effectively reduce the number of dimensions from three (pure fluids) or four (binary mixtures) to two. With hindsight, it is clear that the GpTx surface for binary mixtures is in general constrained, but these constraints take the form of fixed pressure and temperature rather than the more general ones of equality for the fundamental surface. It may be that similar coordinate transformations, which may be advantageous in reducing the number of dimensions, may exist for multicomponent mixtures. Such transformations have been used by Ung and Doherty23 in another context to reduce a constrained problem to an unconstrained one. However, these are substantial considerations and should be the subject of future investigations. One major limitation of the Area method is that it can only be applied to surfaces, where at least two of the variables have different values in coexisting phases. For pure fluids on phase boundaries, for example, it is not possible to apply the maximum area constraint to the GpT surface. For similar reasons, the GpTx surface for binary mixtures cannot be used on systems which exhibit azeotropy since for these systems the compositions are the same in coexisting phases. For these systems, the Area method criterion applied to the AVTx surface, where both A and V are different for coexisting phases, should be more suitable. It is clear, therefore, that, these restrictions apart, the maximum area criterion is of much wider applicability than has been assumed hitherto. In particular, it should now prove possible to analyze the phase behavior of those mixtures, including azeotropes for which the Gibbs energy is unsuited,


with other thermodynamic surfaces and the modified maximum area criterion. The example of case 2 has shown that the Area method is not necessarily limited to continuous functions. In their paper, Eubank et al.11 applied the Area method to a surface for the Gibbs energy of mixing, which was not uniformly functions with continuous over the whole of the composition range. In this example, the discontinuities in the surface and the first derivatives of the surface lay inside the two-phase region and, so, the integral scheme was applied successfully. The Area method is not, therefore, limited only to continuous functions or to continuous first derivatives. Since the Maxwell equal area approach (MEAR), currently under development by IglesiasSilva et al.19, is essentially a derivative of the Area method, it is probable that MEAR may itself be applicable to thermodynamic surfaces other than those for the Gibbs and Helmholtz energy.16,17 Fresh research is required in the field to establish whether or not these expectations are justified. Acknowledgment One of us, R.J.B.C., would like to thank his former colleagues, in particular Dr. A.E. Elhassan and Mr. J.E. Kilner, in the IUPAC Thermodynamic Tables Project Centre, Imperial College, London for many helpful discussions on the above. J.B. would like to acknowledge the support of the Hungarian Reseach Fund (OTKA, projects T 048377 and T 046822). J.B. and R.P.S. would like to acknowledge the support of the HungarianBulgarian Collaboration project No. BGA-20. Appendix In the following, it will be shown that the Area method maximum area criterion as given by Eubank et al.11 can be derived without reference to any thermodynamic criteria. It is not therefore limited to Gibbs energy curves at constant temperature and pressure but has a large number of other applications. From a thermodynamic point of view, the criterion may be applied to other surfaces, and examples of its application to three of these for pure fluids are discussed in the text. Statement. If F(x) is a general function of x with m local minimizers (min1, min2, ...., minm) and m - 1 local maximizers (max1, max2, ..., maxm-1) in the interval [0,1] and L′′(x) is a straight line within the epigraph of conv(F(x)) that crosses (or touches) F(x) at 2m points u1, V1, u2, V2, u3, V3, ..., um, Vm, such that 0 e u1 e V1 e u2 e V2 e ... e um e Vm e 1, then the value of the objective function II (eq 7) within the interval [0,1] is stationary when L′′(x) is a tangent line L′(x) or a common tangent line L′(x) to F(x). Proof. From Figure 16, the area between the line L′′(x) and the curve F(x) is given by m

II ) m

∑ i)1

m-1

II2i-1 +

II2j ) ∫u ∑ j)1

Vm

(F(x) - L′′(x)) dx )

1

m-1

∫u (F(x) - L′′(x)) dx + ∑ ∫V ∑ i)1 i)1 Vi i

ui+1

Figure 16. Positive (+) and negative (-) regions between a nonconvex curve and lines that cut the curve at 2m points: (a) unbroken line; (b) broken line.

Case 1: L′′(x) Is Unbroken in the Interval [u1,Wm]. In this case, L′′(x) is uniquely defined by the points [u1,F(u1)] and [Vm,F(Vm)]. The function II depends therefore on u1 and Vm only and

dII )

( ) ∂II ∂u1

Vm

du1 +

( ) ∂II ∂Vm

dVm

(A-2)

)0

(A-3)

u1

For II to be stationary, dII ) 0, so that

( ) ∂II ∂u1

) 0, and

Vm

( ) ∂II ∂Vm

u1

(F(x) - L′′(x)) dx (A-1)

i

where II2i-1 and II2j are the areas enclosed between F(x) and L′′(x) in the intervals [ui,Vi] and [Vi,ui+1], respectively. As for case 5, the integral is grouped into two sets of terms. In the first set, F(x) e L′′(x) in the interval [ui,Vi] and in the second F(x) g L′′(x) in the interval [Vi,ui+1]. Lines L′′(x) that make II stationary are required. As the value of II depends on the location of the points ui and Vi, then II ) f(ui,Vi), in general. It is helpful to consider cases where the line L′′(x) is unbroken and broken separately.

The equation of the line L′′(x) is given by

L′′(x) - F(u1) F(Vm) - F(u1)

)

x - u1 Vm - u1

(A-4)

since L′′(x) touches F(x) at [u1,F(u1)] and [Vm,F(Vm)]. Therefore at the limits of integration of eq A-1, L′′(u1) ) F(u1) and L′′(Vm) ) F(Vm). Substitution for L′′(x) in terms of eq A-4 into eq A-1 gives the following:

1628

II )


∫uV 1

)

∫uV 1

)

{

m

{

m

[

F(x) -

]}

dx

1 [(V F(u ) - u1F(Vm)) + x(F(Vm) Vm - u1 m 1

}

F(u1))] dx

[

∫uV F(x) dx - Vm -1 u1 (Vm - u1)(VmF(u1) - u1F(Vm)) + m

1

)

F(Vm) - F(u1) (x - u1) Vm - u1

F(x) - F(u1) +

∫uV F(x) dx m

1

(V2m - u21) (F(Vm) - F(u1)) 2

[

VmF(u1) - u1F(Vm) +

]

(Vm + u1) (F(Vm) 2

]

F(u1)) (A-5) Now, consider the partial derivatives of eq A-5 with respect to u1 and Vm. Differentiation with respect to u1 gives the following:

( ) { ∂II ∂u1

[( ) ( )

) -F(u1) - Vm

Vm

(Vm + u1) ∂F(u1) 2 ∂u1

{

∂F(u1) ∂u1 +

Vm

Vm

- F(Vm) -

]}

(F(Vm) - F(u1)) 2

)

( )}

(F(Vm) - F(u1)) (Vm - u1) ∂F(u1) 2 2 ∂u1

(A-6)

Vm

In a similar way, differentiation with respect to Vm gives the following:

( ) { ∂II ∂Vm

u1

( )}

(F(Vm) - F(u1)) (Vm - u1) ∂F(Vm) ) 2 2 ∂Vm

{

( )}

(A-8)

( )}

(A-9)

Vm

F(Vm) - F(u1) ) (Vm - u1)

∂F(Vm) ∂Vm

u1

Equations A-8 and A-9 are both satisfied when u1 ) Vm. When this condition holds, the points [u1,F(u1)] and [Vm,F(Vm)] coincide and the line L′′(x) becomes a tangent line L′(x) to the curve F(x). This situation is discussed under case 1 in the text for maximal values of the objective function II. Rearranging eqs A-8 and A-9 gives the following:

( ) { ( ) { ∂F(u1) ∂u1

∂F(Vm) ∂Vm

(F(Vm) - F(u1))

)

(F(Vm) - F(u1))

Vm

u1

} }

)

(Vm - u1) (Vm - u1)

(A-10) (A-11)

A second possibility for the objective function II to be stationary is apparent from eqs A-10 and A-11. A stationary value is obtained when the straight line L′′(x) is a common tangent line L′(x) to the curve F(x) at the points [u1,F(u1)] and

( ) { ( ) { ∂2II ∂u12

u1

and

{

[Vm,F(Vm)]. The consequences that follow from these two results for maximal values of II are discussed under cases 2-5 in the text. Unlike the arguments based on topology, those given here have only produced conditions for which the objective function II is stationary. The stationary points may be maximal, minimal, or saddle point values. These may be distinguished by considering second derivativessgiven belowsof the objective function, II. However, it is not possible to discriminate further between, for example, various maximal values in this way. The required second derivatives may be found by differentiating eq A-6 with respect to u1 and eq A-7 with respect to Vm. The differentiations yield the following:

(A-7)

Therefore, from eqs A-3, A-6, and A-7, it follows that for stationary conditions

∂F(u1) F(Vm) - F(u1) ) (Vm - u1) ∂u1

Figure 17. Stable, unstable, incipient, and metastable tangent lines for a Gibbs energy curve: IIUN, unstable region (unshaded); IIST, stable region (hatched); IIME, metastable region (vertical - - - -); IIIN, incipient region (horizontal - - - -).

and

∂2II ∂Vm2

) -

Vm

) -

u1

( )}

(A-12)

( )}

(A-13)

(Vm - u1) ∂2F(u1) 2 ∂u12

(Vm - u1) ∂2F(Vm) 2 ∂Vm2

Vm

u1

If both expressions given by eqs A-12 and A-13 are less than zero, then the area between F(x) and the tangent line L′(x) is maximal, and if they are greater than zero, it is minimal. If they alternate in sign, the area is a saddle point value. Since Vm > u1 by definition, the maximal condition will occur when both (∂2F(u1)/∂u12)Vm and (∂2F(Vm)/∂Vm2)u1 are greater than zero or in regions where F(x) is convex. Minimal values will occur when both (∂2F(u1)/∂u12)Vm and (∂2F(Vm)/∂Vm2)u1 are less than zero or in regions where the function F(x) is concave. Typical maximal and minimal regions are illustrated in Figure 17. A special case occurs when the points [u1,F(u1)] and [Vm,F(Vm)] coincide. Here, the second derivatives of the objective function are zero. Therefore at these points, the objective function may be maximal, minimal, or a point of inflection. However, if the surface is convex in the vicinity of these points (i.e., the second derivative of F(x) is positive), then eqs A-12 and A-13 approach zero from the negative side and the objective function II is maximal. This situation corresponds with case 1 in the text. Case 2: L′′(x) Is Broken in the Interval [u1,Wm]. Figure 11, case 4b, is an example where the convex envelope of the function G(x) contains two common tangent lines. This possibility can be treated in part with the derivative scheme by


allowing the line L′′(x) to be broken into segments in the interval [u1,Vm] with the breaks at points V1, u2, V2, u3, V3, ..., um, where the line L′′(x) intersects the function F(x) (See Figure 16b). In this case, L′′(x) is defined by the points [ui,F(ui)] and [Vi,F(Vi)]. The value of the objective function II depends, therefore, on all the ui and Vi so that the total derivative of II becomes m

dII )

()

m

m

∂II

∑ i)1 ∂u

()

∑ i)1 ∂V

dui +

i Vi

dVi +

i ui

∑ j)1

j)1

∂II

∂Vj

( ) ∂II ∂Vi

) 0,

Vi

) 0,

ui

( ) ∂II ∂Vj

) 0,

uj+1

duj+1 (A-14)

( ) ∂II ∂uj+1

∫u

Vi i

{

[

{

j+1

j

[

]}

F(Vi) - F(ui) (x - ui) Vi - ui

dx

F(uj+1) - F(Vj) (x - Vj) uj+1 - Vj

]}

dx (A-17)

Finally, the objective function integrals become, for the II2i-1 terms,

II2i-1(ui,Vi) )

∫uV F(x) dx i

i

[

ViF(ui) - uiF(Vi) +

]

(Vi + ui) (F(Vi) - F(ui)) 2

and, for the II2j terms,

II2j (Vj,uj+1) )

[

∫Vu

j+1

F(x) dx -

j

uj+1F(Vj) - VjF(uj+1) +

ui

]

(uj+1 + Vj) (F(uj+1) - F(Vj)) 2 (A-18)

( )} ( )}

(Vi - ui) ∂2F(ui) 2 ∂ui2

Vi

(Vi - ui) ∂2F(Vi) ) 2 ∂Vi2

ui

) -

)

(A-19)

(

)

∂II2j(Vj,uj+1) ∂uj+1

(

)

{

( )}

uj+1

vj

(F(uj+1) - F(Vj)) (uj+1 - Vj) ∂F(Vj) 2 2 ∂Vj

(

∂Vj2

{

(

)}

(F(uj+1) - F(Vj)) (uj+1 - Vj) ∂F(uj+1) 2 2 ∂uj+1

) { ( )} ) { ( )} ) -

uj+1

∂2II2j(Vj,uj+1) ∂uj+12

uj+1

)

∂2II2j(Vj,uj+1)

F(x) - F(ui) +

F(x) - F(Vj) +

(

∂II2j(Vj,uj+1) ∂Vj

(A-16)

and for the II2j terms

II2j(Vj,uj+1) )

ui

(ii) For the area segment II2j

)0

since the L′′(x,ui,Vi) touch or cross F(x) at [ui,F(ui)] and [Vi,F(Vi)]. Substitution for L′′(x,ui,Vi) in terms of eq A-16 into eq A-1 gives

II2i-1(ui,Vi) )

(F(Vi) - F(ui)) (Vi - ui) ∂F(Vi) 2 2 ∂Vi

Vi

∂Vi2

Vj

For the area segment II2i-1(ui,Vi), the equation of the line L′′(x,ui,Vi) is given by

F(Vi) - F(ui)

ui

( )} ( )} Vi

∂2II2i-1

Vj

x - ui ) Vi - ui

)

∂ui2

(A-15)

L′′(x,ui,Vi) - F(ui)

Vi

(F(Vi) - F(ui)) (Vi - ui) ∂F(ui) 2 2 ∂ui

∂2II2i-1

where the first two summations on the right-hand side are for area segments II2i-1(ui,Vi), for which F(x) < L′′(x,ui,Vi), and the second two are for area segments II2j (Vj,uj+1), for which F(x) > L′′(x,Vj,uj+1). From eq A-1, the II2i-1(ui,Vi) are negative and the II2j(Vj,uj+1) are positive in sign. For II to be stationary, dII ) 0, so that

( )

)

( ) { ( ) {

dVj +

uj+1

∂II

∂uj+1

) { ) {

∂II2i-1(ui,Vi) ∂ui

∂II2i-1(ui,Vi) ∂Vi

() ∑( ) m-1

∂II

m-1

∫Vu

( (

m-1

dII2i-1 + ∑ dII2j ) ∑ i)1 j)1

∂II ∂ui

Now, consider the first and second partial derivatives with respect to ui, Vi, Vj, and uj+1. Comparisons with the derivatives for the unbroken line lead immediately to the following: (i) For the area segment II2i-1

vj

(uj+1 - Vj) ∂2F(Vj) 2 ∂Vj2

(uj+1 - Vj) ∂2F(uj+1) ) 2 ∂uj+12

vj

uj+1

(A-20)

Vj

For stationary values of the objective function II, the first derivatives in eqs A-19 and A-20 should be zero for all i and j. By analogy with eqs A-8 and A-9, stationary values may occur for the following conditions: (i) All the ui ) Vi and all the Vj ) uj+1. In this “trivial” case, all the segments of the broken line are tangents to the curve F(x). The total area vanishes as can be seen from eqs A-17 and A-18. (ii) All the Vj ) uj+1, some of the ui ) Vi, and some of the II2i-1 are such that the conditions (∂F(ui)/∂ui)Vi ) {(F(Vi) - F(ui))/ (Vi - ui)} and (∂F(Vi)/∂Vi)ui ) {(F(Vi) - F(ui))/(Vi - ui)} are satisfied. In this case, all the terms of the objective function II for which F(x) > L′′(x,Vj,uj+1) are zero. The two conditions are equations for the slopes of common tangents L′(x,ui,Vi) to the curve F(x) at the points [ui,F(ui)] and [Vi,F(Vi)] for the various ui,Vi. For these terms, F(x) < L′′(x,ui,Vi), and therefore, the objective function II is negative. (See also eq A-1.) A typical area for this condition is marked IIun in Figure 17. Since the second derivatives of F(x) are less than zero (the function is concave at these points), the second derivatives of the II2i-1 are positive (eq A-19) and the objective function II is minimal for these II2i-1. The number of minimal values will increase with an increasing number of common tangents to the curve F(x),

1630


and subject to the condition that 0 e u1 e V1 e u2 e V2 e ... e um e Vm e 1, the sum of these will be minimal, when all the II2i-1 are areas enclosed between the curve F(x) and the common tangents, L′(x,ui,Vi). (iii) All the ui ) Vi, some of the Vj ) uj+1, and some of the II2j are such that the conditions (∂F(Vj)/∂Vj)uj+1 ) {(F(uj+1) F(Vj))/(uj+1 - Vj)} and (∂F(uj+1)/∂uj+1)Vj ) {(F(uj+1) - F(Vj))/ (uj+1 - Vj)} are satisfied. In this case, all the terms of the objective function II for which F(x) < L′′(x,ui,Vi) and which, in general, make negative contributions to the objective function are zero. The last two equations are slopes of the common tangents L′(x,Vj,uj+1) to the curve F(x) at the points [Vj,F(Vj)] and [uj+1,F(uj+1)] for the various Vj, uj+1. For these terms, F(x) > L′′(x,Vj,uj+1) and, therefore, the objective function II is positive. (See also eq A-1.) A typical area for this condition is marked IIST in Figure 17. Since the second derivatives of F(x) are greater than zero (the function is convex at these points), the second derivatives of the II2j are positive and the objective function is maximal for these II2j. The number of maximal values will increase with an increasing number of common tangents to the curve F(x), and subject to the condition that 0 e u1 e V1 e u2 e V2 e ... e um e Vm e 1, the sum of these will be maximal, when all the II2j are areas enclosed between the curve F(x) and common tangents L′(x,Vj,uj+1). This is essentially the same result as that derived from topological considerations discussed earlier in the text. The foregoing results, like the topological arguments in the main text, are quite general and do not depend on any thermodynamic assumptions. From the viewpoint of thermodynamic stability, the stratagem of splitting the line L′′(x) into broken sections, coupled with the condition 0 e u1 e V1 e u2 e V2 e ... e um e Vm e 1, avoids the possibility of metastable solutions that may result from an unbroken line. A typical metastable region, designated IIME and its associated tangent line is shown in Figure 17. From a computational viewpoint, metastable possibilities can be avoided by maximizing each term of eq A-14 separately, while simultaneously checking that the condition 0 e u1 e V1 e u2 e V2 e ... e um e Vm e 1 holds. If at some point this condition becomes invalid, this implies that an incipient phase is present and the number of terms in eq A-14 should be reduced accordingly. An example of an incipient region, designated IIIN and the associated tangent line is shown in Figure 17. Nomenclature A ) Helmholtz energy C ) convex set conv ) convex envelope of a function cp ) isobaric heat capacity cV ) isochoric (constant volume) heat capacity Dxy ) slope of a line in the x,y plane epi ) epigraph of a function F ) function G ) Gibbs energy G′ ) a tangent line or tangent hyperplane to a Gibbs energy curve or hypersurface H ) enthalpy II ) Area method objective function (geometrically, the area between a curve or a hypersurface and an intersecting line or hyperplane) L′ ) a tangent line or tangent hyperplane to a curve or hypersurface L′′ ) a straight line or hyperplane that intersects a curve or a hypersurface l ) integration path for a line integral of a function

m ) number of local minimizers of a function max ) a local maximizer min ) a local minimizer Nc ) number of components in a system n ) number of moles in a system P ) number of equilibrium phases p ) pressure Pi ) element of a set S ) entropy S ) set T ) thermodynamic temperature t ) tangent point(s) U ) internal energy u ) points of intersection of a straight line with a curve V ) volume V ) points of intersection of a straight line with a curve x ) independent variable of a function x ) mole fraction composition vector of a system y ) independent variable of a function Greek Letters λ ) parameter for which 0 e λ e 1 κT ) isothermal compressibility µ ) chemical potential Superscripts a ) equilibrium phase “a” b ) equilibrium phase “b” Subscripts 1, 2, 3, etc. ) a specific element of a set of numbers or a specific property value i, j, k, m ) a general element of a set of numbers or a generalized property value i, j ) running indices in a summation IN ) incipient phase l ) path of integration lim ) limiting value max ) maximal value ME ) metastable phase stab ) global stability ST ) stable phase UN ) unstable phase Literature Cited (1) Baker, I. E.; Pierce, A. C.; Lucks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982, 22, 731. (2) Michelsen, M. L. The Isothermal Flash Problem. Part 1. Stability. Fluid Phase Equilib. 1982, 9, 1. (3) Stateva, R. P.; Tsvetkov, St. G. A New Method for Thermodynamic Stability Analysis in Multicomponent Systems. Hung. J. Ind. Chem. 1991, 19, 179. (4) Stateva, R. P.; Tsvetkov, St. G. A Rigorous Approach to Stability Analysis as a First Step when solving the Isothermal Multiphase Flash Problem. Technol. Today 1991, 4, 223. (5) McDonald, C. M.; Floudas, C. A. Decomposition Based Branch and Bound Global Optimisation Approaches for the Phase Equilibrium Problem. J. Global Optim. 1994, 5, 205. (6) McDonald, C. M.; Floudas, C. A. Global Optimisation and Analysis for the Gibbs Free Energy Function using UNIFAC, Wilson and ASOG Equations. Ind. Eng. Chem. Res. 1995, 34, 1674. (7) McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase and Chemical Equilibrium Problem: Application to the NRTL Equation. Comput. Chem. Eng. 1995, 19, 1111. (8) McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase Stability Problem. AIChE J. 1995, 41, 1798. (9) McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1997, 21, 1.

Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1631 (10) Harding, S. T.; Floudas, C. A. Phase Stability with CEOS - Global Optimisation Approach. AIChE J. 2000, 46, 1422. (11) 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. (12) 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. (13) Elhassan, A. E.; Tsvetkov, St. G.; Craven, R. J. B.; Stateva, R. P.; Wakeham, W. A. A Rigorous Mathematical Proof of the Area Method for Phase Stability. Ind. Eng. Chem. Res. 1998, 37, 1483. (14) Elhassan, A. E.; Craven, R. J. B.; Stateva, R. P.; Wakeham, W. A. Rebuttal to the Comments by Paul I. Barton and Chyi Hwang in “No Connection between the AREA Criterion and Phase Stability has been Established” on the paper “A Rigorous Mathematical Proof of the Area Method for Phase Stability”. Ind. Eng. Chem. Res. 2000, 39, 3399. (15) Barton, P. I.; Hwang, C. No Connection between the AREA Criterion and Phase Stability has been Established. Ind. Eng. Chem. Res. 2000, 39, 3398. (16) Hanif, N. M.; Shyu, G. S.; Hall, K. R.; Eubank, P. T. Application of the Equal Area Rule to the Calculation of Multiphase Equilibria of Hydrocarbon-Water Mixtures with an EOS Model. Fluid Phase Equilib. 1996, 126, 53.

(17) Hanif, N. M.; Shyu, G. S.; Hall, K. R.; Eubank, P. T. The Area of J. C. Maxell for Pure Component Fluid Phase Equilibria from Equations of State. Ind. Eng. Chem. Res. 1996, 35, 2431. (18) Iglesias-Silva, G. A.; Bonilla-Petriciolet, A.; Eubank, P. T.; Holste, J. C.; Hall, K. R. An Algebraic Method that includes Gibbs Minimization for Performing Phase Equilibrium Calculations for any Number of Components or Phases. Fluid Phase Equilib. 2003, 210, 229. (19) Iglesias-Silva, G. A.; Bonilla-Petriciolet, A.; Hall, K. R. An algebraic Formulation for an Equal Area Rule to determine Phase Compositions in Simple Reactive Systems. Fluid Phase Equilib. 2006, 241, 25. (20) Heidemann, R. A. Non-uniqueness in Phase and Reaction Equilibrium Computations. Chem. Eng. Sci. 1978, 33, 1517. (21) Kadhi, F.; Trad, A. Characterization and Approximation of the Convex Envelope of a Function. JOTA (J. Optim. Theory Appl.) 2001, 110, 457. (22) Brighi, B.; Chipot, M. Approximated Convex Envelope of a Function. SIAM J. Numer. Anal. 1994, 31, 128. (23) Ung, S.; Doherty, M. F. Theory of Phase Equilibria in Multireaction Systems. Chem. Eng. Sci. 1995, 50, 3201.

ReceiVed for reView July 21, 2006 ReVised manuscript receiVed November 15, 2006 Accepted November 20, 2006 IE060949F

The Area Method for Phase Stability Analysis Revisited: Further

Recommend Documents