THF Hydrates - American Chemical Society

Jun 18, 2008 - 15310, Greece, and School of Chemical Engineering, National Technical UniVersity of Athens,. Zografou 15780, Greece. ReceiVed: June 18 ...
3 downloads 0 Views 960KB Size
Download PDF
10294

J. Phys. Chem. C 2008, 112, 10294–10302

Evaluation of the Hydrogen-Storage Capacity of Pure H2 and Binary H2-THF Hydrates with Monte Carlo Simulations N. I. Papadimitriou,*,†,‡ I. N. Tsimpanogiannis,§,† A. Th. Papaioannou,‡ and A. K. Stubos† EnVironmental Research Laboratory, National Center for Scientific Research “Demokritos”, Agia ParaskeVi 15310, Greece, and School of Chemical Engineering, National Technical UniVersity of Athens, Zografou 15780, Greece ReceiVed: June 18, 2007; ReVised Manuscript ReceiVed: April 2, 2008

Grand Canonical Monte Carlo simulations are employed to investigate the hydrogen-storage capacity of pure H2 and binary H2-THF hydrates (both of them are of sII type) at various temperatures, pressures, and THF concentrations. It is found that the storage capacity of pure H2 hydrate could reach 3.0 wt % only at pressures above 380 MPa (at 274 K). Depending on the pressure, the large cavities of this hydrate can accommodate up to four H2 molecules, whereas the small ones are singly occupied even at pressures as high as 450 MPa. For the binary H2-THF hydrate, all large cavities are found to be occupied by a single THF molecule, independent of the THF concentration in the initial solution in a wide range of conditions, with no hydrogen molecules entering the large cavities. Likewise, in the pure H2 hydrate, small cavities are found to be singly occupied over the entire pressure range considered (1-350 MPa). The H2 storage capacity of the binary H2-THF hydrate, at temperatures close to ambient, is estimated to be lower than 1.1 wt %. 1. Introduction Gas hydrates are ice-like, crystalline, nonstoichiometric materials that belong to the class of clathrates (i.e., inclusion compounds). They are composed of a framework of hydrogenbonded water molecules that form cavities (cages) with specific geometries and sizes, inside which small molecules can be encaged (enclathrated).1 Depending on their crystal structure, hydrates are divided in three types: sI, sII, and sH. Each type of hydrate consists of two or three different types of cavities, and the number of cavities of each type per unit cell is fixed.2 Gas hydrates attribute their stability to the interactions between the lattice of water molecules and the enclathrated gas. This is not a chemical bond but a weak, Van der Waals-type, physical interaction. However, in the absence of the encaged gas, hydrates are not stable. Alternatively, they may exhibit stability under conditions where pure ice is not thermodynamically stable (e.g., above 0 °C at ambient pressure). Gas hydrates have been investigated for many years because they are involved in a variety of processes 3–5 of either industrial or scientific/academic interest. Gas hydrates started attracting industrial interest when they were identified as the cause of natural-gas pipeline blocking,6 having significant safety and economical impact. Large amounts of methane hydrates have also been discovered inside oceanic sediments and at permafrost regions, rendering them as possible future energy sources. The energy resources stored in the aforementioned hydrate deposits are estimated to be larger than the combined resources from all other fossil fuels,7 with the issue still under debate in the literature.8 Other potential applications of hydrates include CO2 sequestration in oceanic waters9 and water desalination.10 * Corresponding authors. Fax: +302106525004. E-mail: nikpap@ ipta.demokritos.gr (N.I. Papadimitriou); [email protected] (I.N. Tsimpanogiannis). † Environmental Research Laboratory, National Center for Scientific Research “Demokritos”. ‡ National Technical University of Athens. § Visiting scientist.

Hydrates have been considered as an alternative material for storing and transporting gases such as methane,11 and recently for hydrogen.12,13 Hydrogen storage is an important issue to be resolved before the establishment of the so-called “Hydrogen Economy”. Until recently, the hydrogen molecule was considered too small to stabilize the hydrate cavities and, therefore, incapable of forming hydrate by itself. This picture changed dramatically with the synthesis of pure hydrogen hydrate, initially by Dyadin et al. in 1999,14 and subsequently by Mao et al. 15 in 2002. Hydrogen hydrate was found to be of cubic sII structure, with H2 content up to approximately 5 wt %. The unit cell of the sII hydrate (cubic Fd3m space group) consists of 136 water molecules that form two types of cavities: the small, a pentagonal dodecahedron (512) and the large that is formed by 12 pentagons and 4 hexagons (51264). There are 16 small and 8 large cavities per unit cell, and it was estimated by Mao et al.15 that more than one H2 molecule can enter the same cavity (up to two molecules per small cavity and up to four per large). Patchkovski and Tse16 used Density Functional Theory (DFT) in conjunction with firstprinciples quantum chemistry calculations and confirmed the double occupancy of the small cavities and the quadruple occupancy of the large cavities. Alavi et al.17 used Molecular Dynamics simulations and reported single occupancy of the small cavities and quadruple occupancy of the large cavities. Lokshin et al.18 conducted neutron diffraction studies of the deuterium (D2) clathrate hydrate and while confirming that up to four D2 molecules enter the large cavities they reported single D2 occupancy in the small cavities. The number of hydrogen molecules inside the cavities is a crucial factor that would determine the total storage capacity of the hydrate (e.g., 5.0 wt % H2 with double occupancy of the small cavities, and 3.9 wt % H2 with single occupancy, assuming quadruply occupied large cavities). If hydrates are to be used as hydrogen-storage materials, then specific storage capacity targets should be met.19 In addition, the storing conditions should

10.1021/jp074706s CCC: $40.75  2008 American Chemical Society Published on Web 06/18/2008

Hydrogen-Storage Capacity of Pure H2 and Binary H2-THF Hydrates be moderate, especially in case hydrates are to be used in everyday applications (e.g., hydrogen storage for automotive applications). However, for the stability of pure H2 hydrates, very high pressures are required (200 MPa at 280 K).15 Alternatively, hydrogen hydrate can be stable at ambient pressure (0.1 MPa) for temperatures lower than 140 K.13,15 Both of these cases can be regarded as extreme for practical applications. Florusse et al. 20 managed to produce hydrogen hydrates at low pressures (down to 5 MPa) by adding a hydrate promoter. Hydrate promoters are substances that assist in stabilizing the hydrate structure at moderate conditions by occupying some of the available cavities, allowing H2 molecules to enter the remaining. Obviously, the presence of the promoter reduces the storage capacity of the material. Tetrahydrofuran (THF) was found to be the most effective promoter. If all eight large cavities of the sII structure are occupied by THF and assuming that two H2 molecules enter each small cavity, then the maximum hydrogen content of the hydrate would be 2.1 wt %. However, the maximum hydrogen content would be 1.05 wt % with single hydrogen occupancy of the small cavities. The occupancy of the small cavities of the H2-THF hydrate remains an open issue because there is no consensus among the experimental data.20–23 Furthermore, Lee at al.21 and Kim et al.22 suggested that by adjusting the THF concentration in the formation solution one can “tune” how many of the large cavities would be occupied by THF, with the rest of them remaining available to H2 molecules. They found a critical value for the THF concentration of about 0.15 mol %, which led to a binary H2-THF hydrate having 4.1 wt % H2, stable at 12 MPa. A similar effect was later observed in other binary hydrates.22,24 Alavi et al.25 used Molecular Dynamics to study binary H2-THF hydrates with different H2/THF ratios and reported that the substitution of H2 molecules by THF molecules in the large cavities results in lower configurational energies and enhances the stability of the hydrate. Recent experimental studies 26–29 that attempted to confirm the feasibility of the tuning process of H2-THF hydrate produced results refuting the suggestions by Lee at al.21 and Kim et al.22 The new experimental studies, based on various experimental techniques (gas release measurements,26,27 NMR spectroscopy,26 and Raman spectroscopy28,29), showed a fixed THF content in the hydrate, independent of its concentration in the initial solution. The reported hydrogen content values vary from 0.3 up to 1.05 wt %. These values are far from the requirements for the practical use of this material in mobile applications (hydrogen storage).30 The observed discrepancies highlight the necessity of clarifying the actual uptake of H2-THF hydrates. The main objective of this work is to computationally determine the storage capacity of this material at various conditions using Grand Canonical Monte Carlo (GCMC) simulations. Initially, pure H2 and pure THF hydrates are examined (both known to form cubic sII structure hydrates) and their H2 or THF content is calculated as a function of pressure. These results can be considered as the “adsorption isotherms” keeping in mind that the formation of the hydrate is simulated as an adsorption process. Then, similar simulations are performed for the binary H2-THF hydrate regarding the THF concentration as an extra parameter. Next, we determine the H2 uptake as a function of the THF concentration for comparison with the existing experimental results. Apart from the overall H2 content, the occupancy ratio (i.e., the fraction of cavities of each type occupied by a specific number of H2 molecules) is also investigated.

J. Phys. Chem. C, Vol. 112, No. 27, 2008 10295

2. Background 2.1. Classical Approach. Because of their wide range of applications, hydrates have been investigated experimentally and theoretically for many years. The main approach for the theoretical-computational description of the stability and properties of hydrates was originally developed by Van der Waals and Platteeuw31 (VDWP) and later improved by other researchers.32,33 VDWP theory is probably the most successful application of Statistical Mechanics to industrial problems. In this theory, the calculation of the stability region of the hydrate is achieved by the equation of chemical potential of water in both phases (bulk water or ice phase and hydrate phase), where the chemical potential in the hydrate phase is given by

µΗ ) µ0 + RT

∑ νj ln(1 - ∑ θij) j

(1)

i

where µH and µ0 are the chemical potential of water in the hydrate and in a hypothetical empty hydrate, respectively, νj is the number of water molecules per cavity of type j, R is the gas constant, T is the temperature, and θij is the fraction of cavities of type j occupied by a molecule of component i. The fraction of occupied cavities is a Langmuir-type function of the gas fugacity

θij )

Cij fi 1+

(2)

∑ Cij fi i

where Cij is the Langmuir constant of the adsorption of the component i in the cavity type j, and fi is the fugacity of gas component i. But the most crucial part of this theory is to develop an accurate method to estimate the Langmuir constants, which involves a reliable representation of the interactions between gas molecules and the water lattice. The traditional approach is the use of common potential functions, such as the Kihara spherical core potential.32 Although VDWP theory has been highly successful for many years, it is based on some approximations whose validity has been questioned, especially under conditions involved in novel applications of hydrates (e.g., multiple occupancy of cavities) and a number of modifications for this theory have been proposed.32–37 More recently, the increase of computational power has paved the way for the application of alternative approaches such as Molecular Dynamics (MD) and Monte Carlo (MC) simulations. Such computational methods have given extremely interesting information 17,25,38–41 that is complementary to the classical thermodynamic approach. 2.2. Monte Carlo Simulations. When applying the Grand Canonical Monte Carlo method,42 gas molecules are randomly created, destructed, or displaced using a random-number generator and for each of these events a new configuration (microstate) of the system is constructed. The total potential energy of the new configuration is calculated. This energy consists of the interaction energy between gas and solid molecules as well as between gas molecules themselves. According to the classical approach proposed by Metropolis et al.,43 the probability of accepting the new configuration depends on the total potential energy difference between the new and the old configuration and on the gas fugacity. Specifically, when a new molecule is created, this probability is given by the equation

Fmfn )

fV ∆U exp kT(Nm - 1) kT

(

)

(3)

where, Fmfn is the probability of moving from the old configuration m to the new configuration n, f is the gas fugacity, V is

10296 J. Phys. Chem. C, Vol. 112, No. 27, 2008

Papadimitriou et al.

the volume of the simulation box, k is the Boltzmann constant, T is the temperature, Nm is the number of gas molecules in the old configuration, and ∆U ) Un - Um is the total potential energy difference between the new and old configuration. Similarly, in the case of the destruction of a molecule, the acceptance probability is

Fmfn )

kTNm ∆U exp fV kT

(

)

(4)

When a molecule is displaced within the volume of the simulation box, the acceptance probability depends only on the energy difference:

(

Fmfn ) exp -

∆U kT

)

(5)

In all three cases above, when the value of the acceptance probability is less than one, a random number in the range [0-1] is generated and compared to the probability value. If this number is smaller than the probability value, then the new configuration is accepted; otherwise, it is rejected. If the probability value is more than one, then the new configuration is directly accepted. In the case of rejection of the new configuration, the old configuration is accepted as the new microstate of the system and the next trial step is carried out. After a large number of such steps, the average properties of all of the accepted configurations are calculated (e.g., average number of gas molecules in the simulation box). During a GCMC simulation, the Grand Canonical (µVT) statistical ensemble is sampled, which means that the temperature T, the volume of the simulation box V, and the fugacity f (which is an alternative formulation of chemical potential µ) are fixed, while the number of adsorbed molecules of each one of the adsorbates (if multiple species are adsorbed) can vary during the calculations. GCMC simulations are used widely in the investigation of gas adsorption in solid materials. Gas hydrates have a fixed and welldefined geometry, and this enables the simulation of their formation as a process of “adsorption” of a gas in a solid material where an “adsorption site” can be positioned anywhere within a cavity. In fact, the same concept was implicitly employed by Van der Waals and Platteeuw,31 who described the fraction of filled cavities by a Langmuir-type function. Tester et al.44 were the first to incorporate the Monte Carlo approach in the VDWP theory, in order to accurately calculate the potential within the cavities. Other researchers45,46 utilized GCMC simulations to evaluate the reliability of Lennard-Jones potential function, for several gas molecules, and calculated the Langmuir constants. Tanaka47 used the GCMC method, in the same way as it is applied for gas adsorption, to calculate the occupancy ratios of xenon, ethane, and ethylene in sI hydrates. Recently, Sizov and Piotrovskaya40 applied GCMC to methane hydrates and extended this approach to a nonrigid framework of water molecules where they found that the flexible framework predicts lower methane occupancies than the rigid framework. This behavior was also confirmed by the work of Wierzchowski and Monson48 who combined Monte Carlo simulations with several molecular models in order to calculate the free energies and chemical potentials of methane hydrates. More recently, Katsumasa et al.41 applied the GCMC method to pure H2-hydrates and considered two GCMC approaches: NVT (fixed volume) and NPT (fixed pressure). An important advantage of the GCMC approach is that some of the approximations of the VDWP theory are not a priori necessary and their validity can be evaluated from the results of the simulation. In particular, there is no need to set any limit to the number of gas molecules per cavity. This is a main concern

for the application of the VDWP theory to hydrates of smallsized gases (e.g., H2, N2, Ar) where multiple occupancy of cavities has been identified experimentally49,50 and theoretically.38,51 If multiple occupancy contributes to the stability of the hydrate structure, then such microstates are to be favored by the Monte Carlo simulations (are given higher probabilities), and the final average occupancy of each type of cavity will occur as a result of the simulation. Moreover, interactions between gas molecules, located within the same cavity or at different cavities, can be included in the calculations. This kind of interaction is more important in cases of multiple occupancy of the cavities, where gas molecules are placed close to each other. Finally, with the GCMC approach, each gas molecule interacts with all of the water molecules located within a specified distance and not only with the nearest neighbors. The energy summation is directly carried out by molecule-molecule or atom-atom energy calculations without smearing the interaction with the water molecules around a hypothetical sphere as happens in the VDWP theory. However, the Monte Carlo approach, as applied in this work, cannot provide any information about the stability of the hydrate at the specific conditions. So, the phase diagram cannot be derived from such calculations and the stability of the hydrate must be explicitly assumed before the simulation. 3. Simulation Details 3.1. Water, H2, and THF Models. Simulations are carried out in one unit cell of sII hydrate containing 136 fixed water molecules placed on a rigid lattice. However, in order to ensure that simulations on just one unit cell give sufficiently accurate results, some of the GCMC runs were repeated on a simulation box containing eight unit cells (two along each axis), with the remaining parameters unchanged. The results justify the use of one unit cell as explained further in Section 4. In both cases, three-dimensional periodic boundary conditions are applied. The lattice constant of the unit cell is taken as equal to 17.047 Å.15 The atomic positions of the oxygen atoms are those that have been recently fitted to X-ray diffraction data by Yousuf et al.52 The oxygen atoms of the water molecules are placed exactly on these positions while hydrogen atoms are placed on the line that connects the nearby positions at a distance of 1.000 Å from the oxygen atom. So, the H-O-H angle varies between 105 and 120° depending on the kind of polyhedron the molecule lies on. Among the proton configurations tested, all satisfying the Bernal-Fowler “ice rules”,53 the one with the smallest total dipole moment was selected to be used in the simulations. Water molecules are simulated by the Extended Simple Point Charge (SPC/E) model.54 SPC/E was selected because it gives sufficient correlation of liquid water properties at an acceptable experimental cost38 and has been used widely for molecular simulations on hydrates17,25,38–40,48 with significant success. The efficiency of several water models in predicting hydrate stability through MD simulations is discussed in detail by Frankcombe and Kroes.55 In the SPC/E model, Van der Waals interactions are described with a Lennard-Jones interaction site placed on the oxygen atom. For the electrostatic interactions, a positive partial charge is assigned to each hydrogen atom, and the double negative charge, to the oxygen atom. All parameters connected to the SPC/E model are shown in Table 1. The length of the H-H bond in the H2 molecule is taken as 0.7414 Å.56 The Lennard-Jones interaction site resides in the middle of the H-H bond. Partial charges are assigned to each H atom and to the middle of the H-H bond. Parameters for the H2 molecules are also shown in Table 1. This set of

Hydrogen-Storage Capacity of Pure H2 and Binary H2-THF Hydrates TABLE 1: Interaction Parameters and Partial Charges for Each Molecule molecule

atom

σ (Å)

 (kJ/mol)

charge (e)

H2O

O H H center of mass C O H1a H2a

3.166 0.000 0.000 3.038 3.400 3.000 2.471 2.645

0.6502 0.0000 0.0000 0.2852 0.4577 0.7113 0.0657 0.0657

-0.8476 +0.4238 +0.4932 -0.9864 see Figure 1 see Figure 1 see Figure 1 see Figure 1

H2 THF

a Symbols H1 and H2 refer to the corresponding H atoms in Figure 1.

Figure 1. Partial charges of the optimized THF molecule.

σi + σj and ij ) √i · j 2

isotherms”, we used as the initial configuration for each run the final configuration of the previous run with a lower fugacity value. For the calculation of the occupancy of each type of cavity, a sample configuration every 25 steps is taken and all of the H2 or THF molecules are assigned to their closest center of cavity. At the end of the simulation, the average occupancy of every cavity is calculated and, finally, occupancies are averaged separately for the small and large cavities because it was observed that no cavity presented a different occupancy trend relatively to the other cavities of the same type. 3.3. Fugacity Calculations. The fugacity f of the gas phase is calculated with the Soave-Redlich-Kwong60 equation of state, which was found to give excellent correlation of PVT experimental data of the H2 gas phase61,62 at both lower (77 K) and higher (250-300 K) temperatures. The gas phase is regarded as pure H2 (minor amounts of water or THF vapor are ignored). In case of the binary hydrate, two components are enclathrated (i.e., there are two “adsorbents”): a gas component (H2) and a liquid component (THF) that is dissolved in the aqueous phase. In the GCMC simulations, this binary system should be represented as a mixture of two gases. This approach requires an accurate estimation of the fugacity of THF in the aqueous solution. This is a significant part of the calculation because water and THF are completely miscible and form a highly nonideal solution. The fugacity of the dissolved THF is calculated using the concept of activity coefficient (γ-φ approach):60

( RT1 ∫ V (p)dp)

f ) γ · x · φs · Ps · exp

parameters has been derived from experimental data of the solid and gas phase of H257 and has been used in MD17,25,58 and GCMC41 simulations of hydrogen hydrates. Finally, regarding the THF molecule, a geometry optimization is carried out prior to the GCMC simulation. For this optimization, the AMBER59 forcefield is utilized in the same way as in previous works.25,55 The structural details of the optimized THF geometry are given by Cornell et al.,59 while the interaction parameters and partial charges are presented in Table 1 and Figure 1, respectively. All molecules were considered to be rigid during the simulations. Interaction parameters between different types of atoms are calculated using the Lorentz-Berthelot combining rules60

σij )

J. Phys. Chem. C, Vol. 112, No. 27, 2008 10297

P

Ps

L

(7)

where γ is the activity coefficient, x is the mole fraction of THF in the aqueous solution, Ps is the vapor pressure of pure THF at temperature T, φs is the fugacity coefficient of pure saturated gas THF at temperature T, and VL is the molar volume of liquid THF as a function of pressure. The exponential term is the Poynting correction factor that accounts for the effect of high pressure on the fugacity of liquids. The activity coefficient is calculated from the NRTL equation,60 with parameters derived from VLE equilibrium data and recommended by Gmehling et al.63 For the analytical calculation of the integral in eq 7, THF density data of Back and Woolf64 were fitted using the modified Tait equation that accurately describes the compressibility of liquids at high pressures.65 It has been assumed that THF concentration in the solution remains constant and equal to the initial concentration during any simulation.

(6)

where σ is the distance for zero potential and  is the energy well depth. Potential energy calculation is performed through an atom-atom summation (excluding atoms in the same molecule). A cutoff distance of 20 Å is applied to both Van der Waals and electrostatic interactions. The Ewald summation technique42 is also used for the long-range electrostatic interactions. 3.2. GCMC Runs. Each GCMC run includes 10 000 000 steps, randomly chosen, with equal probabilities for each type of “move”: creation 33.3%, destruction 33.3%, and displacement 33.3%. The maximum displacement distance is set to 1 Å in the case of pure H2 hydrate and to 0.5 Å for both H2 and THF in the case of binary hydrate. These values yield reasonable acceptance ratios42 (40-60% on average). In the case of binary H2-THF hydrate, the two components (H2 and THF) are sampled with equal probabilities. In order to obtain the “adsorption

4. Results and Discussion All of the simulation results reported below have been derived at temperature and pressure conditions where the hydrate has been experimentally confirmed to be stable.13,15,66 In all Figures presented in Sections 4.1-4.4, the empty symbols represent results from runs on a simulation box with eight unit cells. The remaining results (denoted with solid symbols) are from runs on a simulation box with one unit cell. It is apparent that both sizes of simulation box give almost identical results. This fact justifies the selection of only one unit cell for the simulations considering that runs on eight unit cells require significantly more computational power and storage capacity. The large number of steps (107) in each run could partly compensate for the small size of the simulation box and thus yield statistically accurate results. Please also recall that the sII hydrate unit cell has 16 small cages and 8 large ones. A

10298 J. Phys. Chem. C, Vol. 112, No. 27, 2008

Figure 2. Hydrogen content as a function of pressure for H2 hydrate at 77 K. Empty circles denote results from a simulation box with eight unit cells.

Figure 3. Occupancy ratio of large cavities (fraction of cavities occupied by 1, 2, 3, or 4 H2 molecules) as a function of pressure for H2 hydrate at 77 K.

similar effect of the size of the simulation box has been observed in GCMC simulations on argon hydrates.67 4.1. Pure H2 Hydrate at Low Temperature (77 K). The first set of simulations concerns pure H2 hydrate at cryogenic temperatures (77 K). The experimental data of Mao et al.15 show remarkable stability at pressures down to 10 kPa (the quenched hydrate is reported to be stable up to temperature 145 K at ambient pressure). For such low temperatures, the H2 uptake in the hydrate as a function of pressure is shown in Figure 2. In the pressure range examined (0.1-20 MPa), H2 content varies from 1.8 to 2.7 wt %. A more detailed investigation in the occupancy of each type of cavity reveals that small cavities are singly occupied with no occurrences of double occupancy. The fraction of filled small cavities is large even at low pressures. It has a value of 95% at ambient pressure and reaches 99.9% at 5 MPa. The absence of doubly occupied or empty small cavities is in absolute agreement with the experimental data of Lokshin et al.18 However, there is some apparent discrepancy with the same set of experimental data concerning the occupancy of large cavities at low pressures. Lokshin et al. observed an occupancy value of 4.0, almost constant with pressure at 77 K. Our GCMC simulations predict a strong dependence on pressure, as illustrated in Figure 3, and at pressures below 50 MPa, the majority of large cavities seem to be doubly occupied. However, it should be noted that the measurements of Lokshin et al. at ambient pressure were performed after cooling a hydrate that had initially been formed at 220 MPa and 200-270 K.

Papadimitriou et al.

Figure 4. Hydrogen content as a function of pressure for H2 hydrate at 250, 274, and 300 K. Empty symbols denote results from a simulation box with eight unit cells. Solid lines describe the linear regression of the results.

Formation of the hydrate at ambient pressure could possibly lead to reduced occupancy of the cavities. Our results show that quadruply occupied cavities become dominant above 70 MPa and occupancy increases with pressure until it reaches the value of 4.0 at about 150 MPa, consistent with measurements reported by Lokshin et al.18 A further increase in pressure does not affect the occupancy. 4.2. Pure H2 Hydrate at Higher Temperatures (250-300 K). Next, we consider temperatures close to ambient. The pressure dependence of the H2 uptake (in the pressure range 150-450 MPa) at temperatures of 250, 274, and 300 K is shown in Figure 4. H2 content in this pressure range is less than 3.2 wt %. This value is significantly lower than the initial expectations (up to 5.3 wt %),13,15,16 but it still renders H2 hydrate a potential hydrogen-storage material. It seems that there is linear dependence of H2 uptake on pressure, with an average slope of 3.0 × 10-5 MPa-1 (mass of H2 per unit mass of hydrate), which means that an increase in pressure of about 330 MPa is required in order to increase the H2 content by 1.0 wt %. Consequently, raising the pressure to extremely high values does not seem to be an effective way to improve the storage capacity of pure H2 hydrates. As it happens in the case of 77 K and contrary to what has been initially suggested by Mao et al.,15 the small cavities cannot easily accommodate more than one H2 molecule as can bee seen in Figure 5. Doubly occupied cavities are very scarce. Their fraction is less than 0.1% even at pressures as high as 450 MPa. Other recent experimental18 and theoretical17 works confirm the inability of the small cavities to accommodate two H2 molecules. Similar behavior of the filling of the small cavities was predicted by Katsumasa et al.41 Concerning the large cavities (Figure 6), a significant dependence on pressure is again obvious. At pressures up to 340 MPa, the majority of large cavities are doubly occupied. At higher pressures, triple occupancy becomes dominant. Quadruply occupied cavities are present, but their fraction does not exceed 15%. These results are in very good agreement with those of Katsumasa et al.41 who have performed GCMC simulations using the TIP4P model for the water molecules instead of the SPC/E used in this work (results from their NPT ensemble are also shown in Figures 5 and 6). Occupancy values higher than four are virtually not observed. The fraction of cavities occupied by five molecules remained less than 0.05%. Similar behavior is observed at 77 K. This could possibly mean that the large cavity cannot accommodate more than four H2

Hydrogen-Storage Capacity of Pure H2 and Binary H2-THF Hydrates

Figure 5. Fraction of filled small cavities as a function of pressure for H2 hydrate at 274 K. No doubly occupied cavities are observed. Empty symbols denote results from a simulation box with eight unit cells, and the dotted line denotes the NPT results of Katsumasa et al.41

Figure 6. Occupancy ratio of large cavities (fraction of cavities occupied by 1, 2, 3, or 4 H2 molecules) as a function of pressure for H2 hydrate at 274 K. Empty symbols denote results from a simulation box with eight unit cells, and dotted lines denote the NPT results of Katsumasa et al.41

molecules unless extremely high pressures are applied. The insertion of a fifth H2 molecule is, in energy terms, completely unfavorable, and the acceptance probability of such a microstate during a GCMC run is very low. Figure 7a depicts the location of an H2 molecule in a singly occupied small cavity, and Figure 7b presents the location of a cluster of four H2 molecules in a large cavity of a pure H2 hydrate. These images have been taken from a random configuration of the crystal during a GCMC run at 274 K and 200 MPa. The cluster of four H2 molecules seems to present tetrahedral configuration; however, this is just one snapshot over the 107 examined. A more descriptive insight into the distribution of H2 within the hydrate crystal is offered by the snapshots in Figure 8. In particular, Figure 8a shows the probability of finding an H2 molecule inside the cavities of a unit cell of the sII hydrate. Note that red areas correspond to higher probabilities. Figure 8b presents a section of the unit cell that includes the centers of four small and four large cavities. The small cavities present maximum H2 density at their centers because they are singly occupied by H2 molecules, whereas in the large cavities H2 is accumulated in a spherical shell of a radius that is approximately 1/2 of the cavity radius. 4.3. Pure THF Hydrate. We performed a series of simulations that involved pure THF hydrates. Under all of the

J. Phys. Chem. C, Vol. 112, No. 27, 2008 10299

conditions examined (temperatures in the range of 250-300 K, pressures in the range of 0.1-450 MPa, and THF concentrations in the aqueous solution in the range of 0.10-10 mol %), THF occupied all of the large cavities forming a hydrate with a constant THF content of 19.1 wt %. No instance of THF molecules entering the small cavities was observed. Figure 9 illustrates a possible position of a THF molecule in the large cavity of a pure THF hydrate at ambient conditions. 4.4. Binary H2-THF Hydrate. Next, we proceed to the most interesting part of this work, which is the investigation of binary H2-THF hydrates, with either stoichiometric composition (THF · 17H2O) or not. Initially, the stoichiometric H2-THF hydrate is examined. Figure 10 shows the hydrogen content as a function of pressure (up to 200 MPa) at temperatures of 250, 274, and 300 K. As can be concluded from Figure 10, the binary H2-THF hydrate has a rather low H2 storage capacity. Even at pressures as high as 350 MPa (not shown in Figure 10), the H2 uptake is less than 1.1 wt %. In the range below 200 MPa, the H2 content curves look like Langmuir-type adsorption isotherms. The upper limit of 1.05 wt %, shown in the figure, corresponds to single occupancy of all 16 small cavities of the unit cell of the sII hydrate. The fraction of filled small cavities as a function of pressure at 274 K is shown in Figure 11. It can be concluded from Figures 10 and 11 that small cavities are again singly occupied at pressures below 200 MPa, whereas doubly occupied cavities are very rare (less than 0.1%). This finding is consistent with recent experimental data23,28,29 on binary H2-THF hydrates. Therefore, for the binary H2-THF hydrate, and under the temperature, pressure, and THF concentration range considered in this study, the VDWP-theory assumption of single occupancy of all cavities is valid, and the findings of the GCMC simulations can be used to estimate the Langmuir constant for the enclathration of H2 in the small cavities. In this particular case of a binary H2-THF hydrate, eq 2 is written

θHS 2 )

CHS 2 fH2 S 1 + CHS 2 fH2 + CTHF fTHF

(8)

where the subscripts denote the component and the superscripts denote the type of cavity (S ) small, L ) large). After S introducing CTHF ) 0, because THF does not enter the small cavities, eq 8 turns into a simple Langmuir-type equation and the Langmuir constant of H2 for the small cavities can be calculated by regression of the occupancy data derived from the GCMC simulations. We performed this analysis for all of the available “adsorption isotherms” (namely, at 250, 260, 274, 280, 290, and 300 K) to obtain the Langmuir constant as a function of temperature:

CHS 2 )

806.6 exp , ( 0.4519 T ) ( T )

T in K,

CHS 2 in MPa-1

(9)

Although these values are consistent with other estimations of the Langmuir constant of H2,26 their reliability should be considered with caution because they are directly related to the forcefield used for the simulations. Generally, Langmuir constants for hydrate equilibrium are very sensitive to the parameters of the forcefield and an accurate estimation of these parameters is always a crucial issue.68 Figure 12 illustrates the comparison between the GCMCcalculated H2 uptake, at various conditions, with several sets of reported experimental data. The “adsorption isotherm” (at 274 K) calculated in this work seems to sufficiently describe the trend of the experimental data of Strobel et al.26 and Anderson et al.27 Comparatively to the values derived from gas

10300 J. Phys. Chem. C, Vol. 112, No. 27, 2008

Papadimitriou et al.

Figure 7. Possible H2 molecule locations inside the two types of cavities of pure H2 hydrate (at 274 K and 200 MPa): (a) single H2 molecule inside a small cavity, and (b) four H2 molecules inside a large cavity.

Figure 8. Visualization of a GCMC run of a H2 hydrate at 274 K and 200 MPa: (a) Probability distribution of finding a H2 molecule at any point inside the unit cell, and (b) probability distribution of finding an H2 molecule on a section that includes the centers of four small (right) and four large (left) cavities.

evolution measurements, our GCMC results present lower H2 uptake. However, there is complete agreement with the values derived from NMR measurements (Figure 12). If one considers the expected experimental error of such kinds of data and the fact that there are discrepancies among the experimental data themselves depending on the analytical method used (as discussed by the authors themselves26), then it can be concluded that the results of the GCMC simulations represent the behavior of the experimental measurements adequately. The H2 uptake as a function of the THF concentration (including values lower or higher than the stoichiometric value of 5.56 mol %) in the initial solution (at 274 K and 12 MPa) is depicted in Figure 13. The GCMC simulations show that the H2 uptake is completely independent of the THF concentration and has an almost constant value of approximately 0.30 wt %. This result is in contrast to what has been claimed by Lee et al.21 and Kim et al.22 regarding the possibility of tuning the

hydrogen uptake. However, the results of our simulations seem to correlate very reasonably with the measurements of Strobel et al.,26 with our prediction of H2 uptake (0.30 wt %) having a lower value than the experimental (0.42 wt %). By examining the occupancy of cages as derived from the GCMC runs, we observe that over this range of THF concentration (0.10 - 10 mol %) all eight large cavities of the unit cell are occupied by THF molecules and the probability of H2 molecules entering the large cavities is very limited. In fact, the occupancy of these cavities by THF molecules is 98.5% on average. The same behavior was observed at other pressures (5 and 50 MPa). Complete occupancy of the large cavities by THF molecules has also been confirmed by Raman spectroscopy measurements.28,29 In conclusion, the present results offer evidence from simulations that support the statement that the H2 content of a binary H2THF hydrate cannot be tuned by modifying the THF concentration in the initial solution.

Hydrogen-Storage Capacity of Pure H2 and Binary H2-THF Hydrates

J. Phys. Chem. C, Vol. 112, No. 27, 2008 10301

Figure 11. Fraction of filled small cavities as a function of pressure for the stoichiometric (5.56 mol % THF) binary H2-THF hydrate at 274 K. Empty circles denote results from a simulation box with eight unit cells. No doubly occupied cavities are observed. Figure 9. Possible THF molecule location inside the large cavity of pure THF hydrate (at ambient conditions) taken from the lowest energy configuration.

Figure 10. Hydrogen content as a function of pressure for the stoichiometric (5.56 mol % THF) binary H2-THF hydrate at 250, 274, and 300 K. Empty symbols denote results from a simulation box with 8 unit cells. Solid lines have been calculated using values of occupancy of the small cavities derived by the fitted Langmuir (eq 8).

Figure 12. Hydrogen content of H2-THF hydrate as a function of pressure, at 274 K. Comparison between the GCMC results of this work (solid line) and several sets of experimental data: gas release measurements of Anderson et al.27 (crosses); gas release measurements of Strobel et al.26 on 45 µm (squares) or 250 µm (circles) crushed THF hydrate; and NMR measurements (triangles) of Strobel et al.26

5. Conclusions We have performed a series of Grand Canonical Monte Carlo simulations to investigate the hydrogen-storage capacity of pure H2 hydrates and binary H2-THF hydrates. The dependence of H2 uptake on the promoter (THF) concentration in the reacting solution has been a controversial issue with some papers reporting hydrogen-content dependence on THF concentration,21,22 with others reporting a fixed hydrogen content that is independent of THF concentration.26,27 In this work, GCMC simulations have shown that THF molecules occupy all of the large cavities of the hydrate crystal, regardless of its concentration in the initial solution. Moreover, the small cavities cannot accommodate more than one H2 molecule, resulting in an H2 uptake less than 1.1 wt % (at 250-300 K and up to 350 MPa) confirming the recent experimental data.26,27 For pure H2 hydrates at high pressures (above 200 MPa at 274 K), the small cavities are occupied by only one H2 molecule, whereas the large ones can encage from 1 to 4 molecules depending on pressure. The H2 content of a pure H2 hydrate at 274 K and 200 MPa is found to be 2.5 wt %.

Figure 13. Hydrogen content of the binary H2-THF hydrate at 274 K and 12 MPa as a function of THF concentration in the equilibrium solution and comparison with experimental data. Empty circles denote results from a simulation box with eight unit cells.

Acknowledgment. Partial funding by the European Commission DG Research (contract SES6-2006-518271/ NESSHY) is gratefully acknowledged by the authors. The first of the authors wishes to acknowledge the Hellenic State Scholarships Foundation (I.K.Y.) for financial support.

10302 J. Phys. Chem. C, Vol. 112, No. 27, 2008 References and Notes (1) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; Marcel Dekker: New York, 1998. (2) Sloan, E. D. Energy Fuels 1998, 12, 191. (3) Koh, C. A. Chem. Soc. ReV. 2002, 31, 157. (4) Sloan, E. D. Nature 2003, 426, 353. (5) Chatti, I.; Delahaye, A.; Fournaison, L.; Petitet, J. P. Energy ConV. Man. 2005, 46, 1333. (6) Hammerschmidt, E. G. Ind. Eng. Chem. 1934, 26, 851. (7) Kvenvolden, K. A. Org. Geochem. 1995, 23, 997. (8) Milkov, A. V. Earth Sci. ReV. 2004, 66, 183. (9) Brewer, P. G.; Friederich, G.; Peltzer, E. T.; Orr, F. M., Jr Science 1999, 284, 943. (10) Javanmardi, J.; Moshfegian, M. Appl. Therm. Eng. 2003, 23, 845. (11) Thomas, S.; Dawe, R. A. Energy 2003, 28, 1461. (12) Hu, Y. H.; Ruckenstein, E. Angew. Chem., Int. Ed. 2006, 45, 2011. (13) Mao, W. L.; Mao, H. K. Proc. Natl. Acad. Sci. USA 2004, 101, 708. (14) Dyadin, Y. A.; Larionov, E. G.; Manakov, A. Y.; Zhurko, F. V.; Aladko, E. Y.; Mikina, T. V.; Komarov, V. Y. MendeleeV Commun. 1999, 9, 209. (15) Mao, W. L.; Mao, H. K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Science 2002, 297, 2247. (16) Patchkovskii, S.; Tse, J. S. Proc. Natl. Acad. Sci. USA 2003, 100, 14645. (17) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2005, 123, 024507. (18) Lokshin, K. A.; Zhao, Y.; He, D.; Mao, W. L.; Mao, H. K.; Hemley, R. J.; Lobanov, M. V.; Greenblatt, M. Phys. ReV. Lett. 2004, 93, 125503. (19) Schlapbach, L.; Zuttel, A. Nature 2001, 414, 353. (20) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Science 2004, 306, 469. (21) Lee, H.; Lee, J.; Kim, D. Y.; Park, J.; Seo, Y. T.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A. Nature 2005, 434, 743. (22) Kim, D. Y.; Park, Y.; Lee, H. Catal. Today 2007, 120, 257. (23) Hester, K. C.; Strobel, T. A.; Sloan, E. D.; Koh, C. A.; Huq, A.; Schultz, A. J. J. Phys. Chem. B 2006, 110, 14024. (24) Kim, D. Y.; Park, J.; Lee, J.; Ripmeester, J. A.; Lee, H. J. Am. Chem. Soc. 2006, 128, 15360. (25) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 014704. (26) Strobel, T. A.; Taylor, C. J.; Hester, K. C.; Dec, S. F.; Koh, C. A.; Miller, K. T.; Sloan, Jr., E.D. J. Phys. Chem. B 2006, 110, 17121. (27) Anderson, R.; Chapoy, A.; Tohidi, B. Langmuir 2007, 23, 3440. (28) Hashimoto, S.; Murayama, T.; Sugahara, T.; Sato, H.; Ohgaki, K. Chem. Eng. Sci. 2006, 61, 7884. (29) Hashimoto, S.; Sugahara, T.; Sato, H.; Ohgaki, K. J. Chem. Eng. Data 2007, 52, 517. (30) Scuth, F. Nature 2005, 434, 712. (31) Van der Waals, J. H.; Platteeuw, J. C. AdV. Chem. Phys. 1959, 2, 1. (32) Parrish, W. R.; Prausnitz, J. M. Ind. Eng. Chem. Proc. Des. DeV 1972, 11, 26. (33) Holder, G. D.; Corbin, G.; Papadopoulos, K. D. Ind. Eng. Chem. Fund. 1980, 19, 282. (34) John, V. T.; Holder, G. D. J. Phys. Chem. 1981, 85, 1811. (35) Klauda, J. B.; Sandler, S. I. Ind. Eng. Chem. Res. 2000, 39, 3377. (36) Sparks, K. A.; Tester, J. W. J. Phys. Chem. 1992, 96, 11022.

Papadimitriou et al. (37) Anderson, B. J.; Tester, J. W.; Trout, B. L. J. Phys. Chem. B 2004, 108, 18705. (38) Van Klaveren, E. P.; Michels, J. P. J.; Schouten, J. A.; Klug, D. D.; Tse, J. S. J. Chem. Phys. 2001, 114, 5745. (39) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 125, 104501. (40) Sizov, V. V.; Piotrovskaya, E. M. J. Phys. Chem. B 2007, 111, 2886. (41) Katsumasa, K.; Koga, K.; Tanaka, H. J. Chem. Phys. 2007, 127, 044509. (42) Allen, M. P. ; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (43) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (44) Tester, J. W.; Bivins, R. L.; Herick, C. C. AIChE J. 1972, 18, 1220. (45) Natarajan, V.; Bishnoi, P. R. Ind. Eng. Chem. Res. 1995, 34, 1494. (46) Sparks, K. A.; Tester, J. W.; Cao, Z.; Trout, B. L. J. Phys. Chem. B 1999, 103, 6300. (47) Tanaka, H. Fluid Phase Equilib. 1998, 144, 361. (48) Wierzchowski, S. J.; Monson, P. A. J. Phys. Chem. B 2007, 111, 7274. (49) Sasaki, S.; Hori, S.; Kume, T.; Smimizu, H. J. Chem. Phys. 2003, 118, 7892. (50) Kuhs, W. F.; Chazallon, B.; Radaelli, P. G.; Pauer, F. J. Incl. Phenom. Molec. Recognit. Chem. 1997, 29, 65. (51) Itoh, H.; Tse, J. S.; Kawamura, K. J. Chem. Phys. 2001, 115, 9414. (52) Yousuf, M.; Qadri, S. B.; Knies, D. L.; Grabowski, K. S.; Coffin, R. B.; Pohlman, J. W. Appl. Phys. A: Mater. Sci. Process. 2004, 78, 925. (53) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (54) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (55) Frankcombe, T. J.; Kroes, G. J. J. Phys. Chem. C. 2007, 111, 13044. (56) Molecular Structure and Spectroscopy. In CRC Handbook of Chemistry and Physics, 88th ed. (Internet Version 2008); Lide, D. R., Ed.; CRC Press/Taylor and Francis, Boca Raton, FL. http://www.hbcpnetbase.com. (57) Silvera, I. F.; Goldman, V. V. J. Chem. Phys. 1978, 69, 4209. (58) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 204707. (59) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (60) Prausnitz, J. M. ; Lichtenthaler, R. N. ; Avezedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria; Prentice-Hall International Series: New Jersey, 1998. (61) http://webbook.nist.gov. (62) Michels, A.; de Graaff, W.; Wassenaar, T.; Levelt, J. M. H.; Louwerse, P. Physica 1959, 25, 25. (63) Gmehling, J. ; Onken, U. ; Arlt, W. Vapor-Liquid Equilibrium Data Collection; DECHEMA: Frankfurt, Germany, 1978; Vol. 1, part 1a. (64) Back, P. J.; Woolf, L. A. J. Chem. Thermodyn. 1998, 30, 353. (65) Dymond, J. H.; Malhotra, R. Int. J. Thermophys. 1988, 9, 941. (66) Lokshin, K. A.; Zhao, Y. Appl. Phys. Lett. 2006, 88, 131909. (67) Papadimitriou, N. I.; Tsimpanogiannis, I. N.; Papaioannou, A. Th.; Stubos, A. K. Mol. Sim. 2008, In press, DOI: 10.1080/08927020802101734. (68) Papadimitriou, N. I.; Tsimpanogiannis, I. N.; Yiotis, A. G.; Steriotis, T. A.; Stubos, A. K. In Physics and Chemistry of Ice; Kuhs, W. F., Ed.; RSCPublishing: Cambridge, U.K., 2007; pp 475-482.

JP074706S