Abstract
A new method implementing molecular dynamics (MD) simulations for calculating the reference properties of simple gas hydrates has been proposed. The guest molecules affect interaction between adjacent water molecules distorting the hydrate lattice, which requires diverse values of reference properties for different gas hydrates. We performed simulations to validate the experimental data for determining , the chemical potential difference between water and theoretical empty cavity at the reference state, for structure II type gas hydrates. Simulations have also been used to observe the variation of the hydrate unit cell volume with temperature. All simulations were performed using TIP4P water molecules at the reference temperature and pressure conditions. The values were close to the experimental values obtained by the Lee-Holder model, considering lattice distortion.
1. Introduction
Gas hydrates are crystalline solids formed when water forms a complex lattice with gases which occupy the interstices of the hydrogen bonded water molecules [1–4]. These interstices are referred as cages and when empty, they are highly unstable, collapsing into ice structure [4]. Most thermodynamic models consider the theoretical empty cavity as the reference state that is highly unstable and is less likely to be determined experimentally. Computer simulations on the other hand are performed at a molecular level and for very short time intervals close to equilibrium. Simulations play a vital role in relating microscopic details of a system to macroscopic properties. Besides, simulation methods like molecular dynamics (MD) provide averages of the properties at stable states which can be used to determine thermodynamic properties. Theoretically, determining the size, structure, and stability of gas hydrate formation is of importance and computer simulations play a valuable role in building an environment for performing accurate experiments. Simulations can even be applied for the selection of new guest candidates for which limited experimental data are available [2], for example, methane (which normally has structure I) forming structure II in the mixture gas hydrate.
The first equilibrium model for gas hydrate was developed by van der Waals and Platteeuw [5] based on statistical thermodynamics and was generalized by Parrish and Prausnitz [3] to form the basis of all the thermodynamic models even today. These models were further extended by Holder et al. [4] by considering the energy changes due to the restriction of guest molecule movement in the hydrate lattice. These models are the basis of all the gas hydrate equilibrium calculations. Later, Lee and Holder [6] proposed that different guests require different values of reference chemical potential based on the assumption that the lattice can be distorted according to the size of gas molecules in it. They found that a cavity containing a different gas molecule would have a different size to minimize the total energy of the system [7]. This new model was modified to mixture gas hydrates by Lee and Holder [8] and Martin and Peters [1]. We are carrying out molecular-dynamics (MD) simulations using different guest molecules and the NVT ensemble. The purpose is to validate the data on the reference chemical potential difference and the effect of temperature on the size of the unit lattice structure. To apply this ensemble to a lattice distortion assumption, the lattice size was changed simultaneously at a constant temperature, and the total energy with the pressure was calculated at each condition. The MD calculations will be used to determine the reference chemical potential difference and the hydrate unit cell size.
2. Theoretical Background
The resemblance of Langmuir’s theory of gas adsorption with the formation of gas hydrates was first proposed by van der Waals and Platteeuw [4, 5]. The Langmuir constant for both is a function of temperature and the participating components [4, 9]. For the gas hydrate, it is given by where is a cavity in the hydrate lattice occupied by type of guest molecule. Assuming single guest occupancy for every cavity, the fraction of the occupied cavities is given by [9]
The unknown quantity of interest is the Langmuir adsorption coefficient which, assuming spherical symmetry, is given by [4, 10] where is the cell potential of the spherical cavity and is calculated using the Kihara potential model (refer to [3, 4] for the equation form). The classical method was adjusting the Kihara radius (), size (), and energy () parameters for the gas hydrate equilibrium so that the experimental three-phase pressure was in unison with the calculated value [4, 9]. The whole process of determining the Kihara potential parameters depends on the ability to determine the difference of chemical potential, enthalpy, and the volume of water cavity in the empty and occupied hydrate lattice. As a result, the Kihara parameters obtained from experiments are totally different from the Kihara parameters for other systems [7].
At equilibrium, the chemical potential of water in hydrate phase, , is the same as the chemical potential of water in liquid state, [4]. Thus, if , the chemical potential of empty hydrate lattice, is considered as a reference state, the difference of chemical potentials is given by [4–6, 8]
The original model did not account for lattice distortion or the misalignment of the hydrate lattice due to guest molecules [3–5]. Lee and Holder [6–8] were the first to propose this change of geometry of hydrate unit cell. The model accounted for the effect of temperature, pressure, and composition for calculating the chemical potential difference of water
The above equation (4) suggests that at ice point ( K) and 0 pressure, the reference chemical potential . The Langmuir constant in (2) was calculated using the Kihara parameters without any adjustment. In the present study, the value of was obtained from the simulation that is the same as at equilibrium. To calculate these values for different gas hydrates, we calculated the chemical potential for the theoretical empty cavity using MD simulations and also calculated the chemical potential of each gas hydrate using MD simulations. Since has been defined as the difference between these two values at 273.15 K and 0 pressure, we calculated the difference between two values using MD simulation.
Since the primary effect of interaction between adjacent guest molecules is stretching the hydrogen bonds between water molecules [7], hydrate equilibria become more the subject of the interaction between hydrate lattice and the guest molecules. The present study follows the Lee and Holder model [6] which considers different values of the reference chemical potential difference for each gas hydrate.
3. Molecular Dynamics (MD) Simulation
Molecular dynamics describe the application of the Newton’s equations of motion for a set of molecules and consist of integration of forces for short time intervals close to equilibrium. Particle trajectories are generated after several hundred steps from which time-averaged macroscopic properties such as viscosity, thermal conductivity, and diffusivity can be calculated [11]. This technique is a powerful tool to investigate microscopic phenomena and is widely used for simulating water structures. For the MD simulation to generate accurate equations of motion, forces and velocities of the particles have to be within an acceptable range.
To perform the simulations, we used the software MOLDY [12]. It uses the “link cell” to calculate short range forces and Ewald sum technique to calculate long range electrostatic forces [13]. The Gaussian thermostat was used to maintain a constant temperature of 273.15 K and other temperatures. An initial configuration of the gas hydrate was generated from the Jorgensen’s TIP4P model using a method called “skew start” incorporated in the software [14]. The molecules were randomly ordered but guaranteed minimum separation avoiding molecular overlap and provided a rough estimation of the unit cell. It also provided information on the Euler angles in the form of quaternions which were used to develop the empty hydrate lattice. Quaternions lead to equations of motion which are free of singularities resulting in a better numerical stability of the simulation. The highly unsteady empty hydrate lattice was stable for only an instance of almost 0.005 ps. Simulations with various boundary conditions were carried out to obtain the near stable configuration of the empty hydrate associated with quaternions. The instantaneous pressure and total energy were used to observe equilibrium conditions between 0.06 to 0.065 ps or 0.065 to 0.07 ps before the structure collapsed.
The gas hydrate was then simulated to a stable condition with the guest molecules in the empty cavities using the data from the initial stable configuration. Unit cell dimensions were determined at the ice point (273.15 K) for different compounds like propane, isobutene, and cyclopropane. For all the compounds, interactions between pairs of sites within a cutoff radius of 8.5125 Å were included. The Lennard-Jones 12-6 potential was used for representing the interaction between the gas and water molecules as a function of distance between their centers. Interactions between these unlike molecules were approximated using the Lorentz-Berthelot mixing rules [15]. Variation in cell sizes with different guest molecules suggested that each compound distorted the empty hydrate cavity to a different degree, the very foundation of lattice distortion. Using the same unit cell, all the guest molecules were removed and the distorted empty hydrate lattice (empty hydrate unit cell) was simulated at 273.15 K and zero pressure to obtain the reference dissociation energy. This value is used to calculate the reference chemical potential difference. The reference chemical potential difference is obtained as a difference of chemical potential between the occupied and the distorted empty gas hydrate at the ice point.
4. Results and Discussion
Various calculations have been performed over the years to establish the reference properties of gas hydrates through MD simulation techniques. As is shown in Figure 1, we kept the empty hydrate stabilized instantly and got the minimum energy.
A plot of total energy according to time in Figure 2 shows that the curve for empty hydrate crosses the curve for occupied gas hydrate and the pressure became stabilized at the stable structure after the point of intersection. For the calculation of reference chemical potential, the average value of the difference of the total energy of empty hydrate and that of occupied gas hydrate for the proceeding time steps is considered.
Using MD simulations, values for the chemical potential for different gas hydrates at the reference temperature (273.15 K) and pressure (0 kPa) were calculated for near equilibrium conditions. The results obtained from simulation are summarized in Table 1. Results are in close proximity with the experimental values obtained from Lee-Holder equation suggesting a variation in values for different guests, thus supporting the lattice distortion theory. The comparison between Lee-Holder model and our calculation for structure II gas hydrates is in Figure 3.
The unit cell sizes of propane gas hydrate, isobutane gas hydrate, and cyclopropane gas hydrates were observed for a temperature range of 240–278 K as shown in Figure 4. For a given temperature, the unit cell dimension varied till the equilibrium pressure was the same as the experimental pressure. From the equilibrium data obtained from MD simulations, the volume of the unit cell is plotted for different temperatures. It can be seen that the volume of the unit cell increases with temperature. The equilibrium calculations using the optimized unit cell sizes given by MD simulations are in Figure 5, which shows a good agreement between experimental and simulation results.
5. Conclusions
Variations in reference values still exist and the accurate values—the key for better potential parameters [7]—remain ambiguous. A sensitivity analysis [16] has demonstrated the importance of accurate reference values and MD simulation could be one of the techniques to obtain those values. As an attempt to apply MD simulation to calculate the reference state of gas hydrate, we demonstrate lattice distortion of structure II gas hydrates. The reference chemical potential was generally found to increase with the size of the guest molecule. Temperature effect on the unit cell size, which will be used to calculate the enthalpy change due to temperature, has been observed.
Notations
| : | Kihara core radius parameter, pm |
| : | Langmuir constant, kPa−1 |
| : | fugacity coefficient |
| : | Boltzmann constant, erg.K−1 |
| : | pressure, kPa |
| : | gas constant |
| : | temperature, K |
| : | mol fraction of water |
| : | cell potential. |
Greek Letters
| : | molar enthalpy difference between the empty hydrate lattice and pure water at 273.15 K and 0 atm, J/mol |
| : | difference of chemical potential of water and theoretical empty hydrate at 273.15 K and 0 atm, J/mol |
| : | difference of chemical potential of water in the unoccupied hydrate lattice and in the water, J/mol |
| : | difference of chemical potential of water in the unoccupied hydrate lattice and the occupied hydrate, J/mol |
| : | Kihara intermolecular well-depth parameter, J−1 |
| : | chemical potential of water in the unoccupied hydrate lattice, J/mol |
| : | activity coefficient of water |
| : | ratio of small or large cavities to water molecules in a unit cell |
| : | fraction of -type cavities occupied by -type gas molecule |
| : | Kihara core-to-core distance parameter, pm. |