1 ROBUST PHASE STABILITY ANALYSIS USING INTERVAL METHODS
ROBUST PHASE STABILITY ANALYSIS USING
INTERVAL METHODS
Mark A. Stadtherr and Carol A. Schnepper
Department of Chemical Engineering
University of Illinois
600 S. Mathews Ave.
Urbana, IL 61801, USA
Joan F. Brennecke
Department of Chemical Engineering
University of Notre Dame
Notre Dame, IN 46556, USA
Abstract
Conventional equation solving and optimization techniques for solving the phase stability problem may fail to converge or may converge to an incorrect result. A technique for solving the problem with mathematical certainty is needed. One approach to providing such assurance can be found in the use of interval methods. An interval Newton/generalized bisection technique is applied here to solve the phase stability problem. Results for two models of liquid-phase systems, using several different feed compositions, indicate that the technique used is reliable and very efficient.
Keywords
Phase stability, Interval computations, Nonlinear equation solving, Interval Newton method, Tangent plane criterion, NRTL equation
Introduction
The determination of phase stability is a recurrent problem in the computation of phase equilibria, and thus is especially important in the analysis and design of separation operations such as distillation and extraction. The problem is basically to determine whether a phase of given composition, temperature, and pressure will split into multiple phases. This problem is frequently formulated in terms of the tangent plane condition (Baker et al., 1982). Minima in the tangent plane distance are sought, usually by solving a system of nonlinear equations for the stationary points (Michelsen, 1982). If any of these yield a negative tangent plane distance, indicating that the tangent plane intersects (or lies above) the Gibbs energy of mixing surface, the phase is unstable. The difficulty lies in that, in general, given any arbitrary equation of state or activity coefficient model, most computational methods cannot find with complete certainty all the stationary points, and thus no guarantee of stability can be provided.
What is needed are robust techniques that can find all solutions to a system of nonlinear equations, and do so with certainty, or techniques for finding the global optimum of a function, and do so again with certainty. Recent advances in the field of interval mathematics make possible just such techniques.
Interval mathematics involves computation with intervals as opposed to real numbers. Interval Newton methods, when combined with a generalized bisection approach, provide the power to find with confidence all solutions of a system of nonlinear equations (Neumaier, 1990; Kearfott and Novoa, 1990), and to find with total reliability the global minimum of a nonlinear objective function (Hansen, 1992), provided only that upper and lower bounds are available for all variables. Efficient techniques for implementing interval Newton/generalized bisection are a relatively recent development, and thus such methods have not yet been widely applied. Schnepper and
1
2FOCAPD CONFERENCE PROCEEDINGS
Stadtherr (1990) suggested the use of this method for solving chemical process modeling problems, and recently described an implementation (Schnepper and Stadtherr,1994). The technique proved very efficient on a number of small to moderate size problems, including a vapor-liquid equilibrium problem with multiple roots, and was made even more efficient by taking advantage of the fact that the technique is amenable to parallel computing.
For the phase stability problem various approaches have been proposed recently. For example, Sun and Seider (1992,1993) apply a homotopy-continuation method,which will often find all stationary points, but may be initialization dependent and provides no theoretical guarantees that all solutions have been found. Also,McDonald and Floudas (1993,1994) show that for certain activity coefficient models, the problem can be reformulated to make it amenable to solution by global optimization techniques, which do guarantee the correct answer. However, in general there appears to be a need for an efficient general-purpose method that can perform phase stability calculations with complete reliability for any arbitrary equation of state or activity coefficient model.
In this paper, we demonstrate the use of interval methods for phase stability computations. These methods can be applied in connection with any equation of state or activity coefficient model, and when properly implemented are completely reliable.Phase Stability Analysis
The determination of phase stability is often done using tangent plane analysis (Baker et al., 1982;Michelsen, 1982). A phase at specified temperature,pressure, and feed composition z is unstable if the Gibbs energy of mixing versus composition surface m (x ) ?g M (x )/RT ever falls below a plane tangent to the surface at z . That is, if the tangent plane distance
D (x ) = m (x ) - m 0 -
∑i =1
n
(?m /?x i )0 (x i - z i )
(1)
is negative for any composition x , the phase is unstable.The subscript zero indicates evaluation at x = z , and n is the number of components. A common approach for determining if D is ever negative is to minimize D subject to the mole fractions summing to one. It is readily shown that the stationary points in this optimization problem can be found by solving the system of nonlinear equations:()?m ?x i - ()?m ?x i 0 - ()?m ?x n + ()?m
?x n 0
= 0, i = 1:n-1 1 -
∑i =1
n
x i = 0
(2)
The n × n system (2) has a trivial root at x = z , and frequently has multiple nontrivial roots as well. Thus conventional equation solving techniques such as Newton's method may fail by converging to the trivial root or give an incorrect answer to the phase stability problem by converging to a stationary point that is not the global minimum of D . This is aptly demonstrated by the experiments of Green et al. (1993), who show that the pattern of convergence from different initial guesses demonstrates a complex fractal-like response for even the simple liquid phase model used as an example below.
We demonstrate here the use of an interval Newton/generalized bisection method for solving the system (2). The method requires no initial guess, and will find with certainty all the stationary points of D .Interval Computations
A real interval , X , is defined by X = [a ,b ] ={x ∈ ? | a ≤ x ≤ b }, where a ,b ∈ ? and a ≤ b . A real interval vector X = (X i ) = (X 1,X 2,...,X n )T has n real interval components X i , and since it can be interpreted geometrically as an n -dimensional rectangle, is frequently referred to as a box . Several good introductions to computation with intervals are available, including recent monographs by Neumaier (1990) and Hansen (1992).
Of particular interest here are interval Newton methods. Consider the solution of the system of real nonlinear equations f (x ) = 0, where it is desired to find all solutions in an specified initial box X (0). The basic idea in interval Newton methods is to solve the linear interval equation system
F ′(X (k ))(N (k ) - x (k )) = -f (x (k ))
(3)
for the interval N (k ), where k is an iteration counter,
F ′(X (k )) is a suitable interval extension of the real Jacobian J(x ) of f (x ) over the current box X (k ), and x (k ) is a point in the interior of X (k ), usually taken to be the midpoint. It can be shown (Moore, 1966) that any root x * ∈ X (k ) of f (x ) is also contained in N (k ). This suggests the iteration
X (k +1) = X (k ) ∩ N (k ).
(4)
There are various interval Newton methods, which differ in how they determine N (k ) from equation (3) and thus in the tightness with which N (k ) encloses the solution set of (3).
While the iteration scheme given by equations (3) and (4) can be used to tightly enclose a solution (e.g., Shacham and Kehat, 1973), what is of most significance here is the power of (3) to provide an existence and uniqueness test.For several techniques for finding N (k ) from (3), it can be proven (Neumaier, 1990) that if N (k ) ? X (k ), then there is a unique zero of f (x ) in X (k ), and furthermore that Newton's method with real arithmetic will converge to that
Robust Phase Stability Analysis3
solution starting from any point in X(k). Thus, if N(k) is determined using one of these techniques, the computation can be used as a root inclusion test for any interval X(k). If X(k)∩N(k) = ?, then there is no root in X(k); if N(k) ?X(k), then there is exactly one root and Newton's method with real arithmetic will find it; otherwise, no conclusion can be drawn. In the last case, one could then repeat the root inclusion test on the next interval Newton iterate X(k+1), assuming it is sufficiently smaller than X(k), or one could bisect X(k+1) and repeat the root inclusion test on the resulting intervals. This is the basic idea of interval Newton/generalized bisection methods. If f(x) = 0 has a finite number of real solutions in the specified initial box, a properly implemented interval Newton/generalized bisection method can find with mathematical certainty any and all such solutions to a specified tolerance, or can determine with mathematical certainty that there are no solutions in the given box (Kearfott, 1987,1990).
The technique used here for computing N(k) is the preconditioned Gauss-Seidel-like technique developed by Hansen and Sengupta (1981) and Hansen and Greenburg (1983). Neumaier (1985,1990) has proven the existence and uniqueness test for this method of determining N(k).
Since all variables are mole fractions, the initial box X(0) = [0,1] is suitable. In practice the initial lower bound is set to an arbitrarily small positive number ε(10-10 was used) to avoid taking the logarithm of zero in subsequent calculations. Since it is known (Michelson, 1982) that the roots sought lie in the interior of [0,1], this can be done without the loss of reliability providing a sufficiently small value of ε is used.
Our implementation of the interval Newton/generalized bisection method for the phase stability problem is based on appropriately modified routines from the packages INTBIS (Kearfott and Novoa, 1990) and INTLIB (Kearfott et al., 1993). To demonstrate the potential of the technique we use two models of liquid-phase systems, one a simple model of excess Gibbs energy, the other the NRTL model. Simple Model
This model is used by Green et al. (1993) to demonstrate the difficulties in using Newton's method for phase stability problems. It is a three component system with the excess Gibbs energy given by ?g E/RT = 3x1x2. This may be regarded as a special case of the two-suffix Margules equation. We solved the system (2) for this model for the three different feed compositions used by Green et al. The roots found and the corresponding value of D are shown in Table 1, and are the same as those reported by Green et al. The results indicate phase stability for only the third feed composition. Performance data, including the number of root inclusion tests required, the depth reached in the binary tree generated in the bisection, and the computation time on an HP 9000/735 workstation, is given in Table 2.NRTL Model
This example is used by McDonald and Floudas (1993,1994) to demonstrate a global optimization technique for solving the phase stability problem when the NRTL model for excess Gibbs energy is used. There are two components: n-butyl acetate (1) and water (2). The parameters in the NRTL model for this system are G12 = 0.30794, G21 = 0.15904, τ12 = 3.00498, and τ21 = 4.69071. The roots of system (2) found for five feed compositions are shown in Table 3, and corresponding performance data is shown in Table 4. There are as many as five roots for some feed compositions and phase instability is indicated for all five cases.
Table 1. Simple Model: Roots Found
Feed: (z1,z2,z3)Roots: (x1,x2,x3)D
(0.45,0.45,0.10)(0.8509,0.0857,0.0634)-0.0671
(0.0857,0.8509,0.0634)-0.0671
(0.45,0.45,0.10)0.0
(0.60,0.18,0.22)(0.1185,0.6842,0.1973)-0.0284
(0.4675,0.2923,0.2401)0.00152
(0.60,0.18,0.22)0.0
(0.90,0.06,0.04)(0.3447,0.5820,0.0733)0.1664
(0.1193,0.8277,0.0530)0.1479
(0.90,0.06,0.04)0.0
Table 2. Simple Model: Performance Feed: (z1,z2,z3)
Root
Inclusion
Tests
Level
Reached in
Binary Tree
CPU Time:
HP 9000/735
(sec.)
(0.45,0.45,0.10)111140.04
(0.60,0.18,0.22)119150.03
(0.90,0.06,0.04)117140.03
The CPU times indicate performance that is significantly faster than a GAMS implementation of the model-specific procedure of McDonald and Floudas (1993), and that is comparable to a more recent C implementation of their technique (McDonald and Floudas, 1994). For example, for the last feed composition (z1 = 0.65) McDonald and Floudas (1993,1994) report CPU times on a HP 9000/730 of 10.74 s for the GAMS implementation of their technique and 0.11 s for the C implementation of the method, versus 0.06 s on an HP 9000/735 for the interval Newton/generalized bisection method used here. (Note that, based on the well-known LINPACK-100 benchmark, the Model 735 used here is about 70% faster than the Model 730 used by McDonald and Floudas.) In addition to being very efficient, the technique used here has the advantage of being model-independent; it can be easily used to solve the phase stability problem for other models of excess Gibbs energy, and can also be used in connection with equation of state models.
4FOCAPD CONFERENCE PROCEEDINGS Concluding Remarks
The interval computation method demonstrated here
provides an efficient and completely reliable means for
solving the phase stability problem. The method is model-
independent, straightforward to use, and requires no
problem reformulation. Though we have treated the
problem here as one of solving a system of nonlinear
equations for all its roots, the method can readily be
extended to try to take advantage of the fact that this can be
considered a global optimization problem, and that to
determine phase instability all that really needs to be found
is an interval value of D with a negative upper bound.
Since this requires the additional work of doing an interval
function evaluation of D for each box tested, it is not clear
that this will always result in any overall savings, though
for problems in which the phase is unstable, savings are
likely and potentially substantial.
Table 3. NRTL Model: Roots Found
Feed: (z1,z2)Roots: (x1,x2)D
(0.50,0.50)(0.1602,0.8398)0.0278
(0.00421,0.99579)-0.0325
(0.50,0.50)0.0
(0.10,0.90)(0.96346,0.03654)-0.2142
(0.00599,0.99401)-0.0291
(0.10,0.90)0.0
(0.20,0.80)(0.3922,0.6078)-0.00607
(0.00379,0.99621)-0.0743
(0.20,0.80)0.0
(0.15,0.85)(0.5423,0.4577)-0.0388
(0.8549,0.1451)-0.0260
(0.9216,0.0784)-0.0267
(0.00438,0.99562)-0.0557
(0.15,0.85)0.0
(0.65,0.35)(0.7593,0.2407)0.00063
(0.9413,0.0587)-0.00671
(0.1346,0.8654)0.0632
(0.00471,0.99529)0.0150
(0.65,0.35)0.0
Table 4. NRTL Model: Performance
Feed: (z1,z2)
Root
Inclusion
Tests
Level
Reached in
Binary Tree
CPU Time:
HP 9000/735
(sec.)
(0.50,0.50)206130.05 (0.10,0.90)116110.02 (0.20,0.80)204140.04 (0.15,0.85)245150.06 (0.65,0.35)268170.06Acknowledgement
This work was funded in part by the National Science Foundation under grant DDM-9024946.
References
Baker, L. E., A. C. Pierce, and K. D. Luks (1982). Gibbs energy analysis of phase equilibria. Soc. Petrol.
Engrs. J., 22, 731-742.
Green, K. A., S. Zhang, and K. D. Luks (1993). The fractal response of robust solution techniques to the
stationary point problem. Fluid Phase Equilib., 84,
49-78.
Hansen, E. R. (1992). Global Optimization Using Interval Analysis. M. Dekkar, New York.
Hansen, E. R. and R. I. Greenburg (1983). An interval Newton method. Appl. Math. Comput., 12, 89-98. Hansen, E. R. and S. Sengupta (1981). Bounding solutions of systems of equations using interval analysis. BIT,
21, 203-211.
Kearfott, R. B. (1990). Interval arithmetic techniques in the computational solution of nonlinear systems of
equations: Introduction, examples, and comparisons.
Lectures in Applied Mathematics, 26, 337-357. Kearfott, R. B. (1987). Abstract generalized bisection and a cost bound. Math. Comput., 49, 187-202. Kearfott, R. B. and M. Novoa III (1990). INTBIS, A portable interval Newton/bisection package. ACM Trans.
Math. Softw., 16, 152-157.
Kearfott, R. B., M. Dawande, K. Du, and Ch. Hu (1993).
INTLIB: A portable FORTRAN 77 interval standard
function library. Submitted for publication. McDonald, C. and C. A. Floudas (1993). Global optimization for the phase stability problem. Presented at AIChE
Annual Meeting, St. Louis, November. McDonald, C. and C. A. Floudas (1994). Global optimization for the phase stability problem. AIChE J., to appear. Michelsen, M. L. (1982). The isothermal flash problem. Part I: stability. Fluid Phase Equilib., 9, 1-19.
Moore, R. E. (1966). Interval Analysis. Prentice-Hall, Engelwood Cliffs, NJ.
Neumaier, A. (1990). Interval Methods for Systems of Equations. Cambridge University Press, Cambridge. Neumaier, A. (1985). Interval iteration for zeros of systems of equations. BIT, 25, 256-273.
Shacham, M. and E. Kehat (1973). Converging interval methods for the iterative solution of a non-linear
equation. Chem. Eng. Sci., 28, 2187-2193. Schnepper, C. A. and M. A. Stadtherr (1990). On using parallel processing techniques in chemical process
design. Presented at AIChE Annual Meeting,
Chicago, November.
Schnepper, C. A. and M. A. Stadtherr (1994). Robust process simulation using interval methods. Submitted for
publication.
Sun, A. C. and W. D. Seider (1993). Global minimization of the Gibbs free energy. Presented at AIChE Annual
Meeting, St. Louis, November.
Sun, A. C. and W. D. Seider (1992). Homotopy-continuation algorithm for finding globally stable phase
equilibria. Presented at AIChE Annual Meeting,
Miami Beach, November.