Ind. Eng. Chem. Res. 1995,34, 958-966
958
Instability of Successive Substitution Robert A Heidemam* Department of Chemical and Petroleum Engineering, University of Calgary, 2500 University Drive, Calgary, Alberta, Canada T2N 1N4
Michael L. Michelsen Znstitut for Kemiteknik, Bygning 229, Technical University of Denmark, 2800 Lyngby, Denmark
The method of successive substitution for phase equilibrium calculations can only converge to a minimum of the mixture Gibbs free energy, and it is commonly assumed that, given proper initial estimates, convergence will always occur. Unfortunately, although the exemptions are rare, the assumption is incorrect. Divergence, regardless of the quality of the initial estimate, can occur in mixtures that exhibit strong negative deviations from ideal solution behavior, such as polymer solutions and gas hydrates. Unless special precautions are taken, algorithms that implement successive substitution as part of the solution procedure are likely to fail on such mixtures. In this paper we present a n analysis of the thermodynamic conditions that may lead to divergence of successive substitution and present some numerical examples. In addition we briefly discuss some methods capable of solving the appropriate equilibrium equations under conditions where successive substitution fails.
Introduction
yi = YilZy,
In the “successivesubstitution” algorithms for phase equilibrium, the objective is to determine the compositions of one or more phases (typically, xi and yi) using the equilibrium ratios defined by
Ki sz yilxi Two familiar procedures are flash calculations and phase stability calculations, both at fixed temperature and pressure. The flash problem is to find p, the fraction of a mixture with initial compositionzi that is in they phase a t equilibrium. That fraction is updated iteratively by first solving
j
The new mole fractions in the y phase are used to update the Ki values and the process is repeated. The flash calculation process minimizes the Gibbs free energy of the two-phase mixture. The stability computation is also derived from a minimization problem. It locates the point at which the Gibbs free energy surface of they phase is at a minimum distance from the plane drawn tangent to the Gibbs free energy surface of the x phase at the given x composition. As was shown by Michelsen (1982a), the (dimensionless) minimum tangent plane distance is given by the converged value of
The Ki values are updated between iterations in a process that could be written for p, assuming that the Ki are constant. Then the phase mole fractions are evaluated from
(7) in flash calculations or
(3) The thermodynamic models for the two phases are then used to update the Ki values and the process is repeated. The process is known as the Rachford-Rice procedure. In the phase stability computation, the mixture is assumed t o be all in the x phase with known composition. The composition of a y phase is obtained by first calculating
Yi = Kzci
(4)
and then normalizing so that mole fractions sum to unity; i.e.,
* Author
in stability calculations. In either case, convergence is obtained when successiveKi are sufficiently close to each other. The stability calculation is closely related to “saturation point” calculations in which a variable (such as the temperature or pressure) is adjusted in order to satisfy
In general, the process of solving simultaneous equations through successive substitution can be represented by the expression
to whom correspondence should be addressed.
Telephone: (403) 220-8755. FAX: (403) 284-4852. E-mail:
[email protected].
which is presumed to have a solution t* that satisfies
~ a a ~ - ~ a a ~ ~ ~ ~ 0~ 1995 ~ ~American ~ ~ - oChemical ~ ~ ~Society ~ o ~ . o o ~ o
Ind. Eng. Chem. Res., Vol. 34, No. 3,1995 959
t* = f(t*)
(11)
Behavior of the process depends on the eigenvalues of the matrix
E = {afii/at,},,,,
(12)
Convergence to t* is only possible if all the eigenvalues of E are less than unity in absolute value. This property was used by Michelsen (1982a,b) to demonstrate that the successive substitution flash calculation procedure and the tangent plane stability procedure can only converge to minima. Michelsen (1982b) also pointed out that convergence of the Rachford-Rice procedure must necessarily become very slow near a critical point because two of the eigenvalues of the relevant E matrix tend toward +1 at such points. It has been generally presumed that the two successive substitution procedures outlined above will be well behaved, even if slow t o converge in some regions. The observation we wish to report in this paper, however, is that both procedures can be nonconvergent for some systems, far from any points where either of the equilibrium phases can become unstable or approach criticality. The nonconvergent behavior is associated with large negative deviations from the ideal solution model and results when one or more eigenvalues of the E matrix are less than -1. The consequence will be an oscillatory sequence of phase compositions with an error that does not decrease.
Analysis The stability computation of eqs 4, 5, 6, and 8 is the easier t o analyze. In this case, convergence behavior is dictated by the eigenvalues of the matrix
In general, all N eigenvalues of the two matrices are real numbers. An analysis is presented in the Appendix that demonstrates the following properties of the eigenvalues in eq 15 that control the convergence behavior of successive substitution: (a) So long as 0 /3 1,
A&
I A I Amax
(16)
where 1- and Am= represent the smallest and the largest of the nonzero stability eigenvalues, (Ai,A:, ...,
A;, A;, A;, ..., A&).
(b) In the limit, /.? = 0 or p = 1, N - 1 of the flash calculation eigenvalues are identical to the nonzero eigenvalues in the missing phase (i.e., the convergence behavior is dictated by the “small” phase). The last of the eigenvalues depends on the properties of the “large” phase model but is bounded by the eigenvalues in the “small” phase. (c) For a binary mixture with an ideal y phase, there is a single nonzero eigenvalue for all values of /3 and this eigenvalue does not depend on p. The eigenvalues of the E matrix for a binary mixture can be written in a more familiar form. One, of course, is zero. The second can be found through some routine operations using the matrix definition of eq 13 to be
where x2 is regarded as a dependent variable during the differentiation. In terms of the fugacity of the component in the phase
(13) where r$ is a fugacity coefficient or an activity coefficient, depending on the nature of the model used for the missing phase. The eigenvalues of E can be found in a very direct way for some typical binary activity coefficient models. The analysis of the flash calculation procedure involves matrices like E in eq 13 for the two equilibrium phases, i.e., E“ and EY for the x phase and they phase, with eigenvalues (A;, Ai, ..., A;) and (AS;, A;, ..., A&). The analysis requires treating each of these as a product of matrices
E“= -XlW and Ey=-yMy (14) where X and Y are diagonal matrices with elements Xii = xi and Yii = yi. The eigenvalues characterizing the rate of convergence are found from (Michelsen, 1982b):
-[Pcx-’- eeT)+ (1- P)(Y-’ - e e T > l - ’ w+~ (1- P)W]U= AU (15)
In eq 15, the auxiliary vector e is defined as eT= (1,1, ..., 1)and eeTis a matrix with all elements equal to 1. As a consequence of the Gibbs-Duhem equation, both E“ and Ev have a zero eigenvalue, ; 1 = AS; = 0, corresponding to eigenvectors (XI, ..., X N ) and 011, ..., y d T , respectively. For an ideal solution model, matrix M has all zero elements and all the eigenvalues are zero.
When the derivative on the left-hand side of eq 18 is negative, the phase model is “unstable”in the sense that a lower Gibbs free energy can be reached by dividing the phase into two with different compositions. For this reason, A 2 < 1 provides one limit of interest. The second limit is 1 2 > -1, since the successive substitution procedures will become unstable outside this bound. This second limit provides a practical and simple means for determining in advance whether a given system might be computable via successive substitution. If a plot of In f i vs In x1 has a slope greater than +2, successive substitution will not work. Note that the slope of such a plot for an ideal solution would be f l . Another way to state our result is that successive substitution will not work if the nonideal contribution to the phase model exceeds the ideal solution contribution.
Activity Coefficient Examples Some widely used liquid phase activity coefficient models show negative deviations from ideality. For binary systems described by these models, we have obtained the following limits beyond which divergence of successive substitution algorithms will be observed, a t least in some composition ranges. Flory-Huggins Model. For the “athermal”FloryHuggins model, the E matrix will have a real eigenvalue less than -1 in some composition range if the ratio of the molar volumes of the two species is greater than 3
960 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995
+ 6. (This model is incapable of showing liquid instability.) Wilson Model. The Wilson activity coefficient model contains two parameters, 1\12 and 1\21. When the two values are equal, 1\21 = 1112, convergence difficulties will occur when the values exceed 1 &. If both are greater than 1, the limiting sum is around 4.8. The Flory-Huggins limit of 3 & applies for the larger parameter when the product of the parameters is 1. (Like the "athermal" Flory-Huggins model, the Wilson equation cannot show liquid phase instability.) Van Laar Model. A limiting curve for the two A parameters in the van Laar model is given by
+
+
f l ~=) -13.5~/[(2- x)(l + x ) ~ ] A,, =flx),Azl =fll-x)
for 0
< x .c
0.3
sa
0.2
c
0 .c
u
e
LL
v)
0
I
0.1
(19)
1
The limits are (-3.375,0), through (-2,-2), to (0,-3.375). (Note that the same equations, but with positive values for f , define the parameter limits at which liquid-liquid separations occur.) We have located a few binary systems in Vol. I, Part 5, of the Dechema Series (Gmehling et al., 1982) where the correlated van Laar parameters will result in nonconvergent behavior of successive substitution. Examples are; formic acid-N,N-dimethylformamide (p 211, acetic acid-N,N-dimethyl-acetamide (p 1121, and pyridineacetic acid (p 118). It is difficult t o generalize about the behavior of ternary systems, but we discuss in the following a specific example of one system where successive substitution cannot be used. Kang and Sandler (1987) have correlated the liquidliquid equilibrium behavior of a ternary mixture of water with two water-soluble polymers, dextran (DX 17) and poly(ethy1ene glycol) (PEG 6000). The model used was the Flory-Huggins model, including the enthalpic part with binary interaction parameters. The parameters for the model as used by us are as follows: specific volumes (cm3/g),water, 1.0, DX 17, 0.626, PEG 6000, 0.832; molar masses (g/mol),water, 18.0, DX 17,23 000, PEG 6000,6000; interaction parameters, water-DX 17, 0.535, water-PEG 6000, 0.50924, DX 17-PEG 6000, 0.046 83. The behavior of this mixture near its plait point is shown in Figure 1. In Figure 1,the coexistence curve and one tie line are shown as solid lines. We have also presented lines along which one of the eigenvalues of matrix E for the mixture has a constant value. Two sets of lines of constant eigenvalue are shown, i.e., dashed lines for the positive values and dotted lines for the negative values. One line of interest is for il= +l. As was discussed above, this line defines the limit of stability for the homogeneous phase (the so-called spinodal curve) which lies inside the two-phase region. It is, however, the values of the negative eigenvalues that affect the behavior of the successive substitution algorithms. It will be seen from the Figure that at every composition along the coexistence curve the E matrix has a negative eigenvalue less than -1. The consequence of this behavior is that successive substitution algorithms cannot be used to compute the phase behavior of this system. We have calculated the eigenvalues for the flash calculation procedure if applied to mixtures along the tie line shown in Figure 1. The end points of the tie line are at W I = 0.199 769, w 2 = 0.015 230 and w 1 = 0.008 75, w 2 = 0.126 931 where w 1 and w 2 are mass
0.0 0.0
0.1 0.2 Mass Fraction Dextron
0.3
Figure 1. Phase diagram and stability eigenvalues in an aqueous two-polymer mixture of Kang and Sandler. Flory-Huggins Model. 0.6 I
0.2 I 0.0
-10
1
0.0
I
I
I
I
1
I
I
I
I
I
0.2
0.4
0.6
0.8
1 .o
I
I
1
1
0.2 0.4 0.6 0.8 Fraction (Mass Basis) in Phase I1
1 .o
Figure 2. Flash calculation eigenvalues along a tie line in an aqueous two-polymer mixture of Kang and Sandler. FloryHuggins Model.
fractions of dextran and PEG, respectively. The eigenvalues for the stability calculation at the end-point compositions are (-8.7190, 0.4953, 0.0) and (-2.8866, 0.3708, 0.0). Eigenvalues for the flash calculation procedure are shown in Figure 2 where they are plotted against the fraction (in mass units) in the PEG-rich phase. It can be seen from the figure that, whatever the position on the tie line of the feed composition,there is always an eigenvalue less than -1. The variation of the eigenvalue is between the values for the "missing"
Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 961 1.5
I
I
I
I
0.0
0.2
0.4
0.6
0.8
1
I
I
I
I
250
Pyridine-Acetic
Acid
van Loar Liquid Phase A,, =-3.1041
1 .o
n i
Ideal Vapor Phase Vapor Pressures; Pyridine, 242.5 mm Hg Acetic Acid, 207.3 mm Hg
I
E
E
‘
100 0.0
I
I
I
I
0.2
0.4
0.6
0.8
1 .o
Mole Fraction Pyridine
Figure 3. Phase diagram and (a In yla In x ) for pyridine-acetic acid. Van L a a r model.
phase a t the end points. The two positive eigenvalues are bounded at one end by the corresponding eigenvalue for the “missing“phase a t that end of the tie line. There is no zero eigenvalue for the flash calculation process, as is indicated by the analysis in the Appendix. Kang and Sandler (1987) comment on some of the difficulties they experienced in solving the equilibrium equations for this system, an operation required repeatedly in order t o correlate the equilibrium data.
Mixed Model Example One familiar approach to modeling of vapor-liquid equilibrium behavior is to couple an excess free energy model for the liquid phase with an ideal gas model for the vapor. As was pointed out above and as shown in the Appendix, for a binary system the eigenvalues that determine the behavior of flash calculations are the same as the eigenvalues of the matrix in eq 13 for the liquid phase alone. The implications are that the flash calculation will not converge, whatever the feed composition, whenever an eigenvalue for the equilibrium liquid phase is less than -1. The example considered here uses the van Laar coefficients from the Dechema correlation for pyridine (1)-acetic acid (2); i.e., A12 = -3.1041 and A21 = -2.4885. The vapor pressures used are 242.5 and 207.3 mmHg for pyridine and acetic acid, respectively (values that apply at around 80 “C). The bubble-point pressure calculation for this model is noniterative. The calculated equilibrium curves and a plot of (a In y l a In x ) against the pyridine liquid mole fraction are presented in Figure 3. When (a In y l a In x ) exceeds f l at the
composition of the equilibrium liquid phase, oscillatory nonconvergence can be expected in the successive substitution flash calculation procedure. From the curves in Figure 3, it can be concluded that successive substitution will not work on this system whenever the mole fraction of pyridine in the liquid lies in the interval 0.18 Ix I 0.70, or equivalently, at pressures below approximately 158 or 154 mmHg (depending on which side of the azeotrope the feed composition falls). A mixture with 0.85 mole fraction pyridine was flashed at a series of pressures crossing the two-phase region in Figure 3. If for the given K values there was a solution of eq 2 with 0 < p < 1,eq 7 was used to update the K values for the next iteration. Otherwise, if the mixture was found to be all in one phase with either p = 0 or p = 1,then eq 8 was used. The following results were observed: (a)At any pressure above the bubble-point curve, the calculation is noniterative since the Kvalue is independent of the vapor composition. This is a consequence of using the ideal gas model for the vapor phase. (b) At 200 mmHg, convergence to five figures in the mole fractions is obtained in 18 iterations. The vapor fraction is 0.027 40 and the mole fractions of pyridine in the liquid and vapor phases are 0.846 45 and 0.976 11, respectively. (c) At lower pressures, the mole fraction of pyridine in the liquid moves toward the region where (a In y l a In x ) =1, hence toward slower convergence with increasingly oscillatory behavior. At 160 mmHg, the vapor fraction is 0.670 29 and the two pyridine mole fractions are 0.719 43 and 0.914 23. Convergence is reached in 293 iterations. (d) At still lower pressures, the procedure cannot converge. The behavior seen is variable and complex. At 155 mmHg, the process heads toward a “flickering” between two alternating states, one with a vapor fraction of 0.764 and a liquid with pyridine mole fraction 0.683 and a second state with values of 0.709 and 0.724 for the phase fraction and the pyridine mole fraction. Over a range of pressures, the “flickering”is between a state with the whole mixture in the liquid and a state with a finite vapor phase fraction. This sort of behavior continues to pressures below 120 mmHg. Calculations a t 110 and 115 mmHg produce alternating iterations with the feed either all in the liquid or all in the vapor. (e) At 105 mmHg, the “flickering”has evolved into a “limit cycle”with four different states that are repeated in sequence indefinitely. In one of the states, the mixture is all liquid. In the other three, the mixture is all vapor with potential liquid phases computed of differing compositions. (0At 100 mmHg and lower, the computations become chaotic with no repeating pattern of phases. Analysis of the very complex behavior of these equations when they are nonconvergent is perhaps of some theoretical interest. However, in this account we merely wanted to point out that a naive use of very familiar computational methods on this rather simple problem can produce a great variety of unwanted and troublesome results. The source of the trouble can be seen very clearly in the top section of Figure 3 through the presence of a composition range where (a In y l a In x ) is greater than f l .
Equation of State Examples The system HC1-water shows large negative deviations and has been correlated by Stryjek and Vera
962 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 7 ,
5
1
I
Hydrogen Chloride Water at 343.15 and 0.05 MPa
3
\
-1
‘
0.0
d
1
0.2
I
I
0.00 0.0
i I
I
0.1 Mole Fraction HCI
0.10 I
K
I I
0.1 Mole Fraction
Figure 4. Phase diagram and Stryjek-Vera model.
(a
0.2 HCI
In #/a In r) for HC1-water.
(1986) using a relatively simple cubic equation of state. The correlation involves a composition-dependent interaction parameter in the mixing rule for a in the equation of state. In the composition range where HCl is dilute, this interaction parameter reaches rather extreme values (around -1.5). Figure 4 shows the values of (a In $/a In x ) for hydrogen chloride in water liquid at 343.15 K and a pressure of 0.05 MPa that are calculated using the Stryjek and Vera (1986) model. A pressure of 0.05 MPa was chosen for these calculations because it is above the bubble-point pressure of liquid mixtures below the azeotropic composition at the given temperature. However, the results will not be affected in any significant way by the pressure. Bubble point pressures are also shown in Figure 3. As the figure shows, ( 8 In $/a In x ) is greater than +1.0 whenever the liquid HC1 mole fraction is greater than approximately 0.014. The implication is that successive substitution flash calculations will be nonconvergent whenever the equilibrium mole fraction in the liquid phase exceeds about 0.014. (The vapor phase is near-ideal at the low pressures involved, so it is the liquid nonideality that determines the computational behavior.) The calculated bubble-point pressure of a liquid with 0.005 mole fraction HC1 in water a t 343.15 K is 0.031 MPa. The equilibrium vapor is practically free of HC1, and the mole fraction of HC1 in the liquid phase will increase above 0.014 if the mixture is flashed at lower pressures. The following behavior of the successive substitution flash procedure has been observed. (a) P = 0.030 60 MPa, T = 343.15 K. The procedure converges in an oscillatory manner to a vapor fraction
of 0.6118 and with liquid and vapor HC1 mole fractions of 0.012 86 and 0.00001, respectively. Convergence to 4 figures in the mole fractions is obtained in about 100 iterations. The slow oscillatory convergence is due to the negative eigenvalue very close to -1. (b)P = 0.030 55 MPa, T = 343.15 K. The equilibrium liquid mole fraction is 0.013 69, closer to the limit. As a consequence, the iteration count to convergence increases to 270. At P = 0.3050 MPa, the iteration count is much more than the limit set of 300 and convergence was not achieved. (c)P = 0.030 45 MPa, T = 343.15 K. The equilibrium liquid mole fraction exceeds the limit and convergence is no longer possible. The process “gets stuck” and “flickers”between two states in alternate iterations. One shows a vapor fraction of 0.7674 with a HC1 mole fraction of 0.021 44 in the liquid phase. The alternating iteration shows a vapor fraction of 0.3378 with a liquid phase HC1 mole fraction of 0.007 54. (d)P = 0.030 00 MPa, T = 343.15 K. The behavior is similar to the nonconvergence observed in (c) except that one of the two “flickering”states is all liquid with the feed HC1 mole fraction and the alternating iteration shows a vapor fraction of 0.8653 with a liquid HC1 mole fraction of 0.036 96. The behavior observed in these calculations with HCl and water is of the same kind as was noted with the pyridine-acetic acid system and arises from the same source, i.e., a large positive value for (a In $/a In x ) or (a In y/a In x ) . It might be noted that large negative values of ( 8 In $/a In x ) or ( 8 In y/a In x ) are responsible for demixing in EOS or excess free energy models. In equations of state, demixing is induced through use of positive interaction parameters. It is the large negative interaction parameter in the HC1-water model that results in the poor computational behavior. We have carried out a numerical experiment with the nearly ideal mixture of n-pentane and isopentane. The mixture has been made nonideal by introducing an interaction parameter into the quadratic mixing rule for a in the Soave-Redlich-Kwong equation of state. A value of k,j. = 0.17 is sufficient to produce liquid-liquid phase separations in this system. A value of K,j. = -0.17 is sufficient to destabilize the successive substitution algorithms. A number of equations of state have been proposed for polymer-solvent systems which, as can be seen from the Flory-Huggins model, have large negative deviations from ideality in the liquid phase. From the discussion above, it could be expected that such models may not be amenable to solution through the successive substitution algorithms, a t least in some regions of temperature, pressure, and composition. The recent paper by Chen et al. (1993) describes a failure of the successive substitution flash calculation procedure for ethylene-polyethylene mixtures using the SAFT equation. Convergence was obtained near the critical point of the binary mixture examined (and therefore, according t o the analysis by Michelsen (1982b1, where the E matrix has two eigenvalues near +l).We suggest that, away from the critical points, where the negative deviations in the polyethylene-rich phase become dominant, successive substitution would not converge and that this is the source of the difficulty Chen et al. (1993) encountered.
Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 963 One relatively simple procedure to overcome oscillatory nonconvergence of successive substitution in general is to apply a “damping factor” to the computed change in the K values between iterations. The modified process could take the form
Methane Hydrate 50 atm., 280 K
s 10
CI
rg
0.00 4 ,
I
I
I
I
0.05
0.10
0.15
1
I
/
I
I
/I I
Methane Hydrate
in flash calculations or
in stability calculations, where m is the damping factor. It is easy to show that an eigenvalue of the damped process, I * , is related to the corresponding eigenvalue of the original process, I , by
31” = 1 - m ( l - 31)
I
0.00
I
0.05 0.10 Mole Fraction Methane in the Hydrate
0.15
Figure 5. Fugacity of methane in t h e hydrate and (3 Inf l 3 InX I .
Gas Hydrate Example The van der Waals-Platteeuw model for gas hydrates is still another potential source of difficulties with successive substitution equilibrium calculations. Figure 5 shows the fugacity and (a In fla In x) for methane in hydrates of structure I and structure I1 at 280 K and 50 atm. Michelsen’s (1991) procedure for calculating the fugacity and the model parameters of Munck et al. (1988) have been used in generating the figure. The three-phase equilibrium pressure (water-methanemethane hydrate) at 280 K is 52.5 atm, and the stable hydrate structure is structure I with a 0.1436 mole fraction of methane. The full lattice can accommodate a methane mole fraction of 0.148. The derivative (a In f l a In x ) is greater than 2 (hence, (a In +/a In x ) is greater than 1) a t mole fractions in excess of about 0.07 in both structures. Given this behavior of the model, successive substitution algorithms for checking the stability of a methane-rich vapor with respect t o formation of the hydrate solid or for an isothermal flash calculation at conditions where the solid is present would show oscillatory divergence.
Remedies In some cases, problem specific procedures have been developed for dealing with computational difficulties that may have been encountered. A n instructive example is the approach used by Bishnoi et al. (1989) for dealing with gas hydrates in multi-phase flash calculations. The gas components are present in a fluid “referencephase” as well as in the hydrate. They equate the fugacity of a gas component in the hydrate phase with its fugacity in the reference phase. The van der Waals-Platteeuw model is employed only to calculate the fugacity of the water component in the hydrate and the mole fractions. This procedure eliminates the poor behavior of the successive substitution algorithm that results from calculating fugacities from the mole fractions.
(22)
The eigenvalues that are less than -1 can all be brought inside the unit circle (or even made positive, thus ensuring monotonic Convergence) if m is made small enough. The problem with introducing a high degree of damping (i.e., a small m value) is that all positive eigenvalues increase in magnitude and are pushed toward +1,thereby increasing the number of iterations to achieve convergence. It is instructive to consider the effect of “damping” on calculations in the pyridine-acetic acid example. To be specific, suppose that a damping factor of m = 0.5 is implemented. Then, according to eq 22, the zero eigenvalue is moved to I* = 0.5 and the negative eigenvalue (which is numerically equal to -(a In y l a In x ) shown in Figure 4) is shifted so that it takes on values ranging from about -0.2 in the midcomposition range (where (a In y l a In x ) = 1.4) t o 0.50in the dilute ranges (where (a In yla In x ) = 0). The convergence behavior is dominated by the larger eigenvalue, I* = 0.5, and the error should decrease roughly by a factor of 0.5 in each iteration as the final result is approached, regardless of the composition of the liquid phase. (Note that the - In $k-1’)211/2.) To “error” is defined as [Xln reduce the error to should require on the order of 20 iterations, in flash calculations and stability calculations, throughout the pressure-composition region shown in Figure 4. Computations confirm this expectation. Liquid-liquid equilibrium computations on the FloryHuggins mixture shown in Figure 1 would require a smaller “damping factor”. To reduce the magnitude of the negative eigenvalue shown in Figure 3 to within -1 requires a damping factor m < 241 - A) which ranges from 0.205 to 0.515 across the tie line. To assure convergence in a reasonable number of iterations, the value must be smaller. A value m = 0.15 permits convergent behavior at any composition on this tie line or between this tie line and the plait point. Convergence is dominated by the larger of the two positive eigenvalues. This value is near I = 0.5 across the tie line and is increased to I* = 1 - 0.15(1 - 0.5)= 0.925 by a damping factor of 0.15. Therefore, error reduction between iterations would be on the order of 0.925. Computations confirm the picture of convergence behavior deduced from the eigenvalues of the process. Starting with Kvalues of 100,0.01, and 1.0 for dextran, PEG, and water, respectively, convergence to in the K values occurs in 170 iterations at all points on the tie line when m = 0.15. When m = 0.20, convergence is obtained at phase fractions between 0.1 and 0.9 with 125 iterations.
gk)
964 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995
The successive substitution process does provide a direction of change in the K values that is toward a reduced Gibbs free energy. It is therefore possible to find an optimal value of the damping factor in the sense that the Gibbs free energy is minimized along the line of motion. Convergence could always be achieved by monitoring the change in the Gibbs free energy and reducing the step length, if necessary. The damped successive substitution process can also be accelerated using, for instance, the Crowe and Nishio (1975) general dominant eigenvalue method (GDEM); however the process might be erratic. The first-order GDEM appears better behaved on the tie line examined in the Flory-Huggins example than the second-order method, perhaps because numerical precision limitations interfere with calculating two eigenvalues. With a damping factor of m = 0.15, convergence is obtained everywhere on the tie line in around 80 iterations using either the first-order GDEM or the second-order GDEM with acceleration every fifth iteration. The iteration count can be reduced below 60 at most points by accelerating every third iteration, rather than every fifth iteration. Some combinations of feed location, damping factor, and acceleration frequency permit convergence within 40 iterations. The Mehra et al. (1983) acceleration scheme, if applied periodically rather than a t every iteration, has the same effect as the first-order GDEM procedure. Since the flash calculation leads to minimization of the Gibbs free energy and the stability calculation leads to a minimum of the Gibbs tangent plane distance, a number of approaches for minimizing functions can be adapted for use in these computations. Ammar and Renon (1987) describe the use of several optimization methods for phase split computations, including the BFGS algorithm, conjugate gradient methods, and steepest descent procedures. The BFGS algorithm is one general-purpose procedure which can be implemented without the need to evaluate derivatives of the fugacity coefficients or activity coefficients. Full Newton-Raphson iteration can be employed effectively whenever a satisfactory initial guess is available. Chen et al. (1993) have described an approach to Newton iteration in flash calculations with polymer-solvent systems described by an equation of state. Michelsen (1982b), Mehra et al. (19821, Nghiem et al. (19831, and Ammar and Renon (1987) all advocate using successive substitution to reach a region in which Newton iteration will converge. Apparently, an appropriately damped successive substitution procedure could be used in this way, even for systems where the undamped procedure would fail. That Newton’s method is likely to show poor behavior far from the solution in these problems could be anticipated from the large variation in the compositional derivatives of the fugacities that are evident in Figures 2-5, for example. Adaptations of Newton’s method are available that can increase the region of convergence, especially when the equations being solved are the necessary conditions for an optimum (see, for example, Chapter 6 of Dennis and Schnabel, (19831.)
Conclusions The successive substitution procedures for flash calculations and stability analysis are generally assumed to be “safe” computational tools, but that is evidently not always a correct assumption. In this paper we have described several kinds of phase models that can
produce nonconvergent behavior and have demonstrated the behavior with some numerical examples. The analysis presented shows how the difficulties arise. We have also described briefly some procedures for modifying successive substitution in order to restore convergent behavior.
Acknowledgment R.A.H. acknowledges the financial support of the Natural Sciences and Engineering Research Council of Canada.
Nomenclature D = tangent plane distance, measure of stability Ex,EY = matrices in the stability analysis fl = fugacity, substance i Ki = yJxi = equilibrium ratio of mole fractions m = damping factor Mx, MY = matrices in the stability analysis ni = molar amount, substance i R = gas constant T = absolute temperature u,ui = an eigenvector wi = mass fraction of substance i xi = mole fraction of substance i in the x phase X = diagonal matrix, elements xi yi = mole fraction of substance i in they phase Kixi = molar amount before normalizing to mole Yi fraction Y = diagonal matrix, elements yi zi = a feed mole fraction, flash calculations Greek Letters
p = fraction of the mixture in the y phase yi = activity coefficient q!~i = fugacity coefficient II = eigenvalue
II* = eigenvalue, damped successive substitution 9 =DIRT = dimensionless“tangentplane distance”,phase stability measure
Appendix Stability Analysis. The eigenvalues and the corresponding eigenvectors of the stability analysis process satisfy
- m u , = AiUi
(A.1)
We define X1” as the diagonal matrix with elements = X’”hi. This yields
xi“ and a set of vectors vi
-x1’2m1/2vi = sv, = A,Vi
(A.2)
Since S is a symmetric matrix, the eigenvalues are all real numbers and the vi eigenvectors are orthogonal. Let the eigenvectors be scaled such that $vi = 1. Then v Tivj = ad, or
u?x-’u.J = 6.. V
(-4.3)
As was discussed in the main text, EXE -XMx has a zero eigenvalue, ill = 0, with eigenvector u1 = (XI,..., X N ) ~ .Inserting this eigenvector in eq A.3, the result is T -1 T u,X uj=euj=O, j = 2 , 3 ,...,N
where eT= (1,1,..., 1).
(A.4)
Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 966 Equations A.3 and A.4 are required in the analysis of the flash calculation process. Isothermal Flash. So long as 0 < p < 1, eq 15 for the eigenvalues of the isothermal flash process can be rearranged to give
[OW+ ( 1 - /"IU
+ A[~(x-' - eeT>+ (1 - P)(Y-l
- eeT)lu= 0
(A.5)
Au+ABu=O
(A.6)
where A and B are symmetric and positive-definite. All eigenvalues are therefore real (Wilkinson, 1977). A further rearrangement, followed by premultiplication by the transpose of the eigenvector, uT,gives the relation;
+
tW + A(x-' - eeT>lu =o
+
- eeT)lu= 0
(A.7)
Consider the first term. Expand u as a linear combination of the eigenvectors of the stability analysis matrix for the z phase:
UT,
"his equation is satisfied by u = A = 2: for i 2 2. It cannot be concluded, however, that A = A; = 0 and u = x is the only additional linearly independent solution. In fact, eq A.14 is satisfied trivially for any A when u = x. It is also satisfied by any A = A:, u = ui qx. The source of the difficulty is that eq A.5 was arrived at by inverting a matrix which is singular at p = 1 (and at p = 0). To obtain the missing eigenvalue, it is necessary to consider the limiting behavior of eq A.5. We rewrite eq A.5 letting p = 1 - A and consider what happens as A 0. The eigenvector of interest is assumed to be of the form u = x Af where f is a (nonzero) vector. In terms of A, eq A.5 is
-
[W
(A.8)
+
+ A,(x-'- eeT)lu+ A[W - W + A1(Y-' - X-')]U = 0 (A.15)
Substituting into this equation for u in terms of A yields
A[W - W
N
u=
(A.14)
+
This matrix eigenvalue problem is of the form
puT[w II(X-' - eeT>lu (1 - /?)uT[W 4A ( Y 1
Limiting Cases: j9 = 0 or /9 = 1. If the value p = 1 is inserted in eq A.5, the result is
+ A1(Y-' - X-')]X + A[W + Al(X-' - eeT)lf+ A2[ ...I = 0 (A.16)
i=l
Neglecting the second-order term in A and premultiplying by xTresults in
Utilizing eqs A . l , A.3, and A.4 yields
[M" + A(X-' - eeT)]u=
c N
ci(A - A;)ui (A.9)
xT[W- W
+ Al(Y-'
- X-')]X = 0
(A.17)
i=2
This scalar equation may be solved for A. A convenient way to write the result is
and multiplying by uTresults in
+
uT[W A(X-l
c N
- eeT)lu=
c?(A - I,;)
(A.lO)
A1=
-
i=2
Similarly, if u is expanded on the eigenvectors of the stability analysis matrix for the y phase, the second term in eq A.7 becomes
uT[W+ 8 Y - l - eeT)lu=
c
d:(A - Ai)
(A.ll)
and eq A.7 reduces to
(A.18)
The vector x-y is the tie-line direction. This vector may be expanded on the eigenvectors of they phase stability matrix to give
N
i=2
N
(x- YIT W(x - y) (x - yITY-Yx - y)
II, =
N
N
i=2
i=2
2 d : A ~ ld:~
(A.19)
and thus it is evident that A1 is located between the largest and the smallest of the y phase eigenvalues. Symmetric equations apply for the case /3 = 0.
N
Literature Cited The terms on the left-hand side of eq A.12 cannot all have the same sign; hence the eigenvalue for the flash calculation procedure, A, must satisfy
Amin 5 A I"A
(A.13)
as given in eq 16 of the main text. For a binary mixture, the two sums in eq A.12 each contain only a single term. Then, if one of two phases is an ideal solution (say the y phase), with both AT = 0, the only solution possible for A in the flash calculation process is A = Ai; i.e., the eigenvalue controlling the behavior of the flash calculation is identical to the eigenvalue in the nonideal phase, regardless of the distribution between the two phases.
Ammar, M. N.; Renon, H. The Isothermal Flash Problem: New Methods for Phase Split Calculations. AIChE J. 1987,33,926939. Bishnoi, P. R.; Gupta, A. K.; Englezos, P.; Kalogerakis, N. Multiphase Equilibrium Flash Calculations For Systems Containing Gas Hydrates. Fluid Phase Equilib. 1989,53,97-104. Chen, C.-IC; Duran, M. A.; Radosz, M. Phase Equilibria in Polymer Solutions. Block-Algebra,Simultaneous Flash Algorithm Coupled with SAFT Equation of State, Applied to Single-Stage Supercritical Antisolvent Fractionation of Polyethylene. Znd. Eng. Chem. Res. 1993,32,3123-3127. Crowe, A. M.; Nishio, M. Convergence Promotion in the Simulation of Chemical Processes. The General Dominant Eigenvalue Method. AIChE J. 1975,21,528-533. Dennis, J. E., Jr.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; PrenticeHall: Englewood Cliffs, NJ, 1983.
966 Ind. Eng. Chem. Res., Vol. 34, No. 3, 1995 Gmehling, J.; Onken, U.; Arlt, W.; Grenzheuser, P.; Weidlich, U., Kolbe, B.; Rarey-Nies, J. DECHEMA Chemistry Data Series, Vol. I. Vapor-Liquid Equilibrium Data Collection. Part 5. Carboxylic Acids, Anhydrides, Esters; DECHEMA: Frankfurti Main, 1982. Kang, C. H.; Sandler, S. 1. Phase Behavior of Aqueous TwoPolymer Systems. Fluid Phase Equilib. 1987, 38, 245-272. Mehra, R. K.; Heidemann, R. A.; Aziz, K. A. Computation of Multiphase Equilibrium for Compositional Simulation. SOC.Pet. Eng. J. 1982,22, 61-68. Mehra, R. K.; Heidemann, R. A.; Aziz, K. A. An Accelerated Successive Substitution Algorithm. Can. J. Chem. Eng. 1983, 61, 590-596. Michelsen, M. L. The Isothermal Flash Problem. I: Stability. Fluid Phase Equilib. l982a, 9, 1-20. Michelsen, M. L. The Isothermal Flash Problem. 11: Phase-Split Calculation. Fluid Phase Equilib. 1982b, 9, 21-40. Michelsen, M. L. Calculation of Hydrate Fugacities. Chem. Eng. Sei. 1991,46, 1192-1193. Munck, J.;Skjold-J@rgensen,S.; Rasmussen, P. Computations of the Formation of Gas Hydrates. Chem. Eng. Sci. 1988, 43, 2662-2672.
Nghiem, L. X.; Aziz, K. A,; Li, Y. K. A Robust Iterative Method for Flash Calculations Using the Soave-Redlich-Kwong or the Peng-Robinson Equation of State. SOC.Pet. Eng. J. 1983, 23, 51-61. Stryjek, R.; Vera, J. H. Vapor-Liquid Equilibrium of Hydrochloric Acid Solutions with the PRSV Equation of State. Fluid Phase Equilib. 1986,25, 279-290. Wilkinson, J. H. The Algebraic Eigenvalue Problem: Clarendon Press: Oxford 1977. Received for review July 11, 1994 Accepted October 11, 1994 @
IE940430W
Abstract published in Advance ACS Abstracts, January 1, 1995. @
