Ind. Eng. Chem. Res. 1999, 38, 3549-3556
3549
Calculation of Liquid-Liquid Equilibrium Based on the Global Stability Analysis for Ternary Mixtures by Using a Novel Branch and Bound Algorithm: Application to UNIQUAC Equation Zhu Yushan* and Xu Zhihong Institute of Chemical Metallurgy, Chinese Academy of Sciences, Beijing 353# Beijing 100080, People’s Republic of China
In this paper, a novel branch and bound algorithm based on the idea of Horst (Math Program. 1976, 10, 312) was developed to minimize the tangent plane distance function (TPDF) of the stability analysis problem where nonideal liquid phases are described by the UNIQUAC activity coefficient equation. Some simple transformations were used to construct the DC (i.e., difference of two convex functions) programming of TPDF that does not rely on any newly introduced parameters except the inner ones of the UNIQUAC equation. A compact partition of the feasible space and a valid underestimation function of the TPDF depending on the linear characteristics of the stability analysis problem were developed to guarantee the new branch and bound algorithm obtaining the -global convergence to the global solution of the TPDF. In comparison with the algorithms presented by Sun and Seider (Fluid Phase Equilib. 1995, 103, 213) and McDonald and Floudas (AlChE J. 1995, 41, 1798), the novel branch and bound algorithm is slightly slower; however, it can provide complete reliability for solving the phase stability analysis problem. Liquid-liquid equilibrium compositions were then calculated by the Newton-Raphson method (Fluid Phase Equilib. 1982, 9, 21) on the basis of the global minimum of TPDF. The preliminary calculation results for three ternary systems showed that the novel branch and bound algorithm can solve the global stability problem effectively. 1. Introduction The phase and chemical equilibrium problem is crucial for chemical separation and process design. Recently, Seider and Widagdo,5 and Jiang et al.6 have reviewed the analysis and calculation of multiphase equilibria involving chemical reactions. McDonald and Floudas7 introduced their GLOPEQ, a new computational tool for the phase and chemical equilibrium problem. Cisneros et al.8 supplied a new algorithm for the computation of the simultaneous chemical and physical equilibrium involved in simulation of reactive flash operations, calculations of phase diagrams, and the determination of reactive azeotropes. Baker et al.9 proved that a necessary and sufficient condition for an equilibrium solution corresponding to a global minimum of the Gibbs free energy was that the tangent plane corresponding this solution lies on or below the Gibbs free energy surface for all possible values of the composition. Michelsen4,10 first presented different practical implementations of this criterion and used the phase stability results to generate initial points that were applied in the search for an equilibrium solution with a lower value of the Gibbs free energy; these methods have been found to be more reliable and robust than minimizing the Gibbs function directly. Swank and Mullins11 compared various methods for calculating liquid phase-splitting problems. Nagarajan et al.12,13 used the Helmholtz free energy as the proper thermodynamic function describing equilibrium. Eubank et al.14 integrated the Gibbs free energy surface for a single hypothetical phase to obtain a maximal area rather than a minimal tangent plane distance. This has widely been used in all kinds of calculations of phase * Present address: Department of Applied Chemistry, Saga University, Saga 840-8502, Japan. E-mail:
[email protected]. Tel: 0081-952-288670. Fax: 0081-952-288591.
equilibria, especially in critical point region. Recently, this method has obtained many further developments.15-21 Smith et al.22 obtained the necessary and sufficient conditions for stability of the system where chemical reaction may or may not take place, using the KarushKuhn-Tucker optimality conditions as the basis of their consideration. Sun and Seider2 used a two-stage method and employed a homotopy algorithm in an attempt to find all stationary points; however, this method could still fail for certain cases. McDonald and Floudas3,7,23-25 showed that for activity coefficient models the problem could be reformulated to make it amenable to be solved by global optimization techniques, for example, the -global optimization method and the branch and bound algorithm, which did mathematically guarantee that the correct answer was found for phase and chemical equilibrium problem. Wasylkiewicz et al.26 introduced a global algorithm for implementation of the Gibbs tangent plane analysis and liquid-liquid phase equilibrium calculation for complex liquid mixtures which involved locating all of the stationary points of the tangent plane distance function by tracking its ridges and valleys. Hua et al.27-29 used an interval Newton/ generalized bisection (IN/GB) method to solve the phase stability problem. We30 used a simulated annealing algorithm to obtain the global solution of the tangent plane distance function (TPDF) to the liquid-phase equilibrium problem. In this paper, a novel branch and bound algorithm was implemented by introducing a simple DC formulation for TPDF without the separable function assumption, developing a general compact partition for the feasible region on the basis of the characteristics of the linear constrained conditions associated with the phase stability problem, and generating a valid and general underestimation function for TPDF. It should be noted that there is no further assumption on the parameters
10.1021/ie990104m CCC: $18.00 © 1999 American Chemical Society Published on Web 08/07/1999
3550 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999
used in our algorithm except for the inner ones of the UNIQUAC equation, so that the novel branch and bound algorithm used in this paper has complete -global convergence for the phase stability problem, and its calculation speed is acceptable, though slightly slow. Three ternary systems with up to three phases were used to analyze the phase stability at constant temperature and pressure, and the calculation results stated that the novel branch and bound algorithm was completely reliable.
N
F(x) )
xi ∑ i)1
{
Z 2
2. Mathematical Formulation of the Gibbs Tangent Plane Criterion Described by UNIQUAC Equation Gibbs tangent plane criterion was first shown by Gibbs31,32 as a necessary and sufficient condition for absolute stability of a mixture at a temperature T, a pressure P, and overall composition z if the Gibbs energy surface is at no point below the plane tangent to the surface at the given composition. A recent explanation and derivation of this condition was given by Baker et al.,9 Michelsen,4,10 Smith et al.,22 and Jiang et al.6 The optimization formulation for the phase stability analysis with the normalized and nonnegative conditions for the component variables is presented as follows:
∑ i)1
xiµ0i (z)
)
N
(1)
N
xi ) 1 ∑ i)1
0 e xi e 1
(2) (3)
where F(x) is the tangent plane distance function, g(x) is the molar Gibbs free energy surface for a single hypothetical phase, T(x) is the tangent plane to the molar Gibbs free energy surface at the charge z, and µi(x) and µ0i (z) represent the chemical potentials at x and z respectively. xN is the Nth component variable in the above constrained conditions that can be obtained by the following equation according to eq 2:
xi ∑ i)1
(4)
Therefore, there are only N - 1 independent variables in the phase stability analysis problem. In addition, the equality constrained condition, eq 2, can be replaced by the following inequality:
qixi
qi ln
- q′i ln
N
qjxj ∑ j)1
[
N
- µ0i (z)
q′jxj ∑ j)1
∑
∑
∑
∑ l)1
where li is given as
(8)
xi denotes the mole fraction of component i, and ∆Gfi represents the Gibbs free energy of the pure component at the system temperature T, constructed so that the Gibbs energy content of the elemental species is zero at T. ri, qi, and q′i are pure component structure parameters, and τij is the binary interaction parameter. They are all positive real numbers according to their definitions in the UNIQUAC equation. Z is the lattice coordination number (always given as 10). Obviously, F(x) as defined in eq 6 is nonconvex. After being introduced to some simple transformations, this function can be expressed as the difference of two convex functions, as follows:
F(x) ) p(x) - q(x) where N
p(x) )
Z 2
xi ∑ i)1
{ } {
Z
∆Gfi - µ0i (z) -
2
xi ∑ i)1
ln
q′i ln
N
q(x) )
xi ∑ i)1
{
Z 2
(9)
qi ln ri + ln ri +
xi
N
qi ln qi +
(5)
Therefore, this inequality of eqs 5 and 3 constitute the composition space, which is a simplex for the N - 1 independent variables. Obviously, it is a convex set. The tangent plane distance function described by the UNIQUAC equation33 is presented as follows:
(6)
]
N-1
xi e 1 ∑ i)1
}
N
q′jτjixj ∑ j)1
Z Z µi(x) ) ∆Gfi + 1 - qi ln φi + qi ln θi + li 2 2 N N φi N τijθ′j ljxj + q′i - ln τjiθ′j - q′i (7) N xi j)1 j)1 j)1 θ′lτlj
N-1
xN ) 1 -
∑ j)1
li ) Z/2‚(ri - qi) - (ri - 1) ∀i ∈ N
xi {µi(x) - µ0i (z)} ∑ i)1 s.t.
)
where µ0i (z) represents the chemical potential of each component at z, defined as
N
min F(x) ) g(x) - T(x) ) g(x) -
(
rixi Z ∆Gfi + 1 - qi ln + N 2 rjxj
N
rjxj ∑ j)1
Z + qi ln 2
qi ln
xi N
rjxj ∑ j)1
+
N
qjxj ∑ j)1 xi
N
q′j xj + q′i ln ∑ j)1
xi
N
}
(10)
q′jτjixj ∑ j)1
+ q′i ln xi
}
(11)
P(x) and q(x) are all convex functions because they are all linear combinations of the simple linear functions
Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3551
and three different convex functions whose convexities were demonstrated by McDonald and Floudas.23 3. The Novel Branch and Bound Algorithm For nonconvex programming problems with nonseparable functions, a branch and bound algorithm was proposed by Horst1 that forms the basis of the algorithm used in this work. The basic idea of the branch and bound algorithm is to successively refine the feasible region, thereby solving convex subproblems in these subregions and generating a nondecreasing subsequence of lower bounds that will converge on a global solution of the original nonconvex optimization problem. In this section, a novel compact partition of the feasible region is developed on the basis of the linear constrained characteristics of the stability analysis problem, and then a valid convex underestimation function that does not rely on separable conditions is constructed in each feasible subregion. 3.1. Compact Partition of the Feasible Region. The initial convex simplex, i.e., the feasible composition space constituted by eqs 5 and 3, has N vertices and can be described by N superplanes with N - 1 dimensions. The kth superplane of this simplex is described as follows: N-1
aki xi ) bk ∑ i)1
aki
k ) 1, 2, ..., N
(12)
the Mth superplane is given as follows:
xM i
)
1
N-1
xji ∑ N-1 j)1
and are the linear equation constrained where, coefficients which determine kth superplane of above simplex. Since k - 1 vertices except for the kth vertex in this simplex locate in this superplane, then the above constrained equation coefficients can be obtained by solving the following linear equation group: N-1 N
aki xji ) bk ∑ ∑ i)1 j)1
Therefore, the initial simplex can be divided into N - 1 subsimplexes by connecting this gravity center with other vertices using superplanes. The linear constrained conditions can be obtained by solving constrained equations similar to eqs 12 and 13. It should be noted that feasible points lie in the inner space of each subsimplex. The newly generated subsimplexes can be divided into more refined subregions by using the same method as above till the algorithm has converged on the global solution. 3.2. Generating Convex Underestimation Function of the F(x). After the compact partitioning of the feasible region is finished, a valid lower bound of the objective function, which does not rely on separable conditions, can be developed depending on the characteristics of the concave function. The -q(x) in eq 7 is a concave function in each subsimplex, so that a linear function, Q(x), which passes through N vertices of each subsimplex, i.e., {x1, -q(x1)}, {x2, -q(x2)}, ..., and {xN, -q(xN)}, can be constructed in this subregion. Its linear constrained equation coefficients can be obtained by solving the following linear equation group:
k ) 1, 2, ..., N
(13)
k Q k Q aQ ∑ i xi + aN Q(x ) ) b i)1
where xji is the ith coordination of the jth vertex in the above superplane, and j * k. Equation 12 can also be written as follows:
∑ i)1
k
[
∑ j)1
b
xi )
N-1
(akj )2]1/2
(14)
N-1
[
∑ j)1
aki
N-1
D )
∑ i)1
N-1
[
∑ j)1
bk
xi -
(akj )2]1/2
N-1
[
∑ j)1
Q(x) )
bQ
N-1
-
aQ N
∑ i)1
aQ i aQ N
xi
L(x) ) p(x) + Q(x) e F(x)
k ) 1, 2, ..., N
(akj )2]1/2
(19)
Obviously, Q(x) e -q(x), according to the characteristics of the concave function in each subsimplex. Then, the following function
(akj )2]1/2
Therefore, the distance between the kth vertex {xki , xk2, k } and its opposite superplane can be computed ..., xN-1 by k
k ) 1, 2, ..., N (18)
Q Q are the linear function coefwhere aQ i , aN, and b ficients of Q(x) in each subsimplex. It should be noted that these coefficients are different from those defined by eqs 12 and 13, because the latter correspond to the constrained subsimplex. Since Q(x) is a linear equation, it can also be expressed as
j*k
aki
(17)
N-1
bk
N-1
i ) 1, 2, ..., N
(20)
is a valid underestimation function of the original objective function in each subregion. Here, we only use the concave function characteristics of -q(x), and the above underestimation function does not rely on the separable function where p(x) is a convex function presented in eq 10 and Q(x) is a linear function. So, the above the underestimation function, i.e., L(x), is convex. The following convex problem will then yield a valid underestimation of the global solution in the current subregion:
(15)
min L(x) ) p(x) + Q(x)
let
(21)
N-1
D ) min{D } k ) 1, 2, ..., N M
k
(16)
where DM is the minimal one among the those distances defined by eq 14. In this paper, the Mth superplane will be divided into N - 1 subregions. The gravity center of
s.t.
aki xi - bk e 0 ∑ i)1
k ) 1, 2, ..., N
(22)
For successive refinements of the feasible region, the lower bound supplied by the above minimization prob-
3552 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 Table 1. Results of Phase Stability Analysis and Equilibrium Calculation for TAME (1)-Methanol (2)-Water (3) at T ) 298 K and P ) 1 atm, Conditions (i) phase equilibrium calculation comp. 1 2 3 1 2 3
solution charge L φL1 ) 1.0 global LL φL1 ) 0.423
ZL1 i 0.293 3 0.268 5 0.438 2 0.032 57 0.248 85 0.718 58
ZL2 i
phase stability analysis ZL3 i
x*i 0.017 15 0.198 60 0.784 25 0.032 57 0.248 85 0.718 58
0.484 36 0.282 90 0.232 75
F*
CPU (s)
N
-0.054 94
71.79
2572 (188)
95.41
3252 (216)
0.0
Table 2. Results of Phase Stability Analysis and Equilibrium Calculation for TAME (1)-Methanol (2)-Water (3) at T ) 298 K and P ) 1 atm, Conditions (ii) phase equilibrium calculation comp. 1 2 3 1 2 3
solution charge L φL1 ) 1.0 global LL φL1 ) 0.564
ZL1 0.19 0.31 0.50 0.051 43 0.292 93 0.655 64
Z12
0.369 00 0.332 05 0.298 95
lem will be tighter than that supplied by any of its predecessor subregion. Horst and Pardalos34 proved that the algorithm used in this paper converges finitely to an -global solution of the original problem if the subdivision rule described is used. It should be noted that only one underestimation function is required in each subregion. An upper bound is given by computing the function value of F(x) on the global solution of eqs 21 and 22 in each subregion, and it is used to delete the unnecessary regions whose lower bound is higher than the lowest upper bound in all subregions left in each calculation step, in order to save machine memory and increase the computing speed. The detailed steps of the branch and bound algorithm are presented in Appendix A. 4. Phase Stability Calculation Results And discussion There are three ternary systems for phase stability analysis by the novel branch and bound algorithm in this section. They all belong to liquid-liquid phase equilibria without reaction taking place in the system. The Gibbs free energy formation term for the vapor phase is removed through the approximation. All computations were performed on a Pentium 586/120 machine. In this paper, CPU seconds reported represent the real time taken to solve the phase stability problem by the novel algorithm. 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, so the stability test for one of the compositions in the equilibrium phase which satisfies the liquid-liquid equilibrium would be sufficient. The global tolerance of the novel branch and bound algorithm is 10-5, where the Rosen Gradient Projection method35 is used to solve each convex subproblem generated in each subregion. 4.1. Calculation of Phase Equilibria Based on the Stability Analysis. If the global minimum of the TPDF for the overall composition is negative, the system is unstable and will separate into two or more phases at the current temperature and pressure. The local and the global equilibrium compositions in the above phases were calculated by the Newton-Raphson method (with the sum of squared residuals less than 10-8) on the basis of the stability analysis results in this paper this method was first implemented by Michelsen4 as a standard method.
phase stability analysis Zl3
x* 0.036 86 0.260 74 0.702 39 0.051 43 0.292 93 0.655 64
F*
CPU (s)
N
-0.011 23
88.16
2929 (167)
106.67
3684 (243)
0.0
4.2. Calculation Systems. In the following tables, the set of liquid phases is denoted P t {k}, where each is designated Lk ⊆ P. For any postulated solution, the mole fractions of each component i in some phase k are given by zki , so that ∑izki ) 1 ∀k ∈ P. The parameter φk corresponds to the fraction of total mass in phase k, so that ∑iφk ) 1. x*i is the mole vector at a global function, with F* being the corresponding tangent plane distance function. N is the iteration the algorithm took to converge to the global solution for each TPDF, and the numbers in parentheses below show the remaining subregions when the algorithm terminates. System 1. TAME + Methanol + Water. The binary interaction parameters and pure component structure parameters of the UNIQUAC equation for tert-amyl methyl ether (TAME), methanol, and water were given by Arce et al.36 at a temperature of 298 K and a pressure of 1 atm. This is a simple two-phase ternary system. Two different charges in two-phase region were chosen to check the reliability of the novel branch and bound algorithm. Tables 1 and 2 represent the calculation results of stability analysis and phase equilibria for this system at two different charges: z ) {0.2933, 0.2685, 0.4382} and {0.19, 0.31, 0.50}. The final phase number, phase fractions, and the global LLE compositions are in good agreement with the experimental data presented by Arce et al.36 System 2. Water + Nitromethane + 1-Hexanol. A three-phase region for the water, nitromethane, and 1-hexanol ternary mixture was detected experimentally at 21 °C, and the experimental data were presented in the DECHEMA Chemistry Data Series (Vol. V/2, p 69).37 Wasylkiewicz et al.26 studied this system and found the three-phase region by the UNIQUAC binary interaction parameters reported by Sorensen and Arlt (Vol. V/3, pp 422 and 423);37 these parameters are also used in this study. The calculation results of stability analysis and phase equilibria for two different charges, z ) {0.4, 0.4, 0.2} and {0.5, 0.4, 0.1}, are presented in Tables 3 and 4. The global LLLE compositions are in good agreement with the phase diagrams presented by Sorensen and Arlt37 and Wasylkiewicz et al.,26 respectively. System 3. Ethylene Glycol + Lauryl Alcohol + Nitromethane. The UNIQUAC equation binary interaction parameters and pure component structure parameters for ethylene glycol, lauryl alcohol, and nitromethane were supplied by Chakravarty et al.38 This example was
Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3553 Table 3. Results of Phase Stability Analysis and Equilibrium Calculation for Water (1)-Nitromethane (2)-1-Hexanol (3) at T ) 294.15 K and P ) 1 atm, Conditions (i) phase equilibrium calculation comp. 1 2 3 1 2 3 1 2 3
solution charge L φL1 ) 1.0 local LL φL1 ) 0.753 global LLL φL1 ) 0.338 φL2 ) 0.260
ZL1
Zl2
0.4 0.4 0.2 0.212 55 0.522 13 0.265 32 0.125 01 0.762 54 0.112 46
0.970 19 0.028 51 0.001 30 0.970 47 0.028 23 0.001 30
phase stability analysis Zl3
x*
0.261 23 0.336 51 0.402 26
0.978 53 0.020 70 0.000 77 0.119 99 0.773 80 0.106 21 0.125 01 0.762 54 0.112 46
F*
CPU (s)
-0.201 50
134.12
7457 (706)
-0.004 90
189.11
6864 (475)
329.17
14 116 (1053)
0.0
N
Table 4. Results of Phase Stability Analysis and Equilibrium Calculation for Water (1)-Nitromethane (2)-1-Hexanol (3) at T ) 294.15 K and P ) 1 atm, Conditions (ii) phase equilibrium calculation comp. 1 2 3 1 2 3 1 2 3
solution charge L φL1 ) 1.0 local LL φL1 ) 0.419 global LLL φL1 ) 0.424 φL2 ) 0.120
ZL1
Zl2
0.5 0.4 0.1 0.970 60 0.028 06 0.001 37 0.125 01 0.762 54 0.112 46
0.160 85 0.668 04 0.171 10 0.970 47 0.028 23 0.001 30
phase stability analysis Zl3
x*
0.261 23 0.336 51 0.402 26
0.976 90 0.022 41 0.000 69 0.261 56 0.282 55 0.455 89 0.125 01 0.762 54 0.112 46
F*
CPU (s)
-0.19256
145.77
8276 (662)
-0.020 26
115.34
3894 (385)
329.17
14 116 (1053)
0.0
N
Table 5. Results of Phase Stability Analysis and Equilibrium Calculation for Ethylene Glycol (1)-Lauryl Alcohol (2)-Nitromethane (3) at T ) 295 K and P ) 1 atm, Conditions (i) phase equilibrium calculation comp. 1 2 3 1 2 3 1 2 3
solution charge L φL1 ) 1.0 local LL φL1 ) 0.630 global LLL φL1 ) 0.608 φL2 ) 0.029
ZL1 0.4 0.3 0.3 0.270 78 0.473 02 0.256 20 0.278 99 0.491 91 0.229 10
Zl2
0.619 86 0.005 62 0.374 52 0.027 76 0.002 06 0.970 18
phase stability analysis Zl3
x*
0.69280 0.00399 0.30321
0.754 15 0.002 28 0.243 57 0.023 38 0.001 79 0.974 83 0.278 99 0.491 91 0.229 10
F*
CPU (s)
N
-0.11396
293.09
15 300 (869)
-0.05876
157.20
7247 (578)
499.88
23218 (1244)
0.0
Table 6. Results of Phase Stability Analysis and Equilibrium Calculation for Ethylene Glycol (1)-Lauryl Alcohol (2)-Nitromethane (3) at T ) 295 K and P ) 1 atm, Conditions (ii) phase equilibrium calculation comp. 1 2 3 1 2 3 1 2 3
solution charge L φL1 ) 1.0 local LL φL1 ) 0.637 global LLL φL1 ) 0.608 φL2 ) 0.363
ZL1 0.2 0.3 0.5 0.296 72 0.469 50 0.233 78 0.278 99 0.491 91 0.229 10
Zl2
0.030 01 0.002 11 0.967 88 0.027 76 0.002 06 0.970 18
phase stability analysis Zl3
x*
0.692 80 0.003 99 0.303 21
0.012 54 0.001 11 0.986 34 0.714 12 0.003 29 0.282 59 0.278 99 0.491 91 0.229 10
F*
CPU (s)
-0.228 27
101.33
4826 (366)
-0.027 00
342.52
17 250 (507)
499.88
23218 (1244)
0.0
N
Table 7. Comparison of Three Different Algorithms for System 3 at the Overall Composition: z ) {0.2, 0.3, 0.5} algorithm used
compiler
computing machine
stability analysis CPU (s)
final phase equilibria solution
compositions ) {0.030 01, 0.002 10, 0.967 89} zL2 ) {0.296 70, 0.469 50, 0.233 80} L1 z ) {0.278 99, 0.491 91, 0.229 10} zL2 ) {0.027 76, 0.002 06, 0.970 18} zL3 ) {0.692 80, 0.003 90, 0.303 21} zL1 ) {0.278 99, 0.491 91, 0.229 10} zL2 ) {0.027 76, 0.002 06, 0.970 18} zL3 ) {0.692 80, 0.003 90, 0.303 21}
zL1 HOMPEQ
FORTRAN
DEC Station 3100 workstation
22.82
local LLE
GLOPEQ
C
Hewlett-Packard 9000/730 workstation
29.70
global LLLE
in this paper
C
Pentium 586/120
499.88
global LLLE
studied by Null,39 Gautam and Seider,40,41 Sun and Seider,2 and McDonald and Floudas.3 A three-dimensional plot of ∆Gm at 295 K and 1 atm by the UNIQUAC equation was illustrated by Sun and Seider;2 obviously, a three-phase region exists in this system. The calculation results of stability analysis and phase equilibria for two different charges: z ) {0.4, 0.3, 0.3} and {0.2, 0.3, 0.5}, are presented in Tables 5 and 6. They are in
good agreement with those of Sun and Seider2 and McDonald and Floudas.4 4.3. Efficiency and Reliability. The efficiency of the novel branch and bound algorithm was compared with that of a program for multiphase equilibria by using a homotopy-continuation algorithm, HOMPEQ, developed by Sun and Seider2 and GLOPEQ, a package developed by McDonald and Floudas4 for the phase and chemical
3554 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999
equilibrium problem by global optimization algorithms. For the ethylene glycol, lauryl alcohol, and nitromethane system, the comparison of the CPU seconds consumed on one stability analysis problem is somewhat qualitative, because different algorithms were made by different compilers and run on different calculation machines. It should be noted that the HOMPEQ did not get any negative value of the TPDF at the local LLE solution, and then it stopped at this local equilibrium composition for the system. However, the GLOPEQ and our algorithm all found the global minimum of the TPDF, i.e., -0.027 00, described in Table 6; then, the local LLE composition for the system is unstable. The true threephase LLLE compositions described in Tables 6 and 7 were obtained by Newton-Raphson method based on above stability analysis results. It should be noted that a key assumption of the GLOPEQ in using the UNIQUAC equation is that q(x) must be a separable function in the DC formulation for the TPDF,4 but it is unnecessary in the algorithm used in this paper. The principal advantage of the novel branch and bound algorithm used in this paper lies in its complete reliability in spite of its slightly slow computing speed.
eqs ) equations F(x) ) tangent plane distance function ∆G0i ) Gibbs free energy in standard state L(x) ) the underestimation function of the TPDF M ) the initial feasible region min ) minimum P ) pressure, atm p(x) ) convex part of the TPDF -q(x) ) concave part of the TPDF Q(x) ) the underestimation function of the -q(x) q ) structural parameter in the UNIQUAC equation q′ ) structural parameter in the UNIQUAC equation r ) structural parameter in the UNIQUAC equation T ) temperature, K T(x) ) tangent plane of the Gibbs free energy surface at z xi ) mole fractions of component i in the liquid phase z ) the overall composition of a liquid mixture
5. Conclusion
Subscripts
It has been shown how the global solution of the phase stability problem can be obtained when the nonideal liquid phases are modeled by the UNIQUAC equation. A number of simpler changes are made in the way that the TPDF is constructed so that it could be recast as the difference of two convex functions without the separable assumption, and then a novel algorithm based on the idea of Horst1 is developed. The most important advantage of this algorithm is the use of novel compact partitions instead of rectangular ones and a unique refining rule such that the novel branch and bound algorithm does not rely on the concept of convex envelopes and can handle nonseparable function. The preliminary results for three ternary mixtures show that the novel branch and bound algorithm is reliable to obtain the global solution of the TPDF by the UNIQUAC equation with complete -global convergence. This algorithm can also be applied into the UNIFAC, modified Wilson, and ASOG activity coefficient equations because the TPDF described by them can also be transformed into the similar DC formulation.
i ) component index j ) component index L1 ) the first liquid phase at equilibrium L2 ) the second liquid phase at equilibrium L3 ) the third liquid phase at equilibrium
Acknowledgment The authors are grateful to Prof. C. A. Floudas (Princeton University) and Prof. Yuan Yaxiang (ICMSEC, P.R.C.) for their help in providing us with many precious resource data and some helpful suggestions. The authors also gratefully acknowledge financial support from the National Natural Science Foundation of China under Grant 594340402 and the Laboratory of Computer Chemistry in CAS. Nomenclature Notations aki ) coefficient of the linear constrained equation bk ) coefficient of the linear constrained equation N ) number of components Dk ) distance between the kth vertex and its opposite superplane eq ) equation
Greek Symbols ) global convergence tolerance µ0i (z) ) chemical potential of component i at z µi(x) ) chemical potential of component i at x τij ) binary interaction parameters in the UNIQUAC equation φ ) phase mole fraction
Superscripts k ) the iterative number k - 1,pk-1 ) the pk-1th subregion in step k - 1 m ) component index Q ) the denotation of Q(x) Acronyms TAME ) tert-amyl methyl ether TPDF ) tangent plane distance function LLE ) liquid-liquid equilibrium LLLE ) liquid-liquid-liquid equilibrium
Appendix The complete novel branch and bound algorithm for the minimization of the tangent plane distance function when the nonideality of the liquid phase is modeled by the UNIQUAC equation is now presented. The initial feasible region defined by eqs 3 and 5, i.e. a convex simplex, is denoted by M. Step 1. 1.1. Divide the initial simplex M into N - 1 subsimplexes, i.e. {M1,1, ..., M1,N-1}, according to the method described in section 3.1, where M1,1 represents the first subsimplex produced in step 1. 1.2. Determine N - 1 underestimation functions L1,j(x): M1,j f R (i ) 1, ..., N - 1) on above N - 1 subsimplex according to the method described in 3.2, where L1,j(x) is the ith underestimation function produced on subsimplex M1,j in step 1. 1.3. Compute x1,j from
L1,j(x1,j) ) min L1,j(x) i ) 1, 2, ..., N - 1 (A1) on each subsimplex M1,j by the Rosen gradient projection method,35 where x1,j is the minimum of the ith underestimation function. At the same time, compute each TPDF value at x1,j on N - 1 subsimplexes, i.e.,
Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3555
F(x1,j), and compare them; finally, take the minimum among them as the upperestimation value of the TPDF in step 1. 1.4. Compute x1 from
L1(x1) ) min L1,j(x1,j)
(A2)
for i ) 1, ..., N - 1, where x1 is the minimum among all minimization values of the underestimation functions on all subsimplexes in Step 1. 1.5. Compute F(x1). If F(x1) - L1(x1) e , then stop the iteration. If not, then go to step 2. Step k. (k ) 2, 3, ...). Assume xk-1 ∈ Mk-1,Pk-1, where k-1,P k-1 is the pth subsimplex in step k - 1, and xk-1 M is the minimal underestimation value of TPDF in step k - 1. k.1. Divide subsimplex Mk-1,Pk-1 into N - 1 new subsimplexes, i.e., {Mk-1, ..., Mk,N-1}k-1,Pk-1 according to the same method as that described in step 1, where {Mk-1}k-1Pk-1 represents the 1th new subsimplex produced in Mk-1Pk-1 in step k. k.2. Determine N - 1 new underestimation functions Lk,i(x):Mk,i f R ( i ) 1, ..., N - 1) on above N - 1 subsimplexes according to the method described in step 1, where Lk,i(x) is the ith underestimation function produced on subsimplex Mk,i in step k. Obviously, min Lk,i(x) x ∈ Mk,i is computable for those i ∈ (1, ..., N - 1). k.3. Compute xk,i from
Lk,i(xk,i) ) min Lk,i(x)
(A3)
in each new subsimplex Mk,i produced in step k.1 by the same optimization method described in step 1. Then, add these N - 1 minima into the set of step k - 1 and eliminate the original minimum value on the subsimplex Mk-1,Pk-1. Then compute each F(xk,i), take the minimum as the upperestimation in step k, and erase all regions whose minimal values are higher than this upperestimation. k.4. Compute xk from
Lk(xk) ) min Lk,i(xk,i)
(A4)
for each i belonging to the subsimplexes left after step k.3, where xk is the minimal underestimation of the TPDF in step k. k.5. Compute F(xk). If F(xk) - Lk(xk) e , then stop the iteration. If not, and the iterative number is lower than the maximal, then go to step k + 1. Literature Cited (1) Horst, R. An Algorithm for Nonconvex Programming Problems. Math. Program. 1976, 10, 312. (2) Sun, A. C.; Seider, W. D. Homotopy-Continuation for Stability Analysis in the Global Minimization of the Gibbs Free Energy. Fluid Phase Equilib. 1995, 103, 213. (3) McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase Stability Problem. AIChE J. 1995, 41, 1798. (4) Michelsen, M. L. The Isothermal Flash Problem: II. PhaseSplit Calculations. Fluid Phase Equilib. 1982, 9, 21. (5) Seider, W. D.; Widagdo, S. Multiphase Equilibria of Reactive Systems. Fluid Phase Equilib. 1996, 123, 283. (6) Jiang, Y.; Smith, W. R.; Chapman, G. R. On the Geometry of Chemical Reaction and Phase Equilibria. Fluid Phase Equilib. 1996, 118, 77. (7) McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1997, 21, 1.
(8) Cisneros, E. S. P.; Gani, R.; Michelsen, M. L. Reactive Separation Systems-I. Computation of Physical and Chemical Equilibrium. Chem. Eng. Sci. 1997, 52, 527. (9) Baker, L. E.; Pierec, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Petrol. Engrs. J. 1982, 22, 731. (10) Michelsen, M. L. The Isothermal Flash Problem: I. Stability. Fluid Phase Equilib. 1982, 9, 1. (11) Swank, D. J.; Mullins, J. C. Evaluation of Methods for Calculating Liquid-Liquid Phase-Splitting. Fluid Phase Equilib. 1986, 30, 101. (12) Nagarajan, N. R.; Cullick, A. S.; Griewank, A. New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis: I. Stability Analysis and Flash. Fluid Phase Equilib. 1991, 62, 191. (13) Nagarajan, N. R.; Cullick, A. S.; Griewank, A. New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis: II. Critical Point Calculations. Fluid Phase Equilib. 1991, 62, 211. (14) Eubank, P. T.; Elhassen, A. E.; Barrufet, M. A.; Whiting, W. B. Area Method for Prediction of Fluid-Phase Equilibria. Ind. Eng. Chem. Res. 1992, 31, 942. (15) Eubank, P. T.; Hall, K. R. An Equal Area Rule and Algorithm for Determining Phase Compositions. AIChE J. 1995, 41, 924. (16) Shyu, G. S.; Hanif, N. S. M.; Hall, K. R.; Eubank, P. T. Equal Area Rule for Ternary Systems. Ind. Eng. Chem. Res. 1995, 34, 4562. (17) Shyu, G. S.; Hanif, N. S. M.; Hall, K. R.; Eubank, P. T. Maximum Partial Area Rule for Phase Equilibrium Calculations. Ind. Eng. Chem. Res. 1996, 35, 4348. (18) Shyu, G. S.; Hanif, N. S. M.; Hall, K. R.; Eubank, P. T. An Equation of State Model using the W-S Mixture Combining Rule for Carbon Dioxide-Water Phase Equilibrium Calculations. Fluid Phase Equilib. 1997, 130, 73. (19) Hanif, N. S. M.; Shyu, G. S.; Hall, K. R.; Enbank, P. T. Application of the Equal Area Rule to the Calculation of Multiphase Equilibria of Hydrocarbon-Water Mixtures with an EOS Model. Fluid Phase Equilib. 1996, 126, 53. (20) Hanif, N. S. M.; Shyu, G. S.; Hall, K. R.; Enbank, P. T. The Area of J. C. Maxell for Pure Component Fluid Phase Equilibria from Equations of State. Ind. Eng. Chem. Res. 1996, 35, 2431. (21) Elhassan, A. E.; Tsvetkov, S. G.; Craven, R. J. B.; Stateva, R. P.; Wakeham, W. A. A Rigorous Mathematical Proof of the Area Method for Phase Stability. Ind. Eng. Chem. Res. 1998, 37, 1483. (22) Smith, J. V.; Missen, R. W.; Smith, W. R. General Optimality Criteria for Multiphase Multireaction Chemical Equilibrium. AIChE J. 1993, 39, 707. (23) McDonald, C. M.; Floudas, C. A. Decomposition Based and Branch and Bound Global Optimization Approaches for the Phase Equilibrium Problems. J. Global Optim. 1994, 5, 205. (24) 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. 1995, 34, 1674. (25) McDonald, C. M.; Floudas, C. A. Global optimization for the phase and chemical equilibrium problem: Application to the NRTL equation. Comput. Chem. Eng. 1995, 19, 1111. (26) Wasylkiewicz, S. K.; Sridhar, L. N.; Doherty, M. F.; Malone, M. F. Global Stability Analysis and Calculation of Liquid-Liquid Equilibrium in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1996, 35, 1395. (27) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Reliable Prediction of Phase Stability using an Interval-Newton Method. Fluid Phase Equilib. 1996, 116, 52. (28) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Reliable Phase Stability Analysis for Cubic Equation of State Models. Comput. Chem. Eng. 1996, 20, S395. (29) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Enhanced Interval Analysis for Phase Stability: Cubic Equation of State Models. Ind. Eng. Chem. Res. 1998, 37, 1519. (30) Zhu, Y.; Xu, Z. A Reliable Prediction of the Global Phase Stability for Liquid-Liquid Equilibrium through the Simulated Annealing Algorithm: Application to NRTL and UNIQUAC Equations. Fluid Phase Equilib. 1999, 154, 55. (31) Gibbs, J. W. Graphical Methods in the Thermodynamics of Fluids. Trans. Conn. Acad. Arts Sci. 1987, 2, 311.
3556 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 (32) Gibbs, J. W. A Method of Geometrical Representation of the Thermodynamic Properties of Substances by Means of Surfaces. Trans. Conn. Acad. Arts Sci. 1873, 2, 382. (33) Anderson, T. F.; Prausnitz, J. M. Application of the UNIQUAC Equation to Calculation of Multicomponent Phase Equilibria: 1. Vapor-Liquid Equilibria. Ind. Eng. Chem. Pro. Des. Dev. 1978, 17, 552. (34) Horst, R.; Pardalos, P. M. Handbook of Global Optimization; Kluwer Academic Publishers: New York, 1995. (35) Rosen, J. B. The Gradient Projection Method for Nonlinear Programming, Part I: Linear Constraints. SIAM J. Appl. Math. 1960, 8, 181. (36) Arce, A.; Blanco, A.; Blanco, M.; Soto, A.; Vidal, A. LiquidLiquid Equilibia of Water + Methanol + (MTBE or TAME) Mixtures. Can. J. Cem. Eng. 1994, 72, 935. (37) Sorensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection: Ternary Systems; DECHEMA Chemistry Data Series: Frankfurt, 1980.
(38) Chakravarty, T.; White, C. W., III; Seider, W. D. Computation of Phase Equilibrium: Optimization with Thermodynamics Inconsistency. AIChE J. 1985, 31, 316. (39) Null, H. R. Phase Equilibrium in Process Design; WileyInterscience: New York, 1970. (40) Gautam, R.; Seider, W. D. Computation of Phase and Chemical Equilibrium, Part 1: Local and Constrained Minima in Gibbs Free Energy. AIChE J. 1979, 25, 991. (41) Gautam, R.; Seider, W. D. Computation of Phase and Chemical Equilibrium, Part II: Phase-Splitting. AIChE J. 1979, 25, 999.
Received for review February 10, 1999 Revised manuscript received June 2, 1999 Accepted June 7, 1999 IE990104M
