Global Stability Analysis and Calculation of Liquid−Liquid Equilibrium

Ind. Eng. Chem. Res. 1996, 35, 1395-1408

1395

GENERAL RESEARCH Global Stability Analysis and Calculation of Liquid-Liquid Equilibrium in Multicomponent Mixtures† Stanislaw K. Wasylkiewicz,‡ Lakshmi N. Sridhar, Michael F. Doherty,* and Michael F. Malone Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003

A global algorithm for implementing the Gibbs tangent plane stability test in multicomponent, multiphase liquid mixtures is described. The algorithm is self-starting and significantly improves the reliability and robustness of multiphase equilibrium calculations. The main improvement is a result of a new approach for locating the stationary points of the tangent plane distance function. This relies on the geometric properties of the tangent plane distance function but is independent of the specific fluid model. For ternary mixtures, it is shown that stationary points in the tangent plane distance function satisfy a global topological constraint Nmax + Nmin - Nsad ) 1 where Nmax is the number of maxima, Nmin is the number of minima, and Nsad is the number of saddles. This provides a global consistency check on the numerical results that is independent of the specific model or parameters used to describe the mixture. The method has been implemented and tested on a variety of problems with three and four components and up to three liquid phases. Solutions can be found with sufficient speed for interactive design and modeling. In all of the examples, solutions for the entire phase diagram were found without fail, and certain problems that are difficult or impossible to solve using other methods were modeled successfully. Introduction A common phase equilibrium problem in process modeling is to calculate the boiling temperature and the vapor composition in equilibrium with a liquid phase of specified overall composition, z, at a given pressure, P. If the liquid phase-separates, then the number and compositions of all liquid phases plus the vapor phase at the specified overall liquid composition must be found. Given P and z, the first task is to decide whether the liquid is homogeneous or heterogeneous. The typical method is to estimate the boiling temperature, to perform a liquid-phase stability test, and, based on the outcome, to perform either a vapor-liquid equilibrium (VLE) or a vapor-liquid-liquid equilibrium (VLLE) calculation. The result provides a new estimate of the boiling temperature, and the process is repeated until the equilibrium equations are satisfied. During the course of these iterations, liquid phases may come and go, but the correct number of equilibrium phases predicted by the model will typically be found if a reliable stability check is used. If a poor stability check or no stability check is used, it is possible to reach the incorrect conclusion that the liquid is homogeneous when the model actually predicts that there are multiple phases. In this case, multiple spurious solutions to the phase equilibrium equations often appear (Baker et al., 1982; Van Dongen et al., 1983; Lucia and Taylor, 1992). This is a classic problem * To whom correspondence should be addressed. † Parts of this work were presented in paper 151a, AIChE Annual Meeting, St. Louis, MO, Nov 1993. ‡ On leave from Institute of Chemical Engineering and Heating Equipment, Technical University, Wrocław, Poland.

0888-5885/96/2635-1395$12.00/0

in thermodynamicsshow to distinguish between unstable, metastable, and absolutely stable phases. This problem was solved theoretically, though not computationally, more than 100 years ago by Gibbs (1873, 1876), who formulated the Gibbs tangent plane test. Smith et al. (1993) and Jiang et al. (1995) present a more general criterion (but not its numerical implementation) for global multiphase, multireaction, chemical, and phase equilibrium, by applying the first-order Kuhn-Tucker necessary conditions to a new formulation of the general equilibrium problem. The tangent plane test for phase equilibrium is a special case of their criterion. It is only recently that a practical numerical implementation of the Gibbs tangent plane stability test has been developed (Baker et al., 1982; Michelsen, 1982a,b). These papers have revolutionized the way phase equilibrium calculations are performed and have produced a steady series of papers and algorithms during the last decade refining the basic ideas. For examples in azeotropic distillation see Cairns and Furzer (1990) or Widagdo et al. (1992). Some alternative stability tests require the evaluation of the sign definiteness of the Hessian matrix of the Gibbs free energy function with respect to the c - 1 independent mole fractions, e.g., Van Dongen et al. (1983). Stability to infinitesimal disturbances requires the Hessian matrix to be positive definite but does not distinguish between stable and metastable regions. Furthermore, this approach does not provide initial estimates of individual liquid-phase compositions for use in subsequent liquid-liquid equilibrium calculations. Michelsen’s approach, on the other hand, does distinguish absolutely stable from metastable regions and also © 1996 American Chemical Society

1396

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996

provides initial estimates for the solution of the liquidliquid problem. At a fixed temperature and pressure, the problem is to determine whether a liquid mixture at a given overall composition is globally stable with respect to liquidliquid phase splitting and to calculate the compositions of each individual liquid phase in the event of instability. A major difficulty in the solution of this problem rests in the existence of an undesirable trivial solution, where the composition of each liquid phase is equal to the overall liquid composition. This solution satisfies all of the governing equations irrespective of stability (Michelsen, 1993). Another difficulty with existing methods is the lack of certainty that all minima of the tangent plane distance function have been found. Some authors have tried to check “exhaustively” by starting minimization algorithms or root-finders from several points, e.g., all pure component compositions along with points generated from heuristics (Michelsen, 1982a). This approach can be time consuming and still does not guarantee that all minima of F(x) have been found, especially if these minima are very close in composition, e.g., in the vicinity of a critical point. Homotopycontinuation methods for finding the stationary points in the tangent plane distance function are more reliable. For example, the method described by Sun and Seider (1995) solves a much wider collection of problems. However, this method can still fail for certain cases (Sun and Seider (1995), pp 236-237). The need to develop a more reliable numerical implementation of the Gibbs tangent plane test is particularly pressing in connection with data regression, where parameters of thermodynamic models are adjusted to match experimental data. In fact, failure to converge at “difficult” points often causes some points to be ignored in a regression (Michelsen, 1993) even to the extent that the data for some mixtures cannot be treated at all by existing algorithms, e.g., the mixture water, nitromethane, and hexanol in the DECHEMA Chemistry Data Series (Sorensen and Arlt, (1980), Vol. V/2, p 69) for which only experimental data are presented and a three liquid-phase region has been experimentally detected. Using the UNIQUAC binary interaction parameters reported by Sorensen and Arlt (1980, Vol. V/3, pp 422 and 423) in our new algorithm, we predict the three-phase region for this mixture (see example 5). Chemical process synthesis, design, and simulation are other areas where a more reliable and reasonably fast stability test and LLE calculation method would be valuable. This paper describes an improved computational approach, including a more reliable stability test as well as a robust method to calculate the equilibrium compositions. The method was developed to imbed in our design, synthesis, and simulation software (Mayflower) for multicomponent nonideal mixtures (including azeotropic distillation and liquid extraction). The method has been tested extensively for over 2 years on a wide variety of mixtures. It succeeded at every point in all the ternary diagrams reported in this paper, as well as all the quaternary points, such as the one given in example 6. While no algorithm can be guaranteed to work on every problem without failure, we have not experienced any failures so far. The next section gives a brief description of the stability test, distinct from its implementation, which is discussed later in the paper. New results are presented and the solution methods are described for both the stability test and phase equilib-

rium. We present a summary of the algorithm and apply it to several model problems, some of which show difficult behavior. The final section gives our conclusions and remarks on the limitations of the approach. Tangent Plane Stability Test Gibbs (1876) showed that a necessary and sufficient condition for absolute stability of a mixture at a fixed temperature, pressure, and overall composition is that the Gibbs energy surface be at no point below the plane tangent to the surface at the given overall composition. A more modern explanation and derivation of this condition is given by Baker et al. (1982), Smith et al. (1993), and Jiang et al. (1995). Let the overall composition of a mixture be z at a temperature, T, and pressure, P, and let the chemical potential of the ith component of the mixture be µi(T,P,z) at these conditions. We will denote this as µi(z) for the purpose of this paper; other functions will be denoted similarly. All functions and partial derivatives are understood to be held at fixed T and P. The stability condition can be stated quantitatively as follows. First we define c

F(x) ≡ g(x) - L(x;z) )

xi[µi(x) - µi(z)] ∑ i)1

(1)

where g(x) is the Gibbs energy at composition x, and L(x;z) is the tangent plane to the Gibbs energy surface at composition z. F(x) is the tangent plane distance function, which represents the displacement from the tangent plane to the Gibbs energy surface at composition x. The necessary and sufficient condition for global stability is

F(x) g 0

(2)

for all x; cf. Theorems 3 and 4 of Baker et al. (1982). Figure 1 shows g(x) and F(x) for a mixture of heptane and methanol at P ) 1 atm, T ) 0 °C, and overall compositions of 17, 14, and 14.73 mol % methanol. The mixture is unstable (locally metastable) at an overall composition of 17 mol % methanol since the Gibbs free energy surface lies below the tangent plane for some compositions, where F(x) < 0 (Figure 1a). A mixture containing 14 mol % methanol is absolutely stable since F(x) g 0 everywhere (Figure 1b). The common tangent points correspond to equilibrium points, which are globally stable if F(x) g 0 for all x. This happens for a mixture containing 14.73 mol % methanol (Figure 1c). Implementing the Stability Test In an important development, Michelsen (1982a) noted that all of the minima in F(x) are located in the c interior of the permissible region (xi g 0, ∑i)1 xi ) 1). Consequently, F(x) will be nonnegative if it is nonnegative at all stationary points where its derivatives with respect to all of the independent variables (xi; i ) 1, 2, ..., c - 1) vanish. Differentiation of F(x) with respect to the c - 1 independent mole fractions, along with the use of the Gibbs-Duhem equations, yields the stationary conditions

µi(xˆ ) - µi(z) ) µc(xˆ ) - µc(z)

i, ..., c - 1

(3)

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1397

Figure 1. Gibbs free energy and tangent planes for a mixture of heptane and methanol. The NRTL model is used with interaction parameters taken from Sorensen and Arlt (1980), Vol. V/2, p 133. An overall composition of 17 mol % methanol is unstable (metastable) (a), while 14 mol % is stable (b). The common tangent points (c) at 14.73 and 91.48 mol % methanol correspond to the globally stable liquid-liquid equilibrium for all overall compositions between these two values.

This is a system of c - 1 equations in c - 1 unknown mole fractions xˆ , which is the liquid composition at a stationary point. Note that µc(xˆ ) - µc(z) on the righthand side of eq 3 depends on z and xˆ but is independent of the component index i. The corresponding stationary value of the tangent plane distance function is

activity coefficients. The stability test in eq 5 can be rewritten as c

F(xˆ )/RT )

xˆ i[ln xˆ i + ln γi(xˆ ) - ln zi - ln γi(z)] g 0 ∑ i)1 (6)

c

F(xˆ ) )

xˆ i[µc(xˆ ) - µc(z)] ) µc(xˆ ) - µc(z) ∑ i)1

(4)

There is always at least one solution of eq 3 at the trivial stationary point xˆ ) z. When F(x) has additional stationary points, the mixture may phase-separate. The mixture is stable if and only if

F(xˆ ) g 0

(5)

for all stationary points, i.e., solutions of eq 3. Consequently, the mixture is stable if and only if µc(xˆ ) - µc(z) is nonnegative at all stationary points. For liquid mixtures at low to moderate pressures it is common to write the chemical potentials in terms of

where γi(xˆ ) and γi(z) are the activity coefficients evaluated at compositions xˆ and z, respectively. Note that ln zi + ln γi(z) in eq 6 is a constant for a fixed overall composition. Equation 3 can be rewritten as

ln xˆ i + ln γi(xˆ ) - ln zi - ln γi(z) ) ln xˆ c + ln γc(xˆ ) i ) 1, ..., c - 1 (7) ln zc - ln γc(z) Thus, another statement of the stationary conditions (eq 3) is

xˆ iγi(xˆ ) xˆ cγc(xˆ )

ziγi(z) )

zcγc(z)

i ) 1, ..., c - 1

(8)

1398

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996

Equation 8 is similar to the liquid-liquid equilibrium equations except that the activities, xˆ iγi(xˆ ), at the stationary point are not equal but only proportional to the activities, ziγi(z), at the overall composition, by the factor xˆ cγc(xˆ )/zcγc(z). If this factor is unity, the activities of component c are identical at xˆ and z, the stationary point also satisfies the equilibrium equations, and the points xˆ and z are in phase equilibrium (though they are not necessarily at the globally stable phase equilibrium). So if the overall composition z happens to be on the binodal curve, the stationary points from the stability test also provide an exact LLE solution. Otherwise, they are useful as starting points for a subsequent LLE calculation. From inequality (6), it follows that, if the activity of any component, ziγi(z), i ) 1, ..., c, is greater than unity, the phase of overall composition z is unstable (Heidemann, 1983; Cairns and Furzer, 1990). This is because the tangent plane distance function calculated at the point corresponding to pure component i, F(x)ei,z), is less than zero if ziγi(z) > 1. If this happens, there is at least one point in the composition space where the tangent plane distance function is negative. This simple stability test is very fast but not absolute. The phase of composition z is certainly unstable if inequality (6) is violated for any pure component. However, it may still be unstable even if the inequality is satisfied at all points corresponding to the pure components. It should also be noted that the Hessian matrices of F(x) and g(x) are identical and equal to the Jacobian of the gradient of the tangent plane distance function

H[F(x)] ) H[g(x)] ) J[∇F(x)]

(9)

This is because the term L(x;z) in eq 1 is linear in each xi and has no influence on the second derivatives of F(x) with respect to the c - 1 independent mole fractions. Thus, the Hessian matrix of g(x) is positive definite at stationary points corresponding to minima in the tangent plane distance function, and these stationary points are stable or metastable phases. Conversely, all maxima and saddles of the tangent plane distance function lie in the region of unstable phases. Therefore, if g(z) is positive definite, the mixture will be either metastable or globally stable, and the F-test is needed to decide which. If g(z) is not positive definite, the mixture is unstable, and the F-test is needed to generate initial conditions for subsequent liquid-liquid calculations, as discussed below. In either case, the local stability test applied at the overall liquid composition should be augmented by the tangent plane test. The major difficulty in implementing the stability test is that the number of minima of F(x) (or the number of points satisfying Kuhn-Tucker conditions) is unknown a prior, e.g., Michelsen (1982a,b) or Smith et al. (1993). The algorithm presented here locates all the stationary points of F(x) (maxima, minima, and saddles) by tracking the ridges and valleys in F(x) and continuing the exploration of the function domain until the results are globally consistent with a certain topological relationship developed below. This approach has the advantage that the topological constraint sometimes reveals that one or more stationary points remain to be found. This situation often occurs in difficult problems. We conjecture that the ridges and valleys connect every stationary point of F(x) and that they can be completely explored using a line search no matter how many components are present in the mixture.

Solution Methods Michelsen (1982a) introduced new variables Yi ≡ xi exp(-β) where β ) ln xˆ c + ln γc(xˆ ) - ln zc - ln γc(z), and rewrote eq 7 as

ln Yi + ln γi - ln zi - ln γi(z) ) 0

i ) 1, ..., c (10)

The new independent variables, Yi, are related to the c mole fractions by xi ) Yi/∑j)1 Yj. Stationary points can be located as solutions to eq 10, and stability is verified c if β g 0 at all stationary points, which implies ∑i)1 Yi e c 1. Conversely, if ∑i)1Yi > 1 at any of the stationary points, the mixture is unstable. One of the methods suggested by Michelsen (1982a) is to solve eq 10 by the following successive substitution scheme:

) exp(ln zi + ln γi(z) - ln γti) Yt+1 i

(11)

In order to increase the chances of finding all minima of F(x), Michelsen (1982a) also suggested a number of starting points for solving eq 11 (in addition to the pure component compositions). There are also several other numerical methods suggested, but no single method seems to be consistently successful in every case. Cairns and Furzer (1990) describe a Newton-Raphson method to obtain solutions of eq 10. One of the key ideas is to initialize a search for the stationary point using the component with the largest activity, ziγi(z). To select the second starting point, Cairns and Furzer (1990) use the component with the second largest activity. This initialization of the phase compositions is similar to Gautam and Seider (1979), Walraven and Van Rompay (1988), Wasylkiewicz (1992), and Pajak et al. (1993). The reliability of the Michelsen (1982a) implementation of the stability test can be significantly improved by a homotopy-continuation method which locates multiple minima of F(x) along a single path (Sun and Seider, 1995). A robust method requires multiple starting points as described by Sun and Seider (1995). An interesting alternative approach has been proposed by McDonald and Floudas (1994, 1995a,b, 1996). They use a branch and bound method in a global optimization algorithm and show how it is possible to find solutions for the phase and chemical equilibrium problem, when the liquid phase is modeled using NRTL or UNIQUAC equations (McDonald and Floudas, 1994, 1995a). They have also applied the method to UNIFAC, Wilson, Modified Wilson, and ASOG equations (McDonald and Floudas, 1995b, 1996). A different approach to solving phase stability problems has recently been proposed by Stadtherr et al. (1994). They describe the application of interval methods which offer the promise of more robust algorithms for liquid-liquid equilibrium calculations. Another approach, which has been applied successfully to twophase flash problems, is based on homotopy continuation and path following (DeGance, 1993). Further development, testing, and comparison of these methods with other approaches will be interesting. We now develop a new approach for locating the stationary points that will exploit the fact that the minima in F(x) are located in the interior of the c permissible region (xi g 0, ∑i)1 xi ) 1).

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1399 Table 1. Properties of the Gradient System (Eq 12) 1. The singular points of eq 12 occur at stationary points in F(x) and vice versa. 2. The Jacobian is a symmetric matrix and has all real eigenvalues. This is because the Jacobian of the gradient ∇F(x) of the tangent plane distance function is equal to the Hessian matrix of F(x) (eq 9). Consequently, the singular points of eq 12 can be only nodes or saddles (or a combination of them if the singular point is nonelementary, i.e., a bifurcation point). 3. Inspection of the eigenvalues of J[∇F(x)] reveals the following: 3.1. Unstable nodes occur at maxima in F(x). 3.2. Stable nodes occur at minima in F(x). 3.3. Saddle points occur at saddles in F(x). 4. Ridges in F(x) correspond to the unstable manifolds of eq 12, and valleys correspond to the stable manifolds. These manifolds are (c - 2)-dimensional surfaces that are tangent to the eigendirections of J[∇F(xˆ )] at each stationary point of F(x). Trajectories of eq 12 converge toward stable manifolds and diverge from unstable manifolds. 5. The flow of trajectories is always downhill (i.e., F(x) decreases along every solution curve). 6. The flow of trajectories is always into the permissible region.

Our approach uses ideas from differential geometry and the theory of differential equations. First, consider the differential equation

dx/dt ) -∇F(x)

(12)

which is a gradient system with some useful properties (see Table 1). Properties 1-5 (Table 1) are common to all gradient systems, but property 6 is a special property of the tangent plane distance function, F(x). It is this special property, combined with the first four, that allows us to develop the topological criterion for n dimensions (see Appendix B). This topological criterion for a ternary mixture is

Nmax + Nmin - Nsad ) 1

(13)

where Nmax is the number of maxima, Nmin is the number of minima, and Nsad is the number of saddles of F(x). By integrating eq 12 along the eigendirections, we are able to track the ridges and valleys of F(x) to the stationary points. This enables a systematic exploration of the function domain rather than repeated calculations starting from various points generated from heuristics. Finding Stationary Points. In order to find all minima of F(x), we search for all stationary points of F(x) by integrating eq 12 and then check the consistency relationship (eq 13), continuing the exploration of the function domain until all ridges and valleys have been tracked. To locate the stationary points of F(x), we track the ridges and valleys in F(x). Because the flow of trajectories is always into the permissible region, the ridges and valleys must be contained completely inside this region; i.e., a valley or ridge cannot leave the domain and then reenter the domain at another point of the boundary. Valleys intersect each other at minima, ridges intersect each other at maxima, and valleys intersect ridges at saddles. Since the stationary points are isolated and the flow of trajectories is always downhill, neither valleys nor ridges can form closed contours, and therefore isola (i.e., disconnected stationary points) cannot occur. We conjecture that the ridges and valleys of F(x) are all connected to each other and that this connected domain contains all the stationary points. Therefore, it is sufficient to explore only this connected domain in order to locate all the stationary points. Such a search has a known stationary point, since the overall composition, z, is a stationary point of F(x) (see Eq 3). In describing the algorithm, we also explain why we believe a one-dimensional search is sufficient to completely explore the connected domain no matter how many

dimensions (components) are present. We describe tests for two, three, and four components in this paper. Valleys and ridges are tangent to the eigendirections of J[∇F(xˆ )] or H[F(xˆ )] at the stationary points, xˆ (see Appendix A and Figure 11). So we can begin the tracking starting from the trivial solution, z, and then continue from each of the newly found stationary points, xˆ . The search moves out from each stationary point along 2(c - 1) directions, +Ei and -Ei, i ) 1, ..., c - 1, determined by the eigenvectors of the Jacobian matrix J[∇F(xˆ )]. Starting from a minimum, where all eigenvalues λi, i ) 1, ..., c - 1, are positive, we select an eigenvector and integrate eq 29. Starting from a maximum (λi < 0, i ) 1, ..., c - 1), we integrate eq 28. If the starting point is a saddle, we integrate eq 29 when the eigenvalue corresponding to the chosen eigenvector is positive or eq 28 if the eigenvalue is negative. The procedure is discussed in more detail in Appendix A. Tracking the ridges and valleys is very efficient in searching for all stationary points of F(x). It is also very useful for LLE calculations since it provides good initial estimates for the equilibrium phase compositions. Combined Phase Equilibrium and Stability Algorithm Liquid-liquid equilibrium calculations can easily fail without good initial estimates for the number of liquid phases and their composition. The stability test generally does not deliver the exact compositions of phases in equilibrium or even necessarily the correct number of liquid phases. Fortunately, it provides good initial estimates for the equilibrium phase compositions, xIi , xII i , etc., from which it is possible to determine the exact number and composition of liquid phases in globally stable equilibrium. Our algorithm is summarized as follows. 1. Specify a temperature, T, pressure, P, and overall liquid composition, z. 2. Find all stationary points of the tangent plane distance function using the algorithm in Table 2 and check the stability of the liquid. If the liquid is globally stable, then end. Otherwise, go to step 3. 3. Save the coordinates of all minima of F(x). 4. Assume that the number of liquid phases is equal to the number of minima of F(x) (up to the limit allowed by the phase rule) and calculate the LLE initializing the phase compositions using the coordinates of the minima of F(x). The Rachford-Rice formulation (1952) solved with a Newton-Raphson method as described by Wasylkiewicz (1992) and Pajak et al. (1993) was used for the LLE calculations. Note that these calculations may converge to an equilibrium solution containing fewer phases than the number initially assumed if one or more of the phase fractions converge to zero.

1400

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996

Table 2. Search for All Stationary Points of the Tangent Plane Distance Function 1. Take the composition of the original mixture, z, as a starting point, xˆ . Note that F(z) ) 0 and ∇F(z) ) 0, so this is a stationary point of F(x). 2. Calculate the Hessian matrix, H[F(xˆ )], its eigenvalues, λi, and eigenvectors, Ei, i ) 1, ..., c - 1. 3. Add xˆ to the set of stored stationary points (SSP). 4. Add xˆ , λi, and (Ei, i ) 1, ..., c - 1, to the set of stored initial conditions (SIC). 5. Take one point from SIC and integrate the set of c - 1 ordinary differential equations

dx/dt ) -∇F(x)

(28)

dx/dt ) +∇F(x)

(29)

dx/dt ) -H[F(x)]T∇F(x)

(30)

or

or

Integrate up to a new stationary point xˆ NEW or to a constraint from the c linear constraints: xi ) 0, i ) 1, ..., c. 6. If the point xˆ NEW is not in SSP, go to step 2; else continue. 7. If SIC is not empty go to step 5; else continue. 8. Check the topological criterion (13). If satisfied, STOP. Otherwise, adjust the tolerances δ and σ in Table 4 and repeat calculations.

5. Check the global stability of one of the liquid phases present at the calculated equilibrium. (Since all phases in equilibrium have a common tangent plane to the Gibbs free energy surface, it necessarily follows that F(xIeq) is identical to F(xII eq), etc. Therefore, checking the stability of only one of the equilibrium phases is sufficient.) If the calculated equilibrium is globally stable, then end. Otherwise, go to step 6. 6. Save the coordinates of all minima of F(x) found in the last stability test. Go to step 4. Examples Example 1: A Simple Two-Phase System. Using the methods described above, we have calculated the isothermal liquid-liquid binodal and tie-lines for a mixture of methanol, hexane, and heptane at 55 °C using the NRTL model with the binary interaction parameters reported by Sorensen and Arlt (1980), Vol. V/2 p 133. The results are shown in Figure 2a, and contours of the Gibbs free energy function (g/RT) are shown in Figure 2b. (Figure 2 can be computed in approximately 1 s on a DEC workstation (AXP 3000/ 300).) Now consider the stability of this mixture at 55 °C with the overall composition of 30 mol % methanol, 20 mol % hexane, and 50 mol % heptane. Contours of the tangent plane distance function for this overall liquid composition are shown in Figure 3a. There are two minima in F(x); one is the trivial solution at the overall composition of the mixture (where F(z) ) 0), and the second is near the point 80 mol % methanol and 7 mol % hexane. At this point, F(x) is less than zero, so the mixture is unstable. Since the trivial solution is a minimum, the overall composition of the liquid must be in the metastable region (see eq 9). Figure 3b illustrates the search for stationary points of F(x). Each small open circle corresponds to one step of integration, and each filled circle to a stationary point. Two minima are shown; one at the point z, where F(z) ) 0, and one additional negative minimum. One saddle is found between these minima, and the topological relationship (eq 13) is satisfied. The stability test provides good initial estimates for xIi and xII i . With these values, the Newton-Raphson method converges easily to the solution for this and all other overall compositions shown in Figure 2a. Example 2: Three Liquid Phases and Three Minima of the Tangent Plane Distance Function.

Figure 2. (a) Isothermal liquid-liquid equilibrium prediction and (b) contours of the Gibbs free energy surface for the mixture methanol, hexane, and heptane at 55 °C. Based on the NRTL model with binary interaction parameters reported by Sorensen and Arlt (1980), Vol. V/2, p 133.

We have calculated the isothermal equilibrium envelope for the mixture n-propanol, water, and toluene at 5 °C, with the results shown in Figure 4a. (Figure 4 can be

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1401

Figure 3. (a) Contours of the tangent plane distance function and (b) the search for stationary points along ridges and valleys for the mixture methanol, hexane, and heptane at 55 °C. The overall composition is 30 mol % methanol, 20 mol % hexane, and 50 mol % heptane.

Figure 4. (a) Isothermal liquid-liquid equilibrium prediction and (b) contours of the Gibbs free energy surface for the mixture n-propanol, water, and toluene at 5 °C. Based on the NRTL model with binary interaction parameters reported by Sorensen and Arlt (1980), Vol. V/2, p 580.

computed in approximately 1 min on a DEC workstation (AXP 3000/300).) The model predicts three distinct twophase regions bordering one three-phase region. Although the experimental data do not reveal any three liquid-phase region, it is, nevertheless, the most stable solution of the model for some compositions as shown in Figure 4a. Three liquid phases are permissible by the phase rule and are known to occur for some mixtures, e.g., water, hexane, and phenol (Heidemann and Mandhane, 1975); lauryl alcohol, glycol, and nitromethane (Heidemann and Mandhane, 1975; Gautam and Seider, 1979); water, nitromethane, and hexanol (Sorensen and Arlt (1980) Vol. V/2, p 69); bis(2-methylpropyl) ether, phosphoric acid, and water (Marcilla et al., 1994). Contours of the Gibbs free energy function (g/RT) for the mixture are shown in Figure 4b. Now consider an overall composition of 38 mol % n-propanol, 20 mol % water, and 42 mol % toluene at 5 °C. Contours of the tangent plane distance function are shown in Figure 5a. There are three minima of F(x); one occurs at the overall composition where F(z) ) 0; the other two have F(x) less than zero, so the mixture

is unstable (locally metastable). Figure 5b illustrates the search for stationary points of F(x). Each small open circle corresponds to one step of integration, and each filled circle, to a stationary point. Usually the integration step is at its maximum value of 0.01, but this is reduced near the saddle or close to the overall composition. There are two saddles between the three minima, and the global relationship (eq 13) is satisfied. The III stability test gave good estimates for xIi , xII i , and xi . With these values, the Newton-Raphson method converged readily to the three liquid-phase solution (Figure 4a). Example 3: Three Liquid Phases and Two Minima of the Tangent Plane Distance Function. Because the tangent plane at a point z is usually not parallel to the plane suspended on points corresponding to compositions of liquids in equilibrium, the minima of F(x) can differ significantly from the equilibrium points and even their number may not be the same as the number of liquid phases. So, if one finds three minima of F(x) and at least one of them is negative, it does not necessarily mean that the point z is in a three-

1402

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996

Figure 5. (a) Contours of the tangent plane distance function and (b) the search for stationary points in the mixture n-propanol, water, and toluene at 5 °C for an overall composition of 38 mol % n-propanol, 20 mol % water, and 42 mol % toluene.

Figure 6. (a) Contours of the tangent plane distance function and (b) the search for stationary points in the mixture n-propanol, water, and toluene at 5 °C for an overall composition of 30 mol % n-propanol, 20 mol % water, and 50 mol % toluene.

phase region. Also, if there are two minima of F(x) and at least one of them is negative, the point z can be in either a two- or three-phase region. The number of liquid phases in equilibrium cannot be determined on the basis of tangent plane analysis alone. The overall liquid composition in Figure 5 is very close to the composition of one of the three liquid phases in equilibrium (see Figure 4a). Now consider the same mixture of n-propanol (1), water (2), and toluene (3) with another overall composition z1 ) 0.3, z2 ) 0.2, and z3 ) 0.5, which is inside the three-phase region in Figure 4a. Contours of the tangent plane distance function are shown in Figure 6a. The overall liquid composition z is marked in the figure by a filled circle. The procedure, which searches for all stationary points of F(x) (Appendix A), finds two minima (both negative) and one saddle between them at the point z, where F(z) ) 0 (Figure 6b). The global relationship (eq 13) is satisfied. Since the tangent plane distance function F(x) is negative at the minima, the mixture is unstable, but we still do not know how many equilibrium liquid phases will be formed.

The Newton-Raphson method was initiated with the compositions of the minima taken from the previous search for all stationary points of the tangent plane distance function. These correspond to both minima of F(x): one at z1 ) 0.0192 and z2 ) 0.9800 and the other at z1 ) 0.1081 and z2 ) 0.0263. The Newton-Raphson routine converged (sum of squared residuals less than 10-12) in five steps. The resulting tie-line is shown in Figure 7a. The stability test reveals that the composition at one end of the tie line (point no. 4 in Figure 7a) is not stable. There are three negative minima of F(x) at points 1, 3, and 5 and two saddles at points 2 and 4 which correspond to the previously calculated tie-line. It should also be noted that there is no need to perform a stability check for the second, third, etc., liquid phase in “equilibrium” since the points representing these liquids must have a common tangent plane (necessary conditions for equilibrium) even for an unstable tie line. So the number, locations, and values of all stationary points of F(x) must be the same for all liquid phases that satisfy the necessary conditions for equilibrium.

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1403

Figure 8. Stability test for the tie line for the mixture n-propanol, water, and toluene at 5 °C. The overall liquid composition corresponds to point 4. There are three minima of the tangent plane distance function at points 1, 3, and 5 and two saddles at points 2 and 4.

Figure 7. (a) Stability test of the tie line and (b) contours of the tangent plane distance function for the mixture n-propanol, water, and toluene at 5 °C. The overall liquid composition corresponds to point 4. There are three minima of the tangent plane distance function at points 1, 3, and 5 and two saddles at points 2 and 4.

Since the tie line is not globally stable, the next step is to make the three-phase LLE calculations. The initial phase compositions have been taken from the previous search for all stationary points of the tangent plane distance function and correspond to three minima of F(x): the first minimum is at z1 ) 0.0370 and z2 ) 0.9617, the second minimum is at z1 ) 0.4072 and z2 ) 0.3427, and the third minimum is at z1 ) 0.2182 and z2 ) 0.0572. The Newton-Raphson scheme converged in five steps, and the resulting three-phase equilibrium is shown in Figure 4a. This equilibrium point is globally stable. Example 4: Multiple Two-Phase Tie Lines with the Same Overall Composition. Now consider the same mixture of n-propanol (1), water (2), and toluene (3) with another overall liquid composition of z1 ) 0.3, z5 ) 0.2, and z2 ) 0.5 which lies inside the two-phase region A in Figure 4a. We find two minima in F(x) (both negative) and one saddle between them at point z, where F(z) ) 0. The global relationship (eq 13) is satisfied. Since the tangent plane distance function F(x) is negative at some points, the liquid mixture of this

overall composition is unstable; in this case it lies inside the spinodal region. The Newton-Raphson scheme was initiated with the compositions of the minima found in the search for all stationary points of the tangent plane distance function. These minima of F(x) occur at z1 ) 0.0587 and z2 ) 0.0114 and the other at z1 ) 0.0181 and z2 ) 0.9809. The Newton-Raphson scheme converged in eight steps, with the resulting tie line shown in Figure 8 connecting points 2 and 4. The stability test reveals that one end of this tie line (point no. 4 in Figure 8) is unstable. There are three negative minima of F(x) at points 1, 3, and 5 and two saddles at points 2 and 4 which correspond to the previously calculated tie-line compositions. Because the tie line is not globally stable and three negative minima of F(x) have been found, three-phase LLE calculations are performed. The initial phase compositions correspond to the three minima of F(x): the first at z1 ) 0.2424 and z2 ) 0.0698 (point no. 1 in Figure 8), the second at z1 ) 0.4089 and z2 ) 0.3397 (point no. 3), and the third at z1 ) 0.0377 and z2 ) 0.9610 (point no. 5). After two iterations, the NewtonRaphson calculations reach a constraint (zero-phase fraction for the first phase) and the number of liquid phases is reduced to two. The subsequent NewtonRaphson calculations converge in two steps. The resulting tie line is shown in Figure 9. The stability test for one of the equilibrium phases on this tie line indicates that it is globally stable. For the overall composition in the two-phase region, we have found two tie lines: one unstable, which connects points 1 and 2 in Figure 9, and another globally stable, which connects points 3 and 4 in Figure 9. Sections of the Gibbs free energy function (g/RT) along both tie lines are shown in Figure 9b. The lowest Gibbs free energy is for the mixture of liquids at points 3 and 4, and this is the stable solution; it could not be found without previous stability testing of the first unstable tie line 1-2. There is no guarantee that a heuristic initialization of the phase compositions will result in convergence to the globally stable tie line and not to the unstable one.

1404

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 Table 3. Four-Component Liquid-Liquid Equilibrium for the Mixture n-Heptane + Toluene + o-Xylene + Propylene Carbonate at 25 °Ca propylene heptane toluene o-xylene carbonate overall composition phase I (calculated) phase I (experimental) phase II (calculated) phase II (experimental) overall composition phase I (calculated) phase I (experimental) phase II (calculated) phase II (experimental)

46.96 2.18 1.61 92.85 92.30 31.27 2.62 3.08 59.24 59.46

1.92 1.10 0.91 2.76 2.93 10.85 5.89 6.84 15.77 14.87

2.81 0.969 1.12 4.36 4.49 15.97 6.65 8.13 23.38 23.81

48.31 95.74 96.36 0.024 0.28 41.91 84.83 81.95 1.59 1.86

a Table entries are in wt %; data and UNIQUAC model parameters are from Bakar et al. (1994).

At a higher temperature (50 °C), the model predicts two separate two-liquid regions (Figure 10b) and no three-phase region. The bifurcation temperature, where the phase diagram changes qualitatively from three twophase and one three-phase regions to two two-phase regions, is predicted to be between 46.45 °C (Figure 10c) and 46.46 °C (Figure 10d). Example 6: A Four-Component Mixture. We consider a mixture of n-heptane + toluene + o-xylene + propylene carbonate which has been studied experimentally by Bakar et al. (1994). The UNIQUAC model was used for calculations at 25 °C with the interaction parameters reported by Bakar et al. (1994). Our method successfully completed liquid-liquid equilibrium calculations for all of the tie lines reported by Bakar et al. (1994). A comparison of calculated vs experimental values for two of the tie lines is shown in Table 3. Use of the model in the complete design of a 10-stage liquid extractor, solving both mass balances and liquidliquid equilibrium, requires 10-15 s on a DEC workstation (AXP 3000/300). Other Examples Figure 9. (a) Isothermal liquid-liquid equilibrium prediction of the NRTL model and (b) the Gibbs free energy of the unstable tie line (points 1 and 2) and the globally stable tie line (points 3 and 4) connected with the overall liquid composition of 30 mol % n-propanol, 50 mol % water, and 20 mol % toluene at 5 °C.

This effect leads to one sort of difficult problem for other methods and such cases are not uncommon. For example, there are many points in regions A and B of Figure 4a for which two tie lines can be found at the same overall composition. In addition, there are many points in the shaded region of Figure 4a that generate a two-phase tie line as well as a three-phase solution from the same overall composition. Example 5: Bifurcation Temperature. Using the UNIQUAC binary interaction parameters reported by Sorensen and Arlt (1980, Vol. V/3, pp 422 and 423) we find the three-phase region for the mixture shown in Figure 10a. The region has been detected experimentally at 21 °C, and the experimental data are presented in the DECHEMA Chemistry Data Series (Sorensen and Arlt (1980), Vol. V/2, p 69). However, no parameters for thermodynamic models were estimated to match this ternary data, and we used UNIQUAC parameters (Sorensen and Arlt (1980), Vol. V/3) estimated from binary miscibility data.

We have tested this method successfully on a number of other problems but cannot present a detailed description of these here due to space limitations. A brief selection is as follows; the first four are described as difficult by Swank and Mullins (1985). 1. Acrylonitrile + ethyl cyanide + water at 25 °C. 2. Heptane + triolein + furfural at 70 °C. 3. Ethyl acetate + 2-butanol + water at 0 °C. 4. 2-Butanol + tert-butanol + water at 60 °C. 5. The mixture of di-sec-butyl ether + sec-butanol + water, with model parameters from Kovach and Seider (1988), Table XI (which we believe to be the same set used by Sun and Seider, 1995). For this mixture, Sun and Seider (1995) report a number of failures for their method and for a method based on Michelsen’s (1982a,b) approach. The method we describe in this paper experienced no failures in computing liquid-liquid equilibria at 90 °C and 1.16996 atm at overall liquid compositions on a uniform 25 × 25 grid covering the composition space (cf. Sun and Seider, 1995). (We note that some of these points are actually above the boiling temperature surface, and, therefore, the liquid-liquid equilibria computed here would represent only part of a solution strategy for vapor-liquid-liquid phase equilibrium. Sun and Seider include the possibility of a vapor phase, but we have not in this work.)

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1405

Figure 10. Isothermal liquid-liquid equilibrium prediction for the mixture water, nitromethane, and hexanol. Based on the UNIQUAC model with binary interaction parameters from Sorensen and Arlt (1980, Vol. V/3, pp 422 and 423). The three liquid phase region (a) has been experimentally detected (Sorensen and Arlt (1980), Vol. V/2, p 69). The bifurcation temperature is between 46.45 (c) and 46.46 °C (d).

Conclusions A robust global algorithm for implementation of the Gibbs tangent plane analysis for complex liquid mixtures involves locating all the stationary points of the tangent plane distance function, F(x), by tracking its ridges and valleys. We do a systematic exploration of the function domain rather than repeated calculations starting from various points generated using heuristics. This is based on integration of differential equations up to steady states (singular points of eq 12), which correspond to the stationary points of F(x). This approach has the great advantage that it is possible to test whether the results are globally consistent. Since the tangent plane at an overall composition z is usually not parallel to the plane suspended on points corresponding to compositions of liquids in equilibrium, the minima of F(x) can differ significantly from equilibrium points and even their number may not be the same as the number of liquid phases in equilibrium. Therefore, the number of liquid phases in equilibrium cannot be determined on the basis of tangent plane analysis alone, and the combined phase equilibrium and stability algorithm should be used.

The Hessian matrices of F(x) and g(x) are identical. Therefore, all minima of the tangent plane distance function lie in the region where the Hessian matrix of g(x) is positive definite and represent stable or metastable phases. Consequently, the results from the stability analysis are very useful as initial estimates in LLE calculation for the compositions of liquid phases in equilibrium. There is no need to check the stability of all liquid phases in equilibrium since the points representing other liquids must have a common tangent plane (necessary conditions for equilibrium). In other words, the number, locations, and values of all stationary points of F(x) must be the same for each liquid phase in equilibrium and a stability test for any one of the compositions satisfying the LLE is sufficient. Acknowledgment We are grateful to E. I. DuPont Co., Eastman Kodak, Co., ICI, and Union Carbide Corp. for financial support and to Michael Minotti for preparing several test problems.

1406

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 * ) trial point I, II, ... ) phase number I, number II, etc.

Appendix A. Method for Finding Stationary Points A stepwise description of the procedure is given in Table 2. Equation 30 is used near saddles because it is practically impossible to find them by integrating eq 28 or eq 29. Near saddles we recast the original problem as a least-squares problem

1 T min [∇F(x)] ∇F(x) 2 x

(14)

The gradient of the new objective function (eq 14) can be expressed as (see eq 9)

Figure 11. Two valleys intersecting at a minimum (xˆ ) of the tangent plane distance function, the eigenvectors at the minimum (E1 and E2), and the predictor and corrector steps of the valley tracking algorithm.

Nomenclature c ) number of components D ) diagonal matrix e ) unit vector F ) tangent plane distance function ∇F ) gradient of F function g ) molar Gibbs free energy H ) Hessian matrix ind{‚} ) index of a point in a vector field J ) Jacobian matrix L ) oriented line in a vector field L(x;z) ) tangent plane to the Gibbs energy surface at composition z N ) number of stationary points n ) dimension (n ) c - 1) R ) gas constant T ) temperature P ) pressure X ) similarity transformation matrix x ) mole fractions of liquid phase xˆ ) mole fractions of liquid phase at a stationary point of F z ) overall composition of a liquid mixture |‚| ) Euclidean distance Greek Symbols γ ) activity δ ) scaling factor E ) eigenvector λ ) eigenvalue µ ) chemical potential σ ) tolerance Subscripts i ) component index max ) maximum min ) minimum sad ) saddle | ) parallel Superscripts even ) even number of negative eigenvalues odd ) odd number of negative eigenvalues t ) iteration index

{

}

1 ∇ [∇F(x)]T∇F(x) ) J[∇F(x)]T∇F(x) ) 2 H[F(x)]T∇F(x) (15) We show in Appendix C that the Jacobian of the gradient of the least-squares function (eq 14), J{H[F(x)]T∇F(x)}, is a symmetric matrix with real, nonnegative eigenvalues at all stationary points of F(x). This means that the transformation of the tangent plane distance function to the least-squares function turns every stationary point of F(x) into a stable node of eq 30, and the scheme described by eq 30 can pick up saddles of F(x) as well as other stationary points. Note that the least-squares function must have additional saddles between the stable nodes, but these saddles are not singular points of the original function F(x). Equation 30 involves more computation than eq 28 or eq 29 because of the more complicated right-hand side, and it is used only near stationary points which are saddles of F(x). Proper switching from eq 28 or eq 29 to eq 30 is important. It cannot be done too early since the integration can be attracted to the stationary point used as an initial condition. In order to switch properly, we have to be sure that we are in the attraction region of a new stable node of (1/2)[∇F(x)]T∇F(x) which corresponds to a new stationary point of F(x). To achieve this, we calculate the norm of ∇F(x) and check if it increases or decreases along the valley floors or ridge-top loci. A local minimum of |∇F(x)| along the integration path indicates that we are close to a saddle of F(x) within the attraction region of a new stable node of (1/2)[∇F(x)]T∇F(x). The second criterion of switching depends on the value of the cosine of the angle between the old local direction and the new gradient of F(x) at each step of integration. If this value is close to zero (the vectors are almost orthogonal) or its sign has changed (new point is behind a new stationary point), we switch to integrating eq 30. By integrating eq 28 or eq 29 near a stable node, we can easily track valley floors or ridge-top loci. However, when starting the calculations from an unstable node or when going to a saddle, it is practically impossible to track ridges or valleys by simple integration without a special correction technique. We have developed one such technique based on the definition of a ridge (valley) as the maximum (minimum) of F(x) taken in the constraint plane normal to the ridge-top locus (valley floor) (also see Van Dongen and Doherty (1984)). The

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996 1407 Table 4. Valley Tracking Algorithm for One Step of Integration 1. Calculate ∇F(xt) and the increment ∆x ) δ∇F(xt)/|∇F(xt)| (see Figure 11) for a point xt on the valley, where δ is a scaling factor (usually in the range 10-6-10-2). 2. Calculate the trial point x* ) xt + ∆x. 3. Calculate the gradient of the function at the trial point ∇F(x*). 4. Calculate the cosine of the angle between the vectors ∆x and ∇F(x*)

cos R )

∆x‚∇F(x*) |∆x||∇F(x*)|

(31)

(a) If cos R < σ (σ is usually in the range 10-4-10-2), then reduce δ and go to step 1. (b) If cos R > 1 - σ, then put xt+1 ) x*. Increase δ if it has been reduced previously. The point xt+1 is the next point on the valley. (c) If σ e cos R e 1 - σ, then calculate the vector ∇F(x*)|

∇F(x*)| ) ∇F(x*) -

∆x‚∇F(x*) ∆x |∆x||∆x|

(32)

which is the orthogonal projection of ∇F(x*) onto the constraint plane normal to the vector ∆x at point x*. 5. Corrector step. Find the minimum of F(x) in the direction -∇F(x*)|. The minimum is the (corrected) next point, xt+1, on the valley.

algorithm used in this work to track valleys is summarized in Table 4. Appendix B. Topological Criterion for n Dimensions

N

For a vector field in the plane (n ) 2), the rotation of the vector field ∇F(x) (see eq 12) along an oriented line L (not containing singular points of the field ∇F(x)) is defined as an increment of the angle (divided by 2π and computed counterclockwise) which forms the vector ∇F(x) with a specified direction when the point x crosses the line L according to its orientation (Okinski (1992), p 207). If we compute the rotation of the vector field ∇F(x) for a given point x in such a way that L is a closed line and in the region contained by L there are no singular points of the field (except possibly the point x), then such a quantity is called the index of the point x (ind{∇F(x)}). The index of a nonsingular point is zero. The rotation of the vector field ∇F(x) computed along a closed line L is equal to the sum of the indices of the stationary points (which are generally isolated points) contained within the region determined by L. The index of a nondegenerate stationary point xˆ can be determined from the eigenvalues, λj, j ) 1, ..., c - 1, of the Jacobian matrix J[∇F(xˆ )]

where N is a number of stationary points of F(x). From eqs 17-21 it follows that even odd Nmax + (-1)nNmin + Nsad - Nsad ) (-1)n (22)

where Nmax is the number of maxima of F(x), Nmin is even the number of minima, Nsad is the number of saddles with an even number of negative eigenvalues, and odd Nsad is the number of saddles with an odd number of negative eigenvalues of the Jacobian matrix J[∇F(xˆ )]. The topological criterion (eq 22) for n dimensions (c - 1 components) reduces to eq 13 for ternary mixtures (n ) 2, all saddles have an odd number of negative eigenvalues) and to the equation

Nmax - Nmin ) -1

(17)

The Jacobian of the gradient (eq 15) of the leastsquares problem function (eq 14) can be calculated as follows

({

1 J ∇ [∇F(x)]T∇F(x) 2

(18)

T

∂ {J[∇F(x)]T}∇F(x) (24) ∂x

At all stationary points of F(x), ∇F(x) ) 0. So at each stationary point xˆ , the Jacobian of the gradient of the least-squares problem function is

({

}) ) J[∇F(xˆ )] J[∇F(xˆ )] )

1 J ∇ [∇F(xˆ )]T∇F(xˆ ) 2

T

(19)

H[F(xˆ )]TH[F(xˆ )] ) H[F(xˆ )]2 (25)

(20)

(see eq 15). Since the Hessian matrix H[F(x)] is a real symmetric matrix and has all real eigenvalues, the Jacobian of the gradient of the least-squares problem

If the number of negative eigenvalues is odd odd ) -1 ind{∇F(xˆ )}sad

}) ) J(J[∇F(x)] ∇F(x)) )

J[∇F(x)]TJ[∇F(x)] +

where n (n ) c - 1) is the dimension of the vector field ∇F(x). Some eigenvalues at a saddle point are positive and some are negative. If the number of negative eigenvalues is even even ind{∇F(xˆ )}sad )1

(23)

Appendix C. Properties of the Jacobian

Since all eigenvalues at a minimum are negative

ind{∇F(xˆ )}min ) (-1)n

(21)

(16)

The singular points of the gradient system (12) can be only maxima, minima, or saddles. Since all eigenvalues at a maximum are positive

ind{∇F(xˆ )}max ) 1

ind{∇F(xˆ i)} ) (-1)n ∑ i)1

for binary mixtures (n ) 1, no saddles).

c-1

λj ∏ j)1

ind{∇F(xˆ )} ) sgn

From the property of the rotation of the vector field ∇F(x) and property (6) of the gradient system (12), it follows that the sum of indices of all stationary points contained within the permissible region is

1408

Ind. Eng. Chem. Res., Vol. 35, No. 4, 1996

function is also a symmetric matrix and has all real eigenvalues at each stationary poing xˆ . Diagonalization of the Hessian matrix H[F(xˆ )]

D ) X-1H[F(xˆ )]X

(26)

gives a diagonal matrix D, with the eigenvalues of H[F(xˆ )] as the entries on the main diagonal. Here X is the matrix with eigenvectors of H[F(xˆ )] as column vectors. Diagonalization of the square of the Hessian matrix H[F(xˆ )]2

D2 ) X-1H[F(xˆ )]2X

(27)

gives a diagonal matrix D2, with the eigenvalues equal to squares of eigenvalues of H[F(xˆ )] as the entries on the main diagonal. So, at each stationary point xˆ , the eigenvalues of the Jacobian of the gradient of the leastsquares problem function are equal to squares of the eigenvalues of the Hessian matrix H[F(xˆ )]. This means that the Jacobian has all nonnegative eigenvalues at all stationary points of F(x). Literature Cited Bakar, A.; Salem, S. H.; Hamad, E. Z.; Al-Naafa, M. A. Quaternary Liquid-Liquid Equilibrium of n-Heptane + Toluene + o-Xylene + Propylene Carbonate. Ind. Eng. Chem. Res. 1994, 33, 689692. Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Petrol. Eng. AIME 1982, 22, 731-742. Cairns, B. P.; Furzer, I. A. Multicomponent Three-Phase Azeotropic Distillation. 2. Phase-Stability and Phase-Splitting Algorithms. Ind. Eng. Chem. Res. 1990, 29, 1364-1382. Chien, H. H. Formulations for Three-Phase Flash Calculations. AIChE J. 1994, 40, 957-965. DeGance, A. E. Ab Initio Equation of State Phase Equilibria Computations via Path Continuation. Fluid Phase Equilib. 1993, 89, 303-304. Gautam, R.; Seider, W. D. Computation of Phase and Chemical Equilibrium: Part I. Local and Constraint Minima in Gibbs Free Energy. Part II. Phase-Splitting. AIChE J. 1979, 25, 991-1006. Gibbs, J. W. A Method of Geometrical Representation of the Thermodynamic Properties of Substances by Means of Surfaces. Trans. Conn. Acad. Arts Sci., II 1873, 382-404. Reprinted in The Scientific Papers of J. Willard Gibbs; Dover: New York, 1961; Vol. 1, pp 33-54. Gibbs, J. W. On the Equilibrium of Heterogeneous Substances. Trans. Conn. Acad. Arts Sci. III 1876, 108-248. Reprinted in The Scientific Papers of J. Willard Gibbs; Dover: New York, 1961; Vol. 1. See pp 55-129, especially the section Geometrical Illustrations, pp 115-129. Heidemann, R. A. Computation of High Pressure Phase Equilibria. Fluid Phase Equilib. 1983, 14, 55-78. Heidemann, R. A.; Mandhane, J. M. Ternary Liquid-Liquid Equilibria: The Van Laar Equation. Chem. Eng. Sci. 1975, 30, 425-434. Jiang, Y.; Smith, W. R.; Chapman, G. R. Global Optimality Conditions and Their Geometric Interpretation for the Chemical and Phase Equilibrium Problem. SIAM J. Control Optim. 1995, 5, 813-834. Kovach, J. W., III; Seider, W. D. Vapor-Liquid and Liquid-Liquid Equilibria for the System sec-Butyl Alcohol-Di-sec-butyl EtherWater. J. Chem. Eng. Data 1988, 33, 16-20. Lucia, A.; Taylor, R. Complex Iterative Solutions to Process Model Equations? Comput. Chem. Eng. 1992, 16S, S387-S394. Marcilla, A. F.; Ruiz, F.; Sabater, M. C. Two-Phase and ThreePhase Liquid-Liquid Equilibrium for Bis(2-methylpropyl) Ether + Phosphoric Acid + Water. J. Chem. Eng. Data 1994, 39, 1418.

Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982a, 9, 1-19. Michelsen, M. L. The Isothermal Flash Problem. Part II. PhaseSplit Calculation. Fluid Phase Equilib. 1982b, 9, 21-40. Michelsen, M. L. Phase Equilibrium Calculations. What is Easy and What is Difficult? Comput. Chem. Eng. 1993, 17, 431439. McDonald, C. M.; Floudas, C. A. Decomposition Based and Branch and Bound Global Optimization Approaches for the Phase Equilibrium Problem. J. Global Optim. 1994, 5, 205-251. McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase Stability Problem. AIChE J. 1995a, 41, 1798-1814. McDonald, C. M.; Floudas, C. A. Global Optimization and Analysis for the Gibbs Free Energy Function using the UNIFAC, Wilson and ASOG Equations. Ind. Eng. Chem. Res. 1995b, 34, 16741687. McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1996, in press. Okinski, A. Catastrophe Theory, Comprehensive Chemical Kinetics; Elsevier-PWN: Warszawa, 1992; Vol. 33. Pajak, M.; Wasylkiewicz, S. K.; Go´rska, L. Mathematical Model of Extraction Equilibrium for Ethanol Extraction by Tributyl Phosphate. Inz. Chem. Procesowa 1993, 14, 375-392. Rachford, H. H., Jr.; Rice, J. D. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. J. Petrol. Technol. 1952, 4 (10), Section 1, p 19, and Section 2, p 3. Smith, J. M.; Van Ness, H. C. Introduction to Chemical Engineering Thermodynamics; McGraw-Hill: New York, 1987. Smith, J. V.; Missen, R. W.; Smith, W. R. General Optimality Criteria for Multiphase Multireaction Chemical Equilibrium. AIChE J. 1993, 39, 707-710. Sorensen, J. M.; Arlt, W. Liquid-liquid Equilibrium Data Collection. Ternary Systems. DECHEMA Chemistry Data Series; DECHEMA: Frankfurt, 1980; Vol. V. Stadtherr, M. A.; Schnepper, C. A.; Brennecke, J. F. Robust Phase Stability Analysis Using Interval Methods. Presented at Fourth International Conference on Foundations of Computer Aided Process Design (FOCAPD’94), Snowmass, CO, July 1994. Sun, A. C.; Seider, W. D. Homotopy-continuation Method for Stability Analysis in the Global Minimization of the Gibbs Free Energy. Fluid Phase Equilib. 1995, 103, 213-249. Swank, D. J.; Mullins, J. C. Evaluation of Methods for Calculating Liquid-Liquid Phase-Splitting. Fluid Phase Equilib. 1985, 30, 101-110. Van Dongen, D. B.; Doherty, M. F. On the Dynamics of Distillation ProcessessV. The Topology of the Boiling Temperature Surface and its Relation to Azeotropic Distillation. Chem. Eng. Sci. 1984, 39, 883-892. Van Dongen, D. B.; Doherty, M. F.; Haight, J. R. Material Stability of Multicomponent Mixtures and the Multiplicity of Solutions to Phase-Equilibrium Equations. I. Nonreacting Mixtures. Ind. Eng. Chem. Fundam. 1983, 22, 472-485. Van Ness, H. C.; Abbot, M. M. Classical Thermodynamics of Nonelectrolyte Solutions with Applications to Phase Equilibria; McGraw-Hill: New York, 1982. Walraven, F. F. Y.; Van Rompay, P. V. An Improved PhaseSplitting Algorithm. Comput. Chem. Eng. 1988, 12, 777-782. Wasylkiewicz, S. Modelling and Simulation of Multistage Separation Processes; Scientific Papers of the Institute of Chemical Engineering and Heating Equipment, Technical University of Wrocław: Wrocław, Poland, 1992; No. 65; Monograph No. 40. Widagdo, S.; Seider, W. D.; Sebastian, D. H. Dynamic Analysis of Heterogeneous Azeotropic Distillation. AIChE J. 1992, 38, 1229-1242.

Received for review December 21, 1994 Accepted December 29, 1995X IE950049R

Abstract published in Advance ACS Abstracts, February 15, 1996. X