Ind. Eng. Chem. Res. 2000, 39, 1497-1504
1497
A Molecular-Based Model for Normal Fluid Mixtures: Perturbed Lennard-Jones Chain Equation of State Y. P. Lee, Y. C. Chiew,* and G. P. Rangaiah Department of Chemical & Environmental Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore
The perturbed Lennard-Jones chain (PLJC) equation of state, previously used for pure fluids, is extended to normal fluid mixtures. The PLJC equation utilizes the hard-sphere chain mixture as the reference system, while the perturbation term is derived based on the first-order variational theory. Binary vapor-liquid equilibrium data of a wide range of fluids were well described by this model through the use of one binary parameter kij. We showed that the PLJC equation of state was able to predict vapor-liquid equilibrium fairly accurately using binary data. Also, the PLJC equation is able to model retrograde condensation and liquid-liquidvapor three-phase equilibrium of ternary mixtures without the need to reevaluate binary parameters. Theoretical phase equilibrium compositions were determined through an isothermalisobaric flash calculation based on global minimization of the Gibbs free energy using a genetic algorithm. 1. Introduction Prediction of fluid-phase equilibria has always played a crucial role in the design and simulation of separation processes, such as distillation and extraction, in the petrochemical and natural gas processing industries. Over the years, different thermodynamic models have been developed and reported in the literature; these include activity coefficient models that describe liquidphase nonideality and equation of state theories that capture fluid behavior over wide pressure and temperature ranges. Equations of state have found success in modeling high-pressure systems including vapor-liquid equilibria (VLE), solute solubility in supercritical solvents, etc. In recent years several theoretically based equation of state theories for chainlike fluids have been developed.1-5 These models were formulated based on statistical mechanical perturbation theories of liquid state and provide analytical expressions for the Helmholtz energy and compressibility factor. These analytical equations of state take into account contributions from repulsive (reference) and dispersion attractive (perturbation) forces. In these theories, the hard-sphere chain reference fluid consists of a series of linearly jointed hard-sphere segments that exert repulsive forces on other segments. Attractive interactions between chain segments are accounted for by introducing a perturbation on the reference system. In this work we propose the perturbed Lennard-Jones chain (PLJC) equation of state for normal nonassociating fluid mixtures. This is an extension of the work of Chiew et al.5 in which the PLJC equation was derived and applied to pure fluids. The theoretical foundation of the PLJC equation lies in the variational theory of Lennard-Jones chains by O’Lenick and Chiew,6 which was later extended to Lennard-Jones chain mixtures by Von Solms et al.7 and Von Solms and Chiew.8 Similar to other molecular-based chain fluid equation of state * To whom correspondence should be addressed. E-mail:
[email protected]. Fax: (65) 779 1936. Tel.: (65) 874 8044.
models, the PLJC equation consists of three independent parameters, namely, m, σ, and which represent the number of chain segments, segment size, and segment energy, respectively. Values of m, σ, and for a wide range of nonassociating molecules and polymers have been regressed from experimental data and tabulated by Chiew et al.5 The PLJC model is defined in terms of the residual Helmholtz energy, and other thermodynamic properties such as compressibility factor and fugacity coefficients can be readily obtained via thermodynamic relationships. In the mixture PLJC model, the molecular segment based parameters ij for dissimilar molecules follow the Lorentz-Berthelot combining rule with one additional parameter kij, that is, ij ) xiijj(1 - kij). The binary parameter kij is obtained by regressing the model against VLE data for binary mixtures. We adopt a rather unconventional but novel approach to solve phase equilibria calculations. Instead of the traditional practice of solving simultaneously a set of governing equations related to mole balances and equality of fugacities, we minimized the total Gibbs free energy of the system with respect to the mole numbers of the constituent components in all phases under isothermal-isobaric condition. The chief benefit of this approach is that the cumbersome derivation of expressions for fugacities can be avoided and consequently results in simpler problem formulation and easier implementation. In the case of the model considered here, the PLJC mixture theory provides an analytical expression for the Helmholtz energy and hence an analytical expression for the Gibbs free energy. From a mathematical viewpoint, the equation-solving approach at most can only satisfy the necessary or equipotential conditions of phase equilibrium, whereas the Gibbs energy minimization method satisfies both the necessary and sufficient conditions of phase equilibrium.9 In section 2, we present analytical expressions of the Helmholtz energy and compressibility factor for multicomponent mixtures based on the PLJC model. Section 3 briefly describes the regression procedure, in particu-
10.1021/ie990744p CCC: $19.00 © 2000 American Chemical Society Published on Web 04/07/2000
1498
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000
lar the optimization technique that was employed to calculate the equilibrium compositions. Numerical values of the parameter kij for different binary mixtures are reported in section 4. To further investigate the validity of the model, the PLJC equation was applied to ternary VLE and complex phase behavior including retrograde condensation and multiphase equilibria. These results are also reported in section 4. Finally, concluding remarks are summarized in section 5.
2. Perturbed Lennard-Jones Chain Equation of State: Extension to Mixtures Consider a mixture consisting of Nc number of fully flexible homonuclear Lennard-Jones chain molecules, whereby each molecule of component i is made up of mi fully flexible identical segments that are characterized by the segment size σii and segment energy ii. The total residual Helmholtz energy Amixt for the mixture at absolute temperature T and in volume V is given as
Amixt(Nc,V,T) AHSC APERT mixt mixt ) + NckT NckT NckT
(1)
where k is the Boltzmann constant. The residual Helmholtz energy for the hard-sphere chain reference fluid mixture is given by Von Solms et al.:7
AHSC mixt
)
NckT
F
[( ) ζ23
Fc ζ ζ 2 0 3
3ζ1ζ2
- 1 ln(1 - ζ3) +
ζ23
]
ζ0ζ3(1 - ζ3)2
+
[
∑k
ζ0(1 - ζ3)
(mk - 1)Fck Fc
ln(1 - ζ3) -
+
∑i Fi,
Fc )
∑i Fci,
π 6
× 3ζ2dkk
2(1 - ζ3)
∑i miFcidjii
)
21/6 -1 aii
3
)
0.005397mi + 0.006354 kT + 9.44 + mi ii 0.01222mi + 0.005102 (4) 9.947 + mi
The perturbation (attractive) term to the Helmholtz energy was derived by Von Solms et al.7 using a variational first-order perturbation theory based on the Gibbs-Bogoliubov inequality:
APERT mixt
)
NckT
2πFc
∑∑ T* i j
]
( )[ ( ) ]
ij σij xixjmimj 11 aij
3
IAij +
1
aij6
- 1 IBij
(5) Here, xi represents the chain fraction of component i chains and T* is a dimensionless temperature defined as T* ) kT/11. The perturbation integrals IAii and IBii were determined and shown to follow the analytical expressions5
IAii ) (-0.22169mi - 1.0755) - (2.2618mi - 1.0385)ζ3 + (1.47mi - 0.55411)ζ32 0.4571 + mi (6a)
(
) (
)
0.4213 0.0987 + 0.66264 + ζ3 + mi mi 1.7866ζ32 (6b)
(2)
Fi ) miFci, ζj )
(
IBii ) 0.03171 +
where
F)
ing density-independent equation:
The following combining rules are used in eq 5:
σij )
σii + σjj 2
aij )
aiiσii + ajjσjj σii + σjj
(3a-d)
Here, Fci and Fi denote the number densities of component i chains and chain segments, respectively, and dii represents the diameter of the hard-sphere segments. HSC /NckT consists of In eq 2, the Helmholtz energy Amixt the nonbonded hard-sphere mixtures and the bonding terms. The Boublik-Mansoori-Carnahan-Starling expression10,11 is used for the former, while the PercusYevick result12 is used for the latter. Equation 2 differs from the work of Von Solms et al.7 in that the hardsphere chain reference term is taken from Malakhov and Brun.13 The hard segment diameter dii of component i chains depends on temperature, density, and chain length and is conveniently represented by the dimensionless parameter aii ≡ dii/σii. Because the density dependence of aii is small over the density range of interest, Chiew et al.5 suggested that aii can be represented by the follow-
IAij )
ij ) xiijj(1 - kij)
IAii + IAjj 2
dij ) IBij )
dii + djj 2
IBii + IBjj 2
(7a,b) (7c,d)
(7e,f)
where kij is an adjustable binary parameter that is obtained by regressing the model against binary VLE experimental data. Given the Helmholtz energy, the compressibility factor Z can be determined in a straightforward manner via the following thermodynamic relation:
Z ) 1 + Fc
[ ( )] ∂ Amixt ∂Fc NckT
(8)
T*,x
Substituting eqs 1, 2, and 5 into eq 8 yields the following expressions for the compressibility factor:
Z ) ZHSC + ZPERT
(9)
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000 1499
where
Table 1. Basic Features of the GA Employed
[ ]
1
∑i mixi) 1 - ζ
ZHSC ) (
3
ζ2 ζ3
ζ0(1 - ζ3)3
Z
PERT
)
+
3
-
∑k
[
∑i ∑j xixjmimj
2
-
ζ0(1 - ζ3)
1
1 - ζ3
()
ij σij
3ζ23
+
ζ0(1 - ζ3)
(mk - 1)xk
2πFc T*
3ζ1ζ2
+
3
3ζ2dkk
]
2(1 - ζ3)2 (10a)
3
aij
×
[ ( )] 11
JAij +
1
aij6
- 1 JBij (10b)
with
JAii ) (-0.22169mi - 1.0755) - (4.5236mi - 2.077)ζ3 + (4.41mi - 1.6623)ζ32 0.4571 + mi (10c)
(
) (
)
0.4213 0.1974 + 1.3253 + ζ3 + mi mi
JBii ) 0.03171 +
5.3598ζ32 (10d) JAij ) JBij )
JAii + JAjj
(10e)
2 JBii + JBjj
(10f)
2
3. Regression Procedure and Optimization Techniques To obtain the kij values shown in eq 7b for binary mixtures of normal fluids, the PLJC equation of state was regressed against experimentally measured binary VLE data at different temperatures and pressures. This was achieved by minimizing the sum of squared error (SSE) between the calculated and experimental phase compositions with kij as the independent variables: Nd
SSE )
∑
k)1
[
]
ycalc.(Tk,Pk,k12) - yexpt.(Tk,Pk) Nd
∑
k)1
2
+ yexpt.(Tk,Pk) xcalc.(Tk,Pk,k12) - xexpt.(Tk,Pk)
[
xexpt.(Tk,Pk)
]
2
(11)
where Nd is the total number of data points used in the regression. The quantities y and x represent the mole fractions of vapor and liquid phases, respectively; subscripts calc. and expt. refer to calculated and experimental values, respectively. Note that for binary mixtures only one adjustable variable k12 is involved in the regression because kij is a symmetric parameter, i.e., k12 ) k21. Equal weight was assigned to both liquid and vapor phase data. A single-variable optimization technique known as the golden section search method14 was selected to perform the minimization of SSE. Preliminary testing indicated that SSE is unimodal and that
component
method/numerical value
encoding scheme selection scheme fitness technique reproduction technique type of crossover type of mutation population size no. of generations crossover probability mutation probability
real-value representation stochastic universal samplinga linear normalizationb elitist selection modified arithmetic crossoverc nonuniform mutationd 200 100 0.8 0.1
a Developed by: Baker, J. E. Proceedings of the 2nd International Conference on GA; Lawrence Erlbaum Associates: Hillsdale, NJ, 1987; p 14. b Developed by: Davis, L. Handbook of Genetic Algorithms; Van Nostrand Reinhold: New York, 1991; p 31. c Developed by: Haupt, R. L.; Haupt, S. E. Practical Genetic Algorithms; Wiley: New York, 1998; p 59. d Developed by: Michalewicz,17 p 128.
the golden section search is satisfactory for the present application. As indicated earlier, phase equilibrium compositions were calculated by minimizing the total Gibbs free energy of the system instead of solving a set of governing equations involving equality of fugacities. The derivation of the total Gibbs free energy from the Helmholtz energy is shown in the appendix. In this work, we employed a genetic algorithm to minimize the Gibbs energy. Genetic algorithm (GA) is a stochastic global optimization technique that simulates natural evolution in the solution space of the optimization problems. Unlike most conventional gradient-based techniques that sequentially evaluate a single point in the search region, GA operates on a population of potential solutions (individuals) in every iteration (generation) of the minimization process. When some individuals of the previous population are combined via “genetic” operations such as crossover and mutation, a new population that contains “fitter” individuals is consequently produced in the next generation. The selection of the individuals for reproduction is based solely on their objective function values, coupled with the implementation of Darwinian principles of survival of the fittest. In other words, those individuals that possess above average fitness values would have higher chances of being chosen for reproduction and thus survive into the next generation. Theoretically proven by the schema theorem, GA is guaranteed to yield better solutions after every subsequent generation.15 Because no information on derivatives is required and also because of its ability to locate the global optimum, GA has been successfully applied in many engineering optimization problems (e.g., Dasgupta and Michalewicz16), particularly those having multiple local optima. On the basis of these reasons, we have chosen GA for the minimization of the Gibbs free energy. Table 1 briefly describes the important features of the GA procedure that was employed in this work. It should be pointed out that GA is known to be efficient in identifying promising regions where the global minimum lies but it is rather ineffective in producing a sufficiently accurate solution.17 To overcome this problem, a local optimizer (modified simplex method of Nelder and Mead18) was incorporated in our procedure to refine the best solution obtained in the last generation of GA minimization. For a more detailed treatment of GA, refer to Holland,15 Goldberg,19 and Michalewicz.17
1500
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000
Table 2. Binary Parameter k12 for Various Binary Mixturesa system
T (K)
P (bar)
Nd
k12
AAD(y1)
AAD(x1)
hexane (1) + heptane (2) methane (1) + n-decane (2) ethane (1) + n-decane (2) propane (1) + n-decane (2) carbon dioxide (1) + n-decane (2) carbon monoxide (1) + hydrogen sulfide (2) hydrogen sulfide (1) + isobutane (2) carbon dioxide (1) + diethyl ether (2) carbon dioxide (1) + benzene (2) carbon dioxide (1) + n-butene (2) sulfur dioxide (1) + benzene (2) carbon dioxide (1) + toluene (2) tetrachloromethane (1) + chlorobenzene (2) trichloromethane (1) + triethylamine (2) dichloromethane (1) + triethylamine (2) carbon disulfide (1) + cyclohexane (2) 2-butanone (1) + toluene (2) naphthalene (1) + tetradecane (2) carbon dioxide (1) + n-butane (2) methane (1) + n-butane (2) methane (1) + carbon dioxide (2) methane (1) + hydrogen sulfide (2) carbon dioxide (1) + hydrogen sulfide (2)
303.15-323.15 410.93-477.59 277.59-377.59 310.93-410.93 344.26-444.26 293.15 277.65-344.37 298.15-313.15 298.15-313.15 273.15 255.37 311.26-477.04 353.15 283.14 283.15 298.15 330.15 420.05-447.05 255.98-283.15 277.59 270.00 277.59 266.48-305.37 302.59-352.59
0.103-0.498 13.789-293.025 3.447-104.73 3.447-60.191 13.789-188.363 17.934-237.239 2.068-51.152 7.017-72.226 8.938-77.499 1.266-34.855 0.273-0.553 3.337-152.924 0.249-1.077 0.041-0.126 0.056-0.286 0.165-0.452 0.201-0.459 0.133 1.048-41.334 6.922-125.621 35.554-70.207 17.237-82.736 27.579 68.947
15 26 26 25 25 19 24 18 17 13 12 24 19 9 9 10 13 11 12 14 6 13 9 11
-0.0053 0.0061 0.0202 -0.0026 0.1563 0.0516 0.0472 0.0959 0.1196 0.1254 0.0815 0.1375 0.0023 -0.0532 -0.0036 0.0121 0.0046 -0.0091 0.1729 -0.0208 0.1480 0.0487 0.1188 0.1516
0.0233 0.0271 0.0088 0.0067 0.0267 0.0522 0.0404 0.0177 0.0072 0.0210 0.0027 0.0532 0.0080 0.0087 0.0472 0.0103 0.0195 0.0397 0.0085 0.0450 0.1228 0.0399 0.0592 0.1057
0.0345 0.0551 0.0346 0.0178 0.0626 0.0672 0.0545 0.0598 0.0179 0.0226 0.0857 0.0718 0.0091 0.0339 0.0305 0.0127 0.0461 0.0750 0.0456 0.0344 0.0343 0.1201 0.0210 0.0902
a Experimental data taken from: Gmehling, J.; Onken, U. Vapor-Liquid Equilibrium Data Collection; Chemistry Data Series; DECHEMA: Frankfurt, Germany, 1977.
4. Results and Discussion The PLJC equation of state is applied to binary and ternary VLE, three-phase vapor-liquid-liquid equilibria (VLLE) of a ternary mixture, and retrograde condensation. These results are reported in the following subsections. 4.1. Binary VLE. Using the regressing procedure described in the previous section, we calculated the optimal values of the binary parameter k12 for a wide variety of nonassociating binary mixtures. Pure-component parameters, i.e., m, σ, and were obtained from Chiew et al.5 Optimal values of k12, together with the ranges of temperature and pressure of the experimental data, are reported in Table 2. In addition, we calculated the average absolute deviations (AAD) for the mole fractions of component 1 in both vapor and liquid phases based on the optimal binary parameter k12. The AADs are given by
Liquid phase: AAD(x1) )
1
Nd
∑ N k)1 d
Vapor phase: AAD(y1) )
1
Nd
∑ N k)1 d
|
|
x1,calc.(Tk,Pk,k12) - x1,expt.(Tk,Pk)
|
x1,expt.(Tk,Pk)
(12a)
|
y1,calc.(Tk,Pk,k12) - y1,expt.(Tk,Pk) y1,expt.(Tk,Pk)
(12b)
Table 2 shows that the PLJC equation provides a reasonably good fit for most of the selected binary pairs. Most AAD values for both liquid and vapor phases are found to be well below 0.1. In general, the PLJC equation predicts the vapor phase more accurately than the liquid-phase based on the AAD values. Figures 1 and 2 compare the predictions of P-x-y by the PLJC equation of state with experimental data
Figure 1. P-x-y vapor-liquid phase diagram for a hexane (1) + heptane (2) binary mixture.
for binary mixtures of two alkane binary mixtures. Included in these figures are predictions from the SAFT equation2 and the PHSC equation.4 Note that the modified version of the PHSC equation4 was used here. Figure 1 plots the P-x-y behavior of hexane (1) + heptane (2) at 303.15 and 323.15 K. It shows very good agreement between the PLJC predictions and the VLE data and that the performances of the PLJC, SAFT, and PHSC equations are similar. Figure 2 shows the P-x-y diagram of ethane (1) + n-decane (2) at 277.59 and 377.59 K. Predictions from the PLJC equation agree with the experimental data fairly well except in the vicinity of the critical points where theoretical predictions begin to deviate from experimental results. For the two binary systems in Figures 1 and 2, it appears that the PHSC equation performs slightly better than the PLJC and SAFT equations particularly at high pressures. Figures 3-5 report VLE for three different binary systems that do not involve linear hydrocarbons. Figure 3 shows the P-x-y diagram for tetrachloromethane (1)
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000 1501
Figure 2. P-x-y vapor-liquid phase diagram for a ethane (1) + n-decane (2) binary mixture.
Figure 5. P-x-y vapor-liquid phase diagram for a carbon dioxide (1) + benzene (2) binary mixture.
Figure 6. Phase diagram for a methane + carbon dioxide + n-butane ternary mixture. Figure 3. P-x-y vapor-liquid phase diagram for a tetrachloromethane (1) + chlorobenzene (2) binary mixture at 353.15 K.
Figure 4. P-x-y vapor-liquid phase diagram for a 2-butanone (1) + toluene (2) binary mixture at 330.15 K.
+ chlorobenzene (2) at 353.15 K. Again predictions from PLJC, SAFT and PHSC are included in the figure, and they are practically indistinguishable from the experimental data. Figure 4 plots the P-x-y diagram of 2-butanone (1) + toluene (2) at 330.15 K. Excellent agreement is found between model predictions and experimental data except when the mole fraction of 2-butanone approaches unity. Figure 5 displays the P-x-y behavior of carbon dioxide (1) + benzene (2) at 298.15 and 313.15 K. We included predictions of SAFT for comparison. Because the pure-component param-
eters of carbon dioxide for the PHSC model are unavailable, calculations based on this model are not possible. In contrast to the systems considered in Figures 3 and 4, the discrepancies between the PLJC predictions and experimental data in Figure 5 are smaller compared with those obtained from the SAFT equation. It is noteworthy that, even when the system is close to the critical point, PLJC is still able to predict the behavior quite accurately. 4.2. Ternary VLE. In addition to the binary systems considered above, we apply the PLJC equation to model VLE in ternary systems. Using the binary parameters reported in Table 2, we predicted the VLE compositions for two different ternary mixtures using the PLJC equation. Three cases were examined and are plotted in Figures 6-8. Figure 6 shows the ternary phase diagram for a mixture of methane (1), carbon dioxide (2) and n-butane (3) at T ) 277.59 K and P ) 27.58 bar. Excellent agreement between the experimental data and PLJC predictions is seen in the figure with slightly larger deviations in the vapor-phase compositions than in the liquid-phase compositions. Note that no adjustable parameter is used in these ternary calculations other than the kij values determined previously from binary data. Figure 7 displays the PLJC calculations and experimental data for the methane (1) + carbon dioxide (2) + hydrogen sulfide (3) ternary mixture at T ) 277.59 K and P ) 27.58 bar. Again very good agreement between model predictions and experimental data is observed. In Figure 8 we show the predicted VLE phase diagram
1502
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000
Figure 7. Phase diagram for a methane + carbon dioxide + hydrogen sulfide ternary mixture at 27.58 bar.
Figure 9. Equilibrium constant for a methane + carbon dioxide + n-butane ternary mixture.
Figure 8. Phase diagram for a methane + carbon dioxide + hydrogen sulfide ternary mixture at 68.95 bar.
for the same ternary mixture at the same temperature (277.59 K) but at a higher pressure P ) 68.95 bar. It is evident from the figure that discrepancy between model calculations and experimental data begins to increase as the pressure increases. Nevertheless, the PLJC equation is able to capture the qualitative characteristic (i.e., the concavity of the vapor-phase composition) of the coexistent curves rather well. Given that the binary parameters kij used in these ternary calculations were regressed from binary VLE data (and not from actual ternary data), the performance of PLJC is considered to be quite satisfactory. The agreement between model calculations and experimental data may be improved by varying the values of the binary parameters kij. As an example, the use of kCO2-H2S ) 0.19 and kCH4-CO2 ) 0.16 (instead of the values in Table 2) results in a smaller discrepancy between theoretical and experimental values, as shown by the dashed curves in Figure 8. However, a reevaluation of kij by directly regressing against ternary data was not performed in this study. We feel that using the kij’s obtained from binary VLE data allows us to identify the merits and weaknesses of the proposed PLJC model and sheds light on how it can be further improved. In addition, the equilibrium constants Ki ) yi/xi were also calculated for the three systems (shown in Figures 6-8) and are displayed in Figures 9-11. Essentially, the calculated values compare quite favorably with the experimental ones, particularly at lower pressures. 4.3. Complex Phase Equilibria of Ternary Systems. We further applied the PLJC equation to model complex phase behavior for the methane, carbon dioxide, and hydrogen sulfide ternary system. This ternary
Figure 10. Equilibrium constant for a methane + carbon dioxide + hydrogen sulfide ternary mixture at 27.58 bar.
Figure 11. Equilibrium constant for a methane + carbon dioxide + hydrogen sulfide ternary mixture at 68.95 bar.
mixture is commonly encountered in enhanced oil recovery, and it has been experimentally found that it exhibits various phase behaviors such as VLE, liquidliquid equilibrium (LLE), and VLLE under different conditions.20 A feed mixture consisting of 0.4989 mol of methane, 0.0988 mol of carbon dioxide, and 0.4023 mol of hydrogen sulfide was chosen for study. Like our earlier calculations, binary parameters kij reported in section 4.1 were used in our calculations. Two different phenomena were examined: (1) retrograde condensation at 303.15 K and (2) VLLE behavior at 211.45 K. Retrograde Condensation. We employ the PLJC equation to calculate the liquid-phase fraction φ (compared with vapor) for the ternary mixture. Figure 12 shows the liquid-phase fraction φ as a function of pressure at two different temperatures (T ) 253.15 and
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000 1503 Table 3. Vapor-Liquid-Liquid Behavior for a CH4-CO2-H2S Ternary Mixture P (bar)
ya
150.0 143.1 140.0 Figure 12. Liquid-phase fraction φ versus pressure P for a methane + carbon dioxide + hydrogen sulfide ternary mixture.
120.0
303.15 K). At T ) 253.15 K, the system exists as a vapor at pressure below 25 bar (i.e., φ ) 0). As the pressure increases, liquid condensation occurs and the liquidphase fraction φ increases until P ≈ 90 bar where the system is completely condensed (φ ) 1). For the case of T ) 303.15 K, as P increases from approximately 70 bar, the mixture begins to condense and φ increases. However, φ does not increase monotonically to unity. Instead it reaches a maximum value of approximately 0.4 and decreases to zero upon further increases in pressure. This phenomenon is characteristic of retrograde condensation. These calculated results are quite consistent with experimental data reported by Ng et al.20 except that the temperature T ) 303.15 K is slightly beyond the cricondentherm point found by Ng et al. Because the binary parameters are obtained at a lower temperature, slight discrepancies between model predictions and data are unavoidable. Multiphase Behavior in the Ternary System. We computed the phase equilibrium compositions of the CH4-CO2-H2S mixture at 211.45 K as a function of pressures using the PLJC equation of state. Table 3 reports the calculated phase compositions. In Table 3 the vapor-phase (y) mole fractions and the compositions of the liquid phases (xI and xII) are given, together with their corresponding phase quantities (V, LI, and LII). At higher pressure (P ) 150 bar), the mixture would exist as a single liquid phase that is rich in CH4. As the pressure is gradually reduced, the mixture enters into a two-phase liquid-liquid region in which the CH4-rich liquid phase coexists with a second H2S-rich liquid phase (52.7 < P < 143.1 bar). In the pressure range 52.55 < P < 52.65 bar, PLJC predicts VLLE. In the three-phase region, the phase quantity of the CH4-rich liquid decreases as the pressure is further lowered. This CH4-rich liquid phase vanishes at low pressures when the system enters into another two-phase region which consists of a vapor and a H2S-rich liquid. This multiphase behavior described by PLJC agrees quite well with the P-T diagram reported by Ng et al.20 The only difference is that the three-phase VLLE region predicted by the PLJC equation occurs at a lower pressure compared with that reported experimentally which lies approximately between 58 and 59 bar. Again, improvements on PLJC predictions can be made by fitting kij to ternary data directly instead of those used here.
100.0 80.0 60.0 52.7 52.65 52.60 52.55 52.5 40.0 20.0
0.9379 0.0273 0.0348 0.9377 0.0276 0.0347 0.9375 0.0278 0.0347 0.9374 0.0280 0.0346 0.9269 0.0367 0.0364 0.8776 0.0637 0.0587
xI a,b 0.4989 0.0988 0.4023 0.5116 0.0980 0.3904 0.5196 0.0976 0.3828 0.5615 0.0950 0.3435 0.5973 0.0925 0.3102 0.6327 0.0896 0.2777 0.6717 0.0857 0.2426 0.6881 0.0838 0.2281 0.6667 0.0937 0.2396 0.6637 0.0951 0.2412 0.6607 0.0965 0.2428
xII a,b
Vc
LI b,c
LII b,c
1.0 0.3930 0.1048 0.5022 0.3863 0.1052 0.5085 0.3533 0.1076 0.5391 0.3284 0.1097 0.5619 0.3069 0.1121 0.5810 0.2871 0.1149 0.5980 0.2800 0.1162 0.6038 0.2964 0.1268 0.5768 0.2988 0.1282 0.5730 0.3013 0.1296 0.5691 0.3022 0.1306 0.5672 0.1577 0.1483 0.6940 0.0575 0.1397 0.8028
0.8932
0.1068
0.8449
0.1551
0.6993
0.3007
0.6340
0.3660
0.5893
0.4107
0.5508
0.4492
0.5363
0.4637
0.2341
0.1412
0.6247
0.2635
0.0870
0.6495
0.2922
0.0325
0.6753
0.3097
0.6903
0.4436
0.5564
0.5382
0.4618
a The three values in each case refer to mole fractions of CH , 4 CO2, and H2S, respectively. b Superscripts I and II refer to a CH4c rich liquid and a H2S-rich liquid, respectively. V and L refer to phase quantities for vapor and liquid, respectively.
each binary pair. We found that the equation is able to represent binary VLE data very well over wide temperature and pressure ranges. It compares favorably with the SAFT and PHSC equations of state. Further, we demonstrated that the PLJC equation is able to model ternary VLE, retrograde condensation, and LLE, and VLLE in ternary systems with the empirical parameters regressed from binary data. All phase equilibrium calculations in this study were successfully carried out by the Gibbs free energy minimization using a GA; this has the advantages of simplicity and ensuring that a global minimum is found. Because the most significant advantage of the molecular-based equation of state is its ability to describe a wide range of fluid systems, future works on the PLJC model would include application on polymer-solvent mixtures and extensions to associating fluids.
5. Conclusions Acknowledgment The PLJC equation of state developed by Chiew et al.5 was extended to normal fluid mixtures. The mixture PLJC equation of state consists of one parameter kij for
Y.C.C. acknowledges support from the Academic Research Fund, National University of Singapore.
1504
Ind. Eng. Chem. Res., Vol. 39, No. 5, 2000
Appendix: Derivation of Total Gibbs Free Energy
Minimize
The molecular residual Gibbs free energy gr(N,V,T) of a phase can be defined as follows:
G′total
)
kT
kT )
gr ar ) +Z-1 kT kT
(A1)
where ar(N,V,T) is the molecular residual Helmholtz energy and Z is the compressibility factor. The molecular Gibbs free energy g for the phase is thus given as r
(A2)
where gig(N,V,T) is the ideal molecular Gibbs free energy. Because the ideal partial molecular Gibbs energy G h ig i of component i is expressed as 0 0 G h ig i ) Gi + kT ln Pi ) Gi + kT ln(xiP*)
) G0i + kT ln(kT) + kT ln(xiF)
(A3)
where G0i is the molecular Gibbs energy of pure component i at temperature T, Pi is the partial pressure of component i at ideal gas condition, and P* is defined as
P* ) FkT
(A4)
where F ) N/V is the molecular density; therefore, gig(N,V,T) and g can be redefined as follows:
gig
)
kT g
G0i
∑i xi kT ) ∑i xi kT + ln(kT) + ln(xiF)
)
kT
[
G h ig i
ar kT
+Z-1+
∑i
xi
[
G0i
kT
]
]
(A5)
∑i xi ln(xiF)
+ ln(kT) +
(A6)
Thus, the total Gibbs free energy Gtotal of a multicomponent and multiphase system is given as
Gtotal
)
kT )
∑j
gj Nj kT
(
arj
)
∑j Nj kT + Zj + ∑i xij ln xij + ln Fj
[
G0i
+
]
∑i NiT kT + ln(kT) - 1
(A7)
where index j denotes phase j. Note that Nj refers to the total number of moles in phase j and NiT refers to the total number of moles of component i in the system. From eq A7, to eliminate the constant terms for a particular system, the minimization problem can be formulated as
-
(
[
G0i
]
∑i NiT kT + ln(kT) - 1
arj
)
∑j Nj kT + Zj + ∑i xij ln xij + ln Fj
(A8) NiT )
subject to
∑j xijNj
0 e xij e 1 0 e Nj e
ig
a g g ) +Z-1+ kT kT kT
Gtotal
∑i NiT
(A9) (A10) (A11)
Literature Cited (1) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. (2) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse and Associating Molecules: Extension to Fluid Phase Mixtures. Ind. Eng. Chem. Res. 1991, 30, 1994. (3) Song, Y.; Lambert, S. M.; Prausnitz, J. M. A Perturbed Hard-Sphere Chain Equation of State for Normal Fluids and Polymers. Ind. Eng. Chem. Res. 1994, 33, 1047. (4) Song, Y.; Hino, T.; Lambert, S. M.; Prausnitz, J. M. Liquidliquid Equilibria for Polymer Solutions and Blends, including Copolymers. Fluid Phase Equilibria. 1996, 117, 69. (5) Chiew, Y. C.; Chang, D.; Lai, J.; Wu, G. H. A Molecular Based Equation of State for Simple and Chain-like Fluids. Ind. Eng. Chem. Res. 1999, 38, 4951. (6) O’Lenick, R.; Chiew, Y. C. Variational Theory for LennardJones Chains. Mol. Phys. 1995, 85, 257. (7) Von Solms, N.; O’Lenick, R.; Chiew, Y. C. Lennard-Jones Chain Mixtures: Variational Theory and Monte Carlo Simulation Results. Mol. Phys. 1999, 96, 15. (8) Von Solms, N.; Chiew, Y. C. Lennard-Jones Chain Mixtures: Radial Distribution Functions from Monte Carlo Simulation. Mol. Phys. 1999, 97, 997. (9) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982, 22, 731. (10) Boublik, T. Hard-sphere Equation of State. J. Chem. Phys. 1970, 53, 471. (11) Mansoori, G. A.; Carnahan, N. F.; Starling, K. E.; Leland, T. W., Jr. Equilibrium Thermodynamic Properties of the Mixture of Hard Spheres. J. Chem. Phys. 1971, 54, 1523. (12) Chiew, Y. C. Percus-Yevick Integral Equation Theory for Athermal Hard-Sphere Chains. Part I: Equations of State. Mol. Phys. 1990, 70, 129. (13) Malakhov, A. O.; Brun, E. B. Multicomponent Hard-Sphere Heterochain Fluid: Equations of State in a Continuum Space. Macromolecules 1992, 25, 6262. (14) Press: W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77: The Art of Scientific Computing; Cambridge University Press: New York, 1992. (15) Holland, J. H. Adaptation in Natural and Artificial Systems; The University of Michigan Press: Ann Arbor, MI, 1975. (16) Dasgupta, M.; Michalewicz, Z., Eds. Evolutionary Algorithms in Engineering Applications; Springer-Verlag: Berlin, 1997. (17) Michalewicz, Z. Genetic Algorithms + Data Structures ) Evolution Programs; Springer-Verlag: Berlin, 1992. (18) Nelder, R. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7, 308. (19) Goldberg, D. Genetic Algorithms in Searching, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1989. (20) Ng, H.; Robinson, D. B.; Leu, A. Critical Phenomena in a Mixture of Methane, Carbon Dioxide and Hydrogen Sulfide. Fluid Phase Equilib. 1985, 19, 273.
Received for review October 14, 1999 Revised manuscript received February 1, 2000 Accepted February 1, 2000 IE990744P
