Ind. Eng. Chem. Res. 2009, 48, 6877–6881
6877
Tuning the Interaction Parameter for the Soave-Redlich-Kwong Equation of State Using the Critical End Points Quan Yang,*,† Yigui Li,‡ and Shenlin Zhu‡ Key Laboratory of Energy Resources & Chemical Engineering, Ningxia UniVersity, Yinchuan, 750021 China, and Department of Chemical Engineering, Tsinghua UniVersity, Beijing, 100084 China
In the petroleum industry, to solve problems such as the design of the catalytic cracking of heavy oil, correct knowledge of the phase behavior is essential. As a benchmark where the number of degrees of freedom for binary mixtures is zero, the critical end points can represent the type of phase behavior of the binary mixtures according to van Konynenburg and Scott’s classification scheme (van Konynenburg, P. H.; Scott, R. L.Philos. Trans. R. Soc. London A 1980, 298, 495-540). Calculating the phase behavior of binary mixtures using interaction parameters evaluated by other methods can cause wrong determinations of the type of phase behavior. The critical end points were employed to fit the interaction parameters for binary mixtures in the critical state at high pressure. The molar Gibbs energy of mixing was evaluated to locate the ranges of values for the interaction parameters where exploration for the critical end point would not fail. The phase behavior at the critical end point is rather complex. Based on the Soave-Redlich-Kwong equation of state, the method developed by Heidemann and Khalil for calculating critical properties (Heidemann, R. A.; Khalil, A. M. AIChE J. 1980, 26, 769-780) was used to compute the critical points. A technique for determining the equilibrium phase of the critical point and thus determining the critical end point, which was difficult to evaluate, using the tangent plane criterion was developed. The developed algorithm was found to be rather stable and to converge rather quickly. The critical end points of different nonpolar, polar, and associating binary mixtures were calculated employing this algorithm. According to the calculation results and the experimental data for the critical end points, the interaction parameters were fitted successfully. 1. Introduction For binary mixtures, when a vapor or liquid phase is in equilibrium with a critical phase, the corresponding phase point is a critical end point. The critical end points lie at the end of a critical line. As a benchmark where the number of degrees of freedom for binary mixtures is zero, the critical end points can represent the type of the phase behavior according to van Konynenburg and Scott’s1 classification scheme and are eligible to fit interaction parameters for the Soave-Redlich-Kwong equation of state. The interaction parameter of a cubic equation of state is the parameter for the mixing rules.2 It is adjustable and represents the deviation of the temperature functions for the mixture from the geometric mean of the temperature functions for the corresponding pure components. The interaction parameter also shows the energy difference between different molecules. Many correlation equations have been developed for calculating the interaction parameters for cubic equations of state. Chueh and Prausnitz3 proposed a technique for calculating the interaction parameters of alkane mixtures. The critical volumes of the pure components were employed to compute the interaction parameters. Kato et al.4 developed a method to calculate the binary interaction parameters of CO2 and alkane mixtures. In that work, the interaction parameters were presented as functions of the absolute temperature and the eccentric factors of the alkanes. Valderrama and Molina5 published functions for evaluating the interaction parameters for cubic equations of states of mixtures containing H2. Nishiumi et al.6 presented the interaction parameter for cubic equations of state as a function of the critical * To whom correspondence should be addressed. Tel.: +86-9512062328. Fax: +86-951-2062323. E-mail:
[email protected]. † Ningxia University. ‡ Tsinghua University.
volumes. The function was applied to calculate the interaction parameters of the mixtures composed of hydrocarbons, CO2, N2, and H2S. Schulze7 developed an equation for mixtures containing He. Avlonitis et al.8 correlated the equations for mixtures containing N2. Jaubert and Mutelet9 calculated interaction parameters using the group contribution method. Nevertheless, the parameter values for different groups were needed. The technique was employed to calculate properties of alkane mixtures. These mentioned techniques correlated interaction parameters as functions of the eccentric factors, temperatures, critical properties, or numbers of carbon atoms of the components. These functions could be employed to determine the interaction parameters of only certain components or certain kinds of components. The method of fitting interaction parameters using experimental data directly includes two techniques. One takes the errors in the calculated pressures at the bubble points as the objective function, and the other takes the errors in the computed compositions of components as the objective function. The errors in the computed pressures at the bubble points are usually taken as the objective function because the pressures at the bubble points are more sensitive to the interaction parameters. Zhu et al.10 stated that the interaction parameters of cubic equations of state are functions of the temperature and pressure and that different expressions should be employed in different pressure regions. The interaction parameter value fitted using the computed pressures at the bubble points is suitable for calculating the phase behavior in the region near the bubble points. Here, such an interaction parameter was employed to evaluate the corresponding molar Gibbs energy of mixing for ethane and ethanol, at the temperatures and pressures of different critical points whose pressures were located in the range 5.10-5.60 MPa, where the critical end points lie, according to the experimental data. The corresponding Gibbs energy surfaces
10.1021/ie900069t CCC: $40.75 2009 American Chemical Society Published on Web 05/28/2009
6878
Ind. Eng. Chem. Res., Vol. 48, No. 14, 2009
parameters should be modified according to different pressure and temperature regions. 2.2. Algorithm to Compute the Critical End Points. At a critical point, the following equation obtained from the expansion of the Helmholtz free energy in a Taylor series must be satisfied
[
∆A -
]
∑µ i
1 3!
1 2!
)
i0∆ni0
T0,V0
∑ ∑ (∂ A/∂n ∂n )∆n ∆n 2
j
j
∑ ∑ ∑ (∂ A/∂n ∂n ∂n )∆n ∆n ∆n i
j
j
i
j
+
+ O(∆n4) ) 0
3
i
i
i
k
i
j
k
(4)
k
The quadratic term and the cubic term must each equal zero
∑ ∑ (∂ A/∂n ∂n )∆n ∆n 2
j
Figure 1. Molar Gibbs energy of mixing for pure ethane (1) and ethanol (2) at 314.9 K and 5.34 MPa with the interaction parameter for the binary mixture being 0.035.
j
2. Algorithm 2.1. Soave-Redlich-Kwong Equation of State. The Soave-Redlich-Kwong (SRK) equation of state, as one of the most representative cubic equations of state, is rather wellknown and widely used. In this work, the interaction parameter for the SRK equation of state14 was fitted. The equation can be presented in the following form P)
a(T) RT (V - b) V(V + b)
(1)
where a and b are the energy parameter and the covolume parameter of the Soave-Redlich-Kwong equation of state, respectively. The energy parameter and the covolume parameter for the mixtures were evaluated using the van der Waals mixing rules a)
∑ ∑ x x (1 - K )a i j
i
ij
1/2 i
aj1/2
(2)
∑xb
i i
j
)0
3
i
j
(5)
j
k
i
j
k
)0
(6)
k
With eqs 5 and 6, the critical points at different compositions for different mixtures can be evaluated using the Newton-Raphson method for solving equations. After a critical point was determined, the tangent plane criterion was employed to determine the equilibrium phase in the research. The criterion was first used to determine whether a phase was stable. At a given temperature and pressure, the component mole fractions of an M-component mixture are z1, z2,..., zM. If this mixture is separated into two phases, the number of moles in the second phase is ε, where ε is infinitesimal, and the mole fractions of different components in the second phase are y1, y2,..., yM. Because ε is infinitesimal, the compositions of the first phase approximately remain as z1, z2,..., zM. The variation of Gibbs energy for this process is as follows ∆G )
∑ εy µ (yj) - ∑ εy µ (zj) ) ε ∑ y [µ (yj) - µ (zj)] 0 i i
i i
i
i
i
i
0 i
i
(7) To ensure the stability of the system, the variation of the Gibbs energy for the above process should be larger than zero. If ∆G equals zero at composition jye and ∆G is always equal to or larger than zero at another trial composition jy, then the phase with composition jye is the equilibrium phase for the original mixture. That is, if ∆G equals zero at the stationary point of the ∆G function, then the stationary point represents the equilibrium phase. According to the necessary condition of a stationary point, the following equation is obtained from eq 7 To calculate the composition of the equilibrium phase, assume µi(yj) - µ0i (zj) ) 0,
i ) 1, 2,..., M
(8)
µi(yj) - µ0i (zj) ) K,
i ) 1, 2,..., M
(9)
The expression of chemical potential is µi(yj) - µ0i (zj) ) RT[ln yi + ln φi(yj) - ln zi - ln φi(zj)] (10) With eqs 9 and 10, assuming hi ) ln zi + ln φi(zj), θ ) K/RT, and Yi ) exp(-θ)yi, one obtains
j
b)
i
∑ ∑ ∑ (∂ A/∂n ∂n ∂n )∆n ∆n ∆n i
(Figure 1) had only one trough, exhibiting a completely concave surface throughout the composition range. Because a phase split will occur only if the molar Gibbs energy of mixing plot becomes convex at least once throughout the composition range, the shape of the above calculated plot shows that no equilibrium phase exists and that the original critical phase is stable. The results are wrong according to the experimental data and cause wrong determination of the type of binary mixture.1 In this work, the critical end points were employed to fit the interaction parameters for mixtures near the critical state at high pressure. Plots of the molar Gibbs energy of mixing for mixing the pure components were used to locate the proper value ranges for the interaction parameters. In 1980, Heidemann and Khalil11 developed a technique taking the temperature, volume, and compositions as independent variables. Heidemann and Khalil’s method is rather stable and converges rather rapidly. This method was employed to compute the critical points in this work. After the critical points had been obtained, the tangent plane criterion12,13 was used to determine the equilibrium phase.
i
i
(3)
i
Kij is the interaction parameter for components iand j. Usually, Kij ) Kji. According to Zhu et al.,10 the values of the interaction
ln Yi + ln φi(yj) - hi ) 0,
i ) 1, 2,..., M
(11)
The successive substitution method15,16 was employed to solve equation group 11. After Yi was obtained for a given
Ind. Eng. Chem. Res., Vol. 48, No. 14, 2009
6879
iteration, the composition yi for the next iteration was calculated using the equation yi ) Yi /
∑Y
(12)
j
j
If the correct calculation results satisfy the condition ∑jYj > 1 (θ < 0, K < 0, ∆G ) εK < 0), the original mixture is unstable and has the potential to separate into two different phases from the original phase. If ∑jYj < 1, the mixture is stable. In the above two cases, no equilibrium phases exist for the original phase. If ∑jYj ) 1, the variation of the Gibbs energy equals zero, and the corresponding compositions calculated are the compositions of the found equilibrium phase. 2.3. Fitting the Interaction Parameter. If, for a variation of the interaction parameter of ∆Kij, the corresponding variations of the temperature and pressure for the calculated critical end point are ∆T and ∆P, respectively, then the sensitivity ratio of the temperature and pressure to the interaction parameter can be defined as δ)
(∆T/T)/∆Kij ∆T/T ) (∆P/P)/∆Kij ∆P/P
(13)
The value of δ is close to 1, so in this work, the following equation was employed as the objective function to fit the interaction parameter diff )
(
Tcal - Tex Tex
) ( 2
+
Pcal - Pex Pex
)
2
(14)
3. Results and Discussion The molar Gibbs energies of mixing the pure components of different binary mixtures at the temperature and pressure of different critical points were evaluated. The corresponding values for the molar Gibbs energy of mixing were plotted according to the calculation results. For the ethane-ethanol mixture, in the pressure range of 5.10-5.60 MPa, when the interaction parameter was less than 0.06, the corresponding plot had only one trough and was concave throughout the composition range. For an interaction parameter value of 0.035, the molar Gibbs energy of mixing for the ethane-ethanol mixture at 314.9 K and 5.34 MPa is presented in Figure 1, where the surface is concave throughout the composition range. No phase split exists, and no critical end points can be obtained in this case. To obtain the critical end points for the ethane-ethanol mixture, the interaction parameter should have a value larger than 0.06. Figure 2 shows the molar Gibbs energy of mixing for the ethane-ethanol system at 314.9 K and 5.34 MPa with an interaction parameter value of 0.125. In the figure, the Gibbs energy surface becomes convex once. A phase split will be observed, and a critical end point can be obtained. For the methane-n-hexane mixture, in the critical pressure range of 2.79-3.52 MPa, when the interaction parameter is larger than 0.023, the corresponding Gibbs energy of mixing plot at the temperatures and pressures of different critical points is concave throughout the composition range, so no phase split exists and no critical end point can be obtained. To explore the critical end point for the methane-n-hexane mixture, the interaction parameter thus should be less than 0.023. After the proper ranges of the interaction parameter were determined, the developed algorithm was employed to calculate the critical end points corresponding to different interaction parameters in the ranges. The computation results for the ethane-ethanol mixture and the methane-n-hexane mixture are presented in Tables 1
Figure 2. Molar Gibbs energy of mixing for pure ethane (1) and ethanol (2) at 314.9 K and 5.34 MPa (Kij ) 0.125) Table 1. Computation Results of the Critical End Points for the Ethane (1)-Ethanol (2) System Kij
y1
ye,1
T (K)
P (MPa)
0.060 0.080 0.100 0.125 0.135 0.155
0.9361 0.9552 0.9653 0.9729 0.9752 0.9789
0.4207 0.2891 0.2097 0.1462 0.1277 0.0985
330.5 323.6 319.2 315.4 314.3 312.4
6.379 5.897 5.611 5.383 5.316 5.208
Table 2. Computation Results of the Critical End Points for the Methane (1)-n-Hexane (2) System Kij
y1
ye,1
T (K)
P (MPa)
-0.015 0.000 0.005 0.009 0.015 0.017
0.9560 0.9421 0.9360 0.9305 0.9210 0.9175
0.999970 0.999980 0.999990 0.999993 0.999996 0.999997
185.4 182.5 180.9 179.4 176.3 174.9
3.858 3.505 3.334 31.69 2.848 2.683
Table 3. Experimental Values of Critical End Points for Different Mixtures17-20 system
y1
methane (1)-n-hexane (2) 0.9291 ethane (1)-ethanol (2) -a methane (1)-n-heptane (2) propane (1)-phenanthrene (2) propane (1)-triphenylmethane (2) ethane (1)-propanol (2) a
ye,1 0.99995 -
T (K) P (MPa) 182.5 314.7 191.7 377.3 378.8 315.6
3.415 5.556 4.785 4.670 4.760 5.657
No experimental values available.
and 2, respectively. (In Tables 1-4, y1 and ye,1 represent the mole fraction of component 1 in the critical phase and the corresponding equilibrium phase, respectively.) From the data in Tables 1 and 2, it is observed that the calculation results of the critical end points are sensitive to the value of the interaction parameter. If an interaction parameter fitted using other methods, such as that fitted from the bubble point pressures, is employed to calculate the critical end points, the accuracy of the calculated critical end points cannot be ensured. Fitting the interaction parameters according to the experimental values at the critical end points for the methane-nhexane and ethane-ethanol mixtures in Table 3, the results are 0.001 for the methane-n-hexane binary system and 0.124 for the ethane-ethanol binary system. The objective function for the fitting is presented in eq 14. The developed algorithm was stable and converged rather rapidly. The calculated results for mixtures in Table 3 are
6880
Ind. Eng. Chem. Res., Vol. 48, No. 14, 2009
Table 4. Calculation Results of the Critical End Points for Different Mixtures and Errors in the Calculated Pressures Compared with the Corresponding Experimental Values system
Kij
y1
ye,1
T (K)
P (MPa)
∆P/P (%)
methane (1)-n-hexane (2) ethane (1)-ethanol (2) methane (1)-n-heptane (2) propane (1)-phenanthrene (2) propane (1)-triphenylmethane (2) ethane (1)-propanol (2)
0.001 0.124 -0.030 -0.030 -0.015 0.105
0.9462 0.9721 0.9997 0.9986 0.9989 0.9790
0.99998 0.14820 0.85991 0.84243 0.90056 0.26390
182.2 315.5 191.3 373.6 373.6 316.2
3.470 5.402 4.699 4.469 4.546 5.485
1.61 2.77 1.79 4.28 4.50 3.05
presented in Table 4. The interaction parameters in Table 4 were fitted using the corresponding experimental data in Table 3, and the data for the critical end points in Table 4 were determined with the fitted interaction parameters. From Table 4, it is observed that the largest error in the calculated pressure is 4.50%. This value is small compared with the errors in evaluating the critical pressure of associating mixtures employing equations of state.1,21,22 When the interaction parameters obtained using the fitting technique herein are employed to calculate the phase behavior, better results are obtained not only at the critical end points, as shown above, but also at state points far from the critical end points. In Table 4, the fitted value of the interaction parameter for the methane-n-hexane binary system is 0.001. The equilibria of the two phases, the vapor phase and the liquid phase, at 15.0 MPa with an interaction parameter value of 0.001 are presented in Figure 3. The calculation results with an interaction parameter value of 0.026 show that, at 15.0 MPa, no vapor-liquid equilibrium exists for the methane-n-hexane binary system, which is completely inconsistent with the experimental results,23 whereas with an interaction parameter value of 0.001, the correct results are obtained. 4. Conclusions The molar Gibbs energy of mixing was evaluated to locate the ranges of values for the interaction parameters where exploration for the critical end point would not fail. Based on the Soave-Redlich-Kwong equation of state, Heidemann and Khalil’s method was employed to compute the critical properties at the critical points. A technique for determining the equilibrium phase of the critical point and thus the critical end point using the tangent plane criterion was developed. The algorithm was found to be rather stable and to converge rather rapidly and
Figure 3. Tr,y1 isobar for the methane (1)-n-hexane (2) binary system at 15.0 MPa, where Tr represents the ratio of the system temperature to the critical temperature of methane, y1 represents the composition of methane, and V+L represents vapor- and liquid-phase coexistence. (9) Critical points.
was used to calculate the critical end points of different nonpolar, polar, and associating mixtures. The calculation results of the critical end points are sensitive to the values of the interaction parameter. Based on the proper interaction parameter ranges determined, the errors in both the calculated pressure and temperature at the critical end points were employed in the objective function to fit the interaction parameter. The largest error in the computed pressure at the critical end point using the fitted interaction parameter was 4.50%. With the interaction parameter fitted using the evaluated critical end point, which determines the type of phase behavior of the binary mixture, the accuracy of the calculated critical end point is ensured, and improved results are obtained at state points far away from the critical end points as well. Nomenclature a ) attraction parameter, Pa · m6/mol2 A ) Helmholtz free energy, J b ) covolume (size) parameter, m3/mol fi ) fugacity of component i, Pa G ) Gibbs free energy, J Kij ) interaction parameter M ) number of components n ) mole number, mol qij ) element of matrix Q Q ) matrix for the quadratic term of ∆A P ) pressure, MPa R ) universal gas constant, 8.314 J/(mol · K) T ) temperature, K V ) volume, m3 V ) molar volume, m3/mol x ) mole fraction Y ) variables for calculating equilibrium phase y ) mole fraction jy ) composition vector z ) mole fraction jz ) composition vector Greek Symbols δ ) sensitivity ratio of the temperature and the pressure to Kij ε ) number of moles θ ) variables for calculating equilibrium phase µ ) chemical potential φi ) fugacity coefficient of component i Subscripts cal ) calculation results e ) equilibrium phase ex ) experimental values i,j ) component indices T ) total
Literature Cited (1) van Konynenburg, P. H.; Scott, R. L. Critical lines and phase equilibriu in binary van der Waals mixtures. Philos. Trans. R. Soc. London A 1980, 298, 495–540.
Ind. Eng. Chem. Res., Vol. 48, No. 14, 2009 (2) Zhou, W.; Gong, M. Q.; Hong, H.; Wu, J. F. Development of binary interaction parameters. Cryogenics 2006, 1, 1–5. (3) Chueh, P. L.; Prausnitz, J. M. Vapor-liquid equilibria at high pressures: Calculation of partial molar volumes in nonpolar liquid mixtures. AIChE J. 1967, 13, 1099–1107. (4) Kato, K.; Nagahama, K.; Hirata, M. Generalized interaction parameters for the Peng-Robinson equation of state: Carbon n-paraffin binary systems. Fluid Phase Equilib. 1981, 7, 219–231. (5) Valderrama, J. O.; Molina, E. A. Interaction parameter for hydrogen containing mixtures in the Peng-Robinson equation of state. Fluid Phase Equilib. 1986, 31, 209–219. (6) Nishiumi, H.; Arai, T.; Takeuchi, K. Generalization of the binary interaction parameter of the Peng-Robinson equation of state by component familiy. Fluid Phase Equilib. 1988, 42, 43–62. (7) Schulze, W. A simple generalization of the binary temperaturedependent interaction parameters in the Soave-Redlich-Kwong and PengRobinson equations of state for helium mixtures. Fluid Phase Equilib. 1993, 87, 199–211. (8) Avlonitis, G.; Mourikas, G.; Stamataki, S. A generalized correlation for the interaction coefficients of nitrogen-hydrocarbon binary mixtures. Fluid Phase Equlib. 1994, 101, 53–68. (9) Jaubert, J. N.; Mutelet, F. VLE predictions with the Peng-Robinson equation of state and temperature dependent kij calculated through a group contribution method. Fluid Phase Equilib. 2004, 224, 285–304. (10) Zhu, H. B.; Gong, M. Q.; Zhang, Y.; Wu, J. F. Review of the research on interaction coefficients used in two cubic equations of state for phase equilibria predictions. Cryogenics 2005, 5, 7–11. (11) Heidemann, R. A.; Khalil, A. M. The calculation of critical points. AIChE J. 1980, 26, 769–780. (12) Michelsen, M. L. The Isothermal Flash Problem. Part I: Stability. Fluid Phase Equilib. 1982, 9, 1–19. (13) Michelsen, M. L. The Isothermal Flash Problem. Part II: Phase Split Calculation. Fluid Phase Equilib. 1982, 9, 21–40.
6881
(14) Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, 1197–1203. (15) Plybon, B. F. An Introduction to Applied Numerical Analysis; PWSKENT Publishing Company: Boston, MA, 1992. (16) Mehra, R. K.; Heidemann, R. A.; Aziz, K. An Accelerated Successive Substitution Algorithm. Can. J. Chem. Eng. 1983, 61, 590– 596. (17) Lam, D. H.; Jangkamolkulchai, A.; Luks, K. D. Liquid-LiquidVapor Phase Equilibrium Behavior of Certain Binary Ethane + n-Alkanol Mixtures. Fluid Phase Equilib. 1990, 59, 263–277. (18) Lin, Y.-N.; Chen, R. J. J.; Chappelear, P. S.; Kobayashi, R. VaporLiquid Equilibrium of the Methane-n-Hexane System at Low Temperature. J. Chem. Eng. Data 1977, 22, 402–408. (19) Chang, H. L.; Hurt, L. J.; Kobayashi, R. Vapour-Liquid Equilibria of Light Hydrocarbons at Low Temperatures and High Pressures: The Methane-n-Heptane System. AIChE J. 1966, 12, 1212–1216. (20) Peters, C. J.; Rijkers, M. P. W. M.; Roo, D. J. L.; Arons Swaan, J. D. Phase Equilibria in Binary Mixtures of Near Critical Propane and Poly-Aromatic Hydrocarbons. Fluid Phase Equilib. 1989, 52, 373–387. (21) Chen, Z. H.; Yao, Z.; Huang, Z. M. Comparison of Different Mixing Rules in the Calculation of Critical Properties for Binary Mixtures Based on PR Equation of State. J. Chem. Eng. Chin. UniV. 2008, 22, 365–370. (22) Gubblns, K. E.; Shing, K. S.; Streett, W. B. Fluid Phase Equilibria: Experiment, Computer Simulation and Theory. J. Phys. Chem. 1983, 87, 4573–4585. (23) Clever, H. L.; Young, C. L. Solubility Data Series; Pergamon Press: Oxford, U.K., 1987; Vols. 27/28, Methane
ReceiVed for reView January 15, 2009 ReVised manuscript receiVed May 11, 2009 Accepted May 22, 2009 IE900069T
