Molecular dynamic simulation of the molecular characteristic and mechanical property of methane hydrate/ water/ ice mixture

The microscopic molecular characteristic will impact on the mechanical property of hydrate. Thus, molecular dynamics simulation is employed to investigate the molecular characteristic and mechanical property of methane hydrate/ water/ ice mixture system. The brittle fracture occurred during the tensile deformation of the system. Besides, the maximum stress of the hydrate/ water/ ice mixture system is lower than that of intact hydrate system. The fracture strain of studied system is smaller than that of pure hydrate system. The order parameters F3 and F4 can be used for determining the fracture position of mixture system and the changing of micro configuration on the mixture interface.


Introduction
Natural gas hydrate (hereinafter referred as "hydrate") is a kind of special crystalline compound, which is formed by natural gas molecules with water molecules under high-pressure and low-temperature condition [1]. It has great potential application in the fields of energy sources, gas storage and separation, desalination of sea water, etc [2][3][4][5]. Thus, the study of hydrate has been a hot topic for decades. The mechanical properties of hydrate attract certain interest for its impact on safe drilling, geological hazards, and climatic change [6]. Generally, the reservoir conditions of hydrate are changing with the exploitation process. Typically, the reducing of pressure and increasing of temperature would lead to the decomposition of hydrate [7]. So far, the scientists have performed various experiments on testing the mechanical properties of hydrate. Winters et al. [8] investigated the mechanical characteristics of water-bearing sediments. The results showed that the various pore fillers would affect the strength of the hydrate sediments. Ruppel et al. [9,10] tested the bulk modulus, stress-strain curve, Poisson's ratio, and shear strength of synthesized and natural hydrates. Kuniyuki et al. [11] found that the elastic modulus of the hydrate was sharply decreased with the decomposition of hydrates. Masayuki et al. [12] discussed the relationships of methane hydrate failure with temperature, pressure, and hydrate saturation through experiments. The results suggested that the shear strength and elastic modulus of the hydrate are enhanced with the increase in hydrate saturation.
However, the hydrate samples in the experiments contain impurities (e.g. sand sediments and ice) and defects (e.g. porous, vacancy and polycrystalline). Also, the hydrate experiments require special equipment to keep hydrate at a certain temperature and pressure [13]. The mechanisms of hydrate deformation, defect formation and crack propagation are still lack, especially for the hydrate mixture [14,15]. Therefore, theoretical approaches are proposed to study the mechanical properties of hydrates.
Molecular dynamics (MD) simulation is a powerful tool to investigate the micro-mechanism of materials [16][17][18][19][20]. And MD had been successful used to study the thermal and interfacial properties of hydrate [21,22]. Ning et al. [23] investigated the compressibility of CH4/CO2 hydrate mixtures and mechanical instability of mono/poly-crystalline methane hydrate. Their works revealed the mechanical behavior of hydrates from the molecular insight. Cai et al. [24] studied the microscopic structure and mechanical property of defect-containing sI methane hydrate. The simulations denoted that the presence of certain defects in the methane hydrate could promote its mechanical property. And the system mechanical property would be reduced when the defects had reached a certain degree.
Since the hydrate/ water/ ice mixture usually exist in the hydrate system [25], herein, MD is employed to examine the microscopic structure and mechanical property of methane sI hydrate/ water/ ice mixture. This work is also expected to provide useful insights in the structure of water mixture [26].

Simulation model
A complete sI type methane hydrate model is used to build the calculation model. The cage structure of each water molecule is occupied by a methane molecule, in which the atomic coordinates of methane hydrate are demarcated by X-ray experiments [27]. The system energy is minimized and the whole system is optimized after building the model, and the model is relaxed and stabilized at 200 K in NVT ensemble.
Simulations were performed in the LAMMPS software [28], periodic boundary conditions are applied to the X, Y, Z directions of the system, and the time step is 0.2 fs.
Then, the guest methane molecules on both sides of the complete lattice model are deleted, and the intermediate particles are relaxed at 200 K in the NVT ensemble. The system is heated from 300 K to 1200 K to dissolve the crystal structure on both sides into water， the heating process is carried out in 1,000,000 steps by using NVT ensemble. After that, the both ends of the methane molecule will be deleted, heated and then cooled under the NVT ensemble to make the structures on both sides melt into water. Later, the whole system is relaxed at 300 K, 10 MPa in NPT ensemble, and the relaxation time is 1,000,000 steps. 1 is the simulation model for 2 × 2 × 10 unit cells of sI methane hydrate deleting the methane molecules at both ends, in which the red particle is O atom of water molecule, the white particle is H atom of water molecule, and the brown particle is methane molecule. In the simulation, LJ hard sphere model is used to describe the interaction between methane molecules, and TIP4P model [29] is adopted to describe the physical and chemical properties of water molecules. The TIP4P model consists of four particle points, and the electrons and protons of the middle oxygen atom are not coincided. Lorentz-Berthelot mixing rule [30] is employed to simulate the interaction of different particles in the system.
The potential parameters of methane are listed in Table 1.  The computational parameters of TIP4P model for water are shown in Table 2. In TIP4P model, the potential energy of water molecule consists of two parts: one is LJ potential between oxygen atoms and the other is electrostatic potential between charges: The LJ potential energy function of water molecule is as follows: where, ε is the environmental dielectric constant, q is the quantity of electric charge of the ith molecule at position a, r is the distance between the position a of the ith molecule and the position b of the j th molecule, ε and σ are the parameters of LJ potential function between oxygen atoms, and r is the distance between the oxygen atom proton of the ith molecule and the oxygen atom proton of the jth molecule.
The electrostatic potential of water molecules can be expressed by where, N is the number of molecules in the system, M is the number of atoms in each molecule, r is the distance between atom a and atom b, and r is the truncation radius. The function f r can be expressed as: The relaxed model is stretched along the zth direction in the hydrate region (the engineering rate is 1 × 107 / s), the whole process adopts NPT ensemble at 300 K, and the pressure in X direction and Y direction is 10 MPa. The number of steps in the stretching process is 200,000. The engineering strain rate (erate) is expressed as follows:

Fig.2.
Snapshots of a 2×2×10 cell hydrate/ water/ ice mixture system during local stretching Fig. 2 displays the micro change process of a 2 × 2 × 10 cell hydrate/ water/ ice mixture system with local stretching. It can be found that when the deformation condition is applied to the middle hydrate part, the hydrate part has obvious extension in the Zth direction and the proportion of the hydrate and water mixed part in the whole structure increased. These indicate that the interaction between hydrate and water occurred. The fracture of the whole structure occurred near the junction of hydrate/water mixture, for which the fracture of hydrate structure accounted. Different from the whole stretching, the fracture position is near the hydrate/ water junction in terms of microstructure, which can be clearly determined that the fracture occurred in the hydrate structure. From the microstructure, it can be found that the fracture surface is uniform, which proves that the mechanical properties of hydrate remain unchanged although stretching results in the deformation of the hydrate part of the hydrate/ water mixture. Fig.3 is the strain-stress curves in Z axis of a 2×2×10 cell hydrate/ water/ ice mixture system with local stretching. The black curve is the stress-strain curve of hydrate/ water/ ice mixture, which goes through elastic deformation stage, plastic deformation stage and brittle fracture stage with the maximum stress value of 0.7 GPa. The red curve is the stress-strain curve of hydrate with intact lattice. It can be found that the maximum stress of hydrate/ water/ ice mixture is lower than that of intact hydrate, and that the fracture strain of hydrate / water / ice mixture is smaller than that of hydrate. Fig.3. The strain-stress curves in Z axis of the 2×2×10 hydrate/ water/ ice mixture system with local stretching

Order parameters for hydrate/ water/ ice mixture system with local stretching
Microscopically, the arrangement of water molecules should be based on some rules, which are usually expressed by order parameter. The micro configuration of dissolved water system is discussed by using the order parameters F3 and F4 respectively. F3 [31] is a parameter composed of three-body (three H2O) configuration, represents the degree of tetrahedral arrangement between any oxygen atom and other oxygen atoms around, and can be defined as： where φjik is the angle between the designated oxygen atom i and the surrounding two oxygen atoms (j and k), the surrounding area can be defined as 3.2 Å from the oxygen atom. For the crystal structure (hydrate and ice) formed by water, the oxygen atom forms a saturated hydrogen bond with the water molecule in which the four oxygen atoms are located, and F3 of the whole system is zero. However, F3 is about 1.0 for normal water. Fig.4 is F3 order parameter curves for the 2×2×10 hydrate/ water/ ice mixture system during local stretching. By comparing the position of water in hydrate/ water/ ice mixture, it can be found F3 is around 0, which indicates that the oxygen atoms arrangement of water molecule in hydrate/ water/ ice mixture is similar to the crystal structure. As seen from the output of 1,082 step, the order parameter F3has obvious fluctuations at the fracture of micro configuration, which is similar to that of pure hydrate before. Therefore, the fracture position of the stretched hydrate part of hydrate / water / ice mixture can be predicted by the order parameter F3. The order parameter F4 [32] is used to describe the order degree of torsion angle between two adjacent water molecules H-O…O-H in a microscopic water system: where,is the torsion angle between the oxygen atom in two adjacent water molecules and the hydrogen atom farthest away from another oxygen atom. That is,  is the torsion angle composed of HA-OA…OB-HB four atoms. Theoretically, for the hydrogen bond network of a standard hydrate configuration, F4 is close to 1. As a classical order parameter, F4 is widely used to describe the ordered structure of micro water molecules. In the simulations of hydrate, ice and normal water, the corresponding F4 are about 0.7, -0.3 and 0, respectively.  5 depicts the order parameter curves for the 2×2×10 hydrate/ water/ ice mixture system during local stretching. The order parameter F4 of hydrate molecules in the middle part of hydrate/ water/ ice mixture is still around 0.6-0.7, and the crystal structure is maintained. As shown in Fig. 5, the red curve represents the order parameter F4 when the output step is 1, and the crystal structure of the middle hydrate is still maintained. However, the order parameter F4 of the hydrate/ water/ ice binding part fluctuates in the range of -0.3-0.1, which indicates that the structure at the interface of hydrate and water/ice mixture is a mixed structure of the ice-like crystal structure and liquid water. This also confirms the conclusion that the interactions occurred at the interface, which is observed in the microcosmic configuration diagram of hydrate/ water/ ice mixture. By comparing the results of output step 1 and 1082, it can be found that the interface position of hydrate/ water/ ice mixture moves in the Zth direction with the tensile deformation.

Conclusions
In this paper, the change of microstructure and mechanical properties of hydrate/ water/ ice mixture is studied by comparing the local tensile deformation of hydrate. The tensile deformation of local hydrate in hydrate/ water/ ice mixture was simulated, and the microstructure change, stress-strain curve, order parameters F3 and F4 were analyzed. It can be found that the brittle fracture occurred during the tensile deformation of hydrate/ water/ ice mixture, and the maximum stress of the hydrate/ water/ ice mixture is lower than that of intact hydrate. The fracture strain of hydrate/ water/ ice mixture is smaller than that of hydrate. The order parameters F3 and F4 can respectively determine the fracture position of hydrate/ water/ ice mixture and the change of micro configuration on the hydrate/ water/ ice interface.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.