- Email: [email protected]
Computation of multicomponent phase equilibria—Part III. Multiphase equilibria
Compurers and Chemical Engineering Vol. 8, No. 3/4. pp. 163-170, Printed in the U.S.A.
COMPUTATION EQUILIBRIA-PART
1984
009~1354/84 s3.00 + .oo Pergamon Press Ltd.
OF MULTICOMPONENT PHASE III. MULTIPHASE EQUILIBRIA
M. 0. OHANoMaHt and D. W.
THOMPSON$
Department of Chemical Engineering, The University of British Columbia, 2216 Main Mall, Vancouver, B.C., Canada V6T 1WS (Received 29 March 1983; revision received 7 November
1983)
Abstract-A number of algorithms have been developed for multiphase equilibria. The algorithms have been designed to effect a phase elimination without restarting the computation and they have been implemented to handle solid-liquid-liquid-vapour equilibria or any combination, thereof, of two or three phases including at least one liquid phase. The algorithms have been applied to a number of systems and their performance, based on iteration requirements, have been compared. Scope-This work has been directed towards determining the best method for effecting isothermal multicomponent phase equilibrium computation involving up to one vapour phase, two liquid phases and a solid phase. The algorithms are designed to be able to handle any of the following phase combinations: vapour-liquid, liquid-liquid, liquid-solid, vapour-liquidliquid, liquid-liquid-solid, vapour-liquid-solid and vapour-liquid-liquid-solid. A phase-fraction method, similar in check-function formulation to that proposed by Henley L Rosen [9] for vapour-liquid-liquid equilibria has been developed and variously accelerated by a Newton-Raphson approach, a quasi-Newton approach and a vector-projectionaccelerated successive-substitution method. A geometric programming algorithm has also been developed and subjected to the two vector projection methods outlined in the first paper of this series. A method based on the perturbation theory of geometric programming gave poor performance. In addition to the above-mentioned, two alternative formulations, utilizing stripping factors and also employing vector projection, have also been studied. Three initialization schemes have been studied and a “best” chosen from amongst them. All the algorithms have been programmed such that iteration is not restarted at the elimination of a non-existent phase. The ability of each algorithm to do this successfully has been tested by comparing its performance when redundant phases are assumed for a “known” system to when the right number of phases is initially assigned. In the applications, all liquid-phase and vapour-phase nonidealities are rigorously accounted for. The algorithms have not been applied to most of the possible multiphase combinations due to lack of appropriate equilibrium data. Conclusions and signilIcance-Where the number of phases initially assumed to exist is in excess of the number that actually exists, most of the algorithms encounter problems of convergence. This is the true test of their utility for, in practice, the actual number of phases that exists may not be known at any given condition of temperature and pressure. Everything considered, one of the stripping factor formulations, accelerated by tetrahedral projection, gives the best performance. For 16 applications based on the assumption of the right number of phases, it ranks best along with the tetrahedral-projection-accelerated geometricprogramming algorithm. When a wrong number of phases is assumed, it is clearly the best, having been found, in every application, to converge to the right solution in about the same number of iterations as when the right number of phases is assumed. INTRODUCTION
Some algorithms based on the phase-fraction formulation have been proposed for vapour-liquid-liquid equilibria by previous workers (see, for example, Osborne[lS]; Henley & Rosen[9]; Deam & Maddox [4]; Erbar [6]; Peng & Robinson [ 161; Mauri[l3]). Most of these algorithms use a twodimensional Newton-Raphson method of con-
vergence and they variously employ {L,, &}, {-G/J’, &IF}, {FL L&i} and (F/F, L,I(Li + G)) as check variables (all symbols are defined in a Nomenclature at the end of this report). The last form, due to Henley & Rosen[9], ensures clearly defined limits of 0.0 and 1.0. Three multiphase algorithms employing a similar phase-fraction formulation have been studied here.
tpresent address. The Department of Chemical Engineering, The University of Ife, Ile-Ife, Nigeria. SAuthor to whom correspondence should be addressed. 163
M. 0. OHANOMAHand D. W. THOMPSON
164
Free energy minimization algorithms exist for multiphase chemical equilibria (for example: White et al. [ 171; Clasen [3]; Mah & Shipman [ 121; Dluzniewski & Alder [S]; George et al. [8]; Gautam & Seider [7]). In the light of the work of Gautam & Seider[q, the BAND method as implemented by Dluzniewski & Adler[S] was adapted to the specific case of vapourliquid equilibria (see Part I of this series). In spite of modifications designed to improve its efficiency, it was found not to be competitive. It has therefore been excluded from this study. Outlines of the various methods investigated in this work are presented under Items 1 through 4 below. (1) Phare-fraction approach. From a combination of component- and total-mass balances and equilibrium relationships, it can be shown (Ohanomah[l4] for a vapour-liquid-liquid-solid system that xii = zi/{B,KUi+ 0,(1 - B”)K,i + &(l - e,)(l - e,)K, + (1 - &)( 1 - e,)( 1 - e”)} (i = 1,2, . . . ) N) (1) where the equilibrium
ratios are defined by
determined through an appropriate acceleration technique. (e) Any phase whose phase fraction is nonpositive is eliminated. (f) Convergence is tested for. The iteration is terminated in the event of a positive outcome; otherwise, control is transferred to Step (c). A double-loop version of the algorithm, involving an iterative implementation of Step (c) with the objective of correcting for the composition dependence of the K values, was found to perform poorly. The three acceleration methods studied in connection with the phase-fraction algorithm are outlined below. (a) A Newton-Ruphson approach. For a two-phase system, this reduces to the Newton method. For a three-phase system, a 2 x 2 Jacobian matrix is involved and the inverse is determined directly. A four-phase system involves the inverse of a 3 x 3 Jacobian matrix and a Gauss-Jordan elimination method (La Fara[l 11) is employed. (6) A quasi-Newton approach. This approach involves updating the phase-fraction vector, 8, according to tP+‘=tY+t”H”f(B”)
Yl = X2i =
Six/i
K,x,i
(3)
2
and xsi = K,x,.
V/F,
(5)
0, = L,/(F - V)
(6)
6, = S/(F - V - L2).
(7)
and
Define the check functions
(8)
h(e,, e,ye,)= i
Cxli
-
x2J
(9)
i=l
and
where f is the check-function vector and the parameter t is a real number. For a two-phase system, H” is given by H”
(4)
The phase fractions are defined by
e, =
(11)
(2)
A(e,, eb 0,) = i$, (xri - xsi)
(IO)
H” - “/(
=
1 _
p/f?
-
1’)
(12)
with H” = - l/f’(w). For more than two phases, the recurrence formula of Broyden[ 191 is employed. To determine t, a step-limited search, for a reduced value of f’f, is performed through a quadratic-fit minimization. (c) A partitioning method. This method updates 8 through a successive substitution that matches fI,, 0, and 0, with f,, f2and fJ respectively. Each successivesubstitution step utilizes the mean-value theorem (MVT) technique, ‘as documented in Part I of this series. Tetrahedral projection (see Part I) is then applied to these updated values for an improved value of 8. (2) Geometric programming (GP) formulation. Details of the mathematical development are presented by Ohanomah [ 141. For a vapour-liquid-liquid-solid system, the following 3N Eqs. result (i = 1,2,. . . , N): (1 - K,i)u,‘-[F-
t
(&j+~j)-K”i(Z;
- 12)- .ri)
j=l
Equations (I-10) constitute the working equations. A general algorithm has the following main steps: (a) Information is supplied as to which of the four phases are expected to exist (where there is doubt, the maximum possible number of phases is assumed). (b) 0,, 0, and 0, are appropriately initialized, and mole fractions are initialized for phases that have been assumed to exist. (c) Appropriate K values are- computed and the relevant ones amongst Eqs. (14) are solved. (d) The relevant check functions are computed and new values of the relevant phase fractions are
- (1 - &)&jUi + (1 - KJIzi - [F - f
KJ2j - I2i- Si)l,j =
0
(13)
(Vi+ Si>- KJZi - vi - SJ
j-1
- (1 - KJJ,Jl2i + K,kZi - U[- S&l/i= 0
(14)
and (1
-K,)sF-[F-_ (Uj+l2j)-K~~Zi-Ui-Z23 j=l
- (1 - &)1&i
+ K.+(Zi - ui - t’zi)A, = 0
(15)
M. 0. OHANOMAH and D. W. THOMPSON
166
on the basis of the hypothetical two-phase assumption and by further assuming that the fraction of each component in the hypothetical system relative to the actual system is the same for all components. For example, for Liquid Phase 2, we now assume a hypothetical liquid-liquid system with composition N) and liquid phase fraction fl, c$zi (i=1,2,... resulting in Iii = cclZi/[l + (1 -
e,)/e&]
$ = a&/[1 + (1 -
Ci =aZ/U
e,ye$J,
+ (1 - eJ/e&ilv
The steps necessary to implement this algorithm are listed below: (1) Set the phase fractions for phases that cannot exist to zero. Initialize the other phase fractions according to the schemes presented for two-phase systems in Part I (for 6,) and Part II (for 6, and 6,) of this series of papers. (2) Utilizing Eqs. (5)-(7), compute (where appropriate): B = LJF = (1 - 0,)(1 - 6,)(1 - 0,),
(22)
a, = (V + L,)/F = 0, + fl,
(23)
a, = (L, + LJF
(24)
and as = (S + L,)/F = O,(1 - 0,)( 1 - 0,) + p. (25)
Pressure Atm
(27)
$ = dZi/ [ 1 + (1 - 0&/0&J
(28)
IC = Zi( - Uf + 1ij + 3:)
(29)
(4) Where total phase moles need to be initialized, compute them, by summation, from Eqs. (26)-(29). Where phase mole fractions are also required, they are obtained from Eqs. (26)-(29) by division with the total phase moles. Stopping criterion. The convergence criterion employed was the same for all algorithms. In each case, convergence was said to have been attained if the mean absolute difference of component phase molar contents between two successive iterations, taken over all existent phases, was less than 0.0001. Property estimation. For computations involving only one liquid phase, the Wilson model for liquidphase activity coefficient (Wilson[ 181) was used. Where two liquid phases were involved, the Universal Quasichemical (UNIQUAC) model (Abrams and Prausnitz[l]) as modified by Anderson and Prausnitz[2] was employed. For all other thermodynamic
Table 1. Information on systems employed System
(26)
and
(i = 1,2, . . . N).
= 0X1- 0,) + /3,
(3) Compute (where appropriate):
Composition (mole X)
A
1.0
Water (27.217); Ethanol (22.783); Ethyl acetate (27.217); Acetic acid (22.783).
B
1.0
Ethanol (10.0); Chloroform N-Hexane (20.0).
C
1.0
N-Hexane (18.6); Methylcyclopentane (20.7); Toluene (28.6); Cyclohexane (11.8); Benzene (15.3); Ethanol (5.0).
II
1.0
Cyclohexane (50.0); Eicosane (50.0).
E
1.0
Cyclohexane (17.0); Hexadecane (49.0); Tetrachloromethane (17.0); Benzene (17.0).
F
1.0
Water (40.0); Methanol (20.0); Benzene (40.0).
G
1.0
Butyl alcohol (20.2); Water (74.18); Propyl alcohol (5.62).
H
1.0
Water (23.5); Benzene (57.5); Ethanol (19.0).
I
0.980
Methanol (3.38); Water (77.19); Butanol (19.43).
J
1.009
Ethanol (3.17); Water (76.72); Butanol (20.11).
K
0.221
Nitrogen (43.17); Argon (56.83).
L
0.2105
Nitrogen (50.0); Methane (50.0).
(20.0);
Acetone (50.0);
Computation of multicomponent phase equilibria-Part dynamic properties, the methods outlined in Part I (for vapour and liquid phase properties) and Part II (for solid-phase properties) were employed. Applications. In implementing the algorithms, allowance was made for the possibility of L\ reducing to zero while & is greater than zero. In such an event, Liquid Phase 1 is eliminated and Liquid Phase 2 is redefined as Liquid Phase 1, the base liquid phase. Due to lack of information, applications have been restricted to vapour-liquid-liquid multiphase combination. Table 1 contains information on the systems employed. The thermodynamic data are the same as those used in Parts I and II. The applications are in two forms. One form involves assuming the right number of phases known to exist at equilibrium. The other form is based on assuming more phases to be present than actually exist. The results are presented in Tables 2 and 3. (The symbols “S’, “L” and “V” under “Phases” and “vapour” respectively.) denote “solid”, “ hquid” . COMMENTS
III
167
projection and the algorithms using the alternative (stripping factor) formulations give the best results. The quasi-Newton phase-fraction algorithm is unreliable. The partitioned phase-fraction algorithm also performs poorly for multiphase systems. This must be the effect of nonidealities, for when this algorithm (and the same applies to the sensitivity algorithm which has been accorded only a passing mention in this paper) as well as the Newton-Raphson phase-fraction algorithm and the GP algorithm were tested with the sample quarternary liquid-liquid-vapour problem which was solved by Osborne[lS] and by Deam & Maddox[4]-the problem is based on constant K values-all the algorithms converged speedily to the right solution. For applications involving the assumption of a wrong number of phases, the results (Table 3) show that for liquid-liquid systems where it is known with absolute certainty that solid formation cannot occur but where vaporization could take place, the approach whereby the solution procedure is not restarted after a phase elimination seems to be alright with most of the methods. But where, in addition to vapour and liquid phases, the formation of a solid phase is possible, the algorithm employing “Altemative Formulation 1” appears to be the only reliable one.
ON RESULTS
The results for cases when the right number of phases was assumed (Table 2) show that none of the algorithms excels for all systems. However, on the average, the GP algorithm with tetrahedral
Table 2. Iteration counts when the right number of phases is assumed Phase-fraction algorithms System
Temp (K)
Phases
A
356.0
B
GP algorithms
NewtonRaphson
QuasiNewton
Partitioned Triangular Tetrahedral Projection Projection
LV
7
5
7
6
329.5
LV
9
6
9
C
353.0
LV
5
5
D
515.5
LV
5
5
E
470.0
LV
4
F
298.0
LL
G
303.0
LL
H
328.0
Alternative formulations 1
2
6
7
7
9
LO
8
7
5
6
7
7
7
4
2
3
4
4
4
3
3
5
4
5
21
fc
22
9
7
7
8
33
13
32
13
7
LO
8
LL
54
f
54
8
10
9
11
I
330.0
LL
36
fc
35
14
8
9
8
J
330.0
LL
38
fc
38
19
9
11
9
K
78.7
SL
5
6
6
5
5
6
SL
14
fC
14
13
6
a
7
L
92.66
F
335.50
LLV
16
f
f
17
16
13
21
G
364.60
LLV
42
f
f
58
66
44
35
I
364.0
LLV
36
f
42
38
21
29
33
J
365.0
LLV
41
fc
f
38
22
36
32
f = failure fc = faulty convereence
M. 0.
168
and D. W. THOMPSON
OHANOMAH
Table 3. Iteration counts when redundant phases are assumed Phase-fraction algorithms
Phases System
Temp (K)
Expected
Assumed
NeWtOllRaphson
QuasiNewton
Partitioned
GP algorithms
Triangular Projection
Tetrahedral Projection
Alternative formulations 1
2
F
298.0
LL
LLV
28
fc
23
9
8
8
9
G
303.0
LL
LLV
33
fc
32
10
9
9
8
J
330.0
LL
LLV
38
fc
38
14
10
10
10
A
356.0
LV
SLV
9
8
8
9
7
8
9
fc
6
9
9
8
8
5
6
7
5
5
fc
fc
11
fc
C
353.0
LV
SLV
8
D
515.5
LV
SLV
10
fc
F
298.0
LL
SLLV
31
f
PC
J
330.0
LL
SLLV
42
fc
39
16
11
11
11
F
335.5
LLV
SLLV
fc
f
fc
fc
fc
26
28
J
365.0
LLV
SLLV
fc
f
fc
47
41
36
35
Average iteration increase for assuming LLV for LL systems.
2.3
0.3
-2.7
1.3
-0.3
0.7
Average iteration increase for assuming SLV for LV systems
3.3
1.0
3.3
2.3
1.0
1.3
f = failure fc - faulty convergence
NOMENCLATURE
stripping factor check function vector total number of moles in the system negative of the estimate of inverse of Newton-Raphson Jacobian matrix equilibrium ratio moles of a given component in a given liquid phase total moles in a given liquid phase number of components in the system moles of a given component in the solid phase total moles in the solid phase unidimensional search parameter [Eq. (11 )I moles of a given component in the vapour phase total moles in the vapour phase mole fraction of a given component in a given nonvapour phase mole fraction of a given component in the vapour phase mole fraction of a given component in the system total moles of a given component in the system Greek symbols 01 parameters defined by Eqs. (23)-(25)
B parameters defined by Eqs. (21) and (22) 1 parameters defined by Eqs. (16) and (20) 0 phase fraction for a given phase 0 vector of phase fractions Subscripts i component
I s u 1 2
liquid phase solid phase vapour phase liquid phase 1 liquid phase 2
Superscripts n iteration
0 T
standard state ; initial point uanspose derivative
REFERENCES 1. D. S. Abrams, & J. M. Prausnitz, ‘Statistical thermodynamics of liquid mixtures. A new expression for the excess Gibbs energy of partly or completely miscible systems, AZChE J. 21(l), 116 (1975). 2. T. F. Anderson & J. M. Prausnitx, Application of the UNIQUAC equation to calculation of multicomponent phase equilibria-Part I: Vapour-liquid equilibria, Ind. Engng Chem. Process Des. Develop. 17(4), 552 (1978). 3. R. J. Clasen, The numerical solution of the chemical equilibrium problem, in The RAND Corp. RM-4345 PR, January 1965). 4. J. R. Deam & R. N. Maddox, How to figure three-phase flash, Hydrocarbon Processing July, 163 ( 1969). 5. J. H. Dlumiewski t S. B. Adler, Calculation of complex reaction and/or phase equilibria problems, I. Chem. E. Symp. Series 35, 4 (1972). 6. J. H. Erbar, Three phase equilibrium calculation, Natural Gas Processors Assoc. 52, 62 (1973). 7. R. Gautam & W. D. Seider, Computation of phase and chemical equilibria-Part I: Local and constrained minima in Gibbs free energy, AZChE J. u(6), 991 (1979). 8. B. George, L. T. Brown, C. H. Farmer, P. Buthod & F. S. Menning, Computation of multicomponent, multiphase equilibrium, Znd. Engng Chem. Process Des. Develop., lS(3). 372 (1976). 9. E. J. Henley Sr E. M. Rosen, Mareriai and Energy Balance Computations, Wiley, New York (1969). 10. C. D. Holland Fundamentals of Multicomponent Distillation. McGraw-Hill, New York (198 1).
It phase equilibria-Part 11. R. L. La Fara, Computer Methods for Science and Engineering. Hayden Book Co. New Jersey (1973). 12. Y. H. Mah SCC. W. Shipman, On the computation of complex equilibria, AIChE .I. M(2), 299 (1972). 13. C. Mauri, Unified procedure for solving multiphasemulticomponent vapour-liquid equilibrium calculation, Ind. Engng Chem. Process Des. Develop. 19(3), 482 (1980). 14. M. 0. Ohanomah, Computational algorithms for multicomwnent uhase eouihbria and distillation, Ph.D. disset&ion, U&v. of British Columbia, Vancouver (1981). IS. A. Osborne, How to calculate three-phase flash vaporization, Chem. Engng. 21, 97 (1964).
III
169
16. Ding-Yu Peng & D. B. Robinson, Two and three-phase equiiibxium calculations for systems containing water, Ind. Enana Chem. Funabnentals. l!Wh 59 11976). 17. W. B. Whyte, S. M. Johnson & Gl B. ‘D&zig, Chemical equilibrium in complex mixtures, J. Chem. Phys. B(5), 751 (1958). 18. G. M. Wilson, A new expression for the excess free energy of mixing, J. Amer. Chem. Sot. 86, (1964). Reference added in proof 19. C. G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Computat. 19, 577 (1965).
COMPUTATION EQUILIBRIA-PART
1984
009~1354/84 s3.00 + .oo Pergamon Press Ltd.
OF MULTICOMPONENT PHASE III. MULTIPHASE EQUILIBRIA
M. 0. OHANoMaHt and D. W.
THOMPSON$
Department of Chemical Engineering, The University of British Columbia, 2216 Main Mall, Vancouver, B.C., Canada V6T 1WS (Received 29 March 1983; revision received 7 November
1983)
Abstract-A number of algorithms have been developed for multiphase equilibria. The algorithms have been designed to effect a phase elimination without restarting the computation and they have been implemented to handle solid-liquid-liquid-vapour equilibria or any combination, thereof, of two or three phases including at least one liquid phase. The algorithms have been applied to a number of systems and their performance, based on iteration requirements, have been compared. Scope-This work has been directed towards determining the best method for effecting isothermal multicomponent phase equilibrium computation involving up to one vapour phase, two liquid phases and a solid phase. The algorithms are designed to be able to handle any of the following phase combinations: vapour-liquid, liquid-liquid, liquid-solid, vapour-liquidliquid, liquid-liquid-solid, vapour-liquid-solid and vapour-liquid-liquid-solid. A phase-fraction method, similar in check-function formulation to that proposed by Henley L Rosen [9] for vapour-liquid-liquid equilibria has been developed and variously accelerated by a Newton-Raphson approach, a quasi-Newton approach and a vector-projectionaccelerated successive-substitution method. A geometric programming algorithm has also been developed and subjected to the two vector projection methods outlined in the first paper of this series. A method based on the perturbation theory of geometric programming gave poor performance. In addition to the above-mentioned, two alternative formulations, utilizing stripping factors and also employing vector projection, have also been studied. Three initialization schemes have been studied and a “best” chosen from amongst them. All the algorithms have been programmed such that iteration is not restarted at the elimination of a non-existent phase. The ability of each algorithm to do this successfully has been tested by comparing its performance when redundant phases are assumed for a “known” system to when the right number of phases is initially assigned. In the applications, all liquid-phase and vapour-phase nonidealities are rigorously accounted for. The algorithms have not been applied to most of the possible multiphase combinations due to lack of appropriate equilibrium data. Conclusions and signilIcance-Where the number of phases initially assumed to exist is in excess of the number that actually exists, most of the algorithms encounter problems of convergence. This is the true test of their utility for, in practice, the actual number of phases that exists may not be known at any given condition of temperature and pressure. Everything considered, one of the stripping factor formulations, accelerated by tetrahedral projection, gives the best performance. For 16 applications based on the assumption of the right number of phases, it ranks best along with the tetrahedral-projection-accelerated geometricprogramming algorithm. When a wrong number of phases is assumed, it is clearly the best, having been found, in every application, to converge to the right solution in about the same number of iterations as when the right number of phases is assumed. INTRODUCTION
Some algorithms based on the phase-fraction formulation have been proposed for vapour-liquid-liquid equilibria by previous workers (see, for example, Osborne[lS]; Henley & Rosen[9]; Deam & Maddox [4]; Erbar [6]; Peng & Robinson [ 161; Mauri[l3]). Most of these algorithms use a twodimensional Newton-Raphson method of con-
vergence and they variously employ {L,, &}, {-G/J’, &IF}, {FL L&i} and (F/F, L,I(Li + G)) as check variables (all symbols are defined in a Nomenclature at the end of this report). The last form, due to Henley & Rosen[9], ensures clearly defined limits of 0.0 and 1.0. Three multiphase algorithms employing a similar phase-fraction formulation have been studied here.
tpresent address. The Department of Chemical Engineering, The University of Ife, Ile-Ife, Nigeria. SAuthor to whom correspondence should be addressed. 163
M. 0. OHANOMAHand D. W. THOMPSON
164
Free energy minimization algorithms exist for multiphase chemical equilibria (for example: White et al. [ 171; Clasen [3]; Mah & Shipman [ 121; Dluzniewski & Alder [S]; George et al. [8]; Gautam & Seider [7]). In the light of the work of Gautam & Seider[q, the BAND method as implemented by Dluzniewski & Adler[S] was adapted to the specific case of vapourliquid equilibria (see Part I of this series). In spite of modifications designed to improve its efficiency, it was found not to be competitive. It has therefore been excluded from this study. Outlines of the various methods investigated in this work are presented under Items 1 through 4 below. (1) Phare-fraction approach. From a combination of component- and total-mass balances and equilibrium relationships, it can be shown (Ohanomah[l4] for a vapour-liquid-liquid-solid system that xii = zi/{B,KUi+ 0,(1 - B”)K,i + &(l - e,)(l - e,)K, + (1 - &)( 1 - e,)( 1 - e”)} (i = 1,2, . . . ) N) (1) where the equilibrium
ratios are defined by
determined through an appropriate acceleration technique. (e) Any phase whose phase fraction is nonpositive is eliminated. (f) Convergence is tested for. The iteration is terminated in the event of a positive outcome; otherwise, control is transferred to Step (c). A double-loop version of the algorithm, involving an iterative implementation of Step (c) with the objective of correcting for the composition dependence of the K values, was found to perform poorly. The three acceleration methods studied in connection with the phase-fraction algorithm are outlined below. (a) A Newton-Ruphson approach. For a two-phase system, this reduces to the Newton method. For a three-phase system, a 2 x 2 Jacobian matrix is involved and the inverse is determined directly. A four-phase system involves the inverse of a 3 x 3 Jacobian matrix and a Gauss-Jordan elimination method (La Fara[l 11) is employed. (6) A quasi-Newton approach. This approach involves updating the phase-fraction vector, 8, according to tP+‘=tY+t”H”f(B”)
Yl = X2i =
Six/i
K,x,i
(3)
2
and xsi = K,x,.
V/F,
(5)
0, = L,/(F - V)
(6)
6, = S/(F - V - L2).
(7)
and
Define the check functions
(8)
h(e,, e,ye,)= i
Cxli
-
x2J
(9)
i=l
and
where f is the check-function vector and the parameter t is a real number. For a two-phase system, H” is given by H”
(4)
The phase fractions are defined by
e, =
(11)
(2)
A(e,, eb 0,) = i$, (xri - xsi)
(IO)
H” - “/(
=
1 _
p/f?
-
1’)
(12)
with H” = - l/f’(w). For more than two phases, the recurrence formula of Broyden[ 191 is employed. To determine t, a step-limited search, for a reduced value of f’f, is performed through a quadratic-fit minimization. (c) A partitioning method. This method updates 8 through a successive substitution that matches fI,, 0, and 0, with f,, f2and fJ respectively. Each successivesubstitution step utilizes the mean-value theorem (MVT) technique, ‘as documented in Part I of this series. Tetrahedral projection (see Part I) is then applied to these updated values for an improved value of 8. (2) Geometric programming (GP) formulation. Details of the mathematical development are presented by Ohanomah [ 141. For a vapour-liquid-liquid-solid system, the following 3N Eqs. result (i = 1,2,. . . , N): (1 - K,i)u,‘-[F-
t
(&j+~j)-K”i(Z;
- 12)- .ri)
j=l
Equations (I-10) constitute the working equations. A general algorithm has the following main steps: (a) Information is supplied as to which of the four phases are expected to exist (where there is doubt, the maximum possible number of phases is assumed). (b) 0,, 0, and 0, are appropriately initialized, and mole fractions are initialized for phases that have been assumed to exist. (c) Appropriate K values are- computed and the relevant ones amongst Eqs. (14) are solved. (d) The relevant check functions are computed and new values of the relevant phase fractions are
- (1 - &)&jUi + (1 - KJIzi - [F - f
KJ2j - I2i- Si)l,j =
0
(13)
(Vi+ Si>- KJZi - vi - SJ
j-1
- (1 - KJJ,Jl2i + K,kZi - U[- S&l/i= 0
(14)
and (1
-K,)sF-[F-_ (Uj+l2j)-K~~Zi-Ui-Z23 j=l
- (1 - &)1&i
+ K.+(Zi - ui - t’zi)A, = 0
(15)
M. 0. OHANOMAH and D. W. THOMPSON
166
on the basis of the hypothetical two-phase assumption and by further assuming that the fraction of each component in the hypothetical system relative to the actual system is the same for all components. For example, for Liquid Phase 2, we now assume a hypothetical liquid-liquid system with composition N) and liquid phase fraction fl, c$zi (i=1,2,... resulting in Iii = cclZi/[l + (1 -
e,)/e&]
$ = a&/[1 + (1 -
Ci =aZ/U
e,ye$J,
+ (1 - eJ/e&ilv
The steps necessary to implement this algorithm are listed below: (1) Set the phase fractions for phases that cannot exist to zero. Initialize the other phase fractions according to the schemes presented for two-phase systems in Part I (for 6,) and Part II (for 6, and 6,) of this series of papers. (2) Utilizing Eqs. (5)-(7), compute (where appropriate): B = LJF = (1 - 0,)(1 - 6,)(1 - 0,),
(22)
a, = (V + L,)/F = 0, + fl,
(23)
a, = (L, + LJF
(24)
and as = (S + L,)/F = O,(1 - 0,)( 1 - 0,) + p. (25)
Pressure Atm
(27)
$ = dZi/ [ 1 + (1 - 0&/0&J
(28)
IC = Zi( - Uf + 1ij + 3:)
(29)
(4) Where total phase moles need to be initialized, compute them, by summation, from Eqs. (26)-(29). Where phase mole fractions are also required, they are obtained from Eqs. (26)-(29) by division with the total phase moles. Stopping criterion. The convergence criterion employed was the same for all algorithms. In each case, convergence was said to have been attained if the mean absolute difference of component phase molar contents between two successive iterations, taken over all existent phases, was less than 0.0001. Property estimation. For computations involving only one liquid phase, the Wilson model for liquidphase activity coefficient (Wilson[ 181) was used. Where two liquid phases were involved, the Universal Quasichemical (UNIQUAC) model (Abrams and Prausnitz[l]) as modified by Anderson and Prausnitz[2] was employed. For all other thermodynamic
Table 1. Information on systems employed System
(26)
and
(i = 1,2, . . . N).
= 0X1- 0,) + /3,
(3) Compute (where appropriate):
Composition (mole X)
A
1.0
Water (27.217); Ethanol (22.783); Ethyl acetate (27.217); Acetic acid (22.783).
B
1.0
Ethanol (10.0); Chloroform N-Hexane (20.0).
C
1.0
N-Hexane (18.6); Methylcyclopentane (20.7); Toluene (28.6); Cyclohexane (11.8); Benzene (15.3); Ethanol (5.0).
II
1.0
Cyclohexane (50.0); Eicosane (50.0).
E
1.0
Cyclohexane (17.0); Hexadecane (49.0); Tetrachloromethane (17.0); Benzene (17.0).
F
1.0
Water (40.0); Methanol (20.0); Benzene (40.0).
G
1.0
Butyl alcohol (20.2); Water (74.18); Propyl alcohol (5.62).
H
1.0
Water (23.5); Benzene (57.5); Ethanol (19.0).
I
0.980
Methanol (3.38); Water (77.19); Butanol (19.43).
J
1.009
Ethanol (3.17); Water (76.72); Butanol (20.11).
K
0.221
Nitrogen (43.17); Argon (56.83).
L
0.2105
Nitrogen (50.0); Methane (50.0).
(20.0);
Acetone (50.0);
Computation of multicomponent phase equilibria-Part dynamic properties, the methods outlined in Part I (for vapour and liquid phase properties) and Part II (for solid-phase properties) were employed. Applications. In implementing the algorithms, allowance was made for the possibility of L\ reducing to zero while & is greater than zero. In such an event, Liquid Phase 1 is eliminated and Liquid Phase 2 is redefined as Liquid Phase 1, the base liquid phase. Due to lack of information, applications have been restricted to vapour-liquid-liquid multiphase combination. Table 1 contains information on the systems employed. The thermodynamic data are the same as those used in Parts I and II. The applications are in two forms. One form involves assuming the right number of phases known to exist at equilibrium. The other form is based on assuming more phases to be present than actually exist. The results are presented in Tables 2 and 3. (The symbols “S’, “L” and “V” under “Phases” and “vapour” respectively.) denote “solid”, “ hquid” . COMMENTS
III
167
projection and the algorithms using the alternative (stripping factor) formulations give the best results. The quasi-Newton phase-fraction algorithm is unreliable. The partitioned phase-fraction algorithm also performs poorly for multiphase systems. This must be the effect of nonidealities, for when this algorithm (and the same applies to the sensitivity algorithm which has been accorded only a passing mention in this paper) as well as the Newton-Raphson phase-fraction algorithm and the GP algorithm were tested with the sample quarternary liquid-liquid-vapour problem which was solved by Osborne[lS] and by Deam & Maddox[4]-the problem is based on constant K values-all the algorithms converged speedily to the right solution. For applications involving the assumption of a wrong number of phases, the results (Table 3) show that for liquid-liquid systems where it is known with absolute certainty that solid formation cannot occur but where vaporization could take place, the approach whereby the solution procedure is not restarted after a phase elimination seems to be alright with most of the methods. But where, in addition to vapour and liquid phases, the formation of a solid phase is possible, the algorithm employing “Altemative Formulation 1” appears to be the only reliable one.
ON RESULTS
The results for cases when the right number of phases was assumed (Table 2) show that none of the algorithms excels for all systems. However, on the average, the GP algorithm with tetrahedral
Table 2. Iteration counts when the right number of phases is assumed Phase-fraction algorithms System
Temp (K)
Phases
A
356.0
B
GP algorithms
NewtonRaphson
QuasiNewton
Partitioned Triangular Tetrahedral Projection Projection
LV
7
5
7
6
329.5
LV
9
6
9
C
353.0
LV
5
5
D
515.5
LV
5
5
E
470.0
LV
4
F
298.0
LL
G
303.0
LL
H
328.0
Alternative formulations 1
2
6
7
7
9
LO
8
7
5
6
7
7
7
4
2
3
4
4
4
3
3
5
4
5
21
fc
22
9
7
7
8
33
13
32
13
7
LO
8
LL
54
f
54
8
10
9
11
I
330.0
LL
36
fc
35
14
8
9
8
J
330.0
LL
38
fc
38
19
9
11
9
K
78.7
SL
5
6
6
5
5
6
SL
14
fC
14
13
6
a
7
L
92.66
F
335.50
LLV
16
f
f
17
16
13
21
G
364.60
LLV
42
f
f
58
66
44
35
I
364.0
LLV
36
f
42
38
21
29
33
J
365.0
LLV
41
fc
f
38
22
36
32
f = failure fc = faulty convereence
M. 0.
168
and D. W. THOMPSON
OHANOMAH
Table 3. Iteration counts when redundant phases are assumed Phase-fraction algorithms
Phases System
Temp (K)
Expected
Assumed
NeWtOllRaphson
QuasiNewton
Partitioned
GP algorithms
Triangular Projection
Tetrahedral Projection
Alternative formulations 1
2
F
298.0
LL
LLV
28
fc
23
9
8
8
9
G
303.0
LL
LLV
33
fc
32
10
9
9
8
J
330.0
LL
LLV
38
fc
38
14
10
10
10
A
356.0
LV
SLV
9
8
8
9
7
8
9
fc
6
9
9
8
8
5
6
7
5
5
fc
fc
11
fc
C
353.0
LV
SLV
8
D
515.5
LV
SLV
10
fc
F
298.0
LL
SLLV
31
f
PC
J
330.0
LL
SLLV
42
fc
39
16
11
11
11
F
335.5
LLV
SLLV
fc
f
fc
fc
fc
26
28
J
365.0
LLV
SLLV
fc
f
fc
47
41
36
35
Average iteration increase for assuming LLV for LL systems.
2.3
0.3
-2.7
1.3
-0.3
0.7
Average iteration increase for assuming SLV for LV systems
3.3
1.0
3.3
2.3
1.0
1.3
f = failure fc - faulty convergence
NOMENCLATURE
stripping factor check function vector total number of moles in the system negative of the estimate of inverse of Newton-Raphson Jacobian matrix equilibrium ratio moles of a given component in a given liquid phase total moles in a given liquid phase number of components in the system moles of a given component in the solid phase total moles in the solid phase unidimensional search parameter [Eq. (11 )I moles of a given component in the vapour phase total moles in the vapour phase mole fraction of a given component in a given nonvapour phase mole fraction of a given component in the vapour phase mole fraction of a given component in the system total moles of a given component in the system Greek symbols 01 parameters defined by Eqs. (23)-(25)
B parameters defined by Eqs. (21) and (22) 1 parameters defined by Eqs. (16) and (20) 0 phase fraction for a given phase 0 vector of phase fractions Subscripts i component
I s u 1 2
liquid phase solid phase vapour phase liquid phase 1 liquid phase 2
Superscripts n iteration
0 T
standard state ; initial point uanspose derivative
REFERENCES 1. D. S. Abrams, & J. M. Prausnitz, ‘Statistical thermodynamics of liquid mixtures. A new expression for the excess Gibbs energy of partly or completely miscible systems, AZChE J. 21(l), 116 (1975). 2. T. F. Anderson & J. M. Prausnitx, Application of the UNIQUAC equation to calculation of multicomponent phase equilibria-Part I: Vapour-liquid equilibria, Ind. Engng Chem. Process Des. Develop. 17(4), 552 (1978). 3. R. J. Clasen, The numerical solution of the chemical equilibrium problem, in The RAND Corp. RM-4345 PR, January 1965). 4. J. R. Deam & R. N. Maddox, How to figure three-phase flash, Hydrocarbon Processing July, 163 ( 1969). 5. J. H. Dlumiewski t S. B. Adler, Calculation of complex reaction and/or phase equilibria problems, I. Chem. E. Symp. Series 35, 4 (1972). 6. J. H. Erbar, Three phase equilibrium calculation, Natural Gas Processors Assoc. 52, 62 (1973). 7. R. Gautam & W. D. Seider, Computation of phase and chemical equilibria-Part I: Local and constrained minima in Gibbs free energy, AZChE J. u(6), 991 (1979). 8. B. George, L. T. Brown, C. H. Farmer, P. Buthod & F. S. Menning, Computation of multicomponent, multiphase equilibrium, Znd. Engng Chem. Process Des. Develop., lS(3). 372 (1976). 9. E. J. Henley Sr E. M. Rosen, Mareriai and Energy Balance Computations, Wiley, New York (1969). 10. C. D. Holland Fundamentals of Multicomponent Distillation. McGraw-Hill, New York (198 1).
It phase equilibria-Part 11. R. L. La Fara, Computer Methods for Science and Engineering. Hayden Book Co. New Jersey (1973). 12. Y. H. Mah SCC. W. Shipman, On the computation of complex equilibria, AIChE .I. M(2), 299 (1972). 13. C. Mauri, Unified procedure for solving multiphasemulticomponent vapour-liquid equilibrium calculation, Ind. Engng Chem. Process Des. Develop. 19(3), 482 (1980). 14. M. 0. Ohanomah, Computational algorithms for multicomwnent uhase eouihbria and distillation, Ph.D. disset&ion, U&v. of British Columbia, Vancouver (1981). IS. A. Osborne, How to calculate three-phase flash vaporization, Chem. Engng. 21, 97 (1964).
III
169
16. Ding-Yu Peng & D. B. Robinson, Two and three-phase equiiibxium calculations for systems containing water, Ind. Enana Chem. Funabnentals. l!Wh 59 11976). 17. W. B. Whyte, S. M. Johnson & Gl B. ‘D&zig, Chemical equilibrium in complex mixtures, J. Chem. Phys. B(5), 751 (1958). 18. G. M. Wilson, A new expression for the excess free energy of mixing, J. Amer. Chem. Sot. 86, (1964). Reference added in proof 19. C. G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Computat. 19, 577 (1965).