Computation of Multicomponent, Multiphase Equilibrium Boyd George,’ L. P. Brown,* C. H. Farmer, Paul Buthod, and F. S. Manning’ The University of Tulsa, Tulsa, Oklahoma
Equilibrium compositions in multicomponent, multiphase mixtures are computed by minimizing the Gibbs free energy. The Fortran program converts the Gibbs free energy of the mixture into a novel unconstrained form, then minimizes using Powell’s method. Adiabatic equilibrium can be handled and the concentration of trace components for certain cases can be estimated separately. The program handles up to four phases and 15 components. Only 20K storage using double precision and less than 2 min total CPU time on a Sigma 6 are required. Eleven examples illustrate the power, flexibility, and wide applicability of this method.
Introduction The present objective was to develop a general computer program which would calculate the equilibrium compositions in multiphase mixtures. This program should be: (1)usable by nonexperts after a minimum of instruction; (2) versatile in that simultaneous chemical and physical equilibrium must be included. Traditional methods of computing equilibrium compositions involve the simultaneous solution of the appropriate set of equations. In general, three types of constraining equations must be simultaneously satisfied: (1)composition, Z xi = 1; (2) physical equilibrium, such as yi = Kixi; and (3) chemical equilibrium, such as AGO = -RT In K,,,. These constraints usually are easily obtained; however, as the number of equations increases, their simultaneous solution rapidly becomes very tedious, if not prohibitively difficult. When chemical reaction does not occur, the total number of equations is N , N ( N , - l),where N , is the number of phases and N is the number of components. When chemical reaction occurs, one additional restriction must be added for each independent reaction. The total number of components, N , the number of independent components, N I , and the number of independent reactions, N R , are related by
Theory The Gibbs free energy expression is derived first for a single phase and then converted into an unconstrained form. Finally the work is extended to more than two phases. A. Formation of Gibbs Free Energy and Constraints. Calculation of equilibrium conditions involves: (1)expressing the total Gibbs free energy in terms of the desired variables, such as masses and compositions of all phases present; in turn, this requires that the individual component free energies, or chemical potentials, or activities, or fugacities, be related to phase compositions; (2) minimizing the total free energy subject to material-balance restrictions and, for adiabatic processes, energy conservation. The partial free energy of component i in the equilibrium state of a solution at constant temperature and pressure may be expressed as
+
NI=N-NR Accordingly, either N I or N R must be obtained. Perhaps the most direct method of obtaining N I is that described by Amundson (1966). With the advent of the digital computer, much work has been directed toward the computation of multicomponent, simultaneous chemical and physical equilibrium, e.g., Dayhoff et al. (1967);Dluzniewski and Adler (1972); Ma and Shipman (1972); Sanderson and Chien (1973); and White et al. (1958). These methods employ minimization of the Gibbs free energy as the necessary and sufficient condition for thermodynamic equilibrium; however, the minimization algorithms are complicated by the inherent material balance restrictions. Box (1966) has described the advantages of unconstrained optimization very eloquently. “The self-evident conclusion is that the solution of constrained problems is eased considerably by eliminating the constraints, so that one of the more powerful methods for unconstrained optimization, e.g., that of Powell may be employed.” The algorithms required to solve unconstrained problems are simpler and more efficient than constrained-problem algorithms, because optimization can be performed without worrying over “unfeasible space.” Accordingly the present approach recasts the free energy function into a novel unconstrained form. Amoco Production Company, Chicago, Ill. Cities Service Oil Company, Tulsa, Okla.
372 Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 3, 1976
Gi = Gio
fi + R T In a fi
Therefore, the total Gibbs free energy of a mixture is
G=
N
ZiCi =
i=l
N
L
i=l
&Gi0
N
+ R T kx= l ix= l X i k In f i k
bounded by composition and material-balance constraints. For example, the total free energy of a two-phase mixture would be expressed as
G=
N r=l
ZrGro+ RT
N
X I 1lnf,’
1=1
+ RT
N
X I 2In f r 2 (1)
r=l
The composition variables, X L kand 2, are constrained by the following material balances. (1)For physical equilibrium, a material balance on each component present yields
X,’
+
xr2=
z,
(2) When chemical reactions occur, additional balances on each element yield N
xaJ1ZL=B,
6=1,2 , . . . , M )
1=1
where a,, = number of j atoms in the in ith compound, B, = total number of j atoms in mixture, and M = number of atoms in mixture. B. Conversion to Unconstrained Form (Brown, 1971). The above constrained formulation for the Gibbs free energy (eq 1, the material balance and elemental balance, if applicable) is transformed into an unconstrained form by rewriting it in terms of an allocation function, H ( u ) , where: (1)u is an unconstrained variable, Le., - m < u < a; and (2) H ( u )varies over a finite range, in particular 0 < H < 1. Any function
meeting these properties could be used, although good success has been obtained with Brown's (1971) suggestion
H ( u )=
1
+
1 exp(-u)
+
The physical equilibrium material balance for the ith component is xi1
+ xi2 = 2,
where the Zi are considered constants. The moles of i in the first phase, X i 1 , is some fraction of the total moles in the mixture, Zi, and may, therefore, be expressed as
which, through elementary algebraic manipulation yields
Since N
X ; k = total moles in the hth phase
i=l
the mole fractions are readily available for use in evaluating the corresponding fugacity, f i , as discussed above. Incorporating eq 2 and 3 into eq 1 yields the total Gibbs free energy equation as a function of ui. This unconstrained equation is
G
N
=
and 2, is some fraction of Zi,,,, or
The subscript, i N, in the allocation function is intended to show that the U,+N is a different variable from the u, that separates the 2, between phases, but nevertheless can be a part of the u array. By rearranging the matrix, (4),to move the elements (al,/ B,)Z,, i = M 1 to N t o the right-hand side of the equation, the square M X M matrix is left. The Z,, i = M 1 to N are allowed to vary, and the left-hand side is solved for the remaining Z,, i = 1 to M . In order to eliminate the possibility of negative Z,, i = 1 to M , from arising, care should be taken to make sure the first M components will be present in the equilibrium mixture. If one of the first M components is not found to be present in the computer equilibrium mixture, then that component should be replaced by one which will be present and the optimization calculation repeated. C. Extension to More Than Two Phases and the U Array. Additional phases can be either pure or multicomponent phases. The extension of eq 1to three multicomponent phases is now demonstrated. Beginning as before, with the physical equilibrium component balance
+
+
X,'
+ xr2+ x,3= 2,
it is noted that X I 1is some fraction of 2,. Therefore
N
+ R T i= 1 H ( u i ) Z i lnfil N + R T i = l (1 - H ( u i ) ) Z iIn f i 2
Z$i0
i= 1
and (la)
If a chemical reactions occurs, the Zi vary, so in addition to the physical equilibrium component balance, X i 1 + X i 2 = Zi, the elemental balance N
i= 1
XI2 = H(Ur+N) (1 - H(u1))Zr which means
must also be satisfied. Expanded, the matrix appears as
XL3 = (1 - H ( u , + N )(1 - H(u,))Zr
+ a1222 + . . . alNZN = B1 anlZ1 + a2222 + . . . (12NZN = B2 all21
These equations and the corresponding mole fraction relations are incorporated into the following equation
+ a ~ z Z 2+ . . . aMNZN = B M
If each row is divided by its corresponding B element, the matrix
(4)
aM1 21 -
+
+
+ x,3 = (1 - H ( u , ) ) Z ,
from which X I 2is seen to be some fraction of (1- H ( U , ) ) Z , . Introduction of another allocation function allows expression of X I 2as
(j= 1 , 2 , . . . , M )
a,iZi = B;
aMlZ1
x,2
aM2 2 2 . . . aMN Z N = I +BM BM BM is obtained. Because all elements of the above matrix are non-negative, it is self-evident that the maximum value of any (a,i/Rj)Zi is unity. Therefore, for all j , there is a maximum value of aji/Bj which corresponds to the maximum value which Zi may physically assume. Thus
+ RT
N
i= 1
X i 2 In f i 2
N
+ R T i=1 X i 3 In f i 3
(6)
1
as previously shown for eq 1 and la. Under certain circumstances, the existence of a pure phase or phases must also be considered. These are handled very simply by dropping the log terms in equations such as 1 or 6. This method requires the unconstrained optimization of a variable-dimensioned array, u . The size of the u array, or the number of variables to be optimized, can be predetermined. The procedure is: (1)count all possible phases (excluding only the pure phases); (2) if there is no reaction, subtract 1 from the number obtained in (1)and multiply by N; and (3) if there is a reaction, multiply the number obtained in (1) by N and subtract M .
Program Description Figure 1 presents the flow chart of the equilibrium computation program and the entire program listing is included in the microform edition as supplementary material. As the flow sheet (Figure 1) shows, several decisions are made as to the number of variables to be minimized. Minimization then proceeds with Powell's method (Fletcher and Powell, 1963; Powell, 1970) until either convergence or the maximum Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 3, 1976
373
INPUT DATA DECK
Table I
SET NUMBER OF VARIABLES
OF VARIABLD
Free energy of formation at 500 K divided by RT (dimensionless) (AMU) Vapor phase fugacity, atm at 500 K (FV) Liquid phase fugacity, atm at 500 K (FL)
WNVERGENCE FREE ENERGY
a
""&-t-pJ
CALL WHITE
COMPONENT
OUTPUT
SEARCH METHOD
(-) Figure 1. Flow chart of computer program.
number of iterations has occurred. Other decisions are then made to either refine or halt the minimization. Communications occurs between the main or calling program and the listing in the Appendix via both labeled common statements and subroutine dummy variables. The labeled common statements are COMMON/A/Z( 15),FL(15),FV(15),FW(15), AMU(15),ARHS(15)
COMMON/B/IPURE(3),NPURE,IPHASE,K COMMON/C/M,N,B(l5),A(1 5 ~ 5 ) and the subroutine statement is V,ALIQ1,PI,IPRNT). The answers are returned via the dummy arguments in the subroutine list, except for MAX and IPRNT which are user specified. If IPRNT is set equal to 1,the program prints out the Z array after MAX iterations of the first (and most sophisticated) optimization routine, and if a reaction occurs in a "basically" single phase, IPRNT equal to 1will also provide for printing out the Z array for each subsequent iteration. IPRNT equal to zero suppresses all intermediate printout. A "basically" single-phase reaction means either one multicomponent phase, or multiphases in which only one is multicomponent. Also, for the basically single-phase reaction, the P I array is available for the estimation of components which were not part of the maximization calculation. This, of course, is the White, Johnson, and Dantzig method (1958) and is used in the equation = exp(- AMU,
+
M 1=1
a,,PI,)
The variable X in the subroutine dummy list is the array that contains the mole fractions of the first phase, and ALIQ is the total moles of that phase. The same is true of Y and ACAP, 374
Cyclohexane
39.62
34.40
0.0
18.6
16.2
30.0
17.5
15.0
1000.0"
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 3, 1976
Hydrogen
Assumed.
W and ASOL, and V and ALIQ', for the second, third, and fourth phases. The program, as written, can handle up to four phases, but only three multicomponent phases; thus with a four-phase problem, at least one phase must be a pure phase. This is simply a programming limitation, but allows handling of all the examples listed previously. The X, Y, W, V, and PI arrays must be dimensioned in the main or calling program. Powell's (1970) method was selected for minimizing the objective function as a result of extensive comparison of the following alternative optimization procedures: simplex (Nelder and Mead, 1964);rapid descent (Fletcher and Powell, 1963); conjugate gradient (Fletcher and Reeves, 1961); and direct search (Hooke and Jeeves, 1961).When first derivatives of the objective function were required by the optimization method, the derivatives were obtained numerically; therefore analytical differentiation of the objective function is not necessary.
Program Use Use of the program is illustrated by the following example. Example 1. Estimate the composition and quantity of the gas and liquid phases in equilibrium when benzene is hydrogenated at 500 K and 30 atm pressure. Assume only cyclohexane, benzene, and hydrogen are present in the equilibrium mixture. The data listed in Table I are available. The chemical reaction is C6H6 -k 3Hz
SUBROUTINE MF(MAX,X,ALIQ,Y,AVAP,W,ASOL,
Y,
Benzene
F?
C6H12
The hydrogen to benzene ratio is 3.05. Step 1. Adjust the moles of components (or initial reactants) so that the sum of the moles is between 1 and 10. Accordingly, the input data to the program are 1mol of benzene, 3.05 mol of hydrogen, and 0.0 mol of cyclohexane. This will be the feed for the problem and is placed in the Z array. Step 2. Set N equal to the number of components including possible products. Therefore, N = 3. Step 3. Usually for no chemical reaction the AMU array, IPURE array, and NPURE are all set equal to zero. Where there is chemical reaction, step 3 will be as shown after the discussion of step 9. Step 4. Set IPHASE equal to the total number of phases. Here, IPHASE = 2. Step 5. Write a function subprogram in the form FUNCTION ACTCO (K) IMPLICIT REAL*8(A-H, 0 - Z ) where the K refers to the component subscript in the Z array. If activity coefficients are not used, set ACTCO equal to unity in this subprogram. If activity coefficients are used, include a labeled common statement and dimension statement after the IMPLICIT statement. as follows
Table 11. Input Data for Fugacity Information in Two-Phase Equilibrium Equilibrium constant method
Fugacity method Component
FL
1. Benzene 2. Cyclo-
hexane
3. Hydrogen
FV
ACTCO
FL
17.5 15.0
18.6 16.2
1.0 1.0
0.9408 0.9259
1000.0
30.0
1.0
33.3
where the variable “dummy” is the liquid moles of component K and “dummy total” is the reciprocal of the total liquid moles, so that activity coefficients can be related to the concentration of the liquid phase. Step 6. Supply the fugacities for the respective phases. For single phase vapor reactions each term in the FL array is set equal to the total pressure in atmospheres, the FV array is not used, and the activity coefficients may be defined or set to unity as desired. For two-phase equilibrium there are three ways of defining fugacities: (a) Fugacity Method. The fugacity of each component in the liquid phase is read into FL. The fugacity of each component in the vapor phase is read into FV. All values of ACTCO are set to 1. (b) Equilibrium Constant Method. The vapor-liquid equilibrium constant K for each component is put in the FL array. Both FV and ACTCO are set to 1.for all components. (c) Activity Coefficient Method. The vapor pressure of each component is put into FL. The total pressure of the system is put into FV and the activity coefficient for each component is calculated and returned from the user supplied function ACTCO. Example values of each of these three methods are shown in Table 11. All pressures and fugacities must be in atmospheres. For three phases, say two liquids and a vapor, the FW array is introduced. Using the equilibrium constant method, set the FV array to unity and the FL and FW arrays equal to their respective individual equilibrium constant values. S t e p 7. If no reaction occurs, set M equal to zero. If a reaction occurs, see step 9. S t e p 8. Set MAX equal to an integer value of 100 or ten times the number of actual variables to be optimized, whichever is greater. This number is simply a guide when no better value for MAX is available. The above steps complete the input data for an equilibrium problem with no reaction. S t e p 9. If reaction occurs the matrix A must be specified to ensure the correct elemental material balances. For the current example the initial reactants are CsH6 and H2 while the possible products are CsH12. The matrix A is not listed with atoms vertically and components horizontally H2 0 2
FV
ACTCO
1.o 1.0
1.0 1.0
1.o
1.0
FL
FV
ACTCO
20.91 19.45
30.0 30.0
1.350 1.428
30.0
1.0
1000.0
Table 111
FUNCTION ACTCO (K) IMPLICIT REAL*8(A-H, 0 - Z ) DIMENSION dummy (15) COMMON/ALFA/dummy, dummy total
C H
Activity
coefficient method
C6H12
C6H6
6 12
6 6
As previously mentioned, in order to keep the minimization unconstrained, the first M components must be present in the equilibrium mixture. M is always set equal to the number of rows of the A matrix, which in the case is 2. I t is relatively simple to juggle the components so that the first M are present a t equilibrium. In this case H2 is fed in excess of the stoi-
Component
AGf
AGfIRT
CsHs CsHiz
39 240 34 070 0.0
39.62 34.4 0.0
HZ Table IV Component CsHs
Liquid phase mole fraction
Vapor phase mole fraction
CaHiz Hi
3.867 15 X 10-4 0.997 328 0.228 54 X
10-2
10-1
Total moles
0.391 186
0.659 989
3.643 84 X 10-4 0.923 452 0.761 839 X
Total moles 3.917 67 X 10-4 0.999 608 0.051 1745 1.051 175
chiometric minimum and CsHl2 is the only product; therefore the choice of the first M components is obvious. In writing the A matrix, it is not always necessary to use atoms in the row designation; groups of atoms, such as OH, are just as valid. S t e p 3. When a reaction occurs, the following procedure should be used for AMU, IPURE, and NPURE. (1) Set NPURE equal to the number of expected pure phases. In the example, there are no pure phases and NPURE would equal 0. (2) Set the IPURE array. In the example, since NPURE is 0, the IPURE array is also all 0’s. Up to three pure phases can be considered. If pure carbon were present in the present example the A matrix would have to be extended.
C H
H2 0 2
CRHl2 6 12
C6H6 6
6
C 1 0
and NPURE = 1 = number of pure phase; IPURE (I) = 4 = column number of C. (3) Set the AMU array. The AMU array is set equal to the free energies of formation of the components a t the temperature of interest, divided by RT, in consistent units. In the example, the free energies of the components a t 500 K in cal/g-mol are listed in Table 111, in which the last column is used as the AMU array for the corresponding components. I t is imperative t h a t t h e free energies of formation be divided by RT. IPHASE in the example is 2. The computer results are now summarized in Table IV. The computer predictions are valid only when the basic chemical and physical-equilibrium relations are satisfied. In other words, the equilibrium constant must satisfy the free energy data; second, the fugacity of each component must be equal in the two phases. For example 1: 1. K,,, =
fcyc
- e-lGn/RT
fByB (fH.)“l3 where f c = cyclohexane vapor phase fugacity; f B = benzene vapor phase fugacity; f H = hydrogen vapor phase fugacity;y c Ind. Eng. Chern., Process Des. Dev., Vol. 15, No. 3, 1976
375
Table V
Table VI
Component EtOH HAC EtAc HOH Total moles
Liquid mole fraction
Vapor mole fraction
0.048 419 0.230 446 0.145 165 0.575 969 4.065 x 10-lo
0.0753 58 0.0753 58 0.424 642 0.424 642
Component
PYL
XiXcPi
Xi
EXOH HAC EtAc HOH
0.075 358 0.075 358 0.424 642 0.424 642
0.085 584 0.079 694 0.538 780 0.434 367
1.37 1.05 2.84 1.32
10.000 000
Table VI1 = vapor mole fraction of cyclohexane; Y B = vapor mole fraction of benzene; Y H = vapor mole fraction of hydrogen.
2. Yifi = Xifi where xi = liquid mole fraction and f i = liquid phase fugacity. For the first set, K,,, for the equation is K,,, = e5.22= 185. Checking the computer predictions Krxn =
0.923 452 X 16.2 X 18.6 X (30 X 0.076 183 9)3
0.364 384 X
Component EtOH HAC EtAc HOH Total moles (for 10 mol feed)
Liquid mole fraction
Vapor mole fraction
0.183 508 0.401 819 0.136 300 0.278 373 3.699 000
0.253 846 0.125 697 0.351 927 0.268 529 6.301 000
= 185.1
For the second set Benzene Liquid: 0.387 615 X Vapor: 0.364 384 X
X 17.5 = 0.006 768 X 18.6 = 0.006 778
Cyclohexane Liquid: 0.997 328 X 15.0 = 14.96 Vapor: 0.923 452 X 16.2 = 14.96 Hydrogen X 1000 = 2.285 Liquid: 0.228 54 X Vapor: 0.076 183 9 X 30 = 2.286 No other equations need be satisfied since both chemical equilibrium (Set 1)and physical equilibrium (Set 2) are satisfied. The computer predictions would, therefore, be as valid as the data entered for input. Example 2. While the computation of equilibrium compositions is by no means new, it certainly is neither routine nor simple. How easy it is to err may be demonstrated by examining a recent example taken from Sanderson and Chien (1973). “Ethanol is reacted with acetic acid to form ethyl acetate and water. EtOH
+ HAC
EtAc
+ HOH
Five moles each of ethanol and acetic acid are fed to a reactor at 85 “C and 1 atmosphere.” (Relevant thermodynamic data are taken from Sanderson and Chien (1973) and include: (1) vapor phase was assumed ideal; (2) vapor pressures were calculated using Antoine equation; (3) liquid-phase activity coefficients were computed from Wilson parameters. The present program was used with the assumption that both liquid and vapor phases are present and yielded the equilibrium conditions listed in Table V. This reaction should be considered one phase and the equilibrium calculations recast and recomputed on this basis. However, if the two-phase description is continued
zx,= zy, = 1 YLP = XLYlPl
and
376
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 3, 1976
The present program provides the checks listed in Table VI on the vapor-liquid equilibrium. While this vapor-liquid equilibrium agreement is quite approximate, it is satisfactory because of the negligible quantity of liquid present. The reaction rate constant for the vapor phase, K,,,, is computed from the free energy data reported by Stull et al. (1969) and found to be 31.75. (Alternatively, K,,, could be estimated by linearly interpolating K,,, data also supplied by Stull-this yields 32.42) The chemical reaction constraint may now be used to check the computed results of the present program Krxn =
(0.424 642) (0.424 62) = 31.75 (0.075 358) (0.075 358)
This agreement is excellent. The material balance restraints are verified by the total vapor moles because the net change of moles in the reaction is zero. The present equilibrium compositions do not agree with Sanderson and Chien, who reported the data listed in Table VII. While Sanderson and Chien’s results satisfy the material balance and physical equilibrium restraints, they do not satisfy the chemical reaction constraint. Using Sanderson and Chien’s results
The chemical reaction constraint is thus violated by a factor of 10, thus invalidating the computed equilibrium results. A detailed discussion is given by George (1973). Conclusions The present goal was to develop a user-oriented, computer program or “black-box” which would predict multicomponent, multiphase equilibrium compositions. Development of improved optimization techniques lies outside the present scope. The key to the present computational procedure is Brown’s (1971) clever rearrangement of the Gibbs free energy into an unconstrained form which facilitates the minimization calculations. The applicability of the present method has been established to include the following items. Isothermal and adiabatic cases can be solved in addition to chemical and physical equilibrium. The program can handle up to three multicomponent phases, and 15 components. Computer requirements for execution of the program are 20K storage using double precision and less than 2 min CPU time on a Sigma 6. The
abundance of trace components can be estimated following White et al. (1958) whether or not the component is part of the free energy minimization. For example, after the equilibrium compositions in a mixture of C02, H20, and NO have been computed, then the trace concentrations of CO, C302, N20, and NO2 in this exhaust can also be estimated. The program is user-oriented and can be run after a minimum of instruction. The present computational procedure has been applied successfully to the following typical chemical and physical processes. Seven-component vapor-liquid equilibrium, the typical "flash calculation. The two-phase hydrogenation of benzene (Smith and Van Ness, 1959) (a two-phase chemical reaction). The two-phase, nonideal binary mixture carbon tetrachloride and 1,1,2-trichloroethane with activity coefficients (Kissell and Manning, 1962). Three-phase equilibrium flash vaporization posed by Osborn (1964), and solved by Deam and Maddox (1969). Cracking of butane to equilibrium a t 700 K. Zinc smelter problem. Can liquid zinc be formed a t 1100 K and 1atm? Reaction of ethanol and acetic acid (Sanderson and Chien, 1973). Adiabatic reaction of gases from a pyrite burner (Hougen et al., 1959). Adiabatic flame temperature of methane (Smith and Van Ness, 1959). Combustion of methane with two different stoichiometric amounts of air which is an equilibrium check for NO, formation. Combustion of high sulfur Bunker C fuel oil. Notation A = matrix formed with atoms listed vertically, components horizontally a,, = elements of A matrix B, = total amount of j "atoms" in mixture C = coefficient in chemical reaction equation f l = fugacity of component i in solution, pressure units f l o = fugacity of component i in standard state, pressure units f l k =. fugacity of component i in kth phase solution, pressure units F = moles of feed G = total Gibbs free energy, cal GI = partial Gibbs free energy, cal/mol G I o = standard state Gibbs free energy, cal/mol AGO = standard state free energy change, cal/mol AGf = free energy of formation, cal/mol AH = total enthalpy change, cal H ( u , ) = allocation function K , = phase equilibrium constant K,,, = chemical reaction equilibrium constant L = number of phases; total liquid moles M = number of atomic species in mixture N = number of components in mixture N, = numberofphases N I = number of independent components N R = number of independent reactions P = total pressure p , = vapor pressure of component i R = ideal gas constant T = temperature, K u = unconstrained variable x, = mole fraction of component i in liquid phase X I k = moles of component i in k t h phase y, = mole fraction of component i in vapor phase z, = mole fraction of component i in mixture 2, = moles of component i in mixture
Greek Symbols A; = activity coefficient for component i pi =
fugacity coefficient for component i
Computer Notation A = the material balance matrix ALIQ = the total moles of the first phase ALIQl = the total moles of the fourth phase AMU = the Gibbs free energy of formation, at the pressure and temperature of the system, divided by RT ARHS = the portion of the A matrix on the "right hand side" of the equal sign ASOL = the total moles of the third phase AVAP = the total moles of the second phase B = the total atoms or reaction groups in a system FL = the first phase fugacity FV = the second phase fugacity FW = the third phase fugacity IPHASE = the total number of phases present IPRNT = the print out option IPURE = the position of the pure components in the 2 array K = the number of multicomponent phases M = The number of atoms or reaction groups in the A matrix MAX = the maximum number of iterations to be made N = the number of components NPURE = the number of pure component phases PI = the array used for estimation of trace components V = the mole fraction of the fourth phase W = the mole fraction of the third phase X = the mole fractions of the first phase Y = the mole fractions of the second phase Z = the total moles of feed L i t e r a t u r e Cited Amundson. N. R.. "Mathematical Methods in Chemical Engineering: Matrices and Their Application," Prentice-Hall, Englewood Cliffs, N.J.. 1960. Box, M. J.. Comput. J., 9,67-77 (1966). Brown, L. P., Ph.D. Dissertation, The University of Tulsa, 1971. Dayhoff. M. O.,Lippicott, E. R., Eck, R. V., Nagarajan, G., NASA SP-3040, Washington, D.C. (1967). Deam, J. R., Maddox, R. N.. Hydrocarbon Proc., 48 (7), 163-164 (1969). Dluzniewski, J. H., Adler, S. E.. l. Chem.E. Symp. Ser., 4, No. 35, 421-26(1972). Fletcher, R.. Reeves, C. M., Comput. J., 7, 149-54 (1961). Fletcher, R., Powell, M. J. D., Comput. J., 6 (2), 163-168 (1963). George B., Ph.D. Dissertation, The University of Tulsa, 1973. Hooke, R., Jeeves. T. A., J. Assoc. Comp. Mach., 8, (2), 212-229 (1961). Hougen, 0.A., Watson, K . M., Ragatz, R. A,, "Chemical Process Principles, Part II: Thermodynamics," 2d ed. Wiley, New York, N.Y., 1959. Kissell, F. N.. Manning, F. S., J. Chem. Eng. Data, 7,205-206 (1962). Ma, Y. H., Shipman, C. W., AlChE J., 18 (2), 299-304 (1972). Nelder, J. A,, Mead, R., Comput. J., 7,308-13 (1964). Osborn. A,, Chem. Eng., 71 (12), 97-100 (1964). Powell, M. J. D., Atomic Energy Research Establishment Report AERE-R.6469 (1970). Sanderson, R. V., Chien, H. Y.. lnd. Eng. Chem., Process Des. Dev., 12,81-85 (1973). Smith, J. M.. Van Ness, H. C., "Introduction to Chemical Engineering Thermodynamics," 2d ed, McGraw-Hill, New York. N.Y., 1959. Stull, D. R., Westrurn. E. F., Sinke, G. C., "The Chemical Thermodynamics of Organic Compounds," Wiley. New York, N.Y., 1969. White, W. B., Johnson, S. M..Dantzig, G. E., J. Chem. Phys.. 28 (9,751-755, (1958).
Receiued f o r reuiew M a r c h 6, 1975 Accepted F e b r u a r y 17, 1976 Supplementary Material Available: E n t i r e computer program for t h e e q u i l i b r i u m computation (23 pages). Ordering i n f o r m a t i o n is given on a n y c u r r e n t masthead page.
Ind. Eng. Chem., Process Des. Dev., Vol. 15, No. 3, 1976
377
