0 Introduction

Natural gas hydrate, also known as combustible ice, is a special crystalline solid compound composed of water molecules and gas molecules under a low temperature and a high pressure^[1]. It is considered to be one of the new energy sources in the 21st century. In addition, the existence of natural gas hydrate will seriously threaten the safe production of petroleum industry, it is necessary to study the decomposition of natural gas hydrate, and the development process of natural gas hydrate resources includes the decomposition process of solid hydrate. At present, the conventional hydrate decomposition methods mainly include heat injection method, depressurization method and chemical reagent injection method. However, at present, there is no great breakthrough in the development of natural gas hydrate due to the lack of in-depth research on the micro mechanism of hydrate formation and decomposition^[2].

The main component of natural gas hydrate is methane^[3], therefore, natural gas hydrate with more than 99% methane molecular is usually called methane hydrate. Generally, there are three main types of hydrate structures: sI, sII, and sH, depending on the size and properties of gas molecules called guest molecules. Among them, water molecules called host molecules form crystal cavities through hydrogen bond interaction, and guest molecules are filled in the crystal cavities, in which van der Waals force is dominant role between host and guest molecules^[4]. The structural diagrams of these three crystal forms are shown in Fig. 1.

Chemical reagents can be classified into thermodynamic inhibitor, kinetic inhibitor and anti-polymerization agent. In terms of their application, thermodynamic inhibitors cost too much due to large usage^[5], while kinetic inhibitors with a low concentration meet demands under the same conditions^[6]. However, some studies have shown that kinetic inhibitors with low mass fraction cannot inhibit hydrate formation, but accelerate the rate of hydrate formation. Therefore, thermodynamic inhibitors are still widely used at present. The anti-polymerization agent is characterized by preventing the generated gas from agglomerating, not forming gas hydrate in the pipeline^[7]. Alcohol inhibitors, including ethylene glycol and glycerol, are typical thermodynamic inhibitors. Li et al.^[8] studied the decomposition process of propane hydrate with or without methanol additive at different temperatures. The results showed that the hydroxyl group of methanol can form hydrogen bond with partial liquid water and destroy the original structure of hydrate. Zhou et al.^[9] studied the effect of ethanol concentration on methane hydrate decomposition under different conditions based on angular order parameters, radial distribution function, mean square displacement, diffusion coefficient and system energy. The results showed that hydrate decomposition increased with the increase of ethanol concentration. Florianne et al.^[10] studied the decomposition of hydrate under the action of two hydrate inhibitors such as methanol and ethylene glycol, and found that the inhibitor can form hydrogen bond with water molecules, caused in the competition between inhibitor and hydrate water molecules against free liquid water molecules.

Molecular dynamics simulation is a calculation method according to the principle of Newtonian mechanics. Its basic principle is to obtain the trajectory of two atoms in the physical system through Newton's classical mechanics, and then calculate the mechanical, thermodynamic and dynamic properties of the system by using the fixed statistical method^[11]. Sun et al.^[12] studied the decomposition process of methane hydrate in pure water and methanol solution by molecular dynamics, and discussed the effects of temperature and pressure on the decomposition rate. Their results showed that heating and depressurization could significantly improve the decomposition rate of methane hydrate. Dai et al.^[13] studied the effects of alcohol chain length, hydroxyl position and hydroxyl number on methane hydrate decomposition by means of molecular dynamics simulation, obtaining potential energy, radial distribution function and mean square displacement. The results showed that shortening alcohol chain length and increasing alcohol hydroxyl number were conducive to methane hydrate decomposition. Fan et al.^[14] studied the decomposition of hydrate under the action of ethylene glycol by molecular dynamics and found that the concentration of ethylene glycol aqueous solution has an important influence on the decomposition rate of hydrate. Sujith and others^[15] used classical molecular dynamics simulation to study the evolution of dissolved methane in hydrate in the presence of two thermodynamic hydrate inhibitors NaCl and CH₃OH. It was found that the increase of hydrate inhibitor concentration promoted the nucleation of methane nanobubbles in hydrate melt, and CH₃OH molecules gathered around the core, resulting in the decrease of surface tension at the interface between hydrate and water.

In this paper, the decomposition of methane hydrate under the action of ethylene glycol and glycerol aqueous solution will be simulated by molecular dynamics. The final conformations of ethylene glycol-hydrate system and glycerol-hydrate system at different temperatures, pressures, alcohol solution concentrations and ensembles are obtained. The distributions of mean square displacement and radial distribution function are analyzed to explore the effects of temperature, pressure and inhibitor concentration on hydrate decomposition.

1 Simulation System and Method 1.1 Construction of sI Methane Hydrate

In this simulation, sI methane hydrate model is adopted, and the hydrate model is built in Materials Studio software. Referring to the X-ray single crystal diffraction data in Table 1, the initial positions of oxygen atoms in sI methane hydrate and carbon atoms in methane molecules are determined^[16-18], The interaction of water molecules is described by SPC potential energy model. Due to the random distribution of hydrogen atom positions in the cell, the optimization of hydrogen atoms in the cell needs to be completed according to Bernal-Fowler rule^[18].

The formula of potential energy includes two terms, long-range electrostatic interaction and short-range Lennard-Jones interaction, as calculated by

$ U_{i j}\left(r_{i j}\right)=\sum\limits_{A \in i} \sum\limits_{B \in j} \frac{q_A q_B}{r_{A B}}+4 \varepsilon\left[\left(\frac{\sigma}{r_{O O}}\right)^{12}-\left(\frac{\sigma}{r_{O O}}\right)^6\right] $

(1)

where U_ij(r_ij) is the molecular pair interaction potential energy; A and B are respectively represented as A and B atoms cycling on i and j; q_A, q_B represent the charge of A and B atoms respectively; r_AB is the distance between atoms; r_OO is the interaction distance between oxygen atoms of two molecules; σ, ε is the Lennard-Jones action parameter of oxygen atom.

The hydrate structure with minimal intermolecular energy is the most stable, so the optimal structure is established based on the Forcite module. Initially, methane is located in the center of the water cage crystal cavity. Consistent Valence Force Field(CVFF)is selected due to its accuracy in calculating binding energy of the force field, and providing reasonable configuration energy and vibration frequency.

The CVFF potential function is expressed by:

$ \begin{aligned} E_{\mathrm{pot}}= & \sum\limits_b D_b\left[1-e^{-a\left(b-b_0\right)}\right]^2+\sum\limits_\theta H_\theta\left(\theta-\theta_0\right)^2+ \\ & \sum\limits_\phi H_\phi[1-s \cos n \phi]+\sum\limits_\chi H_\chi \chi^2+\\ & \sum\limits_b \sum\limits_{b^{\prime}} F_{b b^{\prime}}\left(b-b_0\right)\left(b^{\prime}-b_0^{\prime}\right)+\sum\limits_\theta \sum\limits_{\theta^{\prime}} F_{\theta \theta^{\prime}}\left(\theta-\theta_0\right) \cdot \\ & \left(\theta^{\prime}-\theta_0^{\prime}\right)+\sum\limits_b \sum\limits_\theta F_{b \theta}\left(b-b_0\right)\left(\theta-\theta_0\right)+ \\ & \sum\limits_\phi F_{\phi \theta \theta^{\prime}} \cos \left(\theta-\theta_0\right)\left(\theta^{\prime}-\theta_0^{\prime}\right)+\sum\limits_\chi \sum\limits_\chi F_{\chi \chi} \chi \chi^{\prime}+ \\ & \sum \varepsilon\left[\left(\frac{r^*}{r}\right)^{12}-2\left(\frac{r^*}{r}\right)^6\right]+\sum \frac{q_i q_j}{\varepsilon r_{i j}} \end{aligned} $

(2)

The first four terms describe the bond length stretching energy, bond angle bending energy, dihedral angle torsion energy and out of plane vibration energy, respectively, and together with the last five cross interaction terms form the expression of intramolecular interaction; the last two terms represent van der Waals interaction and electrostatic Coulomb interaction, indicating intermolecular interaction.

Methane hydrate is composed of water molecules and methane molecules. The intramolecular interactions mainly show the stretching of C-H and O-H bond lengths, the bending of O-H-O and C-H-C bond angles and their related coupling terms; The intermolecular interactions are mainly van der Waals force and electrostatic force between water molecules, between water molecules and methane molecules, and between methane molecules. Table 2 and Table 3 respectively show the interaction parameters of methane hydrate within and between molecules under CVFF force field.

In present study, Ewald method is applied to electrostatic summation. The charged amounts of hydrogen atom and oxygen atom in water molecule are + 0.41e and -0.82e, while the charged amounts of hydrogen atom and carbon atom in methane molecule are +0.1e and -0.4e, respectively. Optimized single crystal cell without hydrogen bond and with hydrogen bond are shown in Figs. 2(a) and (b), where the initial cell size is 22.38Å × 22.38Å × 22.38 Å, 64 CH₄ and 368 H₂O in the optimized 2×2×2 supercell.

There are two difficulties in computer simulation of charged particles. One is that the calculation amount is large because it cannot be cut off. On the other hand, the sum is conditional convergence. Ewald sum gives a solution: Disassemble this long-range force into two parts. One part decays rapidly when r is small, and the other part decays slowly with larger r, this is achieved by skillfully selecting an f(r).

$ \frac{1}{r}=\frac{f(r)}{r}+\frac{1-f(r)}{r} $

(3)

The first part f(r)/r is zero or can be ignored outside the truncation distance r_max, just like the short-range interaction force, it can be directly summed;

In the second part, (1-f(r))/r is a slowly varying function for all r, so its Fourier transform can be represented by several k-vectors of $|\vec{k}| \leqslant \boldsymbol{k}_{\max } $, this means that this part can be calculated efficiently in the reciprocal space.

After the linear superposition of these two parts, the solution of the required problem can be obtained. Ewald^[19] selected complementary error function (erf()) for summation:

$ \operatorname{erf}(x)=\frac{1}{\sqrt{{\rm{ \mathsf{ π} }}}} \int_{-x}^x e^{-t^2} \mathrm{~d} t=\frac{2}{\sqrt{{\rm{ \mathsf{ π} }}}} \int_0^x e^{-t^2} \mathrm{~d} t $

which results in the well-known Ewald formula for the electrostatic energy of the box:

$ E=E^{(r)}+E^{(k)}+E^{(s)}+E^{(d)} $

(4)

where E^(r) is real space energy, E^(k) is reciprocal space energy, E^(s) is self-energy and E^(d) is di-pole correction energy.

1.2 Construction of Simulation System

The application concentration of thermodynamic inhibitor is generally 20wt%-60wt%. Five different concentrations of ethylene glycol aqueous solution and glycerol aqueous solution were selected. The solution system parameters are shown in Tables 4 and 5. The layered structures of methane hydrate and two alcohol solutions are shown in Fig. 3, Fig. 3(a) is ethylene glycol-hydrate system and (b) is glycerol-hydrate system, where the lower layer is the super cell of sI methane hydrate 2 × 2 × 2, containing 64 CH₄ molecules and 368 H₂O molecules, and the upper layer is alcohol solution at different concentrations. The SPC model (water rigid model) is selected for the water molecular force field^[20]. Lennard-Jones potential energy model is used to describe the interaction between methane molecules, water molecules and methane molecules^[21-22]. In the CVFF force field, the charge of oxygen atom in hydrate is -0.82e, and the charge of carbon atom in methane molecule is -0.4e. In the next step, the established model is geometrically optimized by the smart algorithm, namely the steepest descent method combined with conjugate gradient method, and the process of minimizing the energy of the system is done. Fig. 4 represents the energy change and convergence diagram of geometric optimization of ethylene-glycol hydrate and glycerol-hydrate systems. It can be seen from the figure that the energy is basically unchanged after iteration to step 350, the convergence degree is good, and the energy is minimized.

The iterative formula of the steepest descent method is:

(1) Pick x⁽⁰⁾∈Rⁿ;

(2) For k=0, 1, 2, ……;

$ \begin{gathered} r^{(k)}=b-A x^{(k)}, a_k=\frac{\left(r^{(k)}, r^{(k)}\right)}{\left(A r^{(k)}, r^{(k)}\right)} \\ x^{(k+1)}=x^{(k)}+a_k r^{(k)} \end{gathered} $

(5)

(3) When ‖x^(k+1)-x^(k)‖ < ε, the iteration is terminated.

The above algorithm format is simple, and its search direction p^(k) is the residual direction r^(k). The entire iterative process only requires some matrices and vectors for algebraic operation.

The iterative formula of conjugate gradient method is:

(1) Pick x⁽⁰⁾∈Rⁿ, r⁽⁰⁾=b-Ax⁽⁰⁾;

(2) For k=1, 2, ……;

$ \begin{gathered} \beta_{k-1}=-\frac{\left(r^{(k)}, A p^{(k-1)}\right)}{\left(p^{(k-1)}, A p^{(k-1)}\right)} \\ p^{(k)}=r^{(k)}+\beta_{k-1} p^{(k-1)} \\ a_k=\frac{\left(b, p^{(k)}\right)}{\left(A p^{(k)}, p^{(k)}\right)}=\frac{\left(r^{(k)}, p^{(k)}\right)}{\left(A p^{(k)}, p^{(k)}\right)} \\ x^{(k+1)}=x^{(k)}+a_k p^{(k)} \\ r^{(k+1)}=r^{(k)}-a_k A p^{(k)} \end{gathered} $

(6)

(3) When ‖ x^(k+1) - x^(k) ‖ < ε, the iteration is terminated.

1.3 Simulation Method

Through investigating a large number of documents, it is found that the molecular dynamics method is feasible and widely used in this research. Therefore, this method is also applied in this paper. Molecular dynamics simulation is carried out in Materials Studio^[23] software. The time step is set to 1 ps, and the total simulation time is 800 ps. The temperature is set to 250, 273, 298 and 320 K, and the pressure is set to 15, 25, 35, 45 and 50 MPa. To study the effects of temperature and pressure on hydrate decomposition rate, canonical ensemble (NVT)^[24], the particle number and volume of each system are fixed and constant-pressure, constant-temperature (NPT)^[24] ensemble are selected respectively. In addition, NVT ensemble is selected to test the influence of the Alcohol solution concentration on hydrate decomposition at the temperature of 273 K and the pressure of 15 MPa. In previous studies, it was found that the decomposition effect of inhibitors on hydrate was usually analyzed by calculating the radial distribution function and mean square displacement of atoms. Therefore, in this paper, the final simulated conformation of the system was analyzed, and the MSD and RDF of atoms were calculated to discuss the results.

The radial distribution function (RDF) is given by:

$ g_{\alpha \beta}(r)=\frac{V}{N_\alpha N_\beta}\left(\sum\limits_{i=1}^{N_\alpha} \frac{n_{i \beta}(r)}{4 {\rm{ \mathsf{ π} }} r^2 \Delta r}\right) $

(7)

where V represents the volume of lattice volume; N_α and N_β are the number of α and β particles, respectively; n_iβ(r) is the number of β particles in the space between concentric spheres.

The mean square displacement (MSD) is calculated by:

$ \operatorname{MSD}=\left(\left|\bar{r}_i-\bar{r}_{i 0}\right|\right)^2=\frac{1}{N} \sum\limits_{i=1}^N\left(\bar{R}_i(t)-\bar{R}_i\left(t_0\right)\right)^2 $

(8)

where N represents the total number of particles; R_i-(t) represents the position of the particle at the t moment; R_i-(t₀) represents the original position of the particle.

2 Results and Discussions 2.1 Ethylene Glycol

Ethylene glycol (CH₂OH)₂ is abbreviated as EG. It is generally believed that ethylene glycol can decompose hydrate due to that its dihydroxy structure can destroy the hydrogen bond in hydrate. The molecular structure of ethylene glycol is shown in Fig. 5, and 30.4% ethylene glycol is selected.

2.1.1 Effect of temperature on hydrate decomposition rate

After 800 ps kinetic simulation, the conformation of the ethylene glycol system at different temperatures is shown in Fig. 6. It can be seen that ethylene glycol molecules can be attached to the hydrate surface and form the hydrogen bonds with the hydrate surface, which can prevent methane from entering the crystal cavity again to form methane hydrate. In addition, the high temperature facilitates hydrate decomposition, at this moment, methane molecules exist in the form of gas.

The mean square displacement (MSD) is the deviation of a particle's position relative to its reference position which is a measure of the activity level of particles. With regard to a solid system, an upper limit of the mean square displacement is available, while the mean square displacement is linear with time in a liquid system. The MSD diagram of oxygen atom at different temperatures is shown in Fig. 7. The profiles show that the slope becomes great with the increase of temperature, indicating that the methane molecules escape from the hole, and the system changes from solid phase to gas-liquid phase at a higher temperature, thereby the degree of decomposition is increased.

The RDFs of O-O, O-H, H-H, C-O, C-H and C-C at different temperatures are shown in Fig. 8.

In Fig. 8(a), g_O-O has an obvious characteristic peak located near 2.75 Å and the peak decreases with the increase of temperature, indicating that hydrate structures begin to be damaged, and dissociated hydrate crystals are available, after the first peak, subsequently, other small characteristic peaks gradually disappear, and the curves tend to be flat. It means that the hydrate has decomposed on the condition of the hydrate structure.

As shown in Fig. 8(b), there are three obvious characteristic peaks, located at 0.75, 1.75 and 3.25 Å, respectively, and their peak values decline with the increase of temperature. For Fig. 8(c), three obvious characteristic peaks of g_H-H are located near 1.25, 2.25 and 3.75 Å, and the corresponding peak values increase firstly, and then gradually decrease, while the curves of g_H-H tend to be gentle and their peak values are almost the same, meaning that the stability of methane hydrate is destroyed at a high temperature.

In Fig. 8(d), the unique characteristic peak value is around 3.75 Å. After the disappearance of the first peak, the rest of the curve shows an upward trend due to small hole of hydrate being wrapped by methane molecules. For g_C-H, only one characteristic peak locates near 3.75 Å. The curves show an upward trend with the increase of truncation radius after the first peak disappears. The curves in Figs. 8 (d) and (e) show an upward trend, and the reason is that the methane molecules escape from the crystal hole and then gather together.

The g_C-C shown in Fig. 8 (f) has a characteristic peak near 4.25 Å, and the peak increases gradually with the increase of temperature, indicating that methane molecules begin to gather.

2.1.2 Effect of pressure on hydrate decomposition rate

The conformations of ethylene glycol system under different pressures are shown in Fig. 9, and the hydrate decomposition degree becomes greater with the decrease of pressure, and the decomposition rate speeds up. Also, ethylene glycol molecules form hydrogen bonds with the hydrate surface.

The MSD diagram of the oxygen atom is shown in Fig. 10 and the slope of the curve is larger at 25 MPa pressure than at other pressures, indicating that the hydrate is dissociated.

After kinetic simulation, the radial distribution function of O-O and C-C under different pressures are analyzed, as shown in Fig. 11.

In Fig. 11(a), there is an obvious characteristic peak in the g_O-O near 2.75 Å at the pressures of 15 MPa and 25 MPa, their peaks are lower than those at other three pressures, indicating that the hydrate decomposes easily at low pressure, the lowest peak appears at 25 MPa. For Fig. 11(b), there is an obvious characteristic peak in the g_C-C diagram near 4.25 Å. The highest peak appears at the pressure of 25 MPa, demonstrating that methane molecules in the hydrate escape and gather from the crystal cavity, and the hydrate structure is damaged. According to Figs. 9 and 10, the optimal decomposition degree of ethylene glycol system is at the pressure of 25 MPa.

2.1.3 Effect of EG concentration on hydrate decomposition rate

The conformation of the model under different EG concentrations is shown in Fig. 12. It means that methane molecules show a relatively regular aggregation distribution in the solution at a low concentration of ethylene glycol. With the increase of concentration of alcohol solution, methane molecules change from aggregation to dispersion.

The MSD and RDF diagram of oxygen atom is shown in Fig. 13. The slopes of the curves are roughly the same with EG concentrations of 23.7% and 30.4%, and their hydrate decompositions are the same. However, the slope of the curve becomes smaller at the EG concentration greater than 30.4%, and their decomposition rate slows down. For g_O-O, there is an obvious characteristic peak near 2.75 Å, and the low concentration of EG will reduce the peak value, indicating that under this working condition, the two concentrations of EG have better decomposition effect on hydrate.

2.2 Glycerol

Glycerol(C₃H₈O₃), also known as glycerin. Research shows that it is easier to form new hydrogen bonds between hydrate water molecules and glycerol molecules, so as to prevent hydrate regeneration. The molecular structure of glycerol is shown in Fig. 14, and 30.7% glycerol is selected for this section.

2.2.1 Effect of temperature on hydrate decomposition rate

The conformation of the glycerol system at different temperatures is shown in Fig. 15. We can find that glycerol molecules can attach to the hydrate surface to associate hydrogen bonds with the hydrate surface, preventing methane from entering the crystal cavity again, and then forming methane hydrate. In addition, high temperature can promote the decomposition of hydrate, continually confirming that increasing temperature accelerates the decomposition of hydrate.

The MSD of oxygen atom at different temperatures is shown in Fig. 16. The slope of the MSD curve becomes steep with an increase of temperature, indicating that the decomposition rate is faster.

The radial distribution functions of O-O, O-H, H-H, C-O, C-H and C-C at different temperatures are shown in Fig. 17.

In Fig. 17(a), g_O-O has an obvious characteristic peak located near 2.75Å. The peak value becomes lower at a high temperature. After the first characteristic peak, the rest of the curve tends to be gentle, demonstrating that the hydrate has decomposed, the order of water molecules has been broken, and the hydrate has been destroyed.

As shown in Fig. 17(b), there are three obvious characteristic peaks at 0.75, 1.75 and 3.25 Å, respectively. With the increase of temperature, the peak value gradually decreases, the rest of the curve tends to be flat and the peaks gradually disappear. For Fig. 17(c), there are three obvious characteristic peaks at 1.25, 2.25 and 3.75 Å in the curve of g_H-H. When the temperature increases, both peaks decrease, and the other peaks of curve disappear gradually, meaning alcohol molecules form hydrogen bonds with the hydrate surface, resulting in the destruction of hydrate stability.

In Fig. 17(d), the unique characteristic peak value is around 3.75 Å, and the peak becomes lower and lower with the increase of temperature. In addition, except for the first peak, the trend of other curves is generally the same. From Fig. 17(e), the diagram of g_C-O is generally consistent. The curve shows an upward trend along with the distance from 8 Å to 15 Å, which is a sign of methane molecules escaping from the hole and gathering together.

Fig. 17(f) shows the g_C-C diagram, with a characteristic peak near 4.25 Å, which represents that hydrate methane molecules gradually escape from the hole and gather together. In addition, near 6 Å, the curve tends to be consistent, indicating that the arrangement of carbon atoms will not change greatly, and the system is in equilibrium state at this moment.

2.2.2 Effect of pressure on hydrate decomposition rate

The conformations of glycerol system under different pressures are shown in Fig. 18. It also can be seen that, similar to the decomposition mechanism of ethylene glycol, glycerol molecules can form hydrogen bonds with nearby liquid water molecules, destroy the cage structure, and adhere to the hydrate surface to prevent hydrate regeneration.

The MSD diagram of the oxygen atom is shown in Fig. 19. The slope of MSD curve is the steepest at the pressure of 25 MPa, indicating the best decomposition effect, while the slopes of curves at 35 MPa and 45 MPa are the same, so their degrees of decomposition are roughly the same. In addition, the slopes of curves are the lowest at pressures of 15 MPa and 50 MPa.

After kinetic simulation, the g_O-O and g_C-C under different pressures are analyzed, as shown in Fig. 20. In Fig. 20(a), there is an obvious characteristic peak in the g_O-O diagram near 2.75 Å, the peak decreases gradually with the decrease of pressure. In Fig. 20(b) g_C-C diagram, there is also an obvious characteristic peak near 4.25 Å, and the peak is the highest when the pressure is 25 MPa, indicating that under this pressure, the hydrate decomposition degree is the largest, the hydrate crystal is damaged, and the methane molecules escape and then gather together. According to Figs. 19 and 20, at the pressure of 25 MPa, the decomposition effect of hydrate with glycerol is remarkable.

2.2.3 Effect of glycerol concentration on hydrate decomposition rate

The conformation of the glycerol system under different concentrations is shown in Fig. 21. At the concentration of glycerol is 23.5% and 30.7% respectively, methane molecules are aggregated, and when the concentration increases, methane molecules begin to disperse.

The MSD and RDF diagram of oxygen atom is shown in Fig. 22. For Fig. 22(a), at the concentrations of 23.5%, 30.7% and 38.7%, the slopes of the curves are roughly the same, nevertheless, the slope of the curve becomes smaller at the glycerol concentration greater than 38.7%, meaning that their decomposition rate slows down. There is an obvious characteristic peak for g_O-O near 2.75 Å in Fig. 22(b). The peak value increases with the increase of inhibitor concentration, especially, the peak value increases significantly at the concentration exceeding 38.7%. According to Figs. 21 and 22, the hydrate decomposition effect performs excellent at the concentration between 23.5% and 38.7%.

2.3 Comparison of Ethylene Glycol and Glycerol

According to the previous simulation results, glycol has the best effect when the concentration is 23.7%-30.4%, and glycerol has the best effect when the concentration is 23.5%-38.7%. Therefore, when comparing the decomposition effect of the two inhibitors, under the NPT ensemble, the decomposition of hydrate is analyzed at the temperature of 298 K and the pressure of 15 MPa, the concentration of glycol and glycerol is 30.4% and 30.7% respectively. The simulation results are shown in Fig. 23.

In Fig. 24(a), under the same conditions, there is an obvious characteristic peak in g_O-O near 2.75 Å, and the peak value of glycerol system is the lowest. In addition, the MSD diagram of oxygen atoms is shown in Fig. 24(b). Its slope of the curve is almost 0 without inhibitor, indicating that the hydrate exists as a solid phase. While adding alcohol solution as inhibitor, the slope of curve shows an upward trend, indicating that the hydrate is dissociated. According to Figs. 23 and 24, the effect of glycerol shows more effective decomposition of hydrate than that of ethylene glycol, and the effects of the two inhibitors on hydrate decomposition under pressures of 25 MPa and 35 MPa are compared, the results show that the inhibition effect of glycerol is better than that of ethylene glycol, because glycerol owns more hydroxyl groups, the stronger polarity, the greater hydrophilic in water, which are beneficial to the association of hydrogen bonds with its adjacent water molecules to inhibit hydrate formation.

Under the same conditions, the decomposition effects of the two inhibitors at 25 MPa and 35 MPa are compared by calculating the radial distribution function and mean square displacement of the oxygen atom. It can be seen from Fig. 25, the peak value of the g_O-O in the glycerol system is lower than that in the ethylene glycol system, and the MSD slope of oxygen atom in the glycerol system is significantly higher than that in the ethylene glycol system. It is further verified that glycerol inhibitor is better than ethylene glycol in decomposing hydrate.

2.4 Comparison of Decomposition Effects of Hydrate in Two Different Ensembles

Under the same inhibitor concentration, temperature and pressure, the decomposition effects of ethylene glycol solution and glycerol solution under the two ensembles are simulated respectively. In this section, the inhibitor concentration is 38.7%, the temperature is 298 K and the pressure is 15 MPa. The simulation results are shown in Fig. 26, group (a) and group (b) are hydrate decomposition diagrams of ethylene glycol aqueous solution and glycerol aqueous solution under NVT and NPT ensemble.

From Fig. 27, the MSD images of O-O in both ensembles show an upward trend. The slope of the curve under NPT ensemble is significantly greater than that under NVT ensemble with time. The g_O-O diagram in Fig. 28 has roughly the same trend with an obvious characteristic peak near 2.75 Å. However, the peak value for NPT ensemble is significantly smaller than that of NVT ensemble, and the curve gradually flattens, meaning that the hydrate gradually changes from solid phase to gas-liquid phase. According to Figs. 26-28, the decomposition effect of hydrate under NPT ensemble is better than NVT ensemble, the reason is that the NPT ensemble can not only control temperature, but also control pressure, while the NVT ensemble no longer controls pressure after determining the box size.

3 Conclusions

In this paper, the decompositions of sI methane hydrate in ethylene glycol and glycerol aqueous solutions are explored by molecular dynamics simulation. The effects of temperature, pressure and concentration of two alcohol solutions on hydrate structure and decomposition rate are discussed, and the system without inhibitor serves as the control group.

The simulation results show that the decomposition process of hydrate generally includes three stages: the cage structure being destroyed, gas-liquid two-phase coexistence and gas-liquid separation. Adding glycol and glycerol aqueous solution, the hydrate structure is destroyed due to glycol and glycerol molecules attaching to the hydrate surface and forming hydrogen bonds with adjacent water molecules to prevent hydrate formation again.

The analysis on simulated conformation, RDF and MSD demonstrate that high temperature is conducive to hydrate decomposition. The decomposition effect of hydrate under the ethylene glycol and glycerol system is better at the pressure of 25 MPa and the inhibitor concentration ranging 20% - 40%, and the effect of glycerol decomposition hydrate is better than ethylene glycol. The decomposition rate of hydrate with NPT ensemble is better than that of NVT ensemble.

Therefore, the decomposition rate of hydrate can be improved by increasing temperature, reducing pressure, with appropriate concentrations of ethylene glycol and glycerol aqueous solutions.