652
Ind. Eng. Chem. Process Des. Dev. 1985, 2 4 , 652-658
Vapor-LiquM-Liquid
Equilibria Computations
VlJay R. Sampath and Stuart Leiprlger’ Gas Englneering Department, Instltute of Gas Technoiogy/Illinois Instltute of Technobgy, Chlcago, Illinois 606 16
A simplified methodology has been developed for carrying out dew and bubble point calculations and flash calculations for vapor-liquid-liquid equilibrium systems. The methodology applies to systems containing water, and in the methodology, one employs the assumption that the activity of water in the aqueous liquid phase is unity. This assumption can effect a large reductlon in computational time without a significant loss in accuracy. Three-phase equilibrium calculations employing this methodology were carried out for the benzene-water system.
Introduction
zi
The objectives of vapor-liquid-liquid equilibria computations are to determine which of the three phases are present and to determine the composition of the phases. Henley and h e n (1969) presented a rigorous three-phase flash calculation scheme which is an extension of the conventional two-phase flash, bubble point, and dew point calculation methodologies. Erbar (1973) published a three-phase calculation method for water-hydrocarbon systems. His algorithm is similar to that of Henley and Rosen. He calculates K values using a modification of the Chao-Seader method which requires two different values of the solubility parameter for water, one of which is the “true” value for the hydrocarbon liquid phase and the other which is obtained by “inspection” for the aqueous liquid phase. Dluzniewski and Adler (1972), Soares et al. (1982), and Heidemann (1974) proposed a three-phase flash calculation methodology based on the direct minimization of Gibbs free energy. In the absence of chemical reaction, this approach has no particular advantage over the conventional flash calculation method. Heidemann (1974) and Peng and Robinson (1976) used an equation of state to perform three-phase equilibrium calculations. In order to successfully describe three-phase equilibria, an equation of state has to be correlated with temperature-dependent binary interaction parameters. Mauri (1980) published a detailed three-phase calculation algorithm based on the Henley-Rosen method that considers all possible phase conditions. Henley-Rosen Computational Scheme
The method proposed here is a simplification of the Henley-Rosen scheme and is applicable to hydrocarbonwater systems. We define n
CX,H
=1
(1)
r=l
(2)
?yk = r=1
5xrw = 1
(3)
a = V/F
(4)
k=1
p = LH/(LH
+ LW)
(5)
The overall and individual component balances can be written as (6) F=V+LH+LW
* Present address: Chemical Engineering Department, RoseHulman Institute of Technology, Terre Haute, IN 47802. 0196-4305/85/1124-0652$01.50/0
= p(1 - a)x;H + (1 - a)(l - P ) X i W
+ ayi
(7)
The individual phase compositions can be computed from the equation xiH = ~ i / [ @ ( l - a) + (1- a)(l - @)KiH/KiW + y ; = KiHziH
xiw =
(8) (9)
KiHxjH
-
KiW The criteria for vapor-liquid-liquid equilibrium can be written in terms of eq 1 through 3 as N N CXiH - c y i=l
i=l
N
N
i =0
c x i w - cy1 = 0 i=l
(12)
i=l
N N CXiH - CXiW i=l
(11)
=0
(13)
i=l
Equations 11through 13 must be satisfied such that 0 I a I1 and 0 5 p I1 for three phases to coexist. In a two-phase calculation, only one of eq 11 through 13 is required and in a three-phase calculation, two equations chosen from among eq 11 through 13 are required. For a flash calculation we solve for phase indicators a and 8; in a bubble point calculation, we solve for the bubble point temperature or pressure and phase indicator p; and in a two-phase dew point calculation, we solve for the dew point temperature or pressure and phase indicator 0. In a three-phase dew point calculation, we solve for the dew point temperature or pressure and the phase indicator, a. Improved Computational Scheme
In our calculations we employ an assumption which effects a reduction in the number of required constraint equations (11through 13) for a three-phase equilibrium calculation from 2 to 1. This reduction in the number of constraint equations allows us to use a one-dimensional Newton-Raphson convergence scheme which has the advantage of significantly reducing computation time. Since many three-phase systems encountered in industry contain water as one of the components and water tends to form largely immiscible liquid phases, we make the assumption that the activity of water in the aqueous liquid phase is unity. In many systems, particularly hydrocarbon-water systems, the mole fraction of water in the aqueous liquid phase is
0 1985 American Chemical Society
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 3, 1985 653
, (-)
close to unity. If one chooses the pure liquid as the reference state for water in the aqueous liquid phase, the activity coefficient of water is close to unity. When a hydrocarbon liquid is added to water, hydrogen bonds are broken and no new bonds are formed. This results in an increase in the entropy of mixing. Hydrocarbon-water systems therefore exhibit positive deviations from Raoult’s law. This means that yWHzO
-
-
I+ as xwHz0
1-
ASSUME a = 0.33 5: flal
ASSUME IDEAL GAS, IDEAL DILUTE LIQUID SOLUTIONS
(15)
Hence the product uWH& = y W ~ @ W stays ~ pclose to unity even as y W ~ zand O xwHzodepart slightly from unity. This assumption is applicable when two liquid phases are present. The assumption in eq 14 breaks down as we approach the three-phase critical solution end point. Substituting eq 14 into eq 8, 9, and 10 yields
P=
Z H ~ O Y ~ H-~ aO( f l ~ , o-
1)-1
(1- a ) ( p H z ~ / K H H 2 0- 1)
1
CALCULATE K-VALUES
4
CALCULATE PHASE COMPOSITIONS CALCULATE 1 x 7 z y ,
-
EX7 - Zx,w zxw
-1Y,
1
(16)
Equation 16 is valid for 0 < /3 < 1. Since /3 is known as a function of a,we need only solve the appropriate equation from among eq 11through 13 to obtain a for a flash calculation and temperature or pressure for a bubble point or dew point calculation. Knowledge of a,p, and the temperature and pressure will allow us to compute all phase compositions from eq 8 through 10. Flash Calculations. A procedure for two- and threephase flash calculations is shown in Figure 1. Initially, three equilibrium phases are assumed to exist and a search is carried out for values of a (or p) which satisfy one of eq 11through 13. Equations 11and 12 are used to search for a value for a and eq 13 is used to search for a value for 8. Values for a and P are used to determine which phases are present. The sum of the phase compositions, CxiH, Cxiw,and Cy, are used to determine which one of the three equations is to be employed. When Cyi has the lowest value of the three, we search for the value of p that satisfies eq 13. When CxiHis greater than Cxiw,we search for the value of a that satisfies eq 11. When Cxiw is greater than CxIH, we search for the value of LY that satisfies eq 12. The search is carried out by a NewtonFtaphson convergence scheme. The other phase indicator /3 (or a)is then determined from eq 16. If the calculated value of a or P is outside the bounds of 0 and 1, three equilibrium phases do not exist and we then set the phase indicator at the boundary value and the compositions of the two phases can be computed from the same set of eq 8 through 10. If however cy > 1or a = 0 and is outside the bounds, only one equilibrium phase exists and in that case none of the criteria of equilibrium will be satisfied. A solution is attained when the objective function is minimized. The phase composition can be calculated from eq 8 through 10. For binary systems, both a and p cannot be found within the bounds of 0 and 1except at a three-phase equilibrium point. Therefore when a exists within bounds, we set P = 0 or 1 and when exists within bounds, we set cy = 0. Bubble and Dew Point Calculations. The bubble point and dew point calculation procedure is shown in Figure 2. Initially we assume that three equilibrium phases are present. For two-phase dew point calculations we set cu = 1and set /3 = 0 or 1. For three-phase dew point calculations, we set P = 0 or 1depending on whether we are concerned with the incipient hydrocarbon liquid or the aqueous liquid. We then solve for a from eq 16. For bubble point calculations we set a = 0 and calculate p from eq 16. For binary systems, we set p = 0 or 1. If the calculated value of P falls outside the bounds, three
CONVERGED?
a , 8 CONVERGED?
a
OBJECTIVE fcn T X i ” - ZX;” = 0
BINARY?
t I I PRINT RESULTS
/
/
Figure 1. Algorithm for flash calculations.
equilibrium phases do not exist. We set it at the boundary value. We then search for the bubble or dew point temperature or pressure that satisfies either eq 11 or 12. A Newton-Raphson convergence scheme was employed. Initial Estimates. In most three-phase calculation methods, good initial estimates are necessary to avoid a trivial solution. Therefore, most multiphase calculation methods use elaborate procedures to develop initial estimates for the phase distribution and compositions. Since this algorithm is more like that of a two-phase equilibrium calculation, the convergence of the computations is not strongly dependent on the initial estimates. We arbitrarily start all flash calculations with a = 0.33 and calculate j3 from eq 16. This has been sufficient in all cases tested. For bubble or dew point calculations, we include an initial estimate in the data input for the temperature or pressure. It was noticed that convergence is not sensitive to the initial estimates of temperature and pressure. We initialized the computation by assuming ideal gas behavior for the vapor and ideal dilute liquid solution behavior for the liquid in order to calculate the initial set of K values. Convergence. A Rachford-Rice (1952) objective function that is monotonically increasing was used for flash calculations. A one-dimensional Newton-Raphson con-
654
I d . Eng. Chem. Process Des. Dev., Vol. 24, No. 3, 1985
in the equations are determined from a modified UNIQUAC equation. Values for activity coefficient parameters and the vapor phase binary interaction parameter were determined from data. The vapor phase fugacity coefficient diVand 4' were calculated from the modification of DeSantis et d. (1974)of the Redlich-Kwong equation. Henry's constants were represented by a four-parameter empirical equation fitted to data. Pure component molar volumes were used in place of partial molar volumes. The three-phase calculation algorithm is independent of the K value method given above and any other method can be substituted. Results and Discussisn Equilibrium calculations were earried out for the benzene-water system using values for equation parameters developed froq a data reduction of the liquid-liquid equilibrium data of Thompson and Snyder (1964). The calculations were based on the flash calculation algorithms and dew and bubble point algorithms. As discussed earlier, the algorithms are based on s i m p l i i assumptions which significantly reduce computational time. Figure 3 shows the variation of the three-phase locus of the mole fractions of water in the vapor and hydrocarbon liquid phases with temperature for the benzene-water system. At the upper critical solution temperature, 541.5 K, the composition of the saturated vapor and the saturated liquid become identical. At this temperature, the vapor and hydrocarbon liquid phases become indistinguishable; above this temperature, three phases cannot coexist. A comparison between the calculated curves and the data of Rebert and Kay (1959)is shown in the figure. Figure 4 shows the variation of the three-phase pressure with Fmperature up to the critical solution point. The DeSantis-Breedveld-Prausnitz modification of the Redlich-Kwong equation does not require any binary interaction parameters to calculate vapor phase fugacity coefficients. However, a vapor phase binary interaction parameter may be introduced into the equation of state in the traditional manner. The energy interaction term ai,is modified as y.*
FOR BUBBLE POINT CALC SET a = 0 FOR TWO PHASE W W P C + ~CALC T SET a = I FOR THREE PHASE DEWPOINT C A k , SET B i O/I I
ASSUME IDEAL GAS,
CALCULATE K-VALUES
4
CALCULATE PHASE COMPOSITIONS CALCULATE Ex Zy,
-
zxi" -zY,
I
JECTIVE fcn
~~
Figure 2. Algorithm for bubble and dew point calculations.
vergence scheme based on numerically determined derivatives was used for the flash and dew and bubble point calculations. Flash calculations required an average of 15 iterations with a maximum of 35. Dew and bubble point calculations required an average of 10 iterations with a maximum of twenty iterations to converge. K Values. K values are computed by a procedure described by Sampath (1981). In the thermodynamic framework proposed, separate equations were employed to describe liquid and vapor phase behavior and to compute K values. Furthermore, the pure liquid fugacity was chosen as the reference fugacity for all solvent species and Henry's constant was chosen as the reference fugacity for all solute species. This applies for vapor-liquid, liquidliquid, and vapor-liquid-liquid equilibrium systems.
K, =
Y1
- = u,P?y,
exp[(P? - P ) V ? / R q / & " P (17)
Xl
Yi
Kl = - = Hji-yl* exp[P?7jL/Rr]/4jvP
(18)
XI
The subseript j refers to solvent species and subscript j refers to solute species. The vapor pressures in eq 17 were estimated from a curve-fit of pure component data. The pure liquid fugacity coefficient ui was obtained from the correlation of Lyckman et al. (1965).The asymmetric activity coefficients yi and
a[] =
(uiul)l/ql- k,)
(19)
where ki is a vapor-phase binary interaction parameter. We empioyed kij as a fitting parameter. The dashed line in Figure 4 corresponds to the DeSantis-BreedveldPrausnitz equation without the binary interaction parameter (kij = 0); the unbroken line corresponds to the DeSantis-Breedveld-Prausnitz equation with a value of k, = -0.8. The value of -0.8 is large, but its use results in a significant improvement in the prediction of three-phase pressures. As may be expected, at the lower pressures where ideal gas behavior is expected, the differences between the two cases are small. Since inclusion of a vapor-phase binary interaction parameter significantly improves the prediction of pressure (and also vapor phase composition),a value of kij = -0.8 has been utilized in the predictions to follow on the benzene-water system. Rebert and Kay (1959)have experimentally determined the phase boundaries of the benzene-water system over fairly wide ranges of temperature and pressure. Their experimental resulta are compared with predicted results in Figures 5,6,7,and 8. The predicted phase boundaries shown in these figures were obtained with the dew and bubble point pressure algorithms and the two- and three-phase flash calculation algorithms. Figure 5 is a pressure-composition diagram for the benzene-water system at 498 K. The comparison with experimental data is good. Figure 6 is a similar preasure-composition diagram
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 3, 1985 855
BENZENE-WATER SYSTEM THREE-PHASE LOCUS
560
0
SATURATED VAPOR,EXPERIMENTAL
A
SATURATED LIQUID,EXPERIMENTAL
- PREDICTED 540-
I
X
1
520-
W
5 GK L
BI-
5”
480.
460 -
40
BENZENE-WATER SYSTEM THREE-PHASE LOCUS
___
0 EXPERIMENTAL
440
IO
PREDICTED (kii i -0.81 PREDICTED (kij = 0.0)
I
20
I
30
40
50
60
70
BO
90
100
0
PRESSURE (atm)
Figure 4. Three-phase temperaturepressure diagram for benzene-water system.
at a higher temperature. A comparison of Figures 5 and 6 shows that the vapor-hydrocarbon liquid region has diminished considerably. Figure 7 shows the phase boundaries at the experimental three-phase critical point at 541.5 K. As can be seen, the vapor-hydrocarbon liquid region has diminished further and it is on the verge of disappearance. Figure 8 contains experimental data and calculated values at 568 K which is above the three-phase critical point of the binary and also above the critical temperature of benzene. Both the temperature and the pressure of the data in Figure 8 are high and the comparison with experimental data is poor. It should be kept in mind that all of these predictions are based on binary parameters reduced from liquid-liquid
equilibrium data over a temperature range of 311 to 477 K. Large extrapolation of Henry’s constants can lead to large uncertainties in predicted results, particularly at temperatures where Henry’s constants change very rapidly with temperature. The computational time required in the “shortcut” algorithms which employs the assumption that OWHzO = l, is significantly reduced without a loss in accuracy. Table I shows a comparison of computational results and computational times for three-phase equilibrium calculations for the benzene water system. The first column contains the system temperature; the second column contains the three-phase pressure as calculated by a rigorous threephase computational algorithm; the third column contains the three-phase pressure as calculated by the “shortcut”
658
Ind. Eng. Chem. Process Des. Dev., Vol. 24,No. 3, 1985 BENZENE -WATER SYSTEM
-A
SATURATED LIOUID. EXPERIMENTAL PREDICTED
LW
L ~ L~ -
V
10
0.I
0.0
0.2
0.9 QS Q6 MASS FRACTION WATER
0.3
0.7
0.8
0.9
Figure 5. Phase boundaries of benzenewater system (498.2K). 90
-
L" L~
LH
80 -
LW
BENEE-WATER PHASEBWNMRIES
o A
-
30-
20
00
0.1
I
I
0.2
0.3
S A ~ R A T & W & . EXPERIMENTAL SATURATED UPUID. EXPERIMENTAL PREDICTED
I 0.4
0.5
(16
0.7
0.8
0.9
MASS FRACTION WTER
Figure 6. Phase boundaries of benzene-water system (528.2K). Table I. Comparison of Preclioted Three-Phase Pressures press. calcd temp, by rigorous K method, atm 473.2 30.00 493.2 43.03 503.2 50.86 513.2 59.47 523.2 68.55 528.2 73.05 533.2 77.34 538.2 81.26 541.5 83.56 rme error 4.56 computational time, s 0.54
press. calcd by bubble point algoriothm, atm 30.0 43.0 51.0 60.0 71.0 76.0 82.0 81.0 83.0 4.21 0.30
press. calcd by dew point algorithm, atm 30.0 43.0 51.0 60.0 71.0 73.0 78.0 81.0 83.0 4.61 0.40
exptl. press 29.80 43.07 51.17 60.22 70.63 76.42 82.47 88.94 92.81
Ind. Eng. Chem. Process Des. Dev., Voi. 24, No. 3, 1985
657
IO0 Ln- Lw
90
-
A
6
-
----
3
so -
A M
-
A
0
E
c
0
0 y v- L"
v-LW
0
w
s
0
m
Y
0
0
v
60-
BENZENE-WATER SYSTEM PHASE BOUNDARIES T = 541.5 K (3-PHASE CRITICAL POINT)
50 -
MASS FRACTION WATER
Figure 7. Phase boundaries of benzene-water system (541.3 K). A
B
-
A LW
125 -
E
L
O
Y3
115 -
v) v)
w
a a
V
105-
95
-
85 -
75
ao
0
-A 0.I
BENZENE- WATER STSTEM PHASE BOUNDARIES T=568.2K SATURATED VAPOR, EXPERIYENTAL SATURATED LIQUID. EXPERIMENTAL PREDICTED
0.2
0.3
0.4
a5
0.6
0.7
0.8
0.9
MASS FRACTION WATER
Figure 8. Phase boundaries of benzene-water system (568.2 K).
bubble point pressure algorithm; the fourth column contains values of the pressure calculated by the "shortcut" dew point pressure algorithm; and the fifth column contains experimentallymeasured three-phase pressures. The savings in computational time become particularly significant for multicomponent systems. Computations were carried out on an IBM 370/158 computer and computational times for each of the techniques are presented in Table I.
Conclusions A "shortcut" methodology has been developed for carrying out three-phase equilibrium calculations for hydrocarbon-water systems. This methodology differs from the rigorous three-phase equilibrium calculation methodology
proposed by Henley and Rosen (1969) in the respect that it employs the assumption that the activity of water in the aqueous liquid phase is unity. This assumption permits the reduction of a two-dimensional three-phase calculation scheme to a one-dimensional three-phase calculation scheme. It was demonstrated that this assumption can effect a large reduction in computational time without a significant loss in accuracy. The "shortcut" methodology has been coded into computer programs. These programs can be used for two- and three-phase isothermal flash calculations, two- and three-phase bubble and dew point temperature calculations, and two- and three-phase bubble and dew point pressure calculations; for example, they calculate the temperature or pressure at which an incipient equilibrium
Ind. Eng. Chem. Process Des. Dev. lB85, 2 4 , 658-665
658
vapor phase appears in a system where two liquid phases are in equilibrium. This type of calculation allows the development of the phase boundaries for a multicomponent system and the unique three-phase locus for a binary system. Nomenclature a = activity F = number of moles of mixture H = Henry’s constant K = equilibrium ratio L = number of moles of liquid N = number of components P = pressure R = gas constant T = temperature = number of moles of vapor Vi= partial molar volume x = liquid phase mole fraction y = vapor phase mole fraction z = overall feed mole fraction a = vapor mole fraction and phase indicator p = phase indicator (see eq 5 ) y = activity coefficient v = pure liquid fugacity coefficient 4 = vapor phase fugacity coefficient
j = solute species
Subscripts i = solvent species
Acknowledgement is made to the Institute of Gas Technology for support of this work.
Superscripts
H = hydrocarbon rich phase S = saturated W = water rich phase * = normalized in the asymmetric convention
Literature Cited DeSantk, R.; Breedveld, 0. J. F.; Prausnitz, J. M. Ind. f n g . Chem. Process Des. Dev. 1974, 13, 374. Dluzniewski, J. H.; Adler, S. B. Chem. Eng. Symp. S e r . 1872, 35. Erbar, J. H. 52nd Annual Conventlon, Natlonal Gas Processors Association, Dallas, 1973. Hekiemann, R. A. A I C M J. 1974, 20, 847. Henby, E. J.; Rosen, E. M. “Material and Energy Balance Computations”; Wiley: New York, 1969. Lyckman, E. W.; Eckert, C. A.; Prausnltz, J. M. Chem. f n g . Sci. 1985, 20, 685. Mauri, C. I d . Eng. Chem. Process Des. Dev. 1980. 19, 482. Peng, D.-Y.; Robinson, D. 8. Can. J. Chem. €ng. 1876, 5 4 , 595. Rachford, H. H.; Rlce, J. D. J . Pet. Technoi. 1852, IO, 19. Rebert, C. J.; Kay, W. B. AIChE J. 1959, 5, 285. Sampath, V. R. Ph.D. Thesls, Gas Engineerlng Department, Illinois Institute of Technology, Chicago, IL, 1981. Soares, M. E.; Medlna, A. G.; McDermott, C.; Ashton, N. Chem. fng. Scl. 1982, 37, 521. Thompson, W. H.; Snyder, J. R. J. Chem. Eng. Data 1884, 9 , 516.
Receiued for reuiew November 28, 1983 Accepted August 23, 1984
Computation of Vapor-Liquid Equillbria Behavior of Coal Gasificatlon Systems Srlnlvaran Ramanujam, Stuart Lelprlger,” and Sanford A. Well Gas Engineering Department, Institute of Gas Technobgy/IHinois Institute of Technoiogy, Chicego, Iliinois 606 16
A computational technique has been developed to predict vapor-liquid equilibria properties for coal gasification systems. Gas-phase behavior is described by the Prausnltz-Chueh modification of the Redlich-Kwong equation of state and Ilquld-phase behavior is described by the regular solution theory of HiMebrand and Scatchard with a minor modification. The technique Invokes introduction of a solubility 6,. The technique was successfully applied
to several binary systems and a few ternary and quaternary systems.
Introduction Design and operation of modem coal conversion facilities require appropriate thermodynamic and thermophysical property data, and VLE data are among the information required most. Some of the crucial data needs have been discussed in detail by Weinstein (1977). Good predictive techniques must be developed to use with the available VLE data. One should bear in mind that the technique must not only be applicable for elevated pressures and temperatures but also for liquid mixtures that are primarily aromatic. The nature of liquids obtained from coal conversion proceases is dependent on the process as well as the chemical constituents of coal. Often times, the characterization of coal derived liquids is difficult. Sebastian et al. (1981a, b, c) proposed a methodology to predict the gas solubilities in coal liquids at high
* Chemical Engineering Department, be-Hulman Institute
of Technology, Terre Haute, IN 97802.
0196-4305/85/1124-0658$01.50/0
pressures and temperatures. Correlations were developed for the estimation of the fugacity of noncondensable gases in the liquid phase in terms of solubility parameters. However, the applicability of their methodology has not been established in terms of generalization of the vaporphase interaction parameters for a wide spectrum of components and application of the correlations to systems containing mixtures of noncondensable components should be established. Therefore, we propose a somewhat similar technique that addresses the above questions. This methodology has been developed primarily to verify the VLE data of coal hydrogasification liquids. Description of Vapor Phase Behavior The well-known Prausnitz-Cheuh (1968) modification of the Redlich-Kwong (RK) equation of state, because of ita simplicity and ease of generalization,has been used for the computation of the vapor phase fugacity. The RK equation of state is given by 0 1985 American Chemlcal Society
