Ind. Eng. Chem. Res. 2000, 39, 3377-3386
3377
GENERAL RESEARCH A Fugacity Model for Gas Hydrate Phase Equilibria Jeffery B. Klauda and Stanley I. Sandler* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716
A classical thermodynamic model for the equilibrium pressures of different guests of structure I and II gas hydrates is proposed that removes the need for reference energy parameters used in the van der Waals and Platteeuw (vdWP) type models. The assumption of a constant crystal lattice for different guests within a structure, which is not in agreement with quantum chemistry calculations, is removed. This model uses published Kihara cell potential parameters determined from viscosity and second virial coefficient data, unlike previous models that fit these parameters to hydrate data. Quantum mechanical calculations were used to reduce the number of fitted parameters in the new model. This model greatly improves upon the accuracy of previous models and is applicable over a wide temperature range. Percent absolute average deviations in the predicted equilibrium pressures is 3.27% with two or three guest-hydrate specific parameters compared to more than 11% with three adjustable parameters in the vdWP type models and 8.55% using a recent thermodynamic model that uses five guest-hydrate parameters. Introduction Gas hydrates or clathrates are crystalline solids consisting of a guest(s) component(s) and water. Each guest is trapped in a cagelike water structure that forms a crystalline lattice. Hydrates can form at conditions above the normal freezing point of water. The guests of main interest here are methane, carbon dioxide, and other hydrophobic components that are present in the environment. There has been considerable interest in gas hydrates over the past 100 years, especially in the petroleum industry. Large hydrated masses occurring in natural gas pipelines, for example, in the Arctic regions and in the sea, can slow or completely obstruct flow. Also there has been a recent resurgence in developing methods to harvest the huge amounts of methane present in natural gas hydrates. A recent estimate of the amount of methane trapped in hydrates is as much as 300 times that in conventional U.S. reserves.1 The melting and dissociation of gas hydrates at the bottom of the ocean and in permafrost regions, and therefore release of methane, might increase global warming.2,3 This could lead to a cycle of the release of global-warming gases resulting in even higher global temperatures. There also has been recent interest in using hydrates as a sequestering media for carbon dioxide. Brewer et al.4,5 have demonstrated that it is possible to use hydrates to sequester CO2 at the ocean bottom. Hydrate Structure There are three known hydrate crystal structures. Structures I and II form with relatively small guests, e.g., methane, nitrogen, ethane, etc. Structure I (sI) is * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: (302) 831-2945. Fax: (302) 8314466.
a body-centered-cubic (bcc) structure6 with a lattice parameter of 12.03 Å with ethylene oxide as a guest. The cubic cell contains 46 H2O molecules, 2 dodecahedra (512), and 6 large voids (51262), where 512 is used to indicate that the polyhedron contains 12 five-membered ring faces.1 Structure II (sII) is a bcc structure7 with a lattice parameter of 17.31 Å for a double hydrate with tetrahydrofuran and hydrogen sulfide as guests. Each cubic cell contains 136 water molecules, 8 large voids (51264), and 16 dodecahedra. Structure H is only known to form with at least one small guest (i.e., methane) and one large guest,8 e.g., cyclooctane, methylcylcohexane, etc. The structure H (sH) unit cell contains 34 water molecules, 1 large void (51268), 3 dodecahedra, and 2 435663 cavities. The primary focus of this study is structure I and II hydrate formers. The radii of the different cavities in structures I and II are listed in Table 1. There has been a disagreement in the literature as to which set of cavity radii should be used. The cavity radii used in this study were taken directly from crystallography measurements6,7 by averaging the distance of all oxygens in the water molecules from the mass-weighted center of a cage. The radii used by Sloan1 are close to the crystallography measurements but are adjusted slightly to provide a better representation of the cell potential in his model.9 John et al.10 used a different set of radii which best characterized their simulation results of the guest-host interactions in a specific cavity. Other researchers11 have determined that radii from crystallography measurements should be used, because those best represent the actual cavity size. Van der Waals and Platteeuw Type Models The classic model for determining hydrate equilibrium pressures and temperatures was developed by van der Waals and Platteeuw.12 This model was later extended by Parrish and Prausnitz13 to account for
10.1021/ie000322b CCC: $19.00 © 2000 American Chemical Society Published on Web 08/12/2000
3378
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000
where
Table 1. Cage or Cavity Radii cage coordination number
cage radius (Å) structure I
cage 512 51262
II
512 51264
this study 3.906 4.326 3.902 4.682
Sloan1
John et al.10
this study
Sloan1
John et al.10
3.95 4.33 3.91 4.73
3.875 4.152 3.870 4.703
20 24 20 28
20 24 20 28
20 21 20 28
multiple guests in the hydrate structures. The details of the van der Waals and Platteeuw (vdWP) type models can be found in the book by Sloan.1 The equilibrium condition used in the vdWP model is the equality of the chemical potential of water in the hydrate phase, superscript H, and in the other equilibrium phase(s), π, which might be liquid water, ice, or both, i.e., β H β π π ∆µH w ) µw - µw ) µw - µw ) ∆µw
(1)
where µβw is the chemical potential of the empty hydrate, which can be considered to be a metastable or hypothetical phase because a hydrate requires a guest to form a stable phase. The difference between the chemical potential of water in the hypothetical and real (filled) hydrate phases is given by
∆µH w (T,P) ) -RT
νm ln(1 - ∑Θmj) ∑ m j
(2)
where νm is the number of cages of type m per water molecule in the hydrate lattice. The fraction of cages occupied by a guest is given by a Langmuir adsorption relation
Θml(T,P) )
Cml(T) fl(T,P) 1+
(3)
∑j Cmj(T) fj(T,P)
where fl(T,P) is the fugacity of guest l and the Langmuir constant, Cml, in the model is defined as
Cml(T) )
4π kT
∫0R(cell)-a exp(W(r)/kT)r2 dr
(4)
where W(r) is the spherically symmetric cell potential. The following Kihara cell potential with a spherical core was used in the development of the vdWP model
Γ(r) ) ∞, r e 2a Γ(r) ) 4
[(
σ 12 σ 6 , r > 2a r - 2a r - 2a
)]
) (
(5)
where is the characteristic energy, a is the core radius, and σ + 2a is the collision diameter. Equation 5 is used for the interaction between a guest molecule and one water molecule at the cavity wall. Because the cell potential, W(r), is the interaction between a guest and the entire cage, the guest-host interactions are summed over all guest-host interactions in the cage
W(r) ) 2z
[
(
)
σ12 a δ10 + δ11 11 R(cell) R(cell) r σ6 a δ4 + δ5 R(cell) R(cell)5r
(
)]
(6)
δN )
[(
1-
r a R(cell) R(cell)
(
1+
)
-N
-
r a R(cell) R(cell)
) ] -N
/N (7)
and z and R(cell) are the coordination number and the cell radius of the cavity or shell, respectively. Values for the σ, , and a parameters in eqs 6 and 7 in the model have been fit to measured single-component or multicomponent hydrate phase behavior. Although this model results in good correlations of the equilibrium pressures of hydrates, the description of the guest-host interactions may not be realistic. The deviation of these fitted interaction constants from those derived from viscosity and/or second virial coefficient data can be large.10 In addition, the original vdWP model assumed that only the first shell of water molecules was needed to represent total guest-host interactions. Others have suggested that a significant contribution to the interaction energies may come from the second and third shells.9-11,14 Presumably the original vdWP model is relatively successful because by fitting the Kihara potential parameters one indirectly accounts for the effect of the additional shells and other assumptions in this model. Improvements in determining equilibrium pressures of hydrate formation have been made by changing or removing some of the assumptions of the original vdWP model. Early work by John et al.10 replaced the empirically fit Kihara potential parameters with those obtained from the correlation of experimental data.15 For accurate predictions of hydrate equilibrium phase behavior, a correlated asymmetric factor was also included to compensate for the unevenness of the interaction potential within the cavity. The asymmetric factor was entered into the model as a multiplicative factor on the right-hand side of eq 4. Although this correction has a theoretical basis, the correction varied by as much as 5 orders of magnitude, casting doubt on its validity.14 The vdWP model also assumes that the hydrate lattice size is constant and independent of the guest component; that is, there is no lattice distortion. As a result of this assumption, the chemical potential of the hypothetical empty hydrate lattice is the same for all guests that occupy a specific hydrate structure. The difference in chemical potential between the empty hydrate lattice and the additional water phase (ice or liquid) then is
∆µπw(T,P) ∆µw(T0,P0) ) RT RT0
∫T
T
∆Hπw(T)
dT + RT2 π P ∆Vw(T) dP - ln(γwxw) (8) P0 RT 0
∫
where ∆µw(T0,P0) is referred to as the reference chemical potential at the reference temperature and pressure (T0, P0), usually taken to be the freezing point of water and zero pressure. Other factors on the right-hand side of eq 8 correct for temperature, pressure, and the activity of water. To obtain the equilibrium pressure between the hydrate phase and liquid water or the ice phase at a given temperature, one solves (generally by iteration) eqs 2 and 8 until the equality condition of eq 1 is satisfied.
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000 3379
Previous adaptations of the vdWP model assumed that the reference chemical potential remains constant even though experimental16,17 and computational work14,18 has demonstrated that the crystal lattice size depends on the specific guest. In particular, these studies found that the lattice constant can increase by about 3% from the smallest to the largest guest components that are stable in a specific crystal structure. The size of the empty hydrate lattice determines the strength of the hydrogen bonds between the water molecules; i.e., a more dense hydrate of the same structure has stronger hydrogen bonds than a less dense one. Therefore, the reference chemical potential and other hydrate properties should vary depending on the guest that occupies the lattice. To account for this change Zele et al.14 fitted the reference chemical potentials for different species in sI and sII to better represent experimental equilibrium data; they also used Kihara parameters from virial coefficient and viscosity data. The vdWP model is relatively successful, and recent improvements have made it more accurate, but the results are not completely satisfactory. While the vdWP models perform well near the ice point of water because that is where the reference chemical potential(s) are set, deviations at temperatures far away from the normal ice point can be significant. Although further improvements in the heat capacity used in the vdWP model might make it more accurate, we instead propose a new approach for hydrate equilibrium predictions. Theory and Development of a New Model The classical thermodynamic approach using equal fugacities between the hydrate and water is used here, as in studies of hydrates in equilibrium with vapor phase water19-21 and also in two recent papers by Chen and Guo22,23 for the modeling of three-phase equilibria of hydrates. While some success has been achieved with the Chen and Guo model, the model parameters were either related to energy parameters in the vdWP model22 or empirically fitted by fugacity functions23 based on the Antoine equation. Because the Antoine equation is known to lead to inaccurate vapor pressures beyond the range to which it has been fitted,24 the Chen and Guo model can be inaccurate over extended temperature ranges. Also their concept of the basic hydrate, which assumes complete occupation of the large cavity (51262 or 51264), is not supported by experimental evidence. The starting point here is the equality of fugacity for a hydrate in equilibrium with liquid water or ice π fH w (T,P) ) f w(T,P)
(9)
where the fugacity of water in the hydrate phase is
fH w (T,P) ) exp
(
)
(
)
IG β β IG µH µH w - µw w - µw + µw - µw ) exp (10) RT RT
Table 2. Shell Radii
structure
cage
1st
2nd
3rd
1st
2nd
3rd
I
512 51262 512 51264
3.906 4.326 3.902 4.682
6.593 7.078 6.667 7.464
8.086 8.285 8.079 8.782
20 24 20 28
20 24 20 28
50 50 50 50
II
Table 3. Kihara Cell Parameters for Eqs 6 and 7 compound
/k (K)
σ (Å)
a (Å)
CH4 C2H6 C3H8 i-C4H10 H2S N2 CO2 c-C3H6 H2O
232.2 404.3 493.71 628.6 459.6 142.1 513.85 554.4 102.134
3.505 4.022 4.519 4.746 3.607 3.469 3.335 4.185 3.564
0.28 0.574 0.6502 0.859 0.3508 0.341 0.677 0.653 0.000
form of the cell potential is slightly different from that in eq 6 in that an approach similar to John et al.10 is adopted in which the cell potential, W(r), includes contributions from the first, second, and third shells:
W(r) ) W(r)[1] + W(r)[2] + W(r)[3]
(
β fH w (T,P) ) f w(T,P) exp
)
-∆µH w (T,P) RT
σ)
ag + aw σg + σw ; ) xgw; a ) 2 2
where f βw(T,P) is the fugacity of the hypothetical empty hydrate lattice and ∆µH w (T,P) is defined in eq 2. The
(13)
where the subscript g represents the guest in the cavity and w is water. Gaussian 36-point integration was used to determine the Langmuir constants in eq 4.25 The contributions of the second and third shells to the potential are relatively flat over the integration range, so it was assumed constant at the value at the center of the cavity. The fugacities of water in the empty hydrate lattice, ice, and liquid water respectively are
( (
sat,β Psat,β w (T) φw (T) exp
(11)
(12)
The cell radii and coordination numbers for the corresponding shells are listed in Table 2. The first shell radii and coordination numbers are obtained from crystallography measurements.6,7 The second and third shell parameters used those obtained by John et al.10 To avoid a large number of parameters in our model, the shell radii were kept constant even though there is a different degree of lattice distortion for each guest. This is, however, compensated for elsewhere in the model. Also the values for the Kihara cell potential parameters used, listed in Table 3, were those obtained by Tee et al.15 from the correlation of second virial coefficient and gas phase viscosity data or in the case of methane directly from the second virial coefficient. That is, we have not obtained these parameters by fitting hydrate data. The following standard combining rules were used for the a, σ, and parameters:
f βw(T,P) )
Then rearranging eq 10
shell coordination number
shell radius (Å)
f Rw(T,P) )
Psat,R (T) φsat,R (T) exp w w
) )
Vβw(T,P) (P - Psat,β w (T)) (14) RT
VRw(T,P) (P - Psat,R (T)) w (15) RT
3380
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000
f Lw(T,P) ) xw(T,P) γw(xw,T) Psat,L (T) φsat,L (T) × w w exp
(
Table 4. Henry’s Law Constants for Eq 18 in Units of Pa-1
)
VLw(T,P) (P - Psat,L (T)) w (16) RT
where Psat,i w (T) is the vapor pressure of water in phase i. The fugacity coefficient, φsat,i w (T), is taken to be unity because the vapor pressures of the water phases are low. The volume of water in the Poynting correction in the above equations does depend on pressure, but because the change is small, we have assumed it to be constant. The composition of the guest component in the liquid phase is calculated using the Henry’s law constant.1 For most cases the solubility of hydrocarbons in water is extremely low; however, at very high pressures this solubility cannot be assumed to be negligible. Therefore, to be consistent, the solubilities of all guest compounds in water were calculated from
) xguest(s) i
f Vi Hi(T) exp(Z∞i )
(17)
The infinite dilution compressibility factor, Z∞i , for the gaseous species in liquid water is needed to calculate its solubility; therefore, an accurate equation of state is needed for water. The PRSV equation of state26 (a modification to the Peng-Robinson equation of state) was used because of its relative success with water. The following temperature dependence of the Henry’s law constant was used
-ln
(
)
H[2] i
Hi(T) [4] + H[3] ) H[1] i + i ln(T) + Hi T (18) 101325 T
where the constants in this equation are given in Table 4 and the units for Hi(T) are Pa-1. For compounds with appreciable solubility in water, e.g., CO2 and H2S, and at high pressures, the Henry’s law activity coefficient, γw(xw,T), cannot be assumed to be equal to unity. The modified UNIFAC27 activity coefficient model with the parameter set of Dahl et al.28 was used to calculate this term. The UNIFAC parameters fitted by Dahl et al. are for the light gases and hydrocarbons of interest here and were not included in the original modified UNIFAC model. The advantage of this model is its applicability to both low and high pressures that are of particular interest in hydrates. For compounds with appreciable solubility, there is a slight freezing point depression of water. The extent of the freezing point depression can be predicted using 2
∆dT ) Td - Tnfp )
R(Tnfp) ln[xwγw(xw,T)] ∆Hm
(19)
where Td is the freezing point with the solute(s), Tnfp is the normal freezing point of pure water, ∆Hm is the heat of melting pure ice, and the activity coefficient is calculated by the modified UNIFAC model. The molar volumes used in eqs 14 and 16 are pressure dependent, which is important in the Poynting correction because hydrates of interest here form at pressures in the range of 100-3000 MPa. Therefore, an accurate representation of the volume is needed. However, ice in equilibrium with the hydrate forms at relatively low pressures so that the Poynting correction is negligible,
compound i
H[1] i
H[2] i
H[3] i
H[4] i
CH4 C2H6 C3H8 i-C4H10 H2S N2 CO2 c-C3H6
-183.7860 -268.4410 -316.4900 96.1158 -149.5510 -164.9970 -159.8680 -60.9240
9112.582 13369.400 15922.700 -2472.570 8227.328 8433.619 8742.426 9177.534
25.0405 37.5561 44.3285 -17.3680 20.2327 21.5601 21.6712 0.0000
-0.00015 -0.00230 0.00000 0.00000 0.00129 0.00844 -0.00110 0.07278
Table 5. Vapor Pressure Constants in Eq 23 phase
A
B
C
D × 103
% AAD
T range (K)
ice 4.6056 -5501.1243 2.9446 -8.1431 1.507 150-273.15 water 4.1539 -5500.9332 7.6537 -16.1277 0.171 273.15-343.15
and there is no need for the pressure dependency on the molar volume of ice. The molar volumes for ice and liquid water phases were fitted to experimental data from Perry’s Handbook29 and NIST,30 respectively, as follows:
VRw(T [K]) ) 1.912 × 10-5 + 8.387 × 10-10T + 4.016 × 10-12T 2 (20) ln(VLw(T [K],P [MPa])) ) -10.9241 + 2.5 × 10-4(T - 273.15) - 3.532 × 10-4(P - 0.101325) + 1.559 × 10-7(P - 0.101325)2 (21) The objective function for fitting the constants in eqs 20 and 21 was
X({wi}) )
∑ ∀i
|wexp - wcalc i i | wexp i
(22)
where wi is the volume in units of m3/mol. The volume of ice Ih was fitted to temperatures ranging from 150 to 273.15 K with a percent absolute average deviation (% AAD) of 0.025. The liquid molar volume of water was fitted for temperatures varying from 273.15 to 335 K and absolute pressures from 0.01 to 800 MPa with a % AAD of 0.574. [Objective functions analogous to eq 22 were also used in all other parameter estimations.] The final term that is to be evaluated in eqs 14-16 is the vapor pressure of ice and liquid water. The experimental data from the CRC Handbook31 were fitted to the quasi-polynomial, QL1, form24
ln(Psat w [Pa]) ) A ln(T) +
B + C + DT T
(23)
The advantage of the QL1 form over the Antoine equation is its good predictions beyond the range of fitting.24 The values for the constants in eq 23 for ice and liquid water are given in Table 5 along with the % AAD and the temperature range of fitting. Evaluation of Model Parameters The empty hydrate lattice, β, does not exist in nature so there is no direct method of determining its vapor pressure. Ng and Robinson20 fitted the vapor pressure of the hypothetical empty hydrate lattice with an Antoine-type equation using data for hydrates in equilibrium with only the vapor phase (H-V). They assumed
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000 3381
that the empty lattice of each hydrate structure was independent of the guest but dependent on the structure of the hydrate. Although this procedure was successful for H-V equilibria, attempts to use this same vapor pressure for a large set of compounds of ice [or liquid water]-hydrate-vapor (I [L]-H-V) did not yield good results, undoubtedly because the assumption of a constant hydrate lattice for all guests is invalid. As previously mentioned, it is known from experiment that the size of the hydrate lattice depends on the guest, and therefore it is reasonable to expect that the vapor pressure of the hypothetical empty hydrate lattice will also depend on the guest. The vapor pressure of the empty hydrate lattice was fitted for each guest separately and for each hydrate structure using eq 23. Without any loss of accuracy, the parameter C in eq 23 was taken to be constant for all structures and guests, and only fitted to the methane I-H-V data using the objective function in eq 22. In this equation the absolute differences of the left-hand and right-hand sides of eq 9 for experimental temperatures and pressures were used, divided by the left-hand side of eq 9. The vapor pressures of the empty hydrate lattices were fitted to all I-H-V data and a few data points in the L-H-V region near the ice point of water to determine the A, B, and D parameters in eq 23 for all of the single guest hydrates. There are then two different ways of proceeding. One is to fit the A, B, and D parameters to experimental data. Another, which is what we have used here, is to use computational quantum mechanics as a tool to reduce the number of fitted parameters. The heat of vaporization of the empty hydrate lattice can be determined from
∆Hβvap(T) ) RT 2
ln P (d dT )
sat
(24)
and by substitution of eq 23 into eq 24.
∆Hβvap(T) ) R(-B + AT + DT 2)
(25)
Because the values of energy obtained from computational quantum mechanics correspond to those at absolute zero temperature, the B parameter can be obtained from such calculations. The Gaussian 9832 suite of programs was used to compute these energies. Because the computation time increases with the number and size of the atoms, only the first shell, i.e., a single cavity, was included in our calculations. Lattice positions of the water molecules from crystallography measurements6,7 were used as a first guess for the geometry of the cages. Next a guest was inserted into the center of a cage and the geometry determined by minimizing the energy of the complex. The energy minimization was performed using the Hartree-Fock (HF) level of theory and a double-ζ basis set with polarization functions on non-hydrogen atoms, 6-31g(d). Also calculations were done for an isolated water molecule optimized at the HF/6-31g(d) level to represent the vapor phase. Because the Hartree-Fock level of theory does not accurately predict energies of hydrogen-bonded complexes,33,34 the geometries from those calculations were used, but the energies of the empty cage and the isolated water molecule were calculated at the secondorder Møller-Plesset (MP2) level of theory with the same basis set. Corrections for the known basis set superposition error (BSSE)35 were made for the empty
Figure 1. (a) Stable 51262 clathrate cavity with C2H4O as the guest molecule. (b) Unstable 51262 clathrate cavity with i-C4H10 as the guest molecule.
cavity by including the basis set of each guest but removing the atoms themselves. The above procedure was repeated for all guests in sI in the 512 and 51262 cages, but because of size effects, some guests are not stable in the smaller cavity. The instability of such guests in the small cage was seen during the quantum chemical geometry optimization by the loss of the cage structure; see Figure 1a,b. The heat of vaporization was calculated as the product of the energy of the isolated water molecule and the number of water molecules in a cavity minus the energy of the total cavity. This assumes that the zero-point energies of an isolated water molecule and the cavity cancel. Because all guests of interest fit the large cavity, it is assumed that all of the large cavities are completely occupied. Only CH4 and H2S in sI are small enough to fit into the small 512 cavity; the larger guests form hydrates with an empty small cavity. Assuming complete occupation of the relevant cavities, the heat of vaporization can be approximated by
1 3 small,i ∆Hβ,i + ∆Ularge,i vap = ∆Uw w 4 4
(26)
where ∆Uw is the internal energy of vaporization per
3382
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000
Figure 2. Correlated sI values for the B parameter of QL1 in eq 23.
Figure 3. Correlated sII values for the B parameter of QL1 in eq 23. Table 6. Vapor Pressure Constants for the Empty Hydrate Lattice in Eq 23 structure
guest
A
B
C
D × 103
sI
CH4 C2H6 c-C3H6 CO2 H2S N2 C3H8 i-C4H10 c-C3H6
4.6477 4.6766 4.6652 4.6188 4.6446 5.1511 5.2578 4.6818 5.1449
-5242.9790 -5263.9565 -5424.1108 -5020.8289 -5150.3690 -5595.4346 -5650.5584 -5455.2664 -5544.3272
2.7789 2.7789 2.7789 2.7789 2.7789 2.7789 2.7789 2.7789 2.7789
-8.7156 -9.0154 -8.8658 -8.3455 -8.7553 -16.0445 -16.2021 -8.9678 -14.8446
sII
water molecule in each cage type. The weights for the internal energies of vaporization in eq 26 are in the ratio of the number of that cavity type to the total number of cavities in the crystal lattice. This ignores cavity interconnections. However, our goal was not to get an accurate heat of vaporization from quantum chemical calculations but rather to develop a method to correlate the B parameter in eq 25. Figure 2 shows a linear correlation obtained between the fitted B parameter and the calculated value of that parameter from MP2/6-31g(d). The quantum mechanically calculated parameter and this linear correlation were used in the calculations here to obtain the value of the B parameter. Consequently, we then had to obtain only the A and D parameters from correlation. The final set of vapor pressure parameters for the sI compounds obtained in this way are listed in Table 6. A similar procedure was
Figure 4. Comparison of ∆Hvap of ice to empty hydrate lattices.
used for the sII guest but with weights of 2/3 and 1/3, respectively, in eq 26. The correlation obtained for sII compounds is in Figure 3, and the parameters for the QL1 equation are also listed in Table 6. Interestingly, almost identical results for the correlation of hydrate pressures were obtained using the B parameter from quantum mechanics or using a three-parameter (A, B, and D) correlation of the data. Because the hypothetical empty hydrate lattice has longer and thus weaker hydrogen bonds, its heat of vaporization should be lower than that of ice. In Figure 4 the heat of vaporization is plotted versus temperature using eq 25, and as expected the heat of vaporization is lower. Interestingly, the curvatures of the heats of vaporization of all of the empty hydrates are remarkably close to that of ice. For methane hydrates the model fit to the equilibrium data was not completely satisfactory and was not improved by fitting the vapor pressures for the entire temperature range (150-320 K). The errors became noticeable at pressures above 20 MPa and increased as the pressure increased. Because the pressures for this hydrate are extremely high, the Poynting correction is important and an accurate description of the hydrate phase molar volume that depends on temperature as well as pressure is required. The temperature dependence on the volume of sI and sII hydrates was determined previously,36 but to date there has been no description on pressure-induced compression. Therefore, we propose the following equations for the volume of hydrates: -5 V β,I w (T,P) ) (11.835 + 2.217 × 10 T +
10-30NA
2.242 × 10-6T 2)
N βw
- 8.006 × 10-9P +
[
5.448 × 10-12P2
T in K P in MPa
]
(27)
-4 Vβ,II w (T,P) ) (17.13 + 2.249 × 10 T + 2.013 ×
10-30NA
10-6T 2 + 1.009 × 10-9T3)
N βw
-
8.006 × 10-9P + 5.448 × 10-12P2 (28) where the constants multiplying the temperature factors were correlated by Tse from the experimental thermal expansion of hydrates.36 The compressibility
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000 3383 Table 7. Hydrate Equilibrium Pressure Predictions (% AAD from Experiment) guest CH4 C2H6 CO2 H2S c-C3H6
N2 C3H8 i-C4H10
entire data
equilibria
no. of data points
ref
temp range (K)
Sloan1
Parrish and Prausnitz13
Chen and Guo23
this work
I-H(sI)-V L-H(sI)-V overall I-H(sI)-V L-H(sI)-V overall I-H(sI)-V L-H(sI)-V overall I-H(sI)-V L-H(sI)-V overall I-H(sI)-V L-H(sI)-V I-H(sII)-V L-H(sII)-V overall L-H(sII)-V I-H(sII)-V L-H(sII)-V overall I-H(sII)-V L-H(sII)-V overall
18 107 125 12 51 63 22 115 137 10 31 41 8 13 8 6 35 60 15 51 66 17 36 53
1 1 1 1 1 1 1 1 1 1 1 1 39 39 39 39 39 1 1 1 1 1 1 1
148.8-270.9 273.2-319.9 148.8-319.9 200.8-272.0 273.4-287.4 200.8-287.4 151.6-271.8 271.6-282.8 151.6-282.8 250.5-272.8 272.8-302.7 250.5-302.7 237.2-257.92 274.99-289.36 258.09-272.04 273.15-274.25 237.2-289.36 272.0-305.2 247.9-272.9 273.2-278.4 247.9-278.4 241.4-273.1 273.2-275.1 241.4-275.1
6.52 9.28 8.88 10.41 10.94 10.84 24.40 6.43 9.32 10.84 7.22 8.08 N/A N/A N/A N/A N/A 12.50 13.67 3.32 5.67 27.22 39.32 35.44
12.80 20.69 19.55 12.72 7.44 8.45 19.42 11.48 12.75 1.95 18.18 14.31 0.48 0.61 0.49 0.74 0.57 11.02 5.43 3.91 4.25 14.70 1.67 5.85
35.76 9.06 12.90 8.83 2.73 3.89 40.97 1.71 8.01 12.66 6.80 8.19 N/A N/A N/A N/A N/A 17.72 4.73 2.85 3.28 3.21 2.13 2.47
4.69 3.96 4.06 3.60 3.60 3.60 1.83 3.07 2.87 1.74 5.86 4.88 1.32 1.89 2.34 0.67 1.66 1.79 4.05 3.85 3.90 4.68 1.91 2.80
overall
580
max average
39.32 11.74
20.69 11.36
40.97 8.56
5.86 3.27
constants in eqs 27 and 28 were obtained by correlating the data for methane hydrates at high pressures but are assumed to be applicable to all high-pressure hydrates. The linear compressibility factor is 2 orders of magnitude larger than that of pure ice,37 which is to be expected as a result of the longer and weaker bonds in the hydrates. Because the hydrate is a solid phase, the compressibility of the hydrate should be much less than that of liquid water, which is also the case. Recently, high-pressure methane hydrate (up to 700 MPa) experimental work by Hirai et al.38 calculated the molar volume of the hydrate as a function of equilibrium pressure, and the molar volume predicted using eq 27 was within 2.50% AAD. Results and Discussion The prediction of equilibrium pressures of gas hydrates using the model developed here follows the algorithm outlined in Figure 5. The Appendix contains the set of equations for predicting hydrate equilibrium, and a C++ code for the model is available on request from the authors. Iteration on pressure in the inner loop of Figure 5 is done using a numerical NewtonRhaphson technique until the difference between the left-hand side and the right-hand side of eq 9 is less than an acceptable tolerance (10-9). Once a solution is found for the sI hydrate, the program recalculates equilibrium pressures for the sII hydrate, and the equilibrium pressure is taken to be the lower of the two. We compare the results of the model with those of three other models: two vdWP-type models (Sloan1 and Parrish and Prausnitz13) and a fugacity approach (Chen and Guo23). Table 7 and Figure 6 display the statistics for the equilibrium pressures of hydrate structures I and II. The maximum error for our model occurs for L-H-V with hydrogen sulfide (% AAD of 5.9); the experimental data for this system is slightly more scattered than other data sets. For comparison, the maximum errors of the other models are much greater: % AAD of 39.2, 20.7, and 41.0 for the Sloan,1 Parrish and Prausnitz,13 and Chen and Guo23 models, respectively.
Figure 5. Computational algorithm.
One reason for the inaccuracy of the other models is the assumption of a constant reference chemical potential. For some compounds the vdWP model predicts
3384
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000
Figure 6. % AAD from the experiment for all hydrates in the study.
Figure 7. Methane I-H-V equilibrium predictions.
reasonable pressures when ice, not liquid water, is the equilibrium phase or vice versa. For example, the prediction of carbon dioxide in the L-H-V region is reasonable using the Sloan1 model (about 6% error), but in the I-H-V region, the error is about 24%. Similar behavior also occurs for the Parrish and Prausnitz13 adaptation of the vdWP model. A change in the reference chemical potential will shift the equilibrium curve up or down but not change its slope for a specific phase and therefore will not resolve the inaccuracies. A change in the reference enthalpy and/or heat capacity is needed in the temperature integral for eq 8. In most cases the fugacity approach of Chen and Guo is more accurate than the vdWP-type models. However, for some data points the model is less accurate, e.g., I-H-V for CO2 and CH4. In these cases the model resulted in large deviations just beyond the temperature range used for fitting the model parameters. The temperature dependence of the phase equilibria for methane hydrates is shown in Figures 7 and 8. The predictions of our model are more accurate than those of other models over the entire region. The Sloan1 adaptation of the vdWP model is of comparable accuracy to our model in the I-H-V region but not in the L-H-V region, where the model underpredicts the equilibrium pressures. The Parrish and Prausnitz13 version of the vdWP model performs poorly at high temperatures in the L-H-V region and at low temperatures in the I-H-V region. The parameters of the Chen and Guo23 model were fit in the temperature region of 259-303.8 K. Below this temperature region, the model performs poorly, probably because of poor
Figure 8. Methane L-H-V equilibrium predictions.
Figure 9. Isobutane equilibrium predictions.
extrapolation of the fugacity as a result of using the Antoine equation. Figure 9 demonstrates the improvement of the model here compared with other models. All vdWP models predict the behavior of isobutane hydrates in the I-H-V region poorly, although the Chen and Guo23 model predicts the phase behavior well at higher temperatures. The parameters for that model were fit to temperatures above 253.7 K for this system, and below this temperature the model underpredicts the equilibrium pressures. Predicting the equilibrium pressures of cyclopropane hydrates is a more stringent test of the models than
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000 3385
with other hydrates, because cyclopropane forms both sI and sII hydrates, depending on the temperature of formation.39 The structure II hydrate forms between 258.1 and 273.3 K; at higher and lower temperatures, a structure I hydrate forms. Of the other models considered in this study, only the Parrish and Prausnitz13 model has parameters for cyclopropane. That model results a slightly smaller percentage error than the proposed model, although it results in a larger error for all other systems. Interestingly, both models accurately predict the temperature of the phase change from sI to sII and back to sI (Table 7). Conclusions A classical thermodynamic approach using fugacities rather than reference energy parameters has been shown to result in a model that is more accurate than vdWP-type models for the prediction of hydrate equilibrium pressures. An interesting characteristic of this model is the ability to predict the value of one of its parameters from quantum mechanics, thereby reducing the number of empirically fitted parameters. The model proposed here also uses published values of Kihara cell potential parameters from viscosity and second virial coefficient data, reducing the number of adjustable parameters to only three per guest in a hydrate structure or with a quantum mechanics calculation to two per guest. For comparison, in the vdWP models the Kihara cell parameters are fit and there are three adjustable parameters per guest, while the model developed by Chen and Guo has three to six adjustable parameters per guest. Therefore, the new model provides a more accurate description of hydrates with the same number or fewer adjustable parameters. The assumption of a constant hydrate lattice independent of a guest is a limitation of the vdWP-type models and the major reason for the deficiency of this class of models. Making the reference chemical potential and enthalpy in eq 8 dependent on the guest might improve the vdWP-type models but at the expense of additional parameters. Extension of this model to mixture hydrates is underway. Because the vapor pressure of the empty hydrate lattice is guest dependent, a mixing rule for the vapor pressure parameters will be needed to extend the predictions to mixture hydrates. Also, no longer treating the Kihara potential parameters as adjustable opens a method for obtaining the parameters from quantum mechanical calculations and better represents the true potential as well as improves the prediction of cage occupancies. Acknowledgment Financial support of this research was provided from contract DE-FG02-85ER13436 from the Division of Basic Energy Sciences of the U.S. Department of Energy. The authors also thank University of Delaware and National Center for Supercomputing Applications (NCSA) for the use computers for the quantum mechanical calculations. Appendix The equilibrium condition used is the equality of fugacities between the water phases π fH w (T,P) ) f w(T,P)
(A1)
where the fugacity of water in the hydrate phase is
given by
(
β fH w (T,P) ) exp A g ln(T) +
{
Bβg
[ (
T
V βw(T,P) P - exp Aβg ln(T) +
exp
[ ( ∑ν
exp
m
)
+ 2.7789 + DβgT × Bβg T
+ 2.7789 + DβgT
RT
ln 1 -
m
∑ j
Cml(T) fl(T,P) 1+
∑C
mj(T) fj(T,P)
}
)]
)]
×
(A2)
The constants Aβg, Bβg, and Dβg are given in Table 6, the Langmuir constant is calculated using eq 4, and the molar volume of the hypothetical empty hydrate is calculated using eqs 27 and 28. The fugacities of water in the ice and liquid water are
(
f Rw(T,P) ) exp 4.6056 ln(T) -
5501.1243 + T
)
2.9446 - 8.1431 × 10-3T ×
{[
[
(
exp V Rw(T) P - exp 4.6056 ln(T) -
5501.1243 + T
)]] } (A3)
2.9446 - 8.1431 × 10-3T /RT f Lw(T,P) ) xw(T,P) γw(xw,T) ×
(
exp 4.1539 ln(T) -
5500.9332 + T
)
7.6537 - 16.1277 × 10-3T ×
{[
[
(
exp V Lw(T,P) P - exp 4.1539 ln(T) -
5500.9332 + T
)]] } (A4)
7.6537 - 16.1266 × 10-3T /RT
where the molar volumes of ice and liquid water are calculated using eqs 20 and 21. Nomenclature Roman Characters a ) core radius C(T) ) Langmuir constant f ) fugacity H ) enthalpy H[i] ) Henry’s law constants k ) Boltzmann constant NA ) Avogadro’s number Nw ) number of water molecules r ) radial distance P ) pressure R ) gas constant R(cell) ) cell radius sI ) hydrate structure I sII ) hydrate structure II T ) temperature U ) internal energy vdWP ) van der Waals and Platteeuw V ) molar volume W(r) ) cell potential W(r)[i] ) cell potential for shell i x ) molar composition in the liquid phase z ) coordination number of the cage Z ) compressibility factor
3386
Ind. Eng. Chem. Res., Vol. 39, No. 9, 2000
Greek Characters Χ ) objective function ∆ ) change in a property ) characteristic energy φ ) fugacity coefficient γ ) activity coefficient Γ ) Kihara pair interaction potential µ ) chemical potential Θ ) fraction of cage occupation σ ) core distance at zero potential ν ) number of cages per water molecule in the crystal lattice Superscripts and Subscripts calc ) calculated d ) depression exp ) experiment fit ) fitted g ) guest H ) hydrate phase IG ) ideal gas phase L ) liquid water nfp ) normal freezing point 0 ) reference state sat ) saturated conditions vap ) vaporization V ) vapor phase w ) water R ) ice phase β ) empty hydrate phase π ) ice or liquid water
Literature Cited (1) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; Marcel Decker: New York, 1998. (2) Hatzikiriakos, S. G.; Englezos, P. Permafrost Melting and Stability of Offshore Methane Hydrates Subject to Global Warming. Int. J. Offshore Polar Eng. 1994, 4, 162. (3) Hatzikiriakos, S. G.; Englezos, P. The Relationship between Global Warming and Methane Gas Hydrates in the Earth. Chem. Eng. Sci. 1993, 48, 3963. (4) Brewer, P. G.; Friederich, G.; Peltzer, E. T.; Orr, F. M., Jr. Direct Experiments on the Ocean Disposal of Fossil Fuel CO2. Science 1999, 284, 943. (5) Brewer, P. G.; Orr, F. M., Jr.; Friederich, G.; Kvenvolden, K. A.; Orange, D. L. Gas Hydrate Formation in the Deep Sea: In situ Experiments with Controlled Release of Methane, Natural Gas, and Carbon Dioxide. Energy Fuels 1998, 12, 183. (6) McMullan, R. K.; Jeffrey, G. A. Polyhedral Clathrate Hydrates. IX. Structure of Ethylene Oxide Hydrate. J. Chem. Phys. 1965, 42, 2745. (7) Mak, T. C. W.; McMullan, R. K. Polyhedral Clathrate Hydrates. X. Structure of the Double Hydrate of Tetrahydrofuran and Hydrogen Sulfide. J. Chem. Phys. 1965, 42, 2732. (8) Ripmeester, J. A.; Tse, J. S.; Ratcliffe, C. I.; Powell, B. M. A New Clathrate Hydrate Structure. Nature 1987, 325, 135. (9) John, V. T.; Holder, G. D. Contribution of 2nd and Subsequent Water Shells to the Potential-Energy of Guest-Host Interactions in Clathrate Hydrates. J. Phys. Chem. 1982, 86, 455. (10) John, V. T.; Papadopoulos, K. D.; Holder, G. D. A Generalized-Model for Predicting Equilibrium Conditions for Gas Hydrates. AIChE J. 1985, 31, 252. (11) Sparks, K. A.; Tester, J. W.; Cao, Z.; Trout, B. L. Configurational Properties of Water Clathrates: Monte Carlo and Multidimensional Integration versus the Lennard-Jones and Devonshire Approximation. J. Phys. Chem. B 1999, 103, 6300. (12) van der Waals, J. H.; Platteeuw, J. C. Adv. Chem. Phys. 1959, 2, 1. (13) Parrish, W. R.; Prausnitz, J. M. Dissociation Pressures of Gas Hydrates Formed by Gas Mixtures. Ind. Eng. Chem., Process Des. Dev. 1972, 11, 26. (14) Zele, S. R.; Lee, S. Y.; Holder, G. D. A Theory of Lattice Distortion in Gas Hydrates. J. Phys. Chem. B 1999, 103, 10250. (15) Tee, L. S.; Gotoh, S.; Stewart, W. E. Molecular Parameters of Normal Fluids: The Kihara Potential with Spherical Core. Ind. Eng. Chem. Fundam. 1966, 5, 363.
(16) Berecz, E.; Balla-Achs, M. Gas Hydrates; Studies in Inorganic Chemistry. 411; Elsevier: New York (English Translation), 1983. (17) Von Stachelberg, M.; Mu¨ller, H. R. Z. Elektrochem. 1954, 58, 25. (18) Hwang, M.-J.; Holder, G. D.; Zele, S. R. Lattice Distortion by Guest Molecules in Gas-Hydrates. Fluid Phase Equilib. 1993, 83, 437. (19) Sloan, E. D.; Khoury, F. M.; Kobayashi, R. Water-Content of Methane Gas in Equilibrium with Hydrates. Ind. Eng. Chem. Fundam. 1976, 15, 318. (20) Ng, H.-J.; Robinson, D. B. Method for Predicting the Equilibrium Gas-Phase Water-Content in Gas-Hydrate Equilibrium. Ind. Eng. Chem. Fundam. 1980, 19, 33. (21) Anderson, F. E.; Prausnitz, J. M. Inhibition of Gas Hydrates by Methanol. AIChE J. 1986, 32, 1321. (22) Chen, G.-J.; Guo, T.-M. Thermodynamic Modeling of Hydrate Formation Based on New Concepts. Fluid Phase Equilib. 1996, 122, 43. (23) Chen, G.-J.; Guo, T. M. A New Approach to Gas Hydrate Modeling. Chem. Eng. J. 1998, 71, 145. (24) Ru˘ zˇicˇka, K.; Majer, V. Simple and Controlled Extrapolation of Vapor Pressures Toward the Triple Point. AIChE J. 1996, 46, 1723. (25) Carnahan, B.; Luther, H. A.; Wilkes, J. O. Applied Numerical Methods; Wiley: New York, 1969. (26) Stryjek, R; Vera, J. H. PRSVsan Improved Peng-Robinson Equation of State for Pure Compounds and Mixtures. Can. J. Chem. Eng. 1986, 64, 323. (27) Larsen, B. L.; Rasmussen, P.; Fredenslund, A. A Modified UNIFAC Group-Contribution Model for Prediction of Phase Equilibria and Heats of Mixing. Ind. Eng. Chem. Res. 1987, 26, 2274. (28) Dahl, S.; Fredenslund, A.; Rasmussen, P. The MHV2 Model: A UNIFAC-Based Equation of State Model for Prediction of Gas Solubility and Vapor-Liquid Equilibria at Low and High Pressures. Ind. Eng. Chem. Res. 1991, 30 (8), 1936. (29) Green, D. W.; Maloney, J. O. Perry’s Chemical Engineers' Handbook, 7th ed.; McGraw-Hill: New York, 1997. (30) Chase, M. W., Jr. NIST-JANAF Thermochemical Tables, 4th ed. J. Phys. Chem. Ref. Data, Monogr. 1998, 9, 1. (31) Lide, D. R.; Frederiskse, H. P. R. CRC Handbook, 76th ed.; CRC Press: Boca Raton, FL, 1995. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; HeadGordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.6; Gaussian, Inc.: Pittsburgh, PA, 1998. (33) Del Bene, J. E. Basis Set and Correlation Effects on Computed Hydrogen Bond Energies of the Dimers (AHn)2: AHn ) NH3, OH2, and FH. J. Chem. Phys. 1987, 86, 2110. (34) Scheiner, S. Ab-Initio Studies of Hydrogen-Bondssthe Water Dimer Paradigm. Annu. Rev. Phys. Chem. 1994, 45, 23. (35) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553. (36) Tse, J. S. Thermal-Expansion of the Clathrate Hydrates of Ethylene-Oxide and Tetrahydrofuran. J. Phys. (Paris) C1 1987, 48, 543. (37) Smithsonian Physical Tables, 9th ed.; Smithsonian Institution: Washington, DC 1954; p 283. (38) Hirai, H.; Kondo, T.; Hasegawa, M.; Yagi, T.; Yamamoto, Y.; Komai, T.; Nagashima, K.; Sakashita, M.; Fujihisa, H.; Aoki, K. Methane Hydrate Behavior under High Pressure. J. Phys. Chem. B 2000, 104, 1429. (39) Haffemann, D. R.; Miller, S. L. The Clathrate Hydrates of Cyclopropane. J. Phys. Chem. 1969, 73, 1392.
Received for review March 13, 2000 Revised manuscript received June 13, 2000 Accepted June 19, 2000 IE000322B
