Improving the Rigor and Consistency of the Thermodynamic Theory

Article pubs.acs.org/JPCC

Improving the Rigor and Consistency of the Thermodynamic Theory for Clathrate Hydrates through Incorporation of Movement of Water Molecules of Hydrate Lattice Srikanth Ravipati and Sudeep N. Punnathanam* Department of Chemical Engineering, Indian Institute of Science, Bangalore 560012, Karnataka India ABSTRACT: Current applications of statistical thermodynamic theories for clathrate hydrates do not incorporate the translational and rotational movement of water molecules of the hydrate lattice in a rigorous manner. Previous studies have shown that the movement of water molecules has a significant effect on the properties of clathrate hydrates. In this Article, a method is presented to incorporate the effect of water movement with as much rigor as possible. This method is then used to calculate the Langmuir constant of the guest species in a clathrate hydrate. Unlike previous studies on modeling of clathrate hydrate thermodynamics, the method presented in this paper does not regress either the intermolecular potentials or the properties of the empty hydrate from clathrate phase equilibria data. Also the properties of empty hydrate used in the theory do not depend on the nature and composition of the guest molecules. The predicted phase equilibria from the resulting theory are shown to be highly accurate and thermodynamically consistent by comparing them with the phase equilibria computed directly from molecular simulations.

1. INTRODUCTION Clathrate hydrates are inclusion compounds formed by water and typically small gas molecules, such as methane, ethane, propane, argon, nitrogen, carbon-dioxide, etc. These materials find applications in many scientific and technological fields. Hence it is of great interest to accurately predict the thermodynamic properties of clathrate hydrates. Traditionally, the thermodynamics of clathrate hydrates is described by van der Waals and Platteuw (vdWP) theory1 which uses a cell theory approach to model the clathrate hydrates. As currently used, the van der Waals and Platteuw (vdWP) theory suffers from many drawbacks. In the recent years, molecular simulations have been widely used to compute various properties of clathrate hydrates. Molecular simulations provide exact results for a given molecular model and are able to overcome the limitations of various theories for clathrate hydrates. However, molecular simulations of clathrate hydrates are still computationally intensive and developments leading to improvement of current theories for clathrate hydrates are highly desirable. In this study, we focus on improving the vdWP theory for predicting thermodynamic properties of clathrate hydrates. In the original vdWP theory,1 the hydrogen bonded water lattice is modeled as a rigid adsorbent, and the cavities act as adsorption sites with at most one guest molecule per cavity. The interactions of the guest molecule are limited to the nearest neighbor water molecules only. The guest interactions with the water molecules beyond nearest neighbors and the other guest molecules are neglected. The prediction of clathrate properties using the vdWP theory requires knowledge of the © 2015 American Chemical Society

intermolecular interactions and the equation of state of the empty hydrate. To match the experimental data on phase equilibria, the intermolecular potential parameters and the properties of the empty hydrate are usually regressed from the phase equilibria. When this approach was applied to the original theory, John et al.2 and Cao et al.3 showed that the interaction potential parameters regressed using vdWP theory differ from the independent estimates of the parameters4,5 obtained from viscosity measurements and second virial coefficient data. Also there was a wide variation in the estimated values of the properties of the empty hydrate.6 These problems suggested that the success of the theory is due to the presence of the large number of parameters and the theory acts as a data correlator.7 To overcome these problems, there have been a lot of efforts over the years to improve upon the original version of the vdWP theory. John and Holder,8 Rodger,9 Sparks and Tester,10 Kvamme et al.,11 and Klauda and Sandler12−14 showed that the neglected long-range interactions in the theory are significant enough to influence the phase equilibrium predictions. Accordingly, significant improvement in the performance of the theory was achieved through the inclusion of long-range guest−water interactions 8−10 and guest−guest interactions.11−14 The use of rigorous Monte Carlo integration in the evaluation of the Langmuir constant instead of the spherical shell approximation was suggested by Natarajan and Bishnoi15 and Sparks et al.16 Received: February 18, 2015 Revised: April 17, 2015 Published: May 7, 2015 12365

DOI: 10.1021/acs.jpcc.5b01662 J. Phys. Chem. C 2015, 119, 12365−12377

Article

The Journal of Physical Chemistry C In recent years, the use of ab initio calculations3,17−19 has allowed us to obtain the intermolecular potential parameters from first-principles and thus avoiding its regression from phase equilibria. However, regression from phase equilibria is still needed for determining the properties of the empty hydrate. The reason is that the empty hydrate is a hypothetical material and hence its properties cannot be directly measured in a laboratory. Such an indirect estimation of the empty hydrate properties will be valid only if the theory of the clathrate hydrates is accurate. Studies by Klauda and Sandler12−14 and Garapati and Anderson20 have shown that the predictions of the theory can match the experimental data only if the properties of the empty hydrate are allowed to vary depending on the nature and the composition of the guest molecules in the hydrate phase. Since the empty hydrate is only composed of water molecules, this dependence is an indication of the inaccuracy of the current theory for clathrate hydrates. The modifications proposed in this paper to the current existing theory for clathrate hydrates allows us to remove the guest dependence on the properties of the empty hydrate. In 2007, Wierzchowski and Monson21 performed semigrand isothermal−isobaric simulations of methane hydrate and computed the variation of the occupancy of hydrate with fugacity of methane. For the purpose of discussion, we refer to this occupancy versus fugacity variation as the “hydrate occupancy curve”. This hydrate occupancy curve cannot be obtained from experiments since the clathrate is thermodynamically stable only at high fugacity values of guest species. However, the hydrate occupancy curve is predicted by the vdWP theory and can be compared with the simulation results. Since molecular simulations give the exact solution for a given molecular model, any difference between the theory and the simulations can be fully attributed to the shortcomings in the theory. Comparison of the theory with experiments is less helpful since any mismatch between the two sets of data could be due to either inaccuracy of the theory itself or inaccuracy of the molecular model. An accurate estimation of the hydrate occupancy curve is critical for the success of the theory since the chemical potential of water in the hydrate phase is obtained by integrating this curve. Wierzchowski and Monson21 compared the hydrate occupancy predicted from the theory with those computed directly from simulations. The theory included guest−water interactions beyond the neighboring shell, as well as guest−guest interactions using mean field approximation. In spite of the rigor in the theory, the predicted hydrate occupancy results did not match the simulation results. This was attributed to the neglect of translational and rotational movement of water molecules in the hydrate lattice. This issue was further studied in detail by Pimpalgaonkar et al.,22 where it was shown that the neglect of water movement led to a difference of 10−15 K in the values of dissociation temperature for methane hydrates and even larger for ethane hydrates. The significance of the movement of water molecules toward clathrate properties have been known for sometime. A series of studies published by Tanaka and co-workers23−27 used lattice dynamics to include the effects of water movement on the freeenergy of the clathrate hydrate. Although they reported an improvement in prediction of the vdWP theory, it is not clear from these studies whether the mismatch between the theory and experiment was due to shortcomings in the theory itself or the molecular model used for the water and guest molecules. A different approach to incorporate the effects of movement of water molecules was proposed by Ravipati and Punnatha-

nam,28 where the Langmuir constant was computed from isothermal−isobaric ensemble simulations of empty hydrate. The proposed method was shown to give highly accurate predictions of the hydrate occupancy curve. In this study, the method is applied to predict the phase equilibria of hydrates formed by the hydrocarbons methane, ethane, propane, and isobutane. The properties of the empty hydrate at reference conditions are computed directly from simulations using the Frankel−Ladd technique29−31 for free-energy estimation of solids. These predictions are then compared to phase equilibria calculated from molecular simulations. Such a comparison has been made possible due to recent advances22,28,32,33 in computing clathrate phase equilibria directly from molecular simulations. The hydrate systems considered in this study are (i) single hydrates of methane, ethane, propane, and isobutane, (ii) binary hydrates of methane−ethane and methane− propane, and (iii) quaternary hydrate of methane−ethane− propane−isobutane.

2. THEORY Consider a clathrate hydrate that is in equilibrium with a multicomponent vapor mixture. Extending eq 10 of ref 34 for a multicomponent vapor mixture, the fractional occupancy, θiα, of the guest molecule of type i in the cavity of type α can be written as θiα = Ui̅ α =

Ciαfi ̂ exp( −Ui̅ α /kBT ) 1 + ∑j Cjαfĵ exp( −Uj̅ α /kBT )

(1)

∑ ∑ θjβUiαjβ j

(2)

β

In the above equations, Uiαjβ is the interaction of a guest molecule of type i in cavity of type α with guest molecules of type j occupying all the surrounding cavities of type β; U̅ iα is the estimate of the interaction of a guest molecule of type i inside cavity of type α with rest of the guest molecules using the mean-field approximation; Ciα is called the Langmuir constant which describes the interaction of the guest molecule of type i in cavity of type α with all the water molecules in the hydrate lattice; fî is the fugacity of the guest molecule in the vapor phase; kB is the Boltzmann’s constant; and T is the absolute temperature. The total occupancy of guest molecule of type i is obtained from fractional occupancies as θi =

⎛ Mα ⎞ ⎟θ iα M⎠

∑ ⎜⎝ α

(3)

where Mα is the number of cavities of type α and M is the total number of cavities. The Gibbs−Duhem equation for clathrate hydrates with multiple guest molecules at constant temperature and pressure is given by H

Nw d ln fŵ +

H

∑ Mθj d ln f ĵ j

=0 (4)

where Nw is the number of molecules of water, fŵ is the fugacity of water, and fĵ is the fugacity of guest molecule of type j. The superscript H denotes the hydrate phase. Substitution of eq 3 in the Gibbs−Duhem equation followed by integration from a reference state would give that change in the fugacity of water in the hydrate phase because of changes in the guest occupancy. If the reference state is the empty hydrate lattice and guest− guest interactions were neglected then the integration can be 12366


Article

The Journal of Physical Chemistry C performed analytically giving the well-known equation for fugacity of water in the hydrate phase,35 that is H ln fŵ = ln f wβ −

∑ α

Mα ln(1 + Nw

H

∑ Cjαf ĵ

Table 1. Lennard-Jones Parameters and Partial Charges for Nonbonded Interactions molecule

)

j

(5)

where the superscript β denotes an empty hydrate. If the guest−guest interactions are included as in eq 1, then the integration of the Gibbs−Duhem equation has to be performed numerically.36 The estimation of the hydrate occupancy requires calculation of the Langmuir constant, Ciα, and the guest−guest interaction, Uiαjβ. Traditionally, the water lattice in the clathrate hydrate is considered to be rigid and the guest molecules move inside the clathrate cavity subject to a constant potential field caused by water molecules. Previous studies21−28,34 have demonstrated that this approximation is source of significant error in the theory. To incorporate the movement of the water molecules, Ravipati and Punnathanam34 proposed that the Langmuir constant be calculated from an isothermal−isobaric ensemble of configurations of an empty hydrate. Accordingly, the Langmuir constant is evaluated as an ensemble average34 during the simulation of the empty hydrate in an isothermal−isobaric ensemble, that is Ciα =

1 ⟨ kBT

∫V e−u (r)/k T dr⟩N ,P ,T i

w

a

(6)

∫

int i (ω)/ kBT

∫V ∫ e−u (r,ω)+u i

α

dω dr⟩Nw , P , T

q

3.1668

0.0 0.5897 −1.1794 0.0 0.0 0.0

98.0 46.0 10.0

3.75 3.95 4.68

bend group

θ0 (degrees)

kθ/kB (K)

CHx−CH2−CHy CHx−CH−CHy

114 112

62500 62500

The form of the potential is Ubend(θ) = (kθ/2)(θ − θ0)2.

sites are obtained by Lorentz−Berthelot rules. This force field is chosen because the calculated phase diagrams of methane and ethane hydrates using this model are reasonably close to the experimental data.32,33,41 In all the simulations of sI and sII hydrate, a system size consisting of 2 × 2 × 2 unit cells is used. The initial hydrate proton disordered configurations are taken from Takeuchi et al.,42 which have nearly zero net dipole moment and the lowest potential energy. The sI hydrate consists of 16 small (512) and 48 large (51262) cavities and the sII hydrate consists of 128 small (512) and 64 large (51264) cavities. The electrostatic interactions between water molecules are computed using the Ewald summation technique. All the simulations use a fixed cutoff radius of 11.5 Å for the sI hydrate and 14 Å for the sII hydrate while computing Lennard-Jones and real-space Ewald interactions. The Monte Carlo moves consist of (i) translation and rotation around the center of mass for water; (ii) translation, insertion, and deletion for methane; (iii) translation, rotation around the center of mass, insertion, and deletion for ethane; and (iv) translation, rotation around the center of mass, insertion, deletion, and regrow for propane and isobutane. These moves are chosen with equal probability. Simulations at constant pressure also involve isotropic volume changes of the simulation box. The insertion, deletion and regrow moves for propane and isobutane are performed using the coupleddecoupled configurational bias Monte Carlo (CD-CBMC) technique of Martin and Siepmann.40 The values of the CDCBMC parameters40 are nchLJ = 4, nchtor = 100, and nchbend = 1000. The CD-CBMC code for these simulations have been developed in our group. It has been tested by performing Gibbs-Ensemble simulations of the alkanes and reproducing the results published in Martin and Siepmann.39,40 3.2. Computation of Ciα and Uiαjβ. Ciα is computed by performing test insertions of the guest molecules during isothermal−isobaric simulations of the empty hydrate lattice. The simulations consist of 5 × 104 equilibration and 5 × 104 production cycles. During the production step, 5 × 103 test insertions of the guest molecule are performed for every 10th cycle. The test insertion method is analogous to the Widom’s insertion technique for calculation of excess chemical potential. The position to insert the molecule is chosen at random uniformly within the entire simulation box. For rigid polyatomic molecules, the orientation is chosen randomly

int 1 [ e−ui (ω)/ kBT dω]−1 kBT

⟨

σ (Å)

106.1

Table 2. Intramolecular Potential Parameters for Bond Angle Bending in Alkanesa

where Vα is the volume of the cavity of type α and ui(r) is the interaction of the guest molecule with the water molecules in the hydrate lattice and r is center of mass of the guest molecule. The angular brackets represent an average over positions of the water molecules of an empty hydrate lattice in an isothermal− isobaric ensemble. The above equation is valid for a monatomic guest molecule. In case of polyatomic molecules, the equation for the Langmuir constant becomes (see AppendixA for derivation) Ciα =

ϵ/kB (K)

O H M CH3 CH2 CH

alkanes

B

α

atom/group

water

(7)

uint i

where is the intramolecular potential energy of the guest molecule and ω is a vector that contains the generalized coordinates37 representing all degrees of freedom of the guest molecule excluding the translational ones.

3. METHODOLOGY 3.1. Simulation Details. Water is modeled using four site TIP4P/Ice38 model, which is parametrized to reproduce the phase diagram of ice. It is a four site model where three of the sites lie on the oxygen and two hydrogen atoms and the fourth site, M, lies on the bisector in the HOH plane. The ∠HOH is equal to 104.52°, and the O−H and O−M distances are 0.9572 and 0.1577 Å, respectively. The hydrocarbons are modeled using the TRaPPE-UA force field.39,40 In this force field each CHx group is modeled as a united atom interacting with other sites via Lennard-Jones potential. The bond length between two adjacent CHx groups is fixed and equal to 1.54 Å. These Lennard-Jones interaction parameters are parametrized to reproduce the corresponding VLE data of the hydrocarbons. The details of the force fields used are given in Tables 1 and 2. The Lennard-Jones parameters for interactions between unlike 12367


Article

The Journal of Physical Chemistry C

Table 3. Regressed Values of the Coefficients in the Equation Describing the Temperature and Pressure Dependence of the Langmuir Constant structure

cavity

guest

a

b

c (K) × 10−3

d (K)

sI sI sI sI sII sII sII sII sII sII

12

methane methane ethane ethane methane methane ethane ethane propane isobutane

−24.37 −22.68 −26.49 −24.93 −24.37 −21.18 −25.66 −22.95 −24.77 −25.99

0.99 0.99 0.99 0.99 0.99 1.00 0.98 1.00 0.99 0.98

2.87 2.83 1.87 3.96 2.86 2.55 2.04 4.01 4.89 5.11

−0.37 1.57 −9.37 1.12 0.77 1.64 −7.56 2.31 2.04 −1.98

5 51262 512 51262 512 51264 512 51264 51264 51264

of the Langmuir constant where the positions of the water molecules were kept fixed. The predictions of the phase equilibria from these two sets of Ciα allows us to quantify the significance of movement of water molecules toward properties of clathrate hydrates. The Langmuir constant calculated by keeping the water molecules fixed is just a function of temperature. The dimensions of the simulation boxes are 24.06 Å × 24.06 Å × 24.06 Å for sI hydrate and 34.62 Å × 34.62 Å × 34.62 Å for sII hydrate and the configurations were taken from Takeuchi et al.42 The computed values are fitted to the form

and uniformly. The test insertions in the case of rigid molecules such as methane and ethane are used to compute the average value ⟨e−ui/kBT⟩. Since the test insertion is performed over the entire simulation box, this average is then multiplied by ⟨V⟩/Mα where ⟨V⟩ is the average volume of the simulation box. Accordingly, Ciα =

⟨V ⟩ −ui / kBT ⟨e ⟩ MαkBT

(8) −ui/kBT

For this method to work, the value of e should be insignificant when inserted in regions outside the cavity of interest. This is ensured by placing a hard spheres in all cavities that are not of type i during the simulation of the empty hydrate. The size of the hard sphere is chosen to be large enough such that it blocks the insertion of guest molecules into these cavities, but not too large to distort the lattice structure of the clathrate hydrate. If a guest molecule is inserted into a cavity containing a hard sphere, the value of u is ∞, and hence, the test insertion does not contribute to the average ⟨e−ui/kBT⟩. In these simulations, a value of 2.5 Å is used as the contact distance for the hard sphere interaction. In the case of molecules with internal degrees of freedom, such as propane and isobutane, it is computationally more efficient to bias the insertion of the molecules while calculating the Langmuir constant. The biasing technique uses the algorithm of Rosenbluth and Rosenbluth and is analogous to the method used by Smit and Siepmann43 and Maginn et al.37 for the calculation of Henry’s coefficient for adsorption of linear alkanes in zeolites. Accordingly Ciα =

⟨V ⟩ ⟨Wiα⟩ MαkBT ⟨WiIG⟩

ln Ciα = ar +

where Ciα has units of Pa and T has units of K. The values of ar and br are calculated using least-squares regression and are given in Table 4. Table 4. Regressed Values of the Coefficients in the Equation Describing the Temperature Dependence of the Langmuir Constant Calculated Assuming a Rigid Water Lattice

where Wiα is the Rosenbluth weight for insertion of the guest molecule of type i in a cavity of type α and WIG is the i Rosenbluth weight for insertion of the guest molecule of type i in an ideal gas. The Langmuir constant, Ciα, is a function of temperature and pressure. Hence the simulations for computing Ciα are performed in the temperature range of 190−310 K and the pressure range of 0.1−800 bar. The computed values are then fitted to the following form c d + ln P T T

(11) −1

(9)

ln(CiαP) = a + b ln P +

br T

structure

cavity

guest

ar

br (K)

sI sI sI sI sII sII sII sII sII sII

512 51262 512 51262 512 51264 512 51264 51264 51264

methane methane ethane ethane methane methane ethane ethane propane isobutane

2867.60 2714.02 2175.69 3905.44 2843.53 2403.16 1990.58 3783.07 4705.13 5014.25

−23.06 −21.79 −25.35 −23.92 −23.30 −20.43 −26.24 −21.67 −23.49 −25.28

The long-range guest−guest interaction, Uiαjβ, is computed from NVT simulations of the clathrate hydrate with a molecule i in cavity α and molecule j in cavity β. Uiαjβ is the average guest−guest interaction as computed from the simulation. Although the value of Uiαjβ depends on T and V, this dependence is quite weak. Hence a constant value of Uiαjβ computed at T = 273 K, in an NVT simulation is used for all the computations of the occupancy. The dimensions of the simulation boxes for these simulations are 24.06 Å × 24.06 Å × 24.06 Å for sI hydrate and 34.62 Å × 34.62 Å × 34.62 Å for sII hydrate, and the configurations were taken from Takeuchi et al.42 The computed values of Uiαjβ are given in Tables 5 and 6

(10)

−1

where Ciα has units of Pa , P has units of Pa, and T has units of K. The values of a, b, c, and d are calculated using least-squares regression and are given in Table 3. In addition to these computations, Ciα was calculated using the traditional definition 12368


Article

The Journal of Physical Chemistry C Table 5. Computed Values of Uiαjβ for sI Hydrate i

α

j

β

Uiαjβ/kB (K)

methane methane methane methane methane methane methane methane ethane ethane ethane ethane ethane ethane ethane ethane

512 512 512 512 51262 51262 51262 51262 512 512 512 512 51262 51262 51262 51262

methane methane ethane ethane methane methane ethane ethane methane methane ethane ethane methane methane ethane ethane

512 51262 512 51262 512 51262 512 51262 512 51262 512 51262 512 51262 512 51262

−15.18 −230.92 −26.09 −401.52 −76.97 −172.00 −138.55 −294.24 −26.09 −415.65 −44.69 −708.75 −133.84 −294.25 −236.25 −498.29

integration gives the change in the fugacity of water in the hydrate phase from a reference state which is the empty hydrate phase. The computation of fHŵ then requires the values of f βw at different values of T and P. This is achieved by constructing an equation of state from data obtained through isothermal− isobaric simulations of empty hydrates. The temperature and pressure ranges at which these simulations are performed are given in Table 7. The properties sampled from these Table 7. Range of Temperature and Pressure Where Isothermal−Isobaric Ensemble Simulations Are Performed for Different Phases in This Study

i

α

j

β

Uiαjβ/kB (K)

512 512 512 512 512 512 51264 51264 51264 51264 51264 51264 512 512 512 512 512 512 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264 51264

methane methane ethane ethane propane isobutane methane methane ethane ethane propane isobutane methane methane ethane ethane propane isobutane methane methane ethane ethane propane isobutane methane methane ethane ethane propane isobutane methane methane ethane ethane propane isobutane

512 51264 512 51264 51264 51264 512 51264 512 51264 51264 51264 512 51264 512 51264 51264 51264 512 51264 512 51264 51264 51264 512 51264 512 51264 51264 51264 512 51264 512 51264 51264 51264

−208.04 −86.26 −344.67 −148.87 −211.12 −279.86 −172.51 −52.90 −311.34 −86.90 −127.62 −168.93 −346.08 −155.67 −578.65 −260.06 −372.12 −493.36 −297.75 −86.90 −520.13 −150.90 −215.54 −285.31 −422.23 −127.62 −744.25 −215.54 −307.88 −407.55 −559.72 −168.93 −986.73 −285.32 −407.55 −539.48

liquid water

ice

sI empty hydrate

sII empty hydrate

T (K) P (bar)

270.0−310.0 5−1100

140.0−272.0 0.05−25

140.0−315.0 0.05−700

140.0−315.0 0.05−700

simulations are the molar residual enthalpy, HR, and the molar volume, V. These equations are fitted to the following form

Table 6. Computed Values of Uiαjβ for sII Hydrate methane methane methane methane methane methane methane methane methane methane methane methane ethane ethane ethane ethane ethane ethane ethane ethane ethane ethane ethane ethane propane propane propane propane propane propane isobutane isobutane isobutane isobutane isobutane isobutane

phase

V = a1 + a 2T + a3P + a4TP HR = a1P +

⎛ a3 ⎞ 2 ⎜ ⎟P + a + a T 5 6 ⎝2⎠

These forms are chosen such that it fits the simulation data and also satisfies the Maxwell relationship22 ⎛ ∂H /RT 2 ⎞ ⎛ ∂V /RT ⎞ ⎟ ⎟ =⎜ −⎜ R ⎠T ⎝ ∂T ⎠ P ⎝ ∂P

The fugacity of water is then obtained using thermodynamic integration as follows ⎛ fβ ⎞ ln⎜⎜ βw ⎟⎟ = ⎝ f w,0 ⎠ =

∫T

T

0

⎛ 1 ⎞ ⎟ + HR d⎜ ⎝ RT ⎠

∫P

P

0

⎛V ⎞ ⎜ ⎟ dP ⎝ RT ⎠

(12)

P ⎞ ⎛a ⎞ ⎛ a1 ⎞⎛ P ⎟⎜ (13) − 0 ⎟ + ⎜ 2 ⎟(P − P0) ⎝ R ⎠⎝ T T0 ⎠ ⎝ R ⎠ P2 ⎞ ⎛ a ⎞ ⎛ a ⎞⎛ P 2 + ⎜ 3 ⎟⎜ − 0 ⎟ + ⎜ 4 ⎟(P 2 − P02) ⎝ 2R ⎠⎝ T T0 ⎠ ⎝ 2R ⎠ ⎛ a ⎞⎛ 1 1 ⎞ ⎛a ⎞ ⎛ T ⎞ + ⎜ 5 ⎟⎜ − ⎟ − ⎜ 6 ⎟ln⎜ ⎟ ⎝ R ⎠⎝ T T0 ⎠ ⎝ R ⎠ ⎝ T0 ⎠ ⎜

The value of the constants a1, ..., a6 as regressed from simulation data for both sI and sII hydrates are given in Table 8. The fugacity of the empty hydrate at reference conditions, f βw,0 is computed using the Frenkel-Ladd method.21,29−31 The Table 8. Values of the Coefficients in the Equations of State for sI Empty Hydrate and sII Empty Hydrate coefficients a1 a2 a3 a4 a5 a6

3.3. Computation of ln f βw. The values of Ciα and Uiαjβ computed from the methods described above are substituted in eqs 1 to 3 compute the hydrate occupancy from the theory. Substitution of these computed values into eq 4 followed by 12369

(m3 mol−1) (m3 mol−1 K−1) (m3 mol−1 Pa−1) (m3 mol−1 Pa−1 K−1) (J mol−1) (J mol−1 K−1)

sI empty hydrate

sII empty hydrate

21.580 × 10−6 38.978 × 10−10 −92.712 × 10−17 −52.383 × 10−19 −68.691 × 103 23.454

21.885 × 10−6 38.276 × 10−10 −93.671 × 10−17 −54.191 × 10−19 −68.826 × 103 23.389


Article

The Journal of Physical Chemistry C value of the spring constant is 6.25 × 106 K/Å2 for the translation part and 6.25 × 106 K for the rotation part. The simulation details and the computed value of f βw,0is given in Table 9.

molecular simulations for clathrate of hydrocarbon mixtures is described in detail in Ravipati and Punnathanam.33 In brief, semigrand isothermal−isobaric simulations of clathrate hydrates are used to compute fHŵ for a hydrate in equilibrium with the vapor phase. The equation for fHŵ is as follows (see ref 33 for derivation).

Table 9. Fugacity of Water in Liquid State, Ice, sI Empty Hydrate, and sII Empty Hydrate at Reference Conditionsa liquid water ice sI empty hydrate sII empty hydrate

T0 (K)

ρ0 (kg m−3)

P0 (bar)

f w,0 (Pa)

280 260 280 280

988.47 906.80 794.35 784.31

20.00 2.531 32.18 45.06

24.35 3.219 39.45 38.24

⎛ fH ⎞ ln⎜⎜ Hw ⎟⎟ = ⎝ f w,0 ⎠

liquid water

ice 19.123 × 10−6 28.085 × 10−10 −48.476 × 10−17 −54.115 × 10−19 −69.487 × 103 22.365

⎛ V H − ∑n x HV̅ v ⎞ i i=1 i ⎜⎜ ⎟⎟dP H x w RT ⎝ ⎠

(14)

4. RESULTS AND DISCUSSION The accuracy of the proposed method to compute Langmuir Constants is first tested by comparing the predictions of the hydrate occupancy from eqs 1 to 3 with those computed from isothermal−isobaric semigrand simulations. Figure 1 shows such comparison for clathrate hydrates formed by methane, ethane, propane, and isobutane under conditions that are typical of clathrate formation. The methane and ethane occupancy curves are computed for sI hydrate and those for propane and isobutane are computed for sII hydrate. The agreement between the predictions of the theory and those computed from simulations are excellent. Also shown in Figure 1 are predictions from the theory when the water lattice is kept rigid. The assumption of rigid water lattice leads to overestimation of the occupancy which was first reported in Wierzchowski and Monson.21 Since, the area under the occupancy curve gives the magnitude of the change in the fugacity of water in the hydrate phase due to the presence of the guest molecules, accurate predictions of the phase equilibria from the theory is possible only if the prediction of the hydrate occupancy is accurate. Such a comparison is discussed next. The phase diagrams showing dissociation pressures and hydrate occupancy as computed from molecular simulations and theory are shown in Figures 2−6. Also shown in these figures are the predictions using estimates of the Langmuir constant for a rigid water lattice. Figures 2 and 3 shows the phase diagram of sI hydrates formed from pure methane and pure ethane, respectively. The three phases in equilibrium are liquid water, sI hydrate and hydrocarbon vapor. Figures 4 and 5 shows the phase diagram of sII hydrates formed from pure propane and pure isobutane, respectively. In these systems, the three phases in equilibrium are ice, sII hydrate and hydrocarbon vapor. From these figures, one can clearly observe that the dissociation pressures and cage occupancies are accurate when the effect of the movement of water molecules is included. Figure 8 shows the average absolute relative deviation (AARD) in the estimate of the dissociation pressures calculated as follows:

Table 10. Values of the Coefficients in the Equations of State for Liquid Water and Ice 19.074 × 10−6 −30.365 × 10−10 −23.990 × 10−15 56.020 × 10−18 −74.170 × 103 58.817

∫P

P

where the superscript v represents the vapor phase, n is the number of components in the vapor phase, and xi is the mole fraction of component i. The overline indicates partial molar properties of component i in the vapor phase and can be computed from the SRK equation of state. The values of fHŵ obtained from eq 14 are then used to obtain the exact phase equilibria for the systems under study. The comparison of the two sets of phase equilibria (from theory and simulation) is discussed in the following section.

The computed values of Ciα, Uiαjβ, and f βw is then used to calculate the fugacity of water in the hydrate phase, fHŵ . In order to assess the accuracy of the theory, the prediction of fHŵ from the theory has to be compared to the exact results of a molecular simulation. A more illuminating comparison, however, would be between the three phase equilibria for clathrate hydrates predicted from the theory with those computed directly from molecular simulations. 3.4. Computation of the Three-Phase Equilibria. The calculations of the three-phase hydrate equilibria require fHŵ and the equations of state for the other phases, namely, the hydrocarbon vapor and liquid water (or ice). The hydrocarbon vapor phase is modeled using the Soave−Redlich−Kwong (SRK) equation of state. The critical properties and the acentric factor needed for the SRK equation are obtained from VLE computed using TRaPPE-UA model for n-alkanes.40 The equations of state for liquid water and ice are computed from isothermal−isobaric Monte Carlo molecular simulations similar to the method used to compute equation of state for the empty sI and sII hydrates. The temperature and pressure dependence of the fugacities is similar to the form in eq 13. The temperature and pressure ranges at which these simulations are performed are given in Table 7. The regressed values of a1, ..., a6 are given in Table 10, and the fugacities at reference

coefficients

⎛ HH − ∑n x HH̅ v ⎞ ⎛ 1 ⎞ i,R i=1 i ⎟ ⎜⎜ R ⎟⎟d⎜ H x ⎝ ⎠ ⎝ RT ⎠ w

0

The value of the spring constant of the Einstein crystal in Frankel− Ladd simulations was kept equal to 6.25 × 106 Å−2 for translational part and 6.25 × 106 K for the rotational part.

(m3 mol−1) (m3 mol−1 K−1) (m3 mol−1 Pa−1) (m3 mol−1 Pa−1 K−1) (J mol−1) (J mol−1 K−1)

0

+

a

a1 a2 a3 a4 a5 a6

∫T

T

conditions are given in Table 9. The reference fugacities have been computed using the Hamiltonian-Integration31 and the Frenkel−Ladd method21,29−31 for liquid water and ice, respectively. Application of the phase equilibrium criteria to the equations of state and fHŵ computed using methods of the previous sections gives the theoretical prediction of the clathrate phase equilibria. These results are then compared with the phase equilibria directly computed from molecular simulations. The method for direct computation of the phase equilibria using 12370


Article


Figure 1. Guest−occupancy curves for clathrate hydrates of methane, ethane, propane, and isobutane. The circles represent the occupancy values obtained directly from semigrand isothermal−isobaric simulations and the lines represent the theory predictions using eqs 1−3. Red lines represent predictions when movement of water molecules are included and green lines represent predictions assuming rigid water lattice.

Figure 2. Variation of dissociation pressure and cage occupancy with temperature for a system consisting of liquid water, sI hydrate, and pure methane vapor in equilibrium with each other. The circles represent values computed directly from molecular simulations, the solid line represents values calculated from the theory when the effect of movement of water molecules is included, the dashed line represents values calculated from theory assuming a rigid water lattice and the triangles represent experimental data on dissociation pressures. The experimental data is obtained from Adisasmito et al.44 and McLeod and Campbell.45

AARD (%) =

∑# datapoints |Ptheory − Psimulation| /Psimulation #data points

the predictions of the dissociation pressure if one assumes a rigid water lattice, as the theory is currently used, are 59.32%, 62.66%, 46.89%, and 39.63%, respectively. The strategy to reduce this error is to either regress all the parameters from phase equilibria35 or to vary the properties of the empty hydrate according to the guest molecules12−14,20 if independent estimates of the intermolecular potential parameters are available. Such modifications, however, lead to introduction of thermodynamic inconsistencies in the theory, that is,

× 100 (15)

The errors in the estimation of dissociation pressures are 6.87%, 8.45%, 1.89%. and 13.58% for clathrate hydrates of methane, ethane, propane and isobutane, respectively. Considering the fact that no regression was used while applying the theory, the match between the predictions of the theory and the simulation results is very good. The corresponding errors in 12371


Article


Figure 3. Variation of dissociation pressure and cage occupancy with temperature for a system consisting of liquid water, sI hydrate and pure ethane vapor in equilibrium with each other. The meanings of the symbols and lines are the same as in Figure 2. The experimental data on dissociation pressures is obtained from Avlonitis.46 The simulations showed that a tiny amount of ethane occupies the small cavity (512) but only fractional occupancies in the large cavity (51262) are shown.

Figure 4. Variation of dissociation pressure and large cage (51264) occupancy with temperature for a system consisting of liquid water, sII hydrate and pure propane vapor in equilibrium with each other. The meanings of the symbols and lines are the same as in Figure 2. The experimental data on dissociation pressures is obtained from Holder and Godbole.47

Figure 5. Variation of dissociation pressure and large cage (51264) occupancy with temperature for a system consisting of liquid water, sII hydrate and pure isobutane vapor in equilibrium with each other. The meanings of the symbols and lines are the same as in Figure 2. The experimental data on dissociation pressures is obtained from Holder and Godbole.47

predictions of properties that were not used in the regression such as cage occupancies will be inaccurate. The accuracy and the thermodynamic consistency of the theory is maintained even for clathrate hydrates formed from hydrocarbon mixtures. Figure 6 shows dissociation pressures and hydrate occupancies for three mixtures of hydrocarbon vapor. The errors in the estimation of dissociation pressures are 2.33% for the binary methane−ethane sI hydrate, 8.23% for the binary methane-propane sII hydrate and 7.97% for the

quaternary methane−ethane−propane−isobutane sII hydrate. The corresponding errors for the theory that assumes rigid water lattice are 68.68%, 61.87%, and 63.20%, respectively. In addition to accurate predictions of dissociation pressures, the theory also give accurate predictions of the composition of the hydrate phase. The predictions of occupancies using the rigidlattice approximation are also shown in Figure 6. Although in general, the predictions are better when movement of water molecules are included, there are certain cases where the 12372


Article


Figure 6. Variation of dissociation pressure and cage occupancy with temperature for a system consisting of liquid water, hydrate phase, and hydrocarbon vapor of fixed composition. The meanings of the symbols and lines are the same as in Figure 2. The experimental data on dissociation pressures is obtained from Deaton and Frost48 for hydrates formed from vapor mixtures of methane + ethane and methane + propane and Mei et al.49 for hydrates formed from methane + ethane + propane + isobutane vapor mixture.

simulations are 276.17 and 274.67 K for methane mole fractions of 94.2% and 95.3%, respectively. Since the slopes of the two lines are very close to each other the location of the transition temperature is very sensitive to the accuracy of the triple point lines. Although the theory predictions are quite accurate, small deviations from the simulations result in a prediction of 278.41 and 278.49 K for methane mole fractions of 94.2% and 95.3%, respectively. Figure 7c shows the sI−sII phase boundary predicted using vdWP theory assuming both flexible and rigid water lattice. Also shown are the transitions points computed from simulations. 4.1. Development of Guest−Water Potentials. Figures 2−6 also show the experimental data on the dissociation

predictions from the rigid-lattice approximation are closer to simulations because of fortuitous cancellation of errors. Methane−ethane binary hydrate system shows structural transitions between sI−sII hydrates. These were qualitatively predicted by Hendriks et al.50 using vdWP theory and later confirmed via experiments.51−53 The accuracy of the prediction of the transition points is studies by computing the P−T diagram for methane−ethane mixture at mole fractions of 94.2% methane and 95.3% methane. Figure 7a and 7b shows the triple point lines with both sI and sII hydrates for these two mole fractions, respectively. For both these mole fractions, sII is stable at lower temperatures, and sI is stable at higher temperatures. The transition temperatures obtained from 12373


Article


Figure 7. Predictions of structural transitions for methane−ethane binary hydrate from vdWP theory compared with the simulations. The symbols represent data obtained from simulations, the solid line represents predictions from vdWP theory considering the effects of water movement and the dashed line represents predictions assuming rigid water lattice.

for vapors rich in the higher hydrocarbons. Given the high accuracy of the theory, these differences can then be attributed to the shortcomings in the molecular models for water and hydrocarbons used in the simulations, especially the guest-water interactions. The TIP4P/Ice model for water and the TRaPPEUA model for hydrocarbons have been parametrized to give accurate phase equilibria of water−ice system and hydrocarbon VLE, respectively. While these models do have their shortcomings, in our opinion, the main reason for the difference between the theory and experiments lies in the inaccuracy of the guest-water interactions used in simulations and theory. In this work, the cross-interactions between the oxygen molecule of water and the CHx group of the hydrocarbons have been obtained using the Lorentz−Berthelot rules in the simulations and theory. Studies by Docherty et al.54 and Ashbaugh et al.55 have shown that the Lorentz−Berthelot rules do not model the water−hydrocarbon interaction correctly and leads to systematic deviations from experimental data. Hence, we suspect that the use of the Lorentz−Berthelot rules is the main reason for differences in the values of the dissociation pressure between simulations and experiments. The differences between the theory and experiments can be reduced only with improved molecular models, especially the guest-water interactions. In this regard, the ab initio calculations such as those performed by Cao et al.,3 Klauda and Sandler,13 Anderson et al.,17,18 and Velaga and Anderson19 that estimate the guest−water interactions can bring the predictions of the theory closer to experimental data. Current practice is to test these models

Figure 8. Accuracy in the predictions of the hydrate dissociation pressure from theory of clathrate hydrates. The different types of hydrates are as follows: I, methane (100%); II, ethane (100%); III, propane (100%); IV, isobutane (100%); V, methane(56.4%) + ethane; VI, methane (88.3%) + propane; VII, methane (97.25%) + ethane (1.42%) + propane (1.08%) + isobutane.

pressures. The focus of the work presented so far has been on improving the theory for clathrate hydrates, especially through incorporation of the effects of water movement on the Langmuir constant. While the method developed so far gives very accurate and consistent predictions vis-a-vis simulations, the agreement with experiments is not as good. The agreement is reasonably close for hydrate formed from methane and methane-rich hydrocarbon vapors but get progressively worse 12374


Article

The Journal of Physical Chemistry C using the theory with rigid water lattice approximation. Since the rigid water lattice approximation is a source of large errors, future testing of molecular models should use a theory that include (i) the effect of movement of water molecules, (ii) guest−water interaction beyond the first shell, and (iii) guest− guest interactions.

Equation 16 can be rewritten as

5. CONCLUSIONS In this Article, we have proposed a method to compute the Langmuir constant that accounts for movement of water molecules of the hydrate lattice. We then demonstrated that accurate and thermodynamically consistent properties of clathrate hydrates can be predicted without regressing either the intermolecular potential parameters from the phase equilibria and using a single equation of state for the empty hydrate which is independent of the guest molecules. Both the dissociation pressure, as well as composition of the hydrate phase, and for single, as well as mixed, hydrates are accurately predicted when the proposed method is applied to the theory of clathrate hydrates. The immediate application of our work is in the development of guest−water potentials for clathrate hydrates. On the basis of the results from this study, the theory of clathrate hydrates used in such studies has to include (i) effect of movement of water molecules, (ii) guest−water interaction beyond the first shell, and (iii) guest−guest interactions. The proposed method is applicable not only to single site models of guest molecules (such as the widely used Kihara models for guest molecules) but also to more complex multisite models with or without internal degrees of freedom. Although not shown, it is quite straightforward to extend the theory to include guests that show multiple occupancy of cages, such as nitrogen and hydrogen. However, the accuracy of the theory in such cases remains to be investigated.

The occupancy for guest of type i is them obtained by differentiating Γ with respect to fî , that is,

Γ(μ1 , ..., μn , Nw , P , T ) = Δ(Nw , P , T ) n

∏ [1 + ∑ Cjαfĵ e−U̅ /k T ]M jα

α

θi =

=

jα

α

B

■

AUTHOR INFORMATION

■

ACKNOWLEDGMENTS The financial support for this work has been provided by grants from the Department of Science and Technology, Government of India.

■

REFERENCES

(1) van der Waals, J. H.; Platteeuw, J. C. Clathrate solutions. Adv. Chem. Phys. 1959, 2, 1−57. (2) John, V.; Papadopoulos, K.; Holder, G. A generalized model for predicting equilibrium conditions for gas hydrates. AIChE J. 1985, 31, 252−259. (3) Cao, Z.; Tester, J. W.; Sparks, K. A.; Trout, B. L. Molecular computations using robust hydrocarbon-water potentials for predicting gas hydrate phase equilibria. J. Phys. Chem. B 2001, 105, 10950− 10960. (4) Tee, L. S.; Gotoh, S.; Stewart, W. E. Molecular parameters for normal fluids. Lennard-Jones 12−6 Potential. Ind. Eng. Chem. Fundam. 1966, 5, 356−363. (5) Tee, L. S.; Gotoh, S.; Stewart, W. E. Molecular parameters for normal fluids. Kihara potential with spherical core. Ind. Eng. Chem. Fundam. 1966, 5, 363−367. (6) Cao, Z.; Tester, J. W.; Trout, B. L. Sensitivity analysis of hydrate thermodynamic reference properties using experimental data and ab initio methods. J. Phys. Chem. B 2002, 106, 7681−7687. (7) Chialvo, A. A.; Houssa, M.; Cummings, P. T. Molecular dynamics study of the structure and thermophysical properties of model sI clathrate hydrates. J. Phys. Chem. B 2002, 106, 442−451. (8) John, V. T.; Holder, G. D. Contribution of second and subsequent water shells to the potential energy of guest−host interactions in clathrate hydrates. J. Phys. Chem. 1982, 86, 455−459. (9) Rodger, P. M. Stability of gas hydrates. J. Phys. Chem. 1990, 94, 6080−6089. (10) Sparks, K. A.; Tester, J. W. Intermolecular potential energy of water clathrates: The inadequacy of the nearest-neighbor approximation. J. Phys. Chem. 1992, 96, 11022−11029. (11) Kvamme, B.; Lund, A.; Hertzberg, T. The influence of gas−gas interactions on the Langmuir constants for some natural gas hydrates. Fluid Phase Equilib. 1993, 90, 15−44. (12) Klauda, J. B.; Sandler, S. I. A fugacity model for gas hydrate phase equilibria. Ind. Eng. Chem. Res. 2000, 39, 3377−3386.

(16)

int

(17)

In addition, defining the Langmuir constant, Ciα as int

⟨∫ e−β(uj + uj ) dr dω⟩ Vα

int

kBT ∫ e−βuj dω

(20)

The authors declare no competing financial interest.

fĵ

Cjα =

Ciαfi ̂ exp( −Ui̅ α /kBT ) Mα M 1 + ∑j Cjαfĵ exp( −Uj̅ α /kBT )

Notes

int

kBT ∫ e−βuj dω

f[̂i] , P , T

*E-mail: [email protected].

where Δ is the isothermal-isobaric partition function of the empty clathrate hydrate. The angular brackets indicate that the integral is averaged over all the configurations of the water molecules in the hydrate lattice. For molecules with internal degrees of freedom, the activity zj can be expressed in terms of the fugacity, fĵ , as zj =

∑

(19)

Corresponding Author

e−β(uj + uj ) dr dω⟩]Mα

Vα

j=1

α

̂ denotes all fugacities except that of component i. where f[i]

Γ(μ1 , ..., μn , Nw , P , T ) = Δ(Nw , P , T ) n

⎛ ⎞ 1 ⎜ ∂ ln Γ ⎟ M ⎜⎝ ∂ ln fi ̂ ⎟⎠

α

A. LANGMUIR CONSTANT CALCULATION FOR GUEST MOLECULES WITH INTERNAL DEGREES OF FREEDOM The derivation of the theory starts from eq 7 in Ravipati and Punnathanam34 (The notations have been slightly changed to be consistent with equations presented in section 2). The semigrand isothermal-isobaric partition function, Γ, when extended to multiple guests with internal degrees of freedom is

∏ [1 + ∑ zje−U̅ /k T ⟨∫

B

j=1

(18) 12375


Article

The Journal of Physical Chemistry C (13) Klauda, J. B.; Sandler, S. I. Ab initio intermolecular potentials for gas hydrates and their predictions. J. Phys. Chem. B 2002, 106, 5722− 5732. (14) Klauda, J. B.; Sandler, S. I. Phase behavior of clathrate hydrates: A model for single and multiple gas component hydrates. Chem. Eng. Sci. 2003, 58, 27−41. (15) Natarajan, V.; Bishnoi, P. R. Langmuir constant computations for gas hydrate systems. Ind. Eng. Chem. Res. 1995, 34, 1494−1498. (16) 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−6308. (17) Anderson, B. J.; Tester, J. W.; Trout, B. L. Accurate potentials for argon−water and methane−water interactions via ab initio methods and their application to clathrate hydrates. J. Phys. Chem. B 2004, 108, 18705−18715. (18) Anderson, B. J.; Bazant, M. Z.; Tester, J. W.; Trout, B. L. Application of the cell potential method to predict phase equilibria of multicomponent gas hydrate systems. J. Phys. Chem. B 2005, 109, 8153−8163. (19) Velaga, S. C.; Anderson, B. J. Carbon dioxide hydrate phase equilibrium and cage occupancy calculations using ab initio intermolecular potentials. J. Phys. Chem. B 2014, 118, 577−589. (20) Garapati, N.; Anderson, B. J. Statistical thermodynamics model and empirical correlations for predicting mixed hydrate phase equilibria. Fluid Phase Equilib. 2014, 373, 20−28. (21) Wierzchowski, S. J.; Monson, P. A. Calculation of free energies and chemical potentials for gas hydrates using Monte Carlo simulations. J. Phys. Chem. B 2007, 111, 7274−7282. (22) Pimpalgaonkar, H.; Veesam, S. K.; Punnathanam, S. N. Theory of gas hydrates: Effect of the approximation of rigid water lattice. J. Phys. Chem. B 2011, 115, 10018−10026. (23) Tse, J. S.; Klein, M. L.; McDonald, I. R. Molecular dynamics studies of ice Ic and the structure I clathrate hydrate of methane. J. Phys. Chem. 1983, 87, 4198−4203. (24) Tse, J. S.; Klein, M. L.; McDonald, I. R. Computer simulation studies of the structure I clathrate hydrates of methane, tetrafluoromethane, cyclopropane, and ethylene oxide. J. Chem. Phys. 1984, 81, 6146−6153. (25) Tanaka, H.; Kiyohara, K. On the thermodynamic stability of clathrate hydrate. I. J. Chem. Phys. 1993, 98, 4098−4109. (26) Tanaka, H.; Kiyohara, K. The thermodynamic stability of clathrate hydrate. II. Simultaneous occupation of larger and smaller cages. J. Chem. Phys. 1993, 98, 8110−8118. (27) Tanaka, H. The thermodynamic stability of clathrate hydrate. III. Accommodation of nonspherical propane and ethane molecules. J. Chem. Phys. 1994, 101, 10833−10842. (28) Ravipati, S.; Punnathanam, S. N. Analysis of parameter values in the van der Waals and Platteeuw theory for methane hydrates using Monte Carlo molecular simulations. Ind. Eng. Chem. Res. 2012, 51, 9419−9426. (29) Frenkel, D.; Ladd, A. J. C. New Monte Carlo method to compute the free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres. J. Chem. Phys. 1984, 81, 3188−3193. (30) Polson, J. M.; Trizac, E.; Pronk, S.; Frenkel, D. Finite-size corrections to the free energies of crystalline solids. J. Chem. Phys. 2000, 112, 5339−5342. (31) Vega, C.; Sanz, E.; Abascal, J.; Noya, E. Determination of phase diagrams via computer simulation: Methodology and applications to water, electrolytes and proteins. J. Phys.: Condens. Matter. 2008, 20, No. 153101. (32) Jensen, L.; Thomsen, K.; von Solms, N.; Wierzchowski, S.; Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Calculation of liquid water-hydrate-methane vapor phase equilibria from molecular simulations. J. Phys. Chem. B 2010, 114, 5775−5782. (33) Ravipati, S.; Punnathanam, S. N. Calculation of three-phase methane-ethane binary clathrate hydrate phase equilibrium from Monte Carlo molecular simulations. Fluid Phase Equilib. 2014, 376, 193−201.

(34) Ravipati, S.; Punnathanam, S. N. Calculation of chemical potentials and occupancies in clathrate hydrates through Monte Carlo molecular simulations. J. Phys. Chem. C 2013, 117, 18549−18555. (35) Sloan, E. D., Jr; Koh, C. Clathrate Hydrates of Natural Gases; CRC Press: Boca Raton, FL, 2007. (36) In our earlier work,22 an incorrect analytical expression for the change in the fugacity upon inclusion of guest−guest interations was presented. However, the error introduced in the numerical value of the fugacity because of this expression is very small. (37) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. Sorption thermodynamics, siting and conformation of long n-alkanes in silicalite as predicted by configurational-bias Monte Carlo integration. J. Phys. Chem. 1995, 99, 2057−2079. (38) Abascal, J. L. F.; Sanz, E.; García Fernández, R.; Vega, C. A potential model for the study of ices and amorphous water: TIP4P/ Ice. J. Chem. Phys. 2005, 122, No. 234511. (39) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (40) Martin, M. G.; Siepmann, J. I. Novel configurational-bias Monte Carlo method for branched molecules. Transferable potentials for phase equilibria. 2. United-atom description of branched alkanes. J. Phys. Chem. B 1999, 103, 4508−4517. (41) Conde, M. M.; Vega, C. Determining the three-phase coexistence line in methane hydrates using computer simulations. J. Chem. Phys. 2010, 133, No. 064507. (42) Takeuchi, F.; Hiratsuka, M.; Ohmura, R.; Alavi, S.; Sum, A. K.; Yasuoka, K. Water proton configurations in structures I, II, and H clathrate hydrate unit cells. J. Chem. Phys. 2013, 138, No. 124504. (43) Smit, B.; Siepmann, J. I. Computer simulations of the energetics and siting on n-alkanes in zeolites. J. Phys. Chem. 1994, 98, 8442− 8452. (44) Adisasmito, S.; Frank, R. J., III; Sloan, E. D., Jr Hydrates of carbon dioxide and methane mixtures. J. Chem. Eng. Data 1991, 36, 68−71. (45) McLeod, H. O.; Campbell, J. M. Natural gas hydrates at pressures to 10,000 psia. J. Petrol. Technol. 1961, 222, 590−594. (46) Avlonitis, D. Multiphase equilibria in oil−water hydrate forming systems. MSc Thesis, Heriot-Watt University, Edinburgh, Scotland,1988. (47) Holder, G. D.; Godbole, S. P. Measurement and prediction of dissociation pressures of isobutane and propane hydrates below the ice point. AIChE J. 1982, 28, 930−934. (48) Deaton, W.; Frost, E. Gas Hydrates and Their Relation to the Operation of Natural-Gas Pipe Lines; Monograph; American Gas Association: Washington, DC, 1946. (49) Mei, D.-H.; Liao, J.; Yang, J.-T.; Guo, T.-M. Hydrate formation of a synthetic natural gas mixture in aqueous solutions containing electrolyte, methanol, and (electrolyte + methanol). J. Chem. Eng. Data 1998, 43, 178−182. (50) Hendriks, E. M.; Edmonds, B.; Moorwood, R. A. S.; Szczepanski, R. Hydrate structure stability in simple and mixed hydrates. Fluid Phase Equilib. 1996, 117, 193−200. (51) Subramanian, S.; Kini, R.; Dec, S.; S, E., Jr Evidence of structure {II} hydrate formation from methane + ethane mixtures. Chem. Eng. Sci. 2000, 55, 1981−1999. (52) Subramanian, S.; Ballard, A.; Kini, R.; Dec, S.; Sloan, E. Structural transitions in methane + ethane gas hydrates? Part I: Upper transition point and applications. Chem. Eng. Sci. 2000, 55, 5763− 5771. (53) Ballard, A. L., Jr.; Sloan, E. D., Jr. Structural transitions in methane + ethane gas hydrates? Part II: modeling beyond incipient conditions. Chem. Eng. Sci. 2000, 55, 5773−5782. (54) Docherty, H.; Galindo, A.; Vega, C.; Sanz, E. A potential model for methane in water describing correctly the solubility of the gas and the properties of the methane hydrate. J. Chem. Phys. 2006, 125, No. 074510. 12376


Article

The Journal of Physical Chemistry C (55) Ashbaugh, H. S.; Liu, L.; Surampudi, L. N. Optimization of linear and branched alkane interactions with water to simulate hydrophobic hydration. J. Chem. Phys. 2011, 135, No. 054510.

12377


Improving the Rigor and Consistency of the Thermodynamic Theory

Recommend Documents