Computation of multicomponent phase equilibria—Part I. Vapour-liquid equilibria

Computers and Chemical Printed in the U.S.A.

Engineering

0098-1354/84 $3.00 + .oo Pergamon Press Ltd.

Vol. 8, No. 314. pp. 147-156.1984

COMPUTATION EQUILIBRIA-PART

OF MULTICOMPONENT PHASE I. VAPOUR-LIQUID EQUILIBRIA

M. 0. OHANOMAH~and D. W. THOMP~~N~ Department of Chemical Engineering, The University of British Columbia, 2216 Main Mall,

Vancouver, B. C., Canada V6T IWS (Received 29 March 1983; revision received I November 1983)

Abstract--A number of algorithms have been developed for vapour-liquid equilibrium using either the mass-balance (univariate) approach or the free-energy-minimization (multivariate) method. In the former case, both double-loop and single-loop formulations have been studied; in the latter, the technique of geometric programming has been employed. The best of the methods described in the literature have also been implemented. Each algorithm accounts rigorously for nonidealities in both phases and is programmed to use the best of a set of initialization schemes. The various algorithms have been compared, using computation (CPU) time as the measure of performance. Scope-The main objective of this work has been to determine the best method of solving the isothermal vapour-liquid equilibrium problem, without any restrictions on the degree of nonideality of the two phases involved. To achieve this goal, a number of formulations of the mass-balance check function have been introduced and these, studied along with the standard formulation and the formulations of Rachford & Rice[24] and of Barnes & Flores (King,[l4]). Single-loop algorithms as well as double-loop algorithms, which interactively correct for the dependence of equilibrium ratios on phase compositions in an inner loop, have been investigated. A scheme that updates the equilibrium ratios less frequently than in the single-loop algorithms has also been explored. Drawing on the conclusions of Rohl & Suda11[26], the Newton and Richmond acceleration techniques have been employed. In addition, the mean value theorem of differential calculus, the Wegstein projection method and a quadratic form of the Wegstein method have been experimented with. Employing the free-energy-minimization approach, a geometric programming transformation has been used to develop a set of nonlinear equations which have been subjected to a number of solution methods, including the Newton-Raphson method and the method of successive substitution, accelerated by either the “hyperplane linearization” technique of Barreto & Farina[l] or the vector projection methods proposed by Ohanomah[l9]. An algorithm based on the perturbation theory of geometric programming has also been developed. Guided by the comparative study of Gautam 8c Seider[S], the RAND algorithm has been implemented for the specific case of vapour-liquid equilibria. A modification of the RAND method has also been studied. The various algorithms have been implemented based on well established predictive correlations for thermodynamic properties. They have been compared, with computation (CPU) time as the measure of performance. Conelusious and SignIiIeauee-The following conclusions are based on the application of each algorithm to 304 flash calculations spread over the vapour-liquid region for 16 mixtures: (1) The double-loop univariate algorithms are generally comparatively inefficient. (2) The speed, stability and reliability of the computational algorithms are highly dependent on the mode of initialization. With proper initialization, even highly nonideal mixtures can be handled by a single-loop computational algorithm. The initialization scheme proposed here is quite reliable and has been found to yield average computation-time reduction factors of 1.4-1.9, depending on the algorithm, compared to some conventional schemes. (3) The logarithmic forms of the mass-balance check functions have been found to perform slightly worse than their nonlogarithmic equivalents. (4) With a simple terminal initialization scheme (vapour fraction equal to 0.0 or 1.0, depending on the check-function type), the performances of the single-loop univariate algorithms are found to depend very much on the acceleration technique. However, with the proposed initialization scheme, the differences are insignificant. (5) Computing K-values less frequently than “once every iteration” could have a significant effect on the tPresent address: The Department of Chemical Engineering, The University of Ife, Ile.-Ife, Nigeria. SAuthor to whom correspondence should be addressed.

147

M. 0. OHANOMAH and D. W. THOMFWJN

148

quality of results. Where a relative deviation of 0.001 in the phase fraction is tolerable, then computing K-values on alternate iterations should lead to satisfactory results. This produces an average time-saving of about 12% for the single-loop univariate algorithms, and much less (6) The best of the multivariate algorithms for the multivariate algorithms. (sensitivity + tetrahedral projection) and the best univariate algorithm (Rachford-Rice + Wegstein + MVT) give comparable performances. However, the latter is simpler, and also distinctly better when K-values are computed on alternate iterations. INTRODUCrION Situations arise quite frequently in chemical process design or operation that call for the determination of the equilibrium distribution of the components of a multicomponent mixture between a vapour and a liquid phase at isobaric, isothermal conditions. The solution of the problem lies in either the so-called mass-balance approach or in the minimization of the Gibbs free energy of the system. The mass-balance method has experienced various formulations and has been subjected to a number of acceleration methods (see, for example, Rachford & Rice[24]; Rohl & Sudall[26]; Sanderson Jr Chien[27]; Von Rosenberg[28,29]; Holland[ 1I]; Hirose et al. [ IO]; King[l4]. The free energy minimization algorithms have generally been designed for chemical equilibria and have been applied to physical equilibria merely by virtue of the latter being a special case of the former (see, for example, White et al. [32]; Clasen[2]; Ma & Shipman [ 171; Dluzniewski & Adler [31; George et al.[6]; Gautam & Seider,[S]. The various methods studied in this work are briefly discussed below, under three headings labelled (A) through (C). (A) Double-loop univariate method The univariate methods, deriving from the massbalance formulation, treat the vapour fraction, 6,, as the independent variable and define the liquid-phase mole fraction, x, by

(b) Form 2 (“Rachford-Rice”). Rice[24] proposed the check function

Rachford

CYi -x3. i-1

f(k) = f

(i=l,2,...,N)

which is reportedly superior to that of Rachford & Rice[24]. (d) Form 4 (“Log Standard”). This is of the form

fw

= In 15

J.

It derives from Form 1, drawing on the relationship between Forms 2 and 3. (e) Form 5 (“Alternative Form A”). By manipulating the mass-balance and equilibrium relationships, it can be shown that Forms 1 and 2 are related by

(8)

(1)

where N denotes the number of components, zi is-the mole fraction of component i in the system and, in this work, the equilibrium ratio, Ki, is related to system pressure, P, liquid-phase activity coefficient, ye liquid-phase reference fugacity, fi”, and vapourphase fugacity coefficient, & by Ki = yifialbiP. The vapour-phase

(2)

mole fraction is given by yi = Kix,.

f xir=,

In view of this, a study was undertaken function of the form

1.

[i=li1

(9)

where t is some positive real number. Various schemes for estimating the best value for t were studied and it was found that setting t =0.5 was about as good as any more complicated scheme. Thus, the check function was adopted in the form

[i=li1 (10)

f(e,) =eu-o,5 2 xi-

(f) Form 6 (“Alternative Form B”). This is based on the same reasoning as Form 5 and the check function is of the general form

(4)

i=l $xi-l

f(e,)=[a+(l-a)&-I]

Rohl & Sudall[26] dismissed this form as not being capable of converging. As rightly asserted by Holland[ll], it will always converge if 6, is suitably initialized.

of a check

f(e,) = en-li xj-

(3)

The following six formulations have been subjected to the double-loop treatment: (a) Form 1 (“Stundard”). This employs the check function f(e”)=

(5)

This form was favoured by the results of the study of Rohl Jr Sudall[26]. (c) Form 3 (“Barnes-Flares”). According to King[l4], Barnes and Flores have proposed the check function

.

xi = Zi/{1 + (Ki - 1)e”)

&

(

(11) >

where a is a real number. Based on application to a number of systems, a = 0.5 has been found to be the

Computation of multicomponent phase equilibria-Part most suitable. Thus, we have f(&) = OS(1 + e,-‘) f xi - 1 ( is1 )

(12)

which is simply the arithmetic mean of Forms 1 and 2. A general algorithm for the double-loop univariate methods consist of the following main steps:

(1) The Standard Form. (2) The Rachford-Rice Form. (c) Wegstein-projected methods. This implementation of the method of Wegstein[31] involves obtaining a projected value of 6, from

en+’ =pefs-,-’ +(i .+ ”

all components”).

(3) ;; (6)

(7)

(8)

of K and 6, from Eq. (1). Compute y from Eq. (3). Normalize x and y and update K. Using the new values of K, update x. If the new x is not within some tolerance of its old value, then go to Step 3. Otherwise go to Step 7. Test for convergence. If the outcome is positive then terminate the iteration. Otherwise go to Step 8. Update 0, and go to Step 2.

The implementation of Step 8 calls for an acceleration technique. Guided by the conclusions reached by Rohl & Sudall[26], two techniques have been studied here: the Newton method and the third-order Richmond method. (For a discussion of the two methods see Lapidus, [15]). (B) Single-loop &variate methods Experience in this study showed that in the implementation of the univariate methods, the inner loop which corrects for the dependence of K on composition can be done away with-even for the most nonideal systems-given a suitable initialization scheme, with significant savings in computation effort and time. A general single-loop algorithm differs from that outlined above only in the following respect: Step 5 involves updating both x and y, and Step 6 is omitted. The methods studied in this section are discussed under four subsections, labelled (a) through (d), below. (a) Richmond-accelerated methods. Drawing on the experience gained from applying the double-loop univariate algorithms, the Richmond acceleration technique was chosen and applied to the following formulations: (1) (2) (3) (4)

The Standard Form. The Rachford-Rice Form. Alternative Form A. Alternative Form B.

(b) Metho& utilizing the mean-value theorem (MU). The mean-value theorem of differential calculus (HolIand,[l 11) is employed as the acceleration technique. After the nth iteration, 0, is updated from e ”n+ 1 = e: -f(e:yj-ye,m)

(13)

e,m = et - af(e:)pye:)

(14)

ez8

(15)

where P = g”k” g =

- gn-‘),

(16)

e,, - e“9

(17)

and f?,, is obtained from 0, through a “direct iteration” which employs x and y values based on K values computed using normalized values of x and y obtained from yi = KiZi/[l + (Ki - l)@“]* Ii = zi - Fe&

(18) (19)

and

j=1

xi = Ii/ 5 4

(?O)

(Zi is the total moles of i in the system). This approach, which keeps xi and yi in mass rather than equilibrium balance, is found to be generally better. To initiate the procedure, 0, is set equal to e,, after the first iteration. The following “direct iteration” schemes were studied: Scheme 1

[i=l 1’

(21)

e,, = 0, f yi

where t is some positive parameter. A number of t-fixing methods were tried but were found not to be much better than the simple scheme using t = 1.0. This value was therefore adopted. Scheme 2 This scheme involves determining t&, from 0, by applying Richmond’s method to the standard formulation. A similar approach using Newton’s method gave poor performance. Scheme 3 This scheme is like Scheme 2 except that it employs the Rachford-Rice formulation. Scheme 4 This scheme combines the Rachford-Rice formulation with MVT acceleration technique. Scheme 5 O,, is obtained from 8, by applying the Richmond acceleration technique to “Alternative Form B”. (d) Quadratic form of Wegstein’s projection. This method employs the three most current values of (0,,, g>--g is defined in Eq. (17)-to update 0” from en+, ”

where

149

andf(f?J is the derivative off(&). The optimal value of a was found empirically to be 0.6. The technique was applied to two formulations:

(1) Initialize x, y, K and 0, (Boldface implies “for

(2) Compute x corresponding to the current values

I

_gY’r* -

Aen_* 41424%[ c2 OS ---+g:;l+fie:j g”

1

(22)

M. 0. OHANOMAH and D. W. THOMFSON

150

where &=gn-1-gn-2, &=g”-g’-’ and & = g” - g”- I. A new 0,,, is obtained from 8, by applying the MVT technique to the Rachford-Rice formulation. ms choice was based on the outcome of a preliminary study of the schemes in section (c).] There is only one direct iteration between acceleration steps. To initiate the procedure, the normal Wegstcin technique is employed until after the third iteration. (C) Free energy-minimization (multivariate) methodr In a comparative study of free-energy-minimization algorithms for multiphase chemical equilibria, Gautam and Seider[5] concluded that the RAND algorithm as implemented by Dluzniewski and Adler[3] is the best. The RAND algorithm has been developed here for the specific case of vapourliquid phase equilibria. A more efficient modification of it has also been developed. The problem had also been studied from the point of view of geometric programming, which has hitherto received some application in chemical equilibria (Duflin et al. [4]; Passy and Wilde[20]; Lidor and Wilde, [16]). A method based on the perturbation theory of geometric programming has also been studied. The Gibbs free energy, G, for a vapour-liquid system is given in a dimensionless form as a function of vapour molar composition, v, and liquid molar composition, I, by the relationship

B, is an N element vector whosejth element is - p$RT

- In (y;x;),

A, is an N x N matrix whose elements are given by ,4jk=svj

1 forj = k

$with&=

I Oforj#k

A, is defined in a manner similar to A,, with u and V replaced with I and L respectively. With (Av, Al} obtained from solving Eq. (26)-this was done by a Gauss-Siedel matrix inversion technique (see Lapidus, [IS&-a new value of {v, I} is determined from {V+‘,l”+‘} = {v”,I”} +sI{Av”, Al”}.

(27)

A golden-section search (Himmelblau[9]) is implemented to determine 1*, the value of 1 which minimizes G/RT subject to Eq. (24) and (25). (2)The modijied RAND method. The modification has the effect of eliminating the standard Gibbs free energy, reducing the number of independent variables to N and doing away with matrix inversion. It can readily be shown that = - In f

(& - p$/RT

p.

(28)

Eliminating li between Eqs (23) and (24), combining the result with Eq. (28), rearranging and introducing Ki, we have

+l[j$+ln(x$)]} (23) (29) where p” denotes the standard chemical potential. This function is minimized subject to vi+li-zi=o

(i=1,2,...,N)

where

(24)

g(v) =$

- i

Z&RT.

i=l

and

o
(i = 1,2,. . . , N).

(25)

The various methods employed in solving Eq. (23) are briefly outlined below. (1) The RAND method. Applying the RAND method, as adapted by Dluzniewski & Adler[ 31, to ECq. (23) (it should be noted that the method has been modified slightly here by treating the total phase contents, L and V, as dependent variables), we obtain -I,

0

Al

0

-I, A,

-I,

Av

A,

0

-1,

u

=

Bu B, (26)

- ln(&I Py:), ~,

AAv=B

(30)

where B is an N-element vector whose jth element is In (KYx;/yjll) and A is an N x N matrix given by

B,

where u is a vector of Lagrange multipliers uj 0’=1,2,. . . , N), IN is an N x N identity matrix, 0 is an N x N null matrix, Al is an N-element vector whose jth element is (4 - I;), Av is defined in a way similar to Al, B, is an N-element vector whose jth element is - zj, B, is an N-element vector whose jth element is - &/RT

(Note that the last term is constant at isothermal conditions and does not affect the location of the minimum.) When the RAND approach is applied to Eq. (29) the result is

Taking advantage of the special form of matrix A, Eq. (30) can be manipulated to give

where rj = xjn y;/Zj.

Computation of multicomponent phase equilibria-Part I

where

A new value of v is given by v”+‘=v”++*Av.

g(v” + ldv)
-KJli]ui+

KiZJi=O

(i = 1,2,. . . , N)

(33)

where li = i vY In solving the set of N nonlinear i+i

substitution method equations, ‘. a successive (Lapidus [ 151) was found to perform better than the Newton-Raphson method (Henley & Rosen, [S]) and the method of continuity (Perry & Chilton,[21]). In updating vi by the successive substitution method, various arrangements of Eq. (33) have been studied. Updating schemes that use the Newton & Richmond acceleration techniques have also been tried. The best result was obtained with the Newton method. Thus,

u: +’ =

vi” - F(o;)/Fyv,“)

a = (t” - t”-‘)‘t”/(t” - F-‘)r(p

(32)

In determining A*, two schemes were adopted. Scheme 1. This involves determining the value of I that minimizes g(v) subject to Eq. (25). Scheme 2. This involves a step-limited unidimensional search, the search being terminated as soon as a value of 1 is found such that

(1 -KJu?-[F-Z&-(1

151

(34)

-P-l),

f-V,-v

(38)

and v, is obtained from Eq. (34) based on v. Thus, z& replaces 07+ ’ in Eq. (34). The method requires two successive-substitution iterations to get started and one in between acceleration steps. (c) Tetrahedral projection. This method requires three successive-substitution iterations to get started and only one iteration in between acceleration steps. It updates v from vn+’ = P,v,“-2+ P2v;-’ + P,vs”

.:+l=u:+~~~+~,rilB:l(L

- ,&rj”‘)]

(i = 1,2,. . . , N) (35)

and F’(L+)is the derivative of F(uJ with respect to v,. To speed up the convergence of the iteration, it was necessary to apply an acceleration technique. Three acceleration methods were studied in this work. These are discussed below. (a) Hyperpiane linearization. This is an implementation of the method proposed by Barreto & Farina[l] and described by the authors as the accelerated univariate algorithm (AUA). The basic form of the method gave better performance than the alternatives proposed by the exponents. However, the method has two major shortcomings: (1) There are N successive-substitution iterations between accelerations; thus, the accelerating effect is unimportant for large N. (2) Every acceleration step involves matrix inversion,-a time-consuming process. The two vector-projection methods outlined below (They are treated in detail by Ohanomah, [19]) are without these shortcomings. (b) Triangular projection. The vector v is projected according to v”+‘=av,“-‘+(l-a)v:

(36)

(39)

whereP,=l-a+(a-l)m,P2=a(l-m),P,=lm, m = y[(l - a)6 + 2/?C]/a(4AC - S’), I = 2aA/y, y=A+B-C, 8=B-A-C, fi=(E-H)/C, a=(D-_)/A, A=E+D-2G, B=F+D-2I, C = F + E _ 2H, D = f” - 2jTf” - 3, E = f” - IITf’”- I), F=f”T f”, G=f+“7 p-l”, H=f’“-1” f”, and I = fb-217 f” After every negative convergence test, the following assignments could be made to reduce computation: D=E, E=F, G=H, A=C, and a=/3 (These relationships were not taken advantage of in obtaining the results presented later in this work). (4) Sensitivity method. This method is based on application of the perturbation theory of geometric programming (Duliin et al.[4]) to the Gibbs free energy (see Ohanomah,[19]). When the resulting relationship is manipulated to eliminate matrix inversion and a modification incorporating the mean value theorem of differential calculus introduced, the result is

where f(q) = (1 - Ki) v; - [F - ZiKi - (1 - Ki)li] ui + KiZili

(37)

(9

where rj = xpjJZj,

Bj = In (Kirjly,) and “m” denotes a point such that v” is, given by v.m= I

0.n I

B:+j,r;B;l(l (i = 1, 2,. . . , N).

-fir;)]} (41)

Note the similarity between Eq. (31) and (40). A variation of the sensitivity method which employs tetrahedral projection to speed up its convergence was also implemented. In this method, v”+’ from Eq. (40) becomes v; in Eq. (39). In the basic implementation of the free-energyminimization algorithms, K-values are computed once in very iteration. Reducing the frequency of K computation. The effect of the reduction of the frequency of computation of K-values to less than once in every iteration was studied, using the following algorithms, the choice of which was based on the relative merit of the algo-

152

M. 0. OHANOMAHand D. W. THOMPSON

rithms in each class: (1) The Wegstein-projected MVT-accelerated Rachford-Rice method. (2) The tetrahedral-projection-accelerated sensitivity method. (3) The tetrahedral-projection-accelerated GP method. (4) The sensitivity method. The frequency of K computation was varied from “once every iteration” to “once every four iterations”. With the univariate formulation, this amounts to a double-loop iteration with the 6, loop nested with the K loop but with the number of iterations in the inner loop fixed. It is obvious that the final results will be in error to an extent which depends on how sensitive the K-values are to composition variations as well as how much less frequently K is computed. A deviation parameter, E,,was introduced and defined by

and $J (Ki- l)Zi/(l i=l

Fz=2 Step 3. If

+ Ki).

F > 0 then f$‘= 0.5(1 -F),

else

tI,o = (F3- 0.5Fa)/(F, - F4)

(44)

where +2tz’-, i-1 (1 + KJ and

(45)

Fp&. I Step

where subscript “f’ denotes a frequency of “once every f iterations” and V,is the value of V for f = 1. Where it is assumed that E < 0.001 is a tolerable deviation, results indicate that f = 2 is appropriate for the univariate and the sensitivity methods and that f = 3 is satisfactory for the tetrahedralprojection-accelerated GP method. The tetrahedralprojection-accelerated sensitivity method was found to be good only with f = 1. Initialization scheme. Five schemes for univariate methods and four schemes for multivariate algorithms were studied. The general scheme that gave the best performance works as follows: The initial value of the vapour fraction is determined by regula-f&i interpolation based on K computed using x = y = z. The interpolation is first carried out using values of f(O,)=0 and f(O,) = 0.5 estimated from:

4.To guard against overprojection, 0.05 d e,O< 0.95.

Step 5. Where necessary, compute

u: = zi/]l + (1 -

euoyKie,q(i=

fteJ

= i$, 1 + (K”L l)O ” - l.

The details of the method are outlined through 5 below:

in Steps 1

Step 1. Compute K using x = y = z. Step 2. Compute

F = F2/(& - F,)

1,2,. ..,N)

(4)

where initial values are also required for 1,phase total moles and mole fractions (as in the RAND algorithms and the sensitivity algorithm), they are derived from 9 through component balance, summation and normalization respectively. The values of 0,Oobtained by this initialization scheme are generally within + 0.2 of the solution. Stopping criterion. The same stopping criterion was used for all algorithms to ensure a fair comparison. In every application, the n th iteration was considered to have yielded the desired solution if

ii$, If the resulting value of fI, is less than 0.5 it is accepted as &,O,and experience shows that this provides a good starting point. Wowever if the resulting value of 8, is greater than 0.5 it has been found that a better value for 0,Ocan be estimated by making a second interpolation using f(l.0) and f(O.S), obtained from the relationship:

set

10: - ui”-‘I < lo-‘.

(47)

Property estimation. This study has been mindful of the industrial practitioner who is out to make the most judicious use of available thermodynamic information. In view of this, the choice of a predictive correlation for any property has been based on the quality of the resulting estimate, the relative simplicity of the estimation method (compared to other methods that yield results of about the same quality), and the ease of determination of the parameters involved in the correlation. Saturation pressures, P,", were estimated from a six-parameter equation of the type employed by Prausnitz et al. [23], while vapour molar volumes and compressibility ratios were determined from the Wilson modification of the Redlich-Kwong equation of state-WRK--(see Reid et al.[25]). The same EOS was employed in estimating 4i from the relation (see Prausnitz,(22])

(43)

where

F,= $ KiZi - 1 i=l

where V is the total vapour volume, Z is the compressibility and nj is the number of moles of wmponent _j in the phase.

Computation of multicomponent phase equilibria-Part

Liquid molar volumes were estimated from a polynomial fit based on one, two or three. values supplied as a function of temperature. Values of yi were obtained from the Wilson model (Wilson [34]). For fl”, the relationship is Ji” = c$;P; exp [(P - P,“)v,‘IRT]

(49)

where v,.’denotes liquid molar volume. The saturation

I

fugacity coefficient, $f, was determined from the WRK EOS for T < 0.56 T,; otherwise, it was estimated according to the method of Lyckman and co-workers (see Prausnitz et al. [22]). Applications. Table 1 contain‘s information on the systems on which this comparative study was based, while Table 2 contains a summary of the source references for the relevant thermodynamic data. The systems have been ordered in Table 1 in order of

Table 1. Information on systems employed System

Pressure

Composition (mole %)

A

1.0

Cyclohexane (18.6); Benzene (42.0); Isopropanol (35.2); Methyl-ethyl Ketone (4.2).

B

3.7

acetone(25.7); Methanol

C

1.316

Nitrogen (30.36); Argon (29.48); Oxygen (40.16).

D

1.0

Ethanol (10.0); Chloroform (20.0); Acetone (50.0); N-Hexane (20.0).

E

1.0

Benzene (3.6); Chloroform (82.2); Methanol (6.3); Acetone (4.2); Methyl Acetate (3.7).

F

1.0

Acetone (32.4); Chloroform (32.6); Methanol (12.8); Dimethylbutane (22.2).

G

1.0

N-Hexane (37.1); Ethanol (3.4); Benzene (9.4); Methylcyclopentane (50.1).

H

1.0

Benzene (45.6); Chloroform (5.3); Methanol (5.4); Methyl Acetate (43.7).

I

1.0

N-Hexane (18.6); Nethycyclopentane (25.7); Benzene (15.3); Cyclohexane (11.8); Toluene (28.6).

.I

1.0

water (27.217); Ethanol (22.783); Acetic acid (22.783); Ethyl acetate (27.217).

K

1.0

N-Hexane (18.6); Methylcyclopentane (20.7); Benzene (15.3); Cyclohexane (11.8); Toluene (28.6); Ethanol (5.0).

L

1.5

Cyclohexane (50.0); Hexadecane (50.0).

M

1.0

Cyclohexane (17.0); Hexadecane (49.0); Benzene (17.0); Tetrachloromethane (17.0).

N

1.0

Tetrachloromethane (50.0); Hexadecane (50.0).

0

1.0

Benzene (50.0); Heptadecane (50.0).

P

1.0

Cyclohexane (50.0); Eicosane (50.0).

(64.0); Water (10.3).

-

Table 2. Source references for thermodynamic

Data types

data

Data sources

Acentric factors; constants for vapour-pressure relation; critical constants; dipole moment; molar volume

Perry and Chilton (1973); Prausnitz et al. (1967); Reid et al. (1977).

Association constant

Prausnitz et al. (1967).

Wilson parameters

153

Gmehling et al. (1977); Hudson and Van Winkle (1969); Nagata (1973); Prausnitz et al. (1967); Sanderson and Chien (1973); Weatherford and Van Winkle (1970); Willock and Van Winkle (1970).

M. 0.

154

OHANOMAH

increasing boiling range, which ranges from 0.48”C for System A to 196.89”C for system P. The applications involved applying each algorithm to each system at 19 temperatures which correspond to a uniform distribution of equilibrium vapour

and D. W. THOMPSON

fraction from 0.05 to 0.95. The results, based on a PL/I (F-version 5.5) compiler and an AMDAHL-470 computer, are presented in Tables 3-6. The grouping of the systems, for neatness of presentation, is justified by the fact that any variations in the relative

Table 3. Computation time (CPU sec.) for double-loop univariate methods Acceleration 6 Initialization

Newton with proposed initialization

Richmond with proposed initialization

Formulation'

CPU time (xc) for four-system groups A-D E-H M-P I-L

Total CPU time (set).

Standard

3.8342

4.2971

4.5601

2.7067

15.3981

Rachford-Rice

3.7698

4.3276

4.1502

2.4991

14.7467

Barnes-Flores

3.7616

4.3022

4.1325

2.5574

14.7537

Log Standard

3.8592

4.3881

4.5156

2.5721

15.3350

Standard

3.8359

4.2572

2.4837

14.8688

Rachford-Rice

3.7818

4.2920 . 4.2679

4.1719

2.4794

14.7010

Barnes-Flora

3.7974

4.2614

4.1988

2.4879

14.7455

Log Standard

3.8671

4.3228

4.3004

2.5665

15.0568

Alt. Form A

3.8449

4.2929

4.2094

2.4791

14.8263

Alt. Form B

3.7821

4.2234

4.1351

2.4738

14.6144

Newton with 80 " = 1.0

Standard

7.2629

7.0598

7.8942

5.8007

28.0176

Richmond with e; = 1.0

Standard

6.1953

6.0908

5.8790

2.9226

21.0877

Table 4. Computation time (CPU set) for single-loop univariate methods Acceleration

Formulation

CPU time (sex) for four-system groups A-D E-H I-L M-P

Standard

2.3346

2.6508

2.4872

1.6566

Total CPU time (see).

9.1292

Rachford-Rice

2.3106

2.5794

2.4315

1.6008

8.9223

Alt. Form A

2.3151

2.5778

2.4500

1.6068

8.9497

Alt. Form B

2.2860

2.5732

2.4414

1.5941

8.8947

Standard

2.2304

2.5830

2.4089

1.5746

8.7969

Pachford-Rice

2.2115

2.4958

2.3565

1.5062

8.5700

Direct itera2.1199 tion (scheme 1)

2.3403

2.2428

1.3364

8.0394

Standard

2.1019

2.3340

2.2240

1.4381

8.0980

Rachford-Rice

2.0094

2.3004

2.1799

1.4122

7.9019

Wegstein L En'T

Rachford-Rice

2.0096

2.3140

2.1416

1.3951

7.8603

Wegstein & Richmond

Alt. Form B

2.1108

2.3221

2.2058

1.4410

8.0797

Quadratic

Rachford-Rice

2.0363

2.2922

2.1877

1.3855

Richmond

MVT

Wegstein

wegstein

L

Richmond

wegstein & Richmond

Computation of multicomponent phase equilibria-Part

I

155

Table 5. Computation time (CPU set) for free-energy minimization methods Method

CPU time (set) for four-system groups E - 11 I-L A-D M-P

Total CPU time (set).

Modified RAND scheme 2

2.8407 (5.3Xf)

3.2973 (2.6%f)

3.0267 (2.6Xf)

1.7779

10.9426 (2.6%f)

GP with hyperplane linearization

7.3769

6.7042

5.0607

1.9732

21.1150

GP with triangular projection

3.2732

3.6389

2.8734

1.4606

11.2461

GP with tetrahedral projection

2.9047

3.0917

2.8098

1.4935

10.2997

Sensitivity

2.2142

2.6703

2.5427

1.8571

9.2843

Sensitivity with tetrahedral projection

2.4035

2.2833

2.0860

1.5284

8.3012

f = failure

of the algorithms seem to depend mainly on boiling range. Comment on results. Table 3 reveals that the proposed initialization scheme results in significant reduction in computation timefor example, a reduction factor of about 2 for the Newton-accelerated standard method-when compared with a simple terminal initialization scheme. It also has the effect of reducing the disparity in the performance of the algorithms to a near-insignificant level. The results for the single-loop univariate methods (Table 4) show that the Wegstein-projected methods are generally the fastest, the version combining the MVT technique and the Rachford-Rice formulation being slightly the best. A comparison of the freeperformance

energy-minimization algorithms (Table 5) reveals that the modified RAND method is rather unreliable, and that the vector projection methods are far more efficient accelerators of the GP method than the hyperplane linearization technique. The combination of the sensitivity method and the tetrahedral projection technique is the best here. In Table 6 where the best algorithms from the three classes have been compared, it is evident that the double-loop univariate method is not competitive, and that the single-loop univariate approach has a slight advantage over the multivariate method, this advantage becoming more pronounced where a maximum relative “error” of 0.001 in V is tolerable so that K-values are computed on alternate iterations.

Table 6. Final comparison of algorithms Method

CPU time (set) for four-system groups A-D E-H I-L M-P

Total CPU time (set).

Double-loop, Richmondaccelerated Alt. Form B

.3.7422

4.1476

4.0560

2.4589

14.4047

Single loop (Wegstein L MVT 6 Rachford-Rice)

2.0462

2.3306

2.1379

1.4303

7.9450

Sensitivity 6 tetrahedral projection

2.4152

2.2609

2.0736

1.5068

8.2565

1.7627

2.0875

1.8935

1.3245

7.0682

Sensitivity - K computed on alternate iterations

2.1364

2.3380

2.2255

1.4894

8.1893

GP with tetrahedral projection - K computed every third iteration

2.7280

2.8792

2.4291

1.3739

9.4102

Single-loop (Wegstein & MVT & Rachford-Rice) K computed on alternate iterations

M. 0.

156

OHANOMAH and

NOMENCLATURE

! g($ G K I L N P R T ; x Y

;

function liquid-phase reference fugacity totai number of moles in the system parameter defined by Eq. (17) function defiued bv Ea. (29) Gibbs free energy _ _ . ’ equilibrium ratio moles of a given component in the liquid phase total moles in the liquid phase number of components in the system system pressure universal gas constant temperature moles of a given component in the vapour phase total moles in the vapour phase mole fraction of a given component in the liquid 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 y activity coefficient of a given component in the liquid

phase

0” vapour phase fraction Iz li p I#J

uuidimensional search variable parameter defined in Eq. (33) chemical potential vapour-phase fugacity ciefficient for a given cornponent

Subscript i, j, k component

I liquid phase s successive substitution or direct iteration v vapour phase Superscript n iteration

0

standard state; initial value

REFElwNcEs

1. G. F. Barreto & I. H. Farina, An efficient univariate method to solve systems of nonlinear equations, Chem. Engng Sci. 34, 63 (1979). 2. R. J. Clasen, The Numerical Solution of the chemical equilibrium problem, in The RAND Corp. Rh4-4345PR, January (1965). 3. J. H. Dluzniewski & S. B. Adler, Calculation of complex reaction and/or phase equilibria problems, I. Chem. E. Symp. Series 35, 4 (1972).

4. R. J. Dutlin, E. L. Peterson & C. Zener, Geometric Programming-Theory and Application. Wiley, New York (1967). 5. R. Gautam & W. D. Seider, Computation of phase and chemical equilibria-I: local and constrained minima in Gibbs free energy, AZChE J. 25, 991 (1979). 6. B. George, L. T. Brown, C. H. Farmer, P. Buthod & F. S. Menning, Computation of multicomponent, multiphase equilibrium, Ind. Eng. Chem. Process Des. Deve/op. 15(3), 372 (1976). 7. J. Gmehling & U. Gnken, Vapour-Liquid Equilibrium Data Collection, Vol. 1 (Parts 6A, 6B, & 7). DECHEMA, Frankfurt (1977). 8. E. J. Henley & E. M. Rosen, Material and Energy Balance Computations. Wiley, New York (1969). 9. D. M. Himmelblau, Applied Nonlinear Programming. McGraw-Hill, New York (1972).

D. W. THOMPSON 10. Y. Hirose, Y. Kawase & M. Kudoh, General flash calculation by the Newton-Raphson method, J. Chem. Eng. Japan 11(2), 150 (1978). 11. C. D. Holland, Fumhsmentals and Modeling of Separation Processes-Absorption, Distillation, Evaporation, and Extraction. Prentice-Hall, Englewccd Cliffs, New

Jersey (1975). 12. C. D. Holland, Fumhunentals of Multicomponent Distillation. McGraw-Hill, New York (1981). 13. J. W. Hudson Kc M. Van Winkle, Multicomponent vapour-liquid equilibria in system of mixed positive and negative deviations, J. Chem. Engng Data U(3), 310 (1969). Processes, (2nd Edn). 14. C. J. King, Separation McGraw-Hill, New York (1980). 15. L. Lapidus, Digital Computation for Chemical Engineers. McGraw-Hill. New York (1962). 16. G. Lidor & D. J. Wilde, Transcendeutal geometric programs, J. Optimization Theory and Applications, 26(l), 77 (1978). 17. Y. H. Ma & C. W. Shipman, On the computation of complex equilibria, AIChE J. 18(2), 299 (1972). 18. I. Nagata, Prediction accuracy of multicomponent vapour-liquid equilibrium data from binary parameters, J.Chem. Rngng Japan 6(l), 18 (1973). _ _ 19. M. 0. Ohanomah, Computational algorithms for multicomponent phase equilibria and distillation, Ph.D. dissertation, Univ. of British Columbia, Vancouver (1981). 20. U. Passy & D. J. Wilde, A geometric programming algorithm for solving chemical equilibrium problems, 21. Siam J. Appl. Math. 16(2), 363 (1968). R. H. Perry & C. H. Chilton, (eds), Chemical Engineers’ Handbook, (5th Edn.), pp. 2-54. McGraw-Hill, New York (1973). 22. J. M. Prausnitz, Molecular Thermodynamics of FluidPhase Equilibria. Prentice-Hall, Englewocd Cliffs, New Jersey (1969). 23. J. M. Prausnitz, C. A. Eckert, R. V. Grye & J. P. O’Connell, Computer Calculations for Multieomponent Vapour-Liquid Equilibria. Prentice-Hall, New Jersey (1967). 24. H. H. Rachford Jr. & J. D. Rice, Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium, J. Petrol. Technol., Section 1, October, p. 19 (1952). 25. R. C. Reid, J. M. Prausnitz & T. K. Sherwood, The Properties of Gases and Liquids, (3rd Edn.). McGrawHill, New York (1977). 26. J. S. Rohl & N. Sudall, Convergence problems encountered in flash equilibrium calculations using a digital computer, I. Chem. E. Symp. Series 23, 71 (1967). 21. R. V. Sanderson & H. H. Y. Chien, Simultaneous chemical and phase equilibrium calculation, Znd. Engng. Chem. Process Des. Develop. 12(l), 81 (1973). 28. D. U. Von Rosenberg, Determination of the polynomial defining vapour-liquid split of multicomponent mixtures, Chem. Engng Sci. 18, 219 (1963). 29. D. U. Von Rosenberg, On the polynomial defining vapour-liquid split, Chem. Engng Sci. 32, 1546 (1977).30. R. M. Weatherford & M. Van Winkle. Vaoour-liauid equilibria of the quinary system hexane, methylcyclopentane, cyclohexane, benzene and toluene, J. Chem. Engng data H(3), 386 (1970). 31. J. H. Wegstein, Accelerating convergence of iterative processes,-Commun. A.C.M.9(6), 9 (1958). 32. W. B. White. S. M. Johnson & G. B. Dantzie, Chemical equilibrium in complex mixtures, J. Chem. Whys. zs(5), 751 (1958). 33. J. M. Willock & M. Van Winkle, Binarv and ternarv vapour-liquid equilibria of methanoi-acetone-2,3dimethvlbutane. J. Chem. Enana Data 1x2). 281 (1970). 34. G. M.-Wilson,’ A new expression for the excess free energy of mixing, J. Amer. Chem. Sot. 86, 127 (1964).