Molecular dynamics simulation of sI methane hydrate under compression and tension

Abstract Molecular dynamics (MD) analysis of methane hydrate is important for the application of methane hydrate technology. This study investigated the microstructure changes of sI methane hydrate and the laws of stress–strain evolution under the condition of compression and tension by using MD simulation. This study further explored the mechanical property and stability of sI methane hydrate under different stress states. Results showed that tensile and compressive failures produced an obvious size effect under a certain condition. At low temperature and high pressure, most of the clathrate hydrate maintained a stable structure in the tensile fracture process, during which only a small amount of unstable methane broke the structure, thereby, presenting a free-motion state. The methane hydrate cracked when the system reached the maximum stress in the loading process, in which the maximum compressive stress is larger than the tensile stress under the same experimental condition. This study provides a basis for understanding the microscopic stress characteristics of methane hydrate.


Introduction
Methane hydrate [1,2] (hereinafter referred as "hydrate") is a kind of ice-like clathrate crystalline, which is widely distributed in continental margin, continental slope and permafrost zone. Hydrate is regarded as a promising alternative energy source and has gradually attracted the attention of scholars. The study of hydrate has a potential application in the fields of energy source development [3][4][5], gas storage and separation [6][7][8] and desalination of sea water [9][10][11]. Geological hazards and environmental problems caused by gas exploration are simultaneously discussed [12,13]. Hydrate is formed by methane gas (guest molecule) and water (host molecule) under high pressure and low temperature. When small gas molecules exist, the water molecules bind to each other via hydrogen bond interactions; thereby forming a cage-shaped structure with different shapes and sizes with the gas molecule wrapped in it. Gas and water molecules interact with one another through the van der Waals force, which reduces the energy of the entire structure and achieves a stable state [14,15]. If the polyhedron cage-shaped hydrate lattice formed by the water molecule is not occupied by the guest molecule, then the empty hydrate lattice can be regarded as a special kind of ice, but it is unstable and prone to collapse. The filling of a guest molecule is helpful to the stability of the hydrate lattice, and the hydrate will be more stable when many guest molecules are filled.
Currently, various experimental measurements are performed to explore the mechanical characteristics of hydrate. However, experimental hydrate samples often have defects (e.g. porous, vacancy and polycrystalline) [16][17][18] or contain impurities, such as sand sediments and silicas, which result in errors. Also, keeping the hydrate at a certain temperature and pressure through experiment equipment is difficult. Therefore, the mechanisms of hydrate deformation, defect formation and crack propagation still lack systematic in depth research results, especially for the mechanical properties of pure hydrate from a microcosmic perspective. The micro-mechanisms of the mechanical properties of hydrate is a basic task in natural gas exploration and production. However, studies that have been carried out on it are limited. The properties of materials would vary with different sizes at micro/nanoscale. Thus, in this study, we investigate the size effect by molecular dynamics (MD). The normalized stretching force also affect the size relevant properties at nanoscale because the definition of force (pressure) have the parameter of length.
sI methane hydrate is a common type of hydrate, which is a unit cell of body-centred cubic (bcc) structure that contains 46 water molecules [19][20][21]. In this study, MD simulation is used to perform compression and tension of sI methane hydrate in different scales based on the experimental data. MD simulation is a useful method in investigating the micro-mechanisms of materials [22][23][24][25][26]. This method has been successfully applied in the study of the thermal and interfacial properties of hydrate [27,28]. Ning et al. [29,30] investigated the compressibility of CH 4 /CO 2 hydrate mixtures and mechanical instability of mono/poly-crystalline hydrate, recently, and revealed the mechanical behaviour of hydrate at the molecular level. This study is intended to investigate the changes in the microstructure of pure sI methane hydrate and the laws of stress-strain evolution under the condition of compression and tension by MD simulation. Thus, the mechanical characteristics and the microstructure evolution mechanism of sI methane hydrate under different stress conditions are explored to provide theoretical guidance for the practical exploitation of hydrate.

Simulation Model
Hydrate is the compound in which the hydrogen-bonded cages formed by water molecules hold the methane molecules. So far, the lattice structures of hydrate have three types, including sI (cubic) [31], sII face-centred cubic [32] and sH (hexagonal) [33] lattice structures. The sI hydrate was employed to carry out the compression and tension tests and build the simulation model in this study. Its unit cell is composed of 8 guest molecules (methane molecules) and 46 host molecules (water molecules) arranged in the cubic box of 12×12×12 Å 3 . Combined MD simulation with compression and tension tests, 2×2×4, 2×2×6, 2×2×8, 2×2×10, 2×2×12, 2×2×17 and 2×2×22 unit cells (X×Y×Z axis) of sI hydrate were constructed to investigate their dynamic characteristics under various compression and tension.
The methane molecules fully occupied the cage structures. Thus, the coordinates of the atoms of the sI hydrate unit cell were determined by X-ray experiments [34]. The MD simulations were performed using largescale atomic/molecular massively parallel simulator software [35][36][37]. The Lennard-Jones (LJ) potential and the TIP4P model [38,39] were adopted as the force field for methane and water because of its computational simplicity, respectively. Although TIP4P and SPCE models were widely used in describing water interactions, the LJ potential was usually adopted as the force field for methane. The 12-6 LJ potential equation [40,41] is expressed as 12 6 ( ) 4 where U is the potential energy; rij is the distance between two interplay particles i and j; and ε and σ are the energy and length scales for the interaction, respectively. The equation of the TIP4P model [42,43] is shown as Where q is the partial charge, and ε 0 is the dielectric constant. Figure 1 shows the molecular geometry of the TIP4P model for water.  [44], and the electrostatic interactions of the systems were calculated by using the Ewald sum method [45].

Computational Parameters
Initially, the simulation systems based on energy minimization were established, with the whole equilibrated at 200 K in NVT ensemble for 1,000,000 steps. Notably, in all of the simulations of the present study, periodic boundary conditions were applied in the X, Y and Z axes. Then, the time-step was set as 0.2 fs, and the temperature was controlled by Nose-Hoover algorithm. Next, the systems were equilibrated at 200 K and 10 MPa in NPT ensemble for 1,000,000 steps to produce the test samples, such as the 2×2×4 hydrate system shown in Figure 2, for compressing and stretching. Afterwards, the samples were used for the tension test to verify the simulation.
The stretching test was loading in the Z axis because the stress-strain relationships of sI hydrate were almost same in the X, Y and Z axes [30]. More specifically, the stretching test changed the Z dimension of the simulation box at a constant engineering strain rate of 1×10 7 /s at 200 K and 10 MPa. The tests lasted from 2,000,000 to 4,000,000 steps for different samples to break the hydrate. Meanwhile, the compression tests were performed using the original samples with the same constant engineering strain rate of 1×10 7 /s at 200 K and 10 MPa, which lasted 2,000,000 to 4,000,000 steps for different samples.
Ethical approval: The conducted research is not related to either human or animal use.

Strain-stress curves during stretching
Herein, the simulation systems are denoted as Systems 4, 6, 8, 10, 12, 17 and 22 for 2×2×4, 2×2×6, 2×2×8, 2×2×10, 2×2×12, 2×2×17 and 2×2×22 unit cells of sI methane hydrate, respectively. The stretching of the tension test is homogenous and isotropic based on the study of Wu et al. [30]. The stretching stress causes the cracks in the hydrate. Herein, e.g. 2×2×8 unit cells are shown in Figure 3, in which the strain of hydrate in the Z axis before the stretching is 0. As the stress was gradually loaded, the    The stress of the test region in the Z axis is the main stress because the hydrate deforms along the Z axis, as presented in Figure 4. Stress increased with strain, which was similar to elastic deformation. Then, the hydrate reached its maximum stress and cracked rapidly. This maximum value was less than that in the experiments [13], but close to that in the simulation results presented by Ning et al. [29,30]. Stress is calculated on the basis of the pressure tensor of the atoms in the Z direction. However, it did not correspond with practice that the hydrate system still had residual stress in the Z direction when the hydrate broke. Therefore, the stress of the strain-stress curves in Figure 4 (a) is zero when the length of the cracks exceeded 2 Å, considering the interactions of hydrogen bonds.
The hydrates synthesised by the process of experiments contained impurities; whereas the sI hydrate studied here was pure. Further, the maximum stress decreased with an increase in the size of hydrate in the Z axis from 4 to 10 cells (Figure 4 (b)). This finding indicates that the hydrate follows the size effect when the system size is under 10 cells. However, the maximum stress fluctuated when the size was over 10 cells. The maximum tensile stress of the hydrate was in the Z direction, and the cells at a size of 4, 6, 8 and 10 decreased, thereby demonstrating the existence of a size effect. When the hydrate size exceeded 10 cells, the maximum stress value was smaller than that before. Compared with the cells at a size of 4, 6, 8 and 10, no regularity, but value fluctuation, were found.

Strain-stress curves during compressing
The strain-stress curves in the Z axis during compressing are presented in Figure 5. Similar to the stretching process, the compressive stress increased with strain and reached its maximum with strain larger than that of the maximum stress during the stretching process and the maximum tensile stress in the corresponding hydrate system. The explanation is given as follows. The hydrate crystal is composed of a hydrogen bond network of water molecules, in which the hydrogen bond is attributed to the physical intermolecular interactions. Therefore, the physical intermolecular interactions are weaker than the intramolecular interactions, and thus the hydrogen bond is easy to break when stretching. However, when the hydrogen bond is compressing, the microstructure will be adjusted at a certain angle to increase the strain. The microstructure of hydrate will change to fit the compression. Compressive stress is large when the distance between particles is small.
The structure of the hydrates was destroyed after reaching the maximum compressive stress, and stress decreased sharply and, finally, reached stability when the compression system continued to be compressed. This phenomenon was attributed to the periodic boundary conditions in the MD simulation in this research. After  the hydrate was destroyed under pressure, the destroyed particle structure was still compressing until the system could not be compressed any more. In addition, the maximum stress of the hydrate under compression varied with the system size, which was similar to that of the tension. A certain size effect in the 4, 6, 8, 10-cell system was observed, and no regularity was found when the size exceeded 10 cells (Figure 5 (b)).

Conclusions
The stress characteristics of hydrate particles are vital in the investigation of the microscopic properties of hydrate particles. In terms of the mechanisms of hydrate deformation, defect formation and crack propagation are still inconclusive, especially for the pure hydrate from the perspective of mechanical properties microscopically. Moreover, the micro-mechanisms of hydrate's mechanical properties are scarce. In this study, the mechanical properties of hydrates are investigated by MD simulation, which was built by the classic sI methane hydrate. The compression of sI methane hydrate with different sizes and the relationship of stress-strain between the stretching and compressing process were discussed.
In the stretching and compressing process, similar to elastic deformation, stress increased with strain. Then, the hydrate reached its maximum stress and cracked rapidly. In the case of tensile fracture under low temperature and high pressure in a pure hydrate system, the cage-shaped structure of most hydrates was stabilized. Only a small amount of methane broke through and detached from the cage-shaped structure, which exists in a free-motion state. In the case of compressive fracture of this system, after the hydrate was destroyed under pressure, the destroyed particle structure was still compressing until the system could not be compressed any more. Under the 10-cell size (simulation condition), the tensile and compressive failures of the hydrate showed large correlation with the system size. In other words, obvious size effect was observed under certain conditions. When the hydrate size exceeded 10 cells, the maximum stress was smaller than before. Compared with the cells with the size of 4, 6, 8 and 10, no regularity, but value fluctuation, were found. The compressive stress reached its maximum, which was larger than that of tensile stress in the corresponding hydrate system. Furthermore, the strain value of the maximum compressive stress was also larger than that of the maximum tensile stress.
This study aims to investigate the change in microstructure of sI methane hydrate and the law of stress-strain evolution under the condition of stress and tensile tests at low temperature and high pressure. It not only expounded the mechanical characteristic mechanism of sI methane hydrate under a pressured state, but also microscopically obtained the law of stability of natural gas hydrate structure initially. Hopefully, a deepened study of micro-stress mechanism of hydrate particles would provide theoretical guidance and lay a basis for the efficient exploitation of natural gas hydrate in practical projects.

Acknowledgement:
The authors would like to acknowledge the colleagues from the State Key Laboratory of Coal Mine Disaster Dynamics and Control for their perspectives and suggestions related to experimental study and data analyses. This research is supported by the National Natural Science Foundation of China (Grant No. 51904040